test_GeometryDune.cxx File Reference

#include "larcore/Geometry/Geometry.h"
#include <string>
#include <iostream>
#include <fstream>
#include <sstream>
#include <iomanip>
#include "dunecore/ArtSupport/ArtServiceHelper.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "larcorealg/Geometry/Exceptions.h"

Go to the source code of this file.

Classes
struct	SpacePoint

Macros
#define	test_GeometryDune_CXX

Typedefs
typedef readout::TPCsetID	APAID
typedef unsigned int	Index
template<class T >
using	Vector = vector< T >
template<class T >
using	TwoVector = Vector< Vector< T >>
template<class T >
using	ThreeVector = Vector< TwoVector< T >>

Functions
template<class T , class V >
void	resize (Vector< T > &vec1, Index n1, const V &val)
template<class T , class V >
void	resize (TwoVector< T > &vec2, Index n1, Index n2, const V &val)
template<class T , class V >
void	resize (ThreeVector< T > &vec3, Index n1, Index n2, Index n3, const V &val)
void	setExpectedValues (ExpectedValues &ev)
void	setExpectedValuesSpacePoints (Geometry *)
int	test_GeometryDune (const ExpectedValues &ev, bool dorop, Index maxchanprint, bool useExistingFcl)
int	main (int argc, const char *argv[])

Variables
const Index	InvalidIndex = std::numeric_limits<Index>::max()

Macro Definition Documentation

#define test_GeometryDune_CXX

Definition at line 4 of file test_GeometryDune.cxx.

Typedef Documentation

typedef readout::TPCsetID APAID

Definition at line 48 of file test_GeometryDune.cxx.

typedef unsigned int Index

Definition at line 55 of file test_GeometryDune.cxx.

template<class T >

using ThreeVector = Vector<TwoVector<T>>

Definition at line 65 of file test_GeometryDune.cxx.

template<class T >

using TwoVector = Vector<Vector<T>>

Definition at line 62 of file test_GeometryDune.cxx.

template<class T >

using Vector = vector<T>

Definition at line 59 of file test_GeometryDune.cxx.

Function Documentation

int main	(	int	argc,
		const char *	argv[]
	)

Definition at line 541 of file test_GeometryDune.cxx.

                                       {
  ExpectedValues ev;
  setExpectedValues(ev);
  bool dorop = true;
  Index maxchanprint = 10;
  bool useExistingFcl = false;
  if ( argc > 1 ) {
    string sarg = argv[1];
    if ( sarg == "-h" ) {
      cout << argv[0] << ": . [dorop] [maxchan] [usefcl]" << endl;
      cout << argv[0] << ": [ChannelMapClass] [dorop] [maxchan] [usefcl]" << endl;
      cout << argv[0] << ": [Geometry/ChannelMapClass/sorter] [dorop] [maxchan] [usefcl]" << endl;
      cout << "    dorop: If true ROP mapping methods are tested [true]" << endl;
      cout << "  maxchan: Max # channels displayed for each ROP [10]" << endl;
      cout << "   usefcl: Take top-level fcl from current directory [false]" << endl;
      return 0;
    }
    string::size_type ipos = sarg.find("/");
    if ( ipos == string ::npos ) {
      if ( sarg != "." ) ev.chanmap = sarg;
    } else {
      ev.gname = sarg.substr(0, ipos);
      string::size_type jpos = sarg.find("/", ipos+1);
      if ( jpos == string::npos ) {
        ev.chanmap = sarg.substr(ipos+1);
      } else {
        ev.chanmap = sarg.substr(ipos+1, jpos-ipos-1);
        ev.sorter = sarg.substr(jpos+1);
      }
    }
  }
  if ( argc > 2 ) {
    string sarg = argv[2];
    dorop = sarg == "1" || sarg == "true" || sarg == ".";
  }
  if ( argc > 3 ) {
    istringstream ssarg(argv[3]);
    ssarg >> maxchanprint;
  }
  if ( argc > 4 ) {
    string sarg = argv[4];
    useExistingFcl = sarg == "1" || sarg == "true";
  }
  test_GeometryDune(ev, dorop, maxchanprint, useExistingFcl);
  cout << "Tests concluded." << endl;
  return 0;
}

template<class T , class V >

void resize	(	Vector< T > &	vec1,
		Index	n1,
		const V &	val
	)

Definition at line 68 of file test_GeometryDune.cxx.

                                                     {
  vec1.resize(n1, val);
}

template<class T , class V >

void resize	(	TwoVector< T > &	vec2,
		Index	n1,
		Index	n2,
		const V &	val
	)

Definition at line 73 of file test_GeometryDune.cxx.

                                                                  {
  vec2.resize(n1);
  for ( Vector<T>& vec1 : vec2 ) {
    vec1.resize(n2, val);
  }
}

template<class T , class V >

void resize	(	ThreeVector< T > &	vec3,
		Index	n1,
		Index	n2,
		Index	n3,
		const V &	val
	)

Definition at line 81 of file test_GeometryDune.cxx.

                                                                              {
  vec3.resize(n1);
  for ( TwoVector<T>& vec2 : vec3 ) {
    vec2.resize(n2);
    for ( Vector<T>& vec1 : vec2 ) {
      vec1.resize(n3, val);
    }
  }
}

void setExpectedValues

(

ExpectedValues &

ev

)

Definition at line 14 of file test_GeometryDune10kt.cxx.

                                           {
  ev.gname ="dune10kt_workspace_geo";
  ev.fullname = "dune10kt_v1_workspace";
  // Geometry counts.
  ev.ncry = 1;
  ev.ntpc = 8;
  ev.npla = 3;
  ev.napa = 4;
  ev.nrop = 4;
  // Signal type and view for each TPC plane.
  resize(ev.sigType, ev.ncry, ev.ntpc, ev.npla, geo::kMysteryType);
  resize(ev.view, 1, 8, 3, geo::kUnknown);
  for ( Index icry=0; icry<ev.ncry; ++icry ) {
    for ( Index itpc=0; itpc<ev.ntpc; ++itpc ) {
      for ( Index ipla=0; ipla<ev.npla; ++ipla ) {
        ev.sigType[icry][itpc][ipla] = (ipla < 2) ? geo::kInduction : geo::kCollection;
      }
      ev.view[icry][itpc][0] = geo::kU;
      ev.view[icry][itpc][1] = geo::kV;
      ev.view[icry][itpc][2] = geo::kZ;
    }
  }
  // Wire count for each TPC plane.
  resize(ev.nwirPerPlane, ev.ntpc, ev.npla, 0);
  ev.nwirPerPlane[0][0] = 1149;
  ev.nwirPerPlane[0][1] = 1148;
  ev.nwirPerPlane[0][2] =  480;
  for ( Index itpc=1; itpc<ev.ntpc; ++itpc ) {
    for ( Index ipla=0; ipla<ev.npla; ++ipla ) ev.nwirPerPlane[itpc][ipla] = ev.nwirPerPlane[0][ipla];
  }
  // Channel count per ROP.
  resize(ev.nchaPerRop, ev.nrop, 0);
  ev.nchaPerRop[0] = 800;
  ev.nchaPerRop[1] = 800;
  ev.nchaPerRop[2] = 480;
  ev.nchaPerRop[3] = 480;
  for ( Index irop=0; irop<ev.nrop; ++irop ) ev.nchaPerApa += ev.nchaPerRop[irop];
  ev.nchatot = ev.napa*ev.nchaPerApa;
  resize(ev.chacry, ev.nchatot, InvalidIndex);
  resize(ev.chaapa, ev.nchatot, InvalidIndex);
  resize(ev.charop, ev.nchatot, InvalidIndex);
  Index icha = 0;
  for ( Index icry=0; icry<ev.ncry; ++icry ) {
    for ( Index iapa=0; iapa<ev.napa; ++iapa ) {
      for ( Index irop=0; irop<ev.nrop; ++irop ) {
        for ( Index kcha=0; kcha<ev.nchaPerRop[irop]; ++kcha ) {
          ev.chacry[icha] = icry;
          ev.chaapa[icha] = iapa;
          ev.charop[icha] = irop;
          ++icha;
        }
      }
    }
  }
  resize(ev.firstchan, ev.ncry, ev.napa, ev.nrop, raw::InvalidChannelID);
  Index chan = 0;
  for ( Index icry=0; icry<ev.ncry; ++icry ) {
    for ( Index iapa=0; iapa<ev.napa; ++iapa ) {
      for ( Index irop=0; irop<ev.nrop; ++irop ) {
        ev.firstchan[icry][iapa][irop] = chan;
        chan += ev.nchaPerRop[irop];
      }
    }
  }
  // Space points.
  #include "setWorkspaceSpacePoints.dat"
  // Optical detectors.
  ev.nopdet = 10*ev.napa;
  ev.nopdetcha.resize(ev.nopdet, 1);
  ev.nopcha = ev.nopdet;
  resize(ev.opdetcha, ev.nopdet, 1, 0);
  for ( Index iopt=0; iopt<ev.nopdet; ++iopt ) {
    Index icha = iopt;
    for ( Index ioch=0; ioch<ev.nopdetcha[iopt]; ++ ioch ) {
      ev.opdetcha[iopt][ioch] = icha++;
    }
  }
}

void setExpectedValuesSpacePoints

(

Geometry *

)

Definition at line 95 of file test_GeometryDune10kt.cxx.

                                                  {
  const string myname = "setExpectedValuesSpacePoints: ";
  string ofname = "setExpectedValuesSpacePoints.dat";
  cout << myname << "Writing " << ofname << endl;
  const CryostatGeo& crygeo = pgeo->Cryostat(0);
  double origin[3] = {0.0};
  double crypos[3] = {0.0};
  crygeo.LocalToWorld(origin, crypos);
  double cxlo = crypos[0] - crygeo.HalfWidth();
  double cxhi = crypos[0] + crygeo.HalfWidth();
  double cylo = crypos[1] - crygeo.HalfHeight();
  double cyhi = crypos[1] + crygeo.HalfHeight();
  double czlo = crypos[2] - 0.5*crygeo.Length();
  double czhi = crypos[2] + 0.5*crygeo.Length();
  cout << "Cryostat limits: "
       << "(" << cxlo << ", " << cylo << ", " << czlo << "), "
       << "(" << cxhi << ", " << cyhi << ", " << czhi << ")" << endl;
  vector<double> zfs = {0.2, 0.3, 0.4, 0.5, 0.6, 0.7 };
  vector<double> yfs = {0.2, 0.5, 0.8};
  vector<double> xfs = {0.2, 0.3, 0.6 };
  ofstream fout(ofname.c_str());
  for ( double zf : zfs ) {
    double z = czlo + zf*(czhi-czlo);
    for ( double yf : yfs ) {
      double y = cylo + yf*(cyhi-cylo);
      for ( double xf : xfs ) {
        double x = cxlo + xf*(cxhi-cxlo);
        double xyz[3] = {x, y, z};
        cout << "  (" << x << ", " << y << ", " << z << ")" << endl;
        TPCID tpcid = pgeo->FindTPCAtPosition(xyz);
        unsigned int itpc = tpcid.TPC;
        if ( itpc == TPCID::InvalidID ) continue;
        const TPCGeo& tpcgeo = pgeo->TPC(tpcid);
        unsigned int npla = tpcgeo.Nplanes();
        fout << "  ev.posXyz.push_back(SpacePoint(" << x << ", " << y << ", " << z << "));" << endl;
        for ( unsigned int ipla=0; ipla<npla; ++ipla ) {
          PlaneID plaid(tpcid, ipla);
          double xwire = pgeo->WireCoordinate(x, y, plaid);
          fout << "  ev.posTpc.push_back(" << itpc << ");" << endl;
          fout << "  ev.posPla.push_back(" << ipla << ");" << endl;
          fout << "  ev.posWco.push_back(" << xwire << ");" << endl;
        }  // end loop over planes
      }  // end loop over x
    }  // end loop over y
  }  // end loop over z
}

int test_GeometryDune	(	const ExpectedValues &	ev,
		bool	dorop,
		Index	maxchanprint,
		bool	useExistingFcl
	)

Definition at line 204 of file test_GeometryDune.cxx.

                                           {
  const string myname = "test_GeometryDune: ";
  cout << myname << "Starting test" << endl;
#ifdef NDEBUG
  cout << myname << "NDEBUG must be off." << endl;
  abort();
#endif
  string line = "-----------------------------";

  cout << myname << line << endl;
  cout << myname << "   Geometry: " << ev.gname << endl;
  cout << myname << "Channel map: " << ev.chanmap << endl;
  cout << myname << "     Sorter: " << ev.sorter << endl;
  cout << myname << "     Do ROP: " << dorop << endl;

  cout << myname << line << endl;
  cout << myname << "Create configuration." << endl;
  if (useExistingFcl) {
    ArtServiceHelper::load_services("test_GeometryDune.fcl",
                                    ArtServiceHelper::FileOnPath);
  } else {
    std::stringstream config;
    config << "#include \"geometry_dune.fcl\"" << endl;
    config << "services.Geometry:                   @local::" << ev.gname << endl;
    config << "services.ExptGeoHelperInterface:     @local::dune_geometry_helper" << endl;
    if ( ev.chanmap.size() ) {
      config << "services.ExptGeoHelperInterface.ChannelMapClass: " << ev.chanmap << endl;
    }
    //config << "services.ExptGeoHelperInterface.GeoSorterClass: " << ev.sorter << endl;
    ArtServiceHelper::load_services(config);
  }

  cout << myname << line << endl;
  cout << myname << "Get Geometry service." << endl;
  art::ServiceHandle<geo::Geometry> pgeo;

  cout << myname << line << endl;
  check("Default wiggle", pgeo->DefaultWiggle());
  check("Geometry name", pgeo->DetectorName(), ev.fullname);
  cout << myname << "ROOT name: " << pgeo->ROOTFile() << endl;
  cout << myname << "GDML name: " << pgeo->GDMLFile() << endl;

  cout << myname << line << endl;
  double xlo, ylo, zlo;
  double xhi, yhi, zhi;
  cout << myname << "World box"  << endl;
  pgeo->WorldBox(&xlo, &ylo, &zlo, &xhi, &yhi, &zhi);
  check("  xlo", xlo);
  check("  ylo", ylo);
  check("  zlo", zlo);
  check("  xhi", xhi);
  check("  yhi", yhi);
  check("  zhi", zhi);

  cout << myname << line << endl;
  check("SurfaceY", pgeo->SurfaceY());

  cout << myname << line << endl;
  if ( 0  ) check("GetWorldVolumeName", pgeo->GetWorldVolumeName());
  if ( 0 ) check("TotalMass", pgeo->TotalMass());

  cout << myname << line << endl;
  check("CryostatHalfWidth", pgeo->CryostatHalfWidth());
  check("CryostatHalfHeight", pgeo->CryostatHalfHeight());
  check("CryostatLength", pgeo->CryostatLength());
  //check("", pgeo->());

  cout << myname << line << endl;
  Index ncry = pgeo->Ncryostats();
  check("Ncryostats", ncry, ev.ncry);
  check("MaxTPCs", pgeo->MaxTPCs(), ev.ntpc);
  check("MaxPlanes", pgeo->MaxPlanes(), ev.npla);
  check("TotalNTPC", pgeo->TotalNTPC(), ev.ntpc);
  check("Nviews", pgeo->Nviews(), ev.npla);
  check("Nchannels", pgeo->Nchannels(), ev.nchatot);

  cout << myname << line << endl;
  cout << "Check TPC wire plane counts." << endl;
  for ( Index icry=0; icry<ncry; ++icry ) {
    Index ntpc = pgeo->NTPC(icry);
    cout << "  Cryostat " << icry << " has " << ntpc << " TPCs" << endl;
    for ( Index itpc=0; itpc<ntpc; ++itpc ) {
      Index npla = pgeo->Nplanes(itpc, icry);
      cout << "    TPC " << itpc << " has " << npla << " planes" << endl;
      assert( npla == 3 );
      for ( Index ipla=0; ipla<npla; ++ipla ) {
        Index nwir = pgeo->Nwires(ipla, itpc, icry);
        Index ich1 = pgeo->Nchannels();
        Index ich2 = 0;
        for ( Index iwir=0; iwir<nwir; ++iwir ) {
          Index ich = pgeo->PlaneWireToChannel(ipla, iwir, itpc, icry);
          if ( ich < ich1 ) ich1 = ich;
          if ( ich > ich2 ) ich2 = ich;
        }
        cout << "      Plane " << ipla << " has " << setw(4) << nwir << " wires"
             << " readout on channels [" << ich1 << ", " << ich2 << "]"
             << endl;
        assert( nwir == ev.nwirPerPlane[itpc][ipla] );
      }
    }
  }

  cout << myname << line << endl;
  cout << "Check wire planes." << endl;
  const double piOver2 = 0.5*acos(-1.0);
  for ( PlaneID plaid : pgeo->IteratePlaneIDs() ) {
    const geo::PlaneGeo& gpla = pgeo->Plane(plaid);
    cout << "  Plane " << plaid << endl;
    cout << "    Signal type: " << pgeo->SignalType(plaid) << endl;
    cout << "           View: " << pgeo->View(plaid) << endl;
    cout << "     Wire angle: " << gpla.ThetaZ() - piOver2 << endl;
    assert( pgeo->SignalType(plaid) == ev.sigType[plaid.Cryostat][plaid.TPC][plaid.Plane] );
    assert( pgeo->View(plaid) == ev.view[plaid.Cryostat][plaid.TPC][plaid.Plane] );
  }

  cout << myname << line << endl;
  cout << "Check channel-wire mapping." << endl;
  Index itpc1Last = TPCID::InvalidID;
  Index ipla1Last = WireID::InvalidID;
  Index nprint = 0;
  TwoVector<Index> lastwire;
  resize(lastwire, ev.ntpc, ev.npla, 0);
  for ( Index itpc=0; itpc<ev.ntpc; ++itpc )
    for ( Index ipla=0; ipla<ev.npla; ++ipla )
      lastwire[itpc][ipla] = 0;
  for ( Index icha=0; icha<ev.nchatot; ++icha ) {
    vector<WireID> wirids = pgeo->ChannelToWire(icha);
    assert( wirids.size() > 0 );
    WireID wirid1 = wirids[0];
    Index itpc1 = wirid1.TPC;
    Index iapa1 = itpc1/2;
    Index ipla1 = wirid1.Plane;
    if ( itpc1 != itpc1Last || ipla1 != ipla1Last ) nprint = 0;
    itpc1Last = itpc1;
    ipla1Last = ipla1;
    bool print = nprint < maxchanprint;
    if ( print ) ++nprint;
    if ( print ) cout << "  Channel " << setw(4) << icha << " has " << wirids.size() << " wires:";
    ostringstream sposs;
    sposs.precision(2);
    for ( WireID wirid : wirids ) {
      Index itpc = wirid.TPC;
      Index iapa = itpc/2;
      Index ipla = wirid.Plane;
      Index iwir = wirid.Wire;
      if ( print ) cout << " " << itpc << "-" << ipla << "-" << setw(3) << iwir;
      if ( iwir > lastwire[itpc][ipla] ) lastwire[itpc][ipla] = iwir;
      assert( iapa == iapa1 );
      assert( ipla == ipla1 );
      checkval("\nPlaneWireToChannel", pgeo->PlaneWireToChannel(wirid), icha );
      assert( pgeo->SignalType(icha) == ev.sigType[wirid.Cryostat][wirid.TPC][wirid.Plane] );
      checkval("View", pgeo->View(icha), ev.view[wirid.Cryostat][wirid.TPC][wirid.Plane]);
      const WireGeo* pwg = pgeo->WirePtr(wirid);
      TVector3 p1 = pwg->GetStart();
      TVector3 p2 = pwg->GetEnd();
      if ( sposs.str().size() ) sposs << ", ";
      sposs << "(" << setw(7) << std::fixed << p1.x() << ", "
                   << setw(7) << std::fixed << p1.y() << ", "
                   << setw(7) << std::fixed << p1.z() << ") - ";
      sposs << "(" << setw(7) << std::fixed << p2.x() << ", "
                   << setw(7) << std::fixed << p2.y() << ", "
                   << setw(7) << std::fixed << p2.z() << ")";
    }
    if ( print ) cout << " " << sposs.str() << endl;
  }
  for ( Index itpc=0; itpc<ev.ntpc; ++itpc ) {
    for ( Index ipla=0; ipla<ev.npla; ++ipla ) {
      Index nwir = lastwire[itpc][ipla] + 1;
      cout << "  TPC-plane " << itpc << "-" << ipla << " has " << setw(3) << nwir << " wires" << endl;
      assert( nwir == ev.nwirPerPlane[itpc][ipla] );
    }
  }

  if ( dorop ) {
    cout << myname << line << endl;
    cout << "Check ROP counts and channels." << endl;
    check("MaxROPs", pgeo->MaxROPs(), ev.nrop);
    Index icry = 0;
    for ( CryostatID cryid: pgeo->IterateCryostatIDs() ) {
      Index napa = pgeo->NTPCsets(cryid);
      cout << "  Cryostat " << icry << " has " << napa << " APAs" << endl;
      assert( napa == ev.napa );
      for ( Index iapa=0; iapa<napa; ++iapa ) {
        APAID apaid(cryid, iapa);
        Index nrop = pgeo->NROPs(apaid);
        cout << "    APA " << iapa << " has " << nrop << " ROPs" << endl;
        assert( nrop == ev.nrop );
        for ( Index irop=0; irop<nrop; ++irop ) {
          ROPID ropid(apaid, irop);
          Index ncha = pgeo->Nchannels(ropid);
          Index icha1 = pgeo->FirstChannelInROP(ropid);
          Index icha2 = icha1 + ncha - 1;
          cout << "      ROP " << irop << " has " << ncha << " channels: ["
               << icha1 << ", " << icha2 << "]" << endl;
          assert( ncha == ev.nchaPerRop[irop] );
          assert( icha1 == ev.firstchan[icry][iapa][irop] );
        }
      }
      ++icry;
    }
    assert( icry == ncry );
    cout << myname << line << endl;
    cout << "Check channel-ROP mapping." << endl;
    icry = 0;
    for ( Index icha=0; icha<ev.nchatot; ++icha ) {
      ROPID ropid = pgeo->ChannelToROP(icha);
      Index icry = ropid.Cryostat;
      Index iapa = ropid.TPCset;
      Index irop = ropid.ROP;
      assert( icry == ev.chacry[icha] );
      assert( iapa == ev.chaapa[icha] );
      assert( irop == ev.charop[icha] );
    }
    cout << myname << line << endl;
    cout << "Check ROP-TPC mapping." << endl;
    for ( CryostatID cryid: pgeo->IterateCryostatIDs() ) {
      Index napa = pgeo->NTPCsets(cryid);
      cout << "  Cryostat " << icry << " has " << napa << " APAs" << endl;
      assert( napa == ev.napa );
      for ( Index iapa=0; iapa<napa; ++iapa ) {
        APAID apaid(cryid, iapa);
        Index nrop = pgeo->NROPs(apaid);
        cout << "    APA " << iapa << " has " << nrop << " ROPs" << endl;
        assert( nrop == ev.nrop );
        for ( Index irop=0; irop<nrop; ++irop ) {
          ROPID ropid(apaid, irop);
          std::vector<geo::TPCID> rtpcs = pgeo->ROPtoTPCs(ropid);
          Index nrtpc = rtpcs.size();
          assert( nrtpc > 0 );
          cout << "      ROP " << irop << " TPCs:";
          for ( Index irtpc=0; irtpc<nrtpc; ++irtpc ) {
            assert(rtpcs[irtpc].Cryostat == icry);
            cout << (irtpc > 0  ? "," : "") << " " << rtpcs[irtpc].TPC;
          }
          cout << endl;
        }
      }
    }
  } else {
    cout << myname << line << endl;
    cout << "Skipped APA and ROP tests." << endl;
  }  // end dorop

  cout << myname << line << endl;
  Index nspt = ev.posXyz.size();
  if ( nspt == 0 ) {
    cout << myname << "No space points found." << endl;
    cout << myname << "Creating space points." << endl;
    setExpectedValuesSpacePoints(&*pgeo);
  } else {
    cout << myname << "Check " << nspt << " space points." << endl;
    assert( ev.posPla.size() == ev.posTpc.size() );
    assert( ev.posWco.size() == ev.posTpc.size() );
    Index ires = 0;
    int w = 9;
    for ( Index ispt=0; ispt<nspt; ++ispt ) {
      double x = ev.posXyz[ispt].x;
      double y = ev.posXyz[ispt].y;
      double z = ev.posXyz[ispt].z;
      double xyz[3] = {x, y, z};
      TPCID tpcid = pgeo->FindTPCAtPosition(xyz);
      cout << "  (" << setw(w) << x << "," << setw(w) << y << "," << setw(w) << z << "): " << tpcid << endl;
      assert( tpcid.Cryostat != CryostatID::InvalidID );
      unsigned int itpc = tpcid.TPC;
      assert( itpc != TPCID::InvalidID );
      const TPCGeo& tpcgeo = pgeo->TPC(tpcid);
      unsigned int npla = tpcgeo.Nplanes();
      assert( npla == ev.npla );
      assert( npla > 0 );
      for ( unsigned int ipla=0; ipla<npla; ++ipla ) {
        assert( ires < ev.posTpc.size() );
        PlaneID plaid(tpcid, ipla);
        WireID wirid;
        try { wirid = pgeo->NearestWireID(xyz, plaid); }
        catch (geo::InvalidWireError const& e) {
          if (!e.hasSuggestedWire()) throw;
          cerr << "ERROR (non-fatal):\n" << e;
          wirid = e.suggestedWireID();
        }
        double xwire = pgeo->WireCoordinate(x, y, plaid);
        cout << "    TPC " << setw(2) << itpc << " plane " << ipla
             << " nearest wire is " << setw(3) << wirid.Wire
             << " and coordinate is " << xwire << endl;
        assert( itpc == ev.posTpc[ires] );
        assert( ipla == ev.posPla[ires] );
        checkfloat(xwire, ev.posWco[ires], 2.e-4);
        ++ires;
      }  // end loop over planes
    }  // end loop over space points
    cout << "  # checked planes: " << ires << endl;
  }
  // Optical detectors.
  bool dophot = true;
  bool doAssert = true;   // Flag to turn off assertions while testing the test
  if ( dophot ) {
    Index nopt = pgeo->NOpDets();
    check("# Optical detectors expected", nopt, ev.nopdet, true, doAssert);
    Index nopchaReadout = ev.nopcha;
    Index nopchaHardware = ev.nopchaHardware ? ev.nopchaHardware : nopchaReadout;
    cout << " Expected optical readout channel count: " << nopchaReadout << endl;
    cout << "Expected optical hardware channel count: " << nopchaHardware << endl;
    cout << "         Geometry optical channel count: " << pgeo->NOpChannels() << endl;
    check("# Optical channels", pgeo->NOpChannels(), nopchaHardware, true, doAssert);
    cout << "Per-detector optical channel counts:" << endl;
    assert( ev.nopdetcha.size() == nopt );
    bool doAssert = true;  // True enable assertions in the indvidual channel checks.
    for ( Index iopt=0; iopt<nopt; ++iopt ) {
      ostringstream sslab;
      sslab << "  # channels for optical detector " << iopt;
      //bool print = iopt < maxchanprint;
      check(sslab.str(), pgeo->NOpHardwareChannels(iopt), ev.nopdetcha[iopt], true, doAssert);
    }
    assert(doAssert);
    cout << " Opdet channels:" << endl;
    for ( Index iopt=0; iopt<nopt; ++iopt ) {
      Index noch = pgeo->NOpHardwareChannels(iopt);
      for ( Index ioch=0; ioch<noch; ++ioch ) {
        ostringstream sslab;
        sslab << "  Det " << iopt << ", chan " << ioch;
        Index icha = pgeo->OpChannel(iopt, ioch);
        check(sslab.str() + " (channel)", icha, ev.opdetcha[iopt][ioch], true, doAssert);
        if ( nopchaHardware == nopchaReadout ) {
          check(sslab.str() + " (OpDet)", pgeo->OpDetFromOpChannel(icha), iopt, true, doAssert);
        }
        check(sslab.str() + " (HardwareChannel)", pgeo->HardwareChannelFromOpChannel(icha), ioch, true, doAssert);
      }
    }
  }

  cout << myname << line << endl;
  cout << myname << "Done." << endl;
  return 0;
}

Variable Documentation

const Index InvalidIndex = std::numeric_limits<Index>::max()

Definition at line 56 of file test_GeometryDune.cxx.