Public Member Functions | Static Public Member Functions | Private Member Functions | Private Attributes | List of all members
geo::GeometryTestAlg Class Reference

Performs tests on the geometry as seen by Geometry service. More...

#include <GeometryTestAlg.h>

Public Member Functions

 GeometryTestAlg (fhicl::ParameterSet const &pset)
 
virtual ~GeometryTestAlg ()=default
 Virtual destructor. More...
 
virtual void Setup (geo::GeometryCore const &new_geo)
 Runs the test. More...
 
virtual unsigned int Run ()
 Runs the test, returns a number of errors (very unlikely!) More...
 

Static Public Member Functions

static std::array< double, 3 > GetIncreasingWireDirection (const geo::PlaneGeo &plane)
 Returns the direction on plane orthogonal to wires where wire number increases. More...
 

Private Member Functions

void printDetectorIntro () const
 
void printChannelSummary ()
 
void printVolBounds ()
 
void printDetDim ()
 
void printWirePos ()
 
void printWiresInTPC (const TPCGeo &tpc, std::string indent="") const
 
void printAllGeometry () const
 
void testFindVolumes ()
 
void testCryostat ()
 
void testTPC (geo::CryostatID const &cid)
 
void testPlaneDirections () const
 
void testWireOrientations () const
 
void testChannelToROP () const
 
void testChannelToWire () const
 
void testFindPlaneCenters ()
 
void testProject ()
 
void testPlaneProjectionReference () const
 
void testPlanePointDecompositionFrame () const
 
void testPlaneProjectionOnFrame () const
 
void testPlaneProjection () const
 
void testWireCoordFromPlane () const
 
void testParallelWires () const
 
void testPlanePointDecomposition () const
 
void testWireCoordAngle () const
 
void testWirePitch ()
 
void testInterWireProjectedDistance () const
 
void testPlanePitch ()
 
void testStandardWirePos ()
 
void testAPAWirePos ()
 
void testNearestWire ()
 
void testWireIntersection () const
 
void testThirdPlane () const
 
void testThirdPlane_dTdW () const
 
void testStepping ()
 
void testFindAuxDet () const
 
bool shouldRunTests (std::string test_name) const
 
unsigned int testFindWorldVolumes ()
 
unsigned int testFindCryostatVolumes ()
 
unsigned int testFindTPCvolumePaths ()
 
void printAuxiliaryDetectors () const
 Method to print the auxiliary detectors on screen. More...
 
template<typename Stream >
void printAuxDetGeo (Stream &&out, geo::AuxDetGeo const &auxDet, std::string indent, std::string firstIndent) const
 Prints information of an auxiliary detector into the specified stream. More...
 
template<typename Stream >
void printAuxDetGeo (Stream &&out, geo::AuxDetGeo const &auxDet, std::string indent="") const
 Prints information of an auxiliary detector into the specified stream. More...
 
template<typename Stream >
void printAuxDetSensitiveGeo (Stream &&out, geo::AuxDetSensitiveGeo const &auxDetSens, std::string indent, std::string firstIndent) const
 Prints information of the sensitive auxiliary detector into a stream. More...
 
template<typename Stream >
void printAuxDetSensitiveGeo (Stream &&out, geo::AuxDetSensitiveGeo const &auxDetSens, std::string indent="") const
 Prints information of the sensitive auxiliary detector into a stream. More...
 
bool CheckAuxDetAtPosition (double const pos[3], unsigned int expected) const
 Returns whether the auxiliary detector at pos is the expected one. More...
 
bool CheckAuxDetSensitiveAtPosition (double const pos[3], unsigned int expectedDet, unsigned int expectedSens) const
 Returns whether the auxiliary sensitive detector at pos is expected. More...
 
bool isWireAlignedToPlaneDirections (geo::PlaneGeo const &plane, geo::Vector_t const &wireDir) const
 Helper function for testWireIntersection(). More...
 
unsigned int testWireIntersectionAt (geo::TPCGeo const &TPC, TVector3 const &point) const
 Performs the wire intersection test at a single point. More...
 
std::vector< std::pair< geo::PlaneID, double > > ExpectedPlane_dTdW (std::array< double, 3 > const &A, std::array< double, 3 > const &B, const double driftVelocity=-0.1) const
 Returns dT/dW expected from the specified segment A-to-B. More...
 
unsigned int testThirdPlane_dTdW_at (std::vector< std::pair< geo::PlaneID, double >> const &plane_dTdW) const
 Performs the third plane slope test with a single configuration. More...
 

Private Attributes

geo::GeometryCore const * geom
 pointer to geometry service provider More...
 
bool fDisableValidWireIDcheck
 disable test on out-of-world NearestWire() More...
 
std::set< std::stringfNonFatalExceptions
 
std::vector< double > fExpectedWirePitches
 wire pitch on each plane More...
 
std::vector< double > fExpectedPlanePitches
 plane pitch on each plane More...
 
bool fComputeMass = true
 Whether to print the detector mass. More...
 
testing::NameSelector fRunTests
 test filter More...
 

Detailed Description

Performs tests on the geometry as seen by Geometry service.


Configuration parameters

Definition at line 109 of file GeometryTestAlg.h.

Constructor & Destructor Documentation

geo::GeometryTestAlg::GeometryTestAlg ( fhicl::ParameterSet const &  pset)
explicit

Definition at line 118 of file GeometryTestAlg.cxx.

119  : geom(nullptr)
120  , fDisableValidWireIDcheck( pset.get<bool>("DisableWireBoundaryCheck", false) )
121  , fExpectedWirePitches( pset.get<std::vector<double>>("ExpectedWirePitches", {}) )
122  , fExpectedPlanePitches( pset.get<std::vector<double>>("ExpectedPlanePitches", {}) )
123  , fComputeMass( pset.get("ComputeMass", true) )
124  {
125  // initialize the list of non-fatal exceptions
126  std::vector<std::string> NonFatalErrors(pset.get<std::vector<std::string>>
127  ("ForgiveExceptions", std::vector<std::string>()));
128  std::copy(NonFatalErrors.begin(), NonFatalErrors.end(),
129  std::inserter(fNonFatalExceptions, fNonFatalExceptions.end()));
130 
131  // initialize the list of tests to be run
132  //
133  // our name selector accepts everything by default;
134  // the default set skips the following:
135  fRunTests.AddToDefinition("default",
136  "-CheckOverlaps", "-ThoroughCheck", "-PrintWires"
137  );
138  fRunTests.ParseNames("@default"); // let's start from default
139 
140  // read and parse the test list from the configuration parameters (if any)
141  fRunTests.Parse(pset.get<std::vector<std::string>>("RunTests", {}));
142 
143  std::ostringstream sstr;
145 
146  mf::LogInfo("GeometryTestAlg") << "Test selection:" << sstr.str();
147 
148  } // GeometryTestAlg::GeometryTestAlg()
std::set< std::string > fNonFatalExceptions
std::vector< double > fExpectedWirePitches
wire pitch on each plane
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
void AddToDefinition(std::string set_name, NAMES...names)
Parses a list of names and add them to the specified definition.
Definition: NameSelector.h:127
void ParseNames(NAMES...names)
Parses a list of names and adds them to the selector.
Definition: NameSelector.h:97
void PrintConfiguration(std::ostream &) const
Prints the configuration into a stream.
std::vector< double > fExpectedPlanePitches
plane pitch on each plane
testing::NameSelector fRunTests
test filter
void Parse(LIST const &items)
Parses the names in the list and adds them to the selector.
Definition: NameSelector.h:77
bool fDisableValidWireIDcheck
disable test on out-of-world NearestWire()
T copy(T const &v)
geo::GeometryCore const * geom
pointer to geometry service provider
bool fComputeMass
Whether to print the detector mass.
virtual geo::GeometryTestAlg::~GeometryTestAlg ( )
virtualdefault

Virtual destructor.

Member Function Documentation

bool geo::GeometryTestAlg::CheckAuxDetAtPosition ( double const  pos[3],
unsigned int  expected 
) const
private

Returns whether the auxiliary detector at pos is the expected one.

Definition at line 3490 of file GeometryTestAlg.cxx.

3491  {
3492  unsigned int foundDet = std::numeric_limits<unsigned int>::max();
3493  try {
3494  foundDet = geom->FindAuxDetAtPosition(pos);
3495  }
3496  catch (cet::exception const& e) {
3497  mf::LogProblem("GeometryTestAlg")
3498  << "Caught an exception while looking for aux det around "
3499  << lar::dump::array<3U>(pos) << " (within aux det #" << expected
3500  << "); message:\n" << e.what();
3501  return false;
3502  }
3503  if (foundDet != expected) {
3504  mf::LogProblem("GeometryTestAlg")
3505  << "Auxiliary detector at position " << lar::dump::array<3U>(pos)
3506  << ", expected within aux det #" << expected << ", was returned to be "
3507  << foundDet << " instead";
3508  return false;
3509  }
3510  return true;
3511  } // GeometryTestAlg::CheckAuxDetAtPosition()
unsigned int FindAuxDetAtPosition(double const worldLoc[3], double tolerance=0) const
Returns the index of the auxiliary detector at specified location.
const char expected[]
Definition: Exception_t.cc:22
MaybeLogger_< ELseverityLevel::ELsev_error, true > LogProblem
const double e
static int max(int a, int b)
geo::GeometryCore const * geom
pointer to geometry service provider
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
bool geo::GeometryTestAlg::CheckAuxDetSensitiveAtPosition ( double const  pos[3],
unsigned int  expectedDet,
unsigned int  expectedSens 
) const
private

Returns whether the auxiliary sensitive detector at pos is expected.

Definition at line 3515 of file GeometryTestAlg.cxx.

3517  {
3518  size_t foundDet = std::numeric_limits<unsigned int>::max();
3519  size_t foundSensDet = std::numeric_limits<unsigned int>::max();
3520  try {
3521  geom->FindAuxDetSensitiveAtPosition(pos, foundDet, foundSensDet);
3522  }
3523  catch (cet::exception const& e) {
3524  mf::LogProblem("GeometryTestAlg")
3525  << "Caught an exception while looking for aux det sensitive around "
3526  << lar::dump::array<3U>(pos) << " (within aux det #" << expectedDet
3527  << " sensitive volume #" << expectedSens << "); message:\n" << e.what();
3528  return false;
3529  }
3530  if ((foundDet != expectedDet) || (foundSensDet != expectedSens)) {
3531  mf::LogProblem("GeometryTestAlg")
3532  << "Auxiliary detector at position " << lar::dump::array<3U>(pos)
3533  << ", expected within aux det #" << expectedDet
3534  << ", sensitive volume #" << expectedSens
3535  << ", was returned to be in aux det #"
3536  << foundDet << " sensitive volume #" << foundSensDet << " instead";
3537  return false;
3538  }
3539  return true;
3540  } // GeometryTestAlg::CheckAuxDetAtPosition()
void FindAuxDetSensitiveAtPosition(geo::Point_t const &point, std::size_t &adg, std::size_t &sv, double tolerance=0) const
Fills the indices of the sensitive auxiliary detector at location.
MaybeLogger_< ELseverityLevel::ELsev_error, true > LogProblem
const double e
static int max(int a, int b)
geo::GeometryCore const * geom
pointer to geometry service provider
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
std::vector< std::pair< geo::PlaneID, double > > geo::GeometryTestAlg::ExpectedPlane_dTdW ( std::array< double, 3 > const &  A,
std::array< double, 3 > const &  B,
const double  driftVelocity = -0.1 
) const
private

Returns dT/dW expected from the specified segment A-to-B.

Definition at line 2961 of file GeometryTestAlg.cxx.

2965  {
2966  // This function returns a list of entries, one for each plane:
2967  // - plane ID
2968  // - slope of the projection of AB from the plane, in dt/dw units
2969 
2970  // find which TPC we are taking about
2971  geo::TPCID tpcid = geom->FindTPCAtPosition(A.data());
2972 
2973  if (!tpcid.isValid) {
2974  throw cet::exception("GeometryTestAlg")
2975  << "ExpectedPlane_dTdW(): can't find any TPC containing point A ("
2976  << A[0] << "; " << A[1] << "; " << A[2] << ")";
2977  }
2978 
2979  if (geom->FindTPCAtPosition(B.data()) != tpcid) {
2980  throw cet::exception("GeometryTestAlg")
2981  << "ExpectedPlane_dTdW(): point A ("
2982  << A[0] << "; " << A[1] << "; " << A[2] << ") is in "
2983  << std::string(tpcid)
2984  << " while point B (" << B[0] << "; " << B[1] << "; " << B[2]
2985  << ") is in " << std::string(geom->FindTPCAtPosition(B.data()));
2986  }
2987 
2988  geo::TPCGeo const& TPC = geom->TPC(tpcid);
2989 
2990  // conversion from X coordinate to tick coordinate
2991  double dT_over_dX = 1. / driftVelocity;
2992  switch (TPC.DriftDirection()) {
2993  case geo::kPosX:
2994  // if the drift direction is toward positive x, higher x will reach the
2995  // plane earlier and have smaller t, hence the flip in sign
2996  dT_over_dX = -dT_over_dX;
2997  break;
2998  case geo::kNegX: break;
2999  default:
3000  throw cet::exception("InternalError")
3001  << "GeometryTestAlg::ExpectedPlane_dTdW(): drift direction #"
3002  << ((int) TPC.DriftDirection()) << " of " << std::string(tpcid)
3003  << " not supported.\n";
3004  } // switch drift direction
3005 
3006  const unsigned int nPlanes = TPC.Nplanes();
3007  std::vector<std::pair<geo::PlaneID, double>> slopes(nPlanes);
3008 
3009  for (geo::PlaneID::PlaneID_t iPlane = 0; iPlane < nPlanes; ++iPlane) {
3010  geo::PlaneID pid(tpcid, iPlane);
3011 // const double wA = geom->WireCoordinate(A[1], A[2], pid);
3012 // const double wB = geom->WireCoordinate(B[1], B[2], pid);
3013  const double wA
3014  = geom->WireCoordinate(geo::vect::makeFromCoords<geo::Point_t>(A), pid);
3015  const double wB
3016  = geom->WireCoordinate(geo::vect::makeFromCoords<geo::Point_t>(B), pid);
3017 
3018  slopes[iPlane]
3019  = std::make_pair(pid, ((B[0] - A[0]) * dT_over_dX) / (wB - wA));
3020 
3021  } // for iPlane
3022 
3023  return slopes;
3024  } // GeometryTestAlg::ExpectedPlane_dTdW()
geo::Length_t WireCoordinate(double YPos, double ZPos, geo::PlaneID const &planeid) const
Returns the index of the nearest wire to the specified position.
std::string string
Definition: nybbler.cc:12
unsigned int Nplanes() const
Number of planes in this tpc.
Definition: TPCGeo.h:165
unsigned int PlaneID_t
Type for the ID number.
Definition: geo_types.h:473
The data type to uniquely identify a Plane.
Definition: geo_types.h:472
bool isValid
Whether this ID points to a valid element.
Definition: geo_types.h:211
Geometry information for a single TPC.
Definition: TPCGeo.h:38
Drift towards negative X values.
Definition: geo_types.h:162
geo::TPCID FindTPCAtPosition(double const worldLoc[3]) const
Returns the ID of the TPC at specified location.
The data type to uniquely identify a TPC.
Definition: geo_types.h:386
DriftDirection_t DriftDirection() const
Returns an enumerator value describing the drift direction.
Definition: TPCGeo.h:144
Drift towards positive X values.
Definition: geo_types.h:161
TPCGeo const & TPC(unsigned int const tpc=0, unsigned int const cstat=0) const
Returns the specified TPC.
geo::GeometryCore const * geom
pointer to geometry service provider
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
std::array< double, 3 > geo::GeometryTestAlg::GetIncreasingWireDirection ( const geo::PlaneGeo plane)
inlinestatic

Returns the direction on plane orthogonal to wires where wire number increases.

Definition at line 2089 of file GeometryTestAlg.cxx.

2090  {
2091  TVector3 IncreasingWireDir = plane.GetIncreasingWireDirection();
2092  // BUG the double brace syntax is required to work around clang bug 21629
2093  // (https://bugs.llvm.org/show_bug.cgi?id=21629)
2094  return
2095  {{ IncreasingWireDir.X(), IncreasingWireDir.Y(), IncreasingWireDir.Z() }};
2096  } // GeometryTestAlg::GetIncreasingWireDirection()
Vector GetIncreasingWireDirection() const
Returns the direction of increasing wires.
Definition: PlaneGeo.h:457
bool geo::GeometryTestAlg::isWireAlignedToPlaneDirections ( geo::PlaneGeo const &  plane,
geo::Vector_t const &  wireDir 
) const
private

Helper function for testWireIntersection().

Definition at line 2348 of file GeometryTestAlg.cxx.

2349  {
2350  /*
2351  * Returns `true` if `wireDir` is aligned with plane frame or wire direction
2352  */
2353 
2354  auto const isOrthogonalTo = [&wireDir](geo::Vector_t const& other)
2355  { return std::abs(geo::vect::dot(wireDir, other)) < 1.0e-5; };
2356 
2357  // we dislike wires aligned to plane frame:
2358  if (isOrthogonalTo(plane.WidthDir<geo::Vector_t>()))
2359  return true;
2360  if (isOrthogonalTo(plane.DepthDir<geo::Vector_t>()))
2361  return true;
2362 
2363  if (isOrthogonalTo(plane.GetIncreasingWireDirection<geo::Vector_t>()))
2364  return true;
2365 
2366  return false;
2367 
2368  } // isWireAlignedToPlaneDirections()
constexpr auto dot(Vector const &a, Vector const &b)
Return cross product of two vectors.
ROOT::Math::DisplacementVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Vector_t
Type for representation of momenta in 3D space.
Definition: geo_vectors.h:164
T abs(T value)
void geo::GeometryTestAlg::printAllGeometry ( ) const
private

Definition at line 567 of file GeometryTestAlg.cxx.

567  {
568  const unsigned int nCryostats = geom->Ncryostats();
569  mf::LogVerbatim("GeometryTest") << "Detector " << geom->DetectorName()
570  << " has " << nCryostats << " cryostats:";
571  for(unsigned int c = 0; c < nCryostats; ++c) {
572  const geo::CryostatGeo& cryostat = geom->Cryostat(c);
573  const unsigned int nTPCs = cryostat.NTPC();
574  mf::LogVerbatim("GeometryTest") << " cryostat #" << c << " at "
575  << lar::dump::vector3D(cryostat.GetCenter()) << " cm has "
576  << nTPCs << " TPC(s):";
577  for(unsigned int t = 0; t < nTPCs; ++t) {
578  const geo::TPCGeo& tpc = cryostat.TPC(t);
579  printWiresInTPC(tpc, " ");
580  } // for TPC
581  } // for cryostat
583  mf::LogVerbatim("GeometryTest") << "End of detector "
584  << geom->DetectorName() << " geometry.";
585  } // GeometryTestAlg::printAllGeometry()
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
geo::Point_t GetCenter() const
Returns the geometrical center of the cryostat.
Definition: CryostatGeo.h:124
auto vector3D(Vector3D const &v)
Returns a manipulator which will print the specified vector.
Definition: DumpUtils.h:301
Geometry information for a single TPC.
Definition: TPCGeo.h:38
Geometry information for a single cryostat.
Definition: CryostatGeo.h:43
unsigned int Ncryostats() const
Returns the number of cryostats in the detector.
std::string DetectorName() const
Returns a string with the name of the detector, as configured.
void printAuxiliaryDetectors() const
Method to print the auxiliary detectors on screen.
unsigned int NTPC() const
Number of TPCs in this cryostat.
Definition: CryostatGeo.h:181
CryostatGeo const & Cryostat(geo::CryostatID const &cryoid) const
Returns the specified cryostat.
void printWiresInTPC(const TPCGeo &tpc, std::string indent="") const
const TPCGeo & TPC(unsigned int itpc) const
Return the itpc&#39;th TPC in the cryostat.
Definition: CryostatGeo.cxx:93
geo::GeometryCore const * geom
pointer to geometry service provider
template<typename Stream >
void geo::GeometryTestAlg::printAuxDetGeo ( Stream &&  out,
geo::AuxDetGeo const &  auxDet,
std::string  indent,
std::string  firstIndent 
) const
private

Prints information of an auxiliary detector into the specified stream.

Definition at line 605 of file GeometryTestAlg.cxx.

609  {
610 
612 
613  std::array<double, 3U> center, normal;
614  auxDet.GetCenter(center.data());
615  auxDet.GetNormalVector(normal.data());
616 
617  out << firstIndent << "\"" << auxDet.Name()
618  << "\" centered at " << lar::dump::array<3U>(center)
619  << " cm, size ( " << (2.0 * auxDet.HalfWidth1());
620  if (coordIs.nonEqual(auxDet.HalfWidth1(), auxDet.HalfWidth2()))
621  out << "/" << (2.0 * auxDet.HalfWidth2());
622  out << " x " << (2.0 * auxDet.HalfHeight())
623  << " x " << auxDet.Length() << " ) cm"
624  << ", normal facing " << lar::dump::array<3U>(normal);
625  unsigned int nSensitive = auxDet.NSensitiveVolume();
626  switch (nSensitive) {
627  case 0: break;
628  case 1:
629  out << "\n" << indent << "with sensitive volume ";
631  std::forward<Stream>(out),
632  auxDet.SensitiveVolume(0U), indent + " ", ""
633  );
634  break;
635  default:
636  out << "\n" << indent
637  << "with " << auxDet.NSensitiveVolume() << " sensitive detectors:";
638  for (unsigned int iSens = 0; iSens < nSensitive; ++iSens) {
639  out << "\n" << indent << " [#" << iSens << "] ";
640  printAuxDetSensitiveGeo(std::forward<Stream>(out),
641  auxDet.SensitiveVolume(iSens), indent + " ", "");
642  } // for
643  break;
644  } // if sensitive detectors
645 
646  } // GeometryTestAlg::printAuxDetGeo()
Provides simple real number checks.
const double e
def center(depos, point)
Definition: depos.py:117
void printAuxDetSensitiveGeo(Stream &&out, geo::AuxDetSensitiveGeo const &auxDetSens, std::string indent, std::string firstIndent) const
Prints information of the sensitive auxiliary detector into a stream.
template<typename Stream >
void geo::GeometryTestAlg::printAuxDetGeo ( Stream &&  out,
geo::AuxDetGeo const &  auxDet,
std::string  indent = "" 
) const
inlineprivate

Prints information of an auxiliary detector into the specified stream.

Definition at line 195 of file GeometryTestAlg.h.

197  { printAuxDetGeo(std::forward<Stream>(out), auxDet, indent, indent); }
void printAuxDetGeo(Stream &&out, geo::AuxDetGeo const &auxDet, std::string indent, std::string firstIndent) const
Prints information of an auxiliary detector into the specified stream.
template<typename Stream >
void geo::GeometryTestAlg::printAuxDetSensitiveGeo ( Stream &&  out,
geo::AuxDetSensitiveGeo const &  auxDetSens,
std::string  indent,
std::string  firstIndent 
) const
private

Prints information of the sensitive auxiliary detector into a stream.

Definition at line 651 of file GeometryTestAlg.cxx.

655  {
656 
658 
659  std::array<double, 3U> center, normal;
660  auxDetSens.GetCenter(center.data());
661  auxDetSens.GetNormalVector(normal.data());
662 
663  out << firstIndent << "centered at " << lar::dump::array<3U>(center)
664  << " cm, size ( " << (2.0 * auxDetSens.HalfWidth1());
665  if (coordIs.nonEqual(auxDetSens.HalfWidth1(), auxDetSens.HalfWidth2()))
666  out << "/" << (2.0 * auxDetSens.HalfWidth2());
667  out << " x " << (2.0 * auxDetSens.HalfHeight())
668  << " x " << auxDetSens.Length() << " ) cm"
669  << ", normal facing " << lar::dump::array<3U>(normal);
670 
671  } // GeometryTestAlg::printAuxDetSensitiveGeo()
Provides simple real number checks.
const double e
def center(depos, point)
Definition: depos.py:117
template<typename Stream >
void geo::GeometryTestAlg::printAuxDetSensitiveGeo ( Stream &&  out,
geo::AuxDetSensitiveGeo const &  auxDetSens,
std::string  indent = "" 
) const
inlineprivate

Prints information of the sensitive auxiliary detector into a stream.

Definition at line 208 of file GeometryTestAlg.h.

212  {
214  (std::forward<Stream>(out), auxDetSens, indent, indent);
215  }
void printAuxDetSensitiveGeo(Stream &&out, geo::AuxDetSensitiveGeo const &auxDetSens, std::string indent, std::string firstIndent) const
Prints information of the sensitive auxiliary detector into a stream.
void geo::GeometryTestAlg::printAuxiliaryDetectors ( ) const
private

Method to print the auxiliary detectors on screen.

Definition at line 589 of file GeometryTestAlg.cxx.

589  {
590 
591  mf::LogVerbatim log("GeometryTest");
592 
593  unsigned int const nAuxDets = geom->NAuxDets();
594  log << "There are " << nAuxDets << " auxiliary detectors:";
595  for (unsigned int iDet = 0; iDet < nAuxDets; ++iDet) {
596  log << "\n[#" << iDet << "] ";
597  printAuxDetGeo(log, geom->AuxDet(iDet), " ", "");
598  } // for
599 
600  } // GeometryTestAlg::printAuxiliaryDetectors()
AuxDetGeo const & AuxDet(unsigned int const ad=0) const
Returns the specified auxiliary detector.
void printAuxDetGeo(Stream &&out, geo::AuxDetGeo const &auxDet, std::string indent, std::string firstIndent) const
Prints information of an auxiliary detector into the specified stream.
geo::GeometryCore const * geom
pointer to geometry service provider
unsigned int NAuxDets() const
Returns the number of auxiliary detectors.
void geo::GeometryTestAlg::printChannelSummary ( )
private

Definition at line 416 of file GeometryTestAlg.cxx.

417  {
418  static unsigned int OneSeg = 0;
419  static unsigned int TwoSegs = 0;
420  static unsigned int ThreeSegs = 0;
421  static unsigned int FourSegs = 0;
422  uint32_t channels = geom->Nchannels();
423  if(geom->NTPC() > 1) channels /= geom->NTPC()/2;
424 
425  for(uint32_t c = 0; c < channels; c++){
426 
427  unsigned int ChanSize = geom->ChannelToWire(c).size();
428 
429  if (ChanSize==1) ++OneSeg;
430  else if(ChanSize==2) ++TwoSegs;
431  else if(ChanSize==3) ++ThreeSegs;
432  else if(ChanSize==4) ++FourSegs;
433 
434  }
435 
436  mf::LogVerbatim("GeometryTest") << "OneSeg: " << OneSeg
437  << ", TwoSegs: " << TwoSegs
438  << ", ThreeSegs: " << ThreeSegs
439  << ", FourSegs: " << FourSegs;
440 
441  }
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
std::vector< geo::WireID > ChannelToWire(raw::ChannelID_t const channel) const
Returns a list of wires connected to the specified TPC channel.
unsigned int Nchannels() const
Returns the number of TPC readout channels in the detector.
unsigned int NTPC(unsigned int cstat=0) const
Returns the total number of TPCs in the specified cryostat.
geo::GeometryCore const * geom
pointer to geometry service provider
void geo::GeometryTestAlg::printDetDim ( )
private

Definition at line 478 of file GeometryTestAlg.cxx.

479  {
480  for(unsigned int c = 0; c < geom->Ncryostats(); ++c){
481 
482  mf::LogVerbatim("GeometryTest") << "Cryo " << c;
483  mf::LogVerbatim("GeometryTest") << " width: "
484  << geom->CryostatHalfWidth(c);
485  mf::LogVerbatim("GeometryTest") << " height: "
487  mf::LogVerbatim("GeometryTest") << " length: "
488  << geom->CryostatLength(c);
489 
490  mf::LogVerbatim("GeometryTest") << " TPC 0";
491  mf::LogVerbatim("GeometryTest") << " width: "
492  << geom->DetHalfWidth(0,c);
493  mf::LogVerbatim("GeometryTest") << " height: "
494  << geom->DetHalfHeight(0,c);
495  mf::LogVerbatim("GeometryTest") << " length: "
496  << geom->DetLength(0,c);
497  }
498  }
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
geo::Length_t CryostatHalfHeight(geo::CryostatID const &cid) const
Returns the height of the cryostat (y direction)
geo::Length_t DetHalfWidth(geo::TPCID const &tpcid) const
Returns the half width of the active volume of the specified TPC.
unsigned int Ncryostats() const
Returns the number of cryostats in the detector.
geo::Length_t DetHalfHeight(geo::TPCID const &tpcid) const
Returns the half height of the active volume of the specified TPC.
geo::Length_t DetLength(geo::TPCID const &tpcid) const
Returns the length of the active volume of the specified TPC.
geo::Length_t CryostatHalfWidth(geo::CryostatID const &cid) const
Returns the half width of the cryostat (x direction)
geo::Length_t CryostatLength(geo::CryostatID const &cid) const
Returns the length of the cryostat (z direction)
geo::GeometryCore const * geom
pointer to geometry service provider
void geo::GeometryTestAlg::printDetectorIntro ( ) const
private

Definition at line 386 of file GeometryTestAlg.cxx.

386  {
387 
388  geo::WireGeo const& testWire = geom->Wire(geo::WireID(0, 0, 1, 10));
389  mf::LogVerbatim log("GeometryTest");
390  log
391  << "Wire Rmax " << testWire.RMax()
392  << "\nWire length " << 2.*testWire.HalfL()
393  << "\nWire Rmin " << testWire.RMin()
394  ;
395 
396  if (fComputeMass) {
397  log
398  << "\nTotal mass " << geom->TotalMass();
399  }
400 
401  log
402  << "\nNumber of views " << geom->Nviews()
403  << "\nNumber of channels " << geom->Nchannels()
404  << "\nMaximum number of:"
405  << "\n TPC in a cryostat: " << geom->MaxTPCs()
406  << "\n planes in a TPC: " << geom->MaxPlanes()
407  << "\n wires in a plane: " << geom->MaxWires()
408  << "\nTotal number of TPCs " << geom->TotalNTPC()
409  << "\nAuxiliary detectors " << geom->NAuxDets()
410  ;
411 
412  } // GeometryTestAlg::printDetectorIntro()
Geometry description of a TPC wireThe wire is a single straight segment on a wire plane...
Definition: WireGeo.h:65
unsigned int TotalNTPC() const
Returns the total number of TPCs in the detector.
unsigned int MaxWires() const
Returns the largest number of wires among all planes in this detector.
unsigned int Nchannels() const
Returns the number of TPC readout channels in the detector.
double TotalMass() const
Returns the total mass [kg] of the specified volume (default: world).
unsigned int MaxPlanes() const
Returns the largest number of planes among all TPCs in this detector.
double RMax() const
Returns the outer half-size of the wire [cm].
Definition: WireGeo.cxx:91
double HalfL() const
Returns half the length of the wire [cm].
Definition: WireGeo.h:128
WireGeo const & Wire(geo::WireID const &wireid) const
Returns the specified wire.
unsigned int Nviews() const
Returns the number of views (different wire orientations)
unsigned int MaxTPCs() const
Returns the largest number of TPCs a cryostat in the detector has.
double RMin() const
Returns the inner radius of the wire (usually 0) [cm].
Definition: WireGeo.cxx:95
geo::GeometryCore const * geom
pointer to geometry service provider
unsigned int NAuxDets() const
Returns the number of auxiliary detectors.
bool fComputeMass
Whether to print the detector mass.
void geo::GeometryTestAlg::printVolBounds ( )
private

Definition at line 444 of file GeometryTestAlg.cxx.

445  {
446  double origin[3] = {0.};
447  double world[3] = {0.};
448  for(unsigned int c = 0; c < geom->Ncryostats(); ++c){
449  geom->Cryostat(c).LocalToWorld(origin, world);
450 
451  mf::LogVerbatim("GeometryTest") << "Cryo " << c;
452  mf::LogVerbatim("GeometryTest") << " -x: " << world[0] - geom->Cryostat(c).HalfWidth();
453  mf::LogVerbatim("GeometryTest") << " +x: " << world[0] + geom->Cryostat(c).HalfWidth();
454  mf::LogVerbatim("GeometryTest") << " -y: " << world[1] - geom->Cryostat(c).HalfHeight();
455  mf::LogVerbatim("GeometryTest") << " +y: " << world[1] + geom->Cryostat(c).HalfHeight();
456  mf::LogVerbatim("GeometryTest") << " -z: " << world[2] - geom->Cryostat(c).Length()/2;
457  mf::LogVerbatim("GeometryTest") << " +z: " << world[2] + geom->Cryostat(c).Length()/2;
458 
459  for(unsigned int t = 0; t < geom->NTPC(c); ++t){
460  geom->Cryostat(c).TPC(t).LocalToWorld(origin, world);
461 
462  mf::LogVerbatim("GeometryTest") << " TPC " << t;
463  mf::LogVerbatim("GeometryTest") << " -x: " << world[0] - geom->Cryostat(c).TPC(t).HalfWidth();
464  mf::LogVerbatim("GeometryTest") << " +x: " << world[0] + geom->Cryostat(c).TPC(t).HalfWidth();
465  mf::LogVerbatim("GeometryTest") << " -y: " << world[1] - geom->Cryostat(c).TPC(t).HalfHeight();
466  mf::LogVerbatim("GeometryTest") << " +y: " << world[1] + geom->Cryostat(c).TPC(t).HalfHeight();
467  mf::LogVerbatim("GeometryTest") << " -z: " << world[2] - geom->Cryostat(c).TPC(t).Length()/2;
468  mf::LogVerbatim("GeometryTest") << " +z: " << world[2] + geom->Cryostat(c).TPC(t).Length()/2;
469  }
470  }
471 
472  }
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
unsigned int Ncryostats() const
Returns the number of cryostats in the detector.
double Length() const
Length is associated with z coordinate [cm].
Definition: TPCGeo.h:115
double HalfWidth() const
Half width of the cryostat [cm].
CryostatGeo const & Cryostat(geo::CryostatID const &cryoid) const
Returns the specified cryostat.
unsigned int NTPC(unsigned int cstat=0) const
Returns the total number of TPCs in the specified cryostat.
double HalfHeight() const
Half height of the cryostat [cm].
double HalfHeight() const
Height is associated with y coordinate [cm].
Definition: TPCGeo.h:111
void LocalToWorld(const double *cryo, double *world) const
Transform point from local cryostat frame to world frame.
Definition: CryostatGeo.h:387
const TPCGeo & TPC(unsigned int itpc) const
Return the itpc&#39;th TPC in the cryostat.
Definition: CryostatGeo.cxx:93
void LocalToWorld(const double *tpc, double *world) const
Transform point from local TPC frame to world frame.
Definition: TPCGeo.h:563
double HalfWidth() const
Width is associated with x coordinate [cm].
Definition: TPCGeo.h:107
geo::GeometryCore const * geom
pointer to geometry service provider
constexpr Point origin()
Returns a origin position with a point of the specified type.
Definition: geo_vectors.h:227
double Length() const
Length of the cryostat [cm].
Definition: CryostatGeo.h:107
void geo::GeometryTestAlg::printWirePos ( )
private

Definition at line 502 of file GeometryTestAlg.cxx.

503  {
504  unsigned int cs = 0;
505 
506  for(unsigned int t=0; t<std::floor(geom->NTPC()/12)+1; ++t){
507  for(unsigned int p=0; p<3; ++p){
508  for(unsigned int w=0; w<geom->Cryostat(0).TPC(t).Plane(p).Nwires(); w++){
509 
510  double xyz[3] = {0.};
511  geom->Cryostat(0).TPC(t).Plane(p).Wire(w).GetCenter(xyz);
512 
513  std::cout << "WireID (" << cs << ", " << t << ", " << p << ", " << w
514  << "): x = " << xyz[0]
515  << ", y = " << xyz[1]
516  << ", z = " << xyz[2] << std::endl;
517  }
518  }
519  }
520  }
WireGeo const & Wire(unsigned int iwire) const
Definition: PlaneGeo.cxx:506
p
Definition: test.py:223
CryostatGeo const & Cryostat(geo::CryostatID const &cryoid) const
Returns the specified cryostat.
unsigned int NTPC(unsigned int cstat=0) const
Returns the total number of TPCs in the specified cryostat.
const TPCGeo & TPC(unsigned int itpc) const
Return the itpc&#39;th TPC in the cryostat.
Definition: CryostatGeo.cxx:93
unsigned int Nwires() const
Number of wires in this plane.
Definition: PlaneGeo.h:269
PlaneGeo const & Plane(geo::View_t view) const
Return the plane in the tpc with View_t view.
Definition: TPCGeo.cxx:263
void GetCenter(double *xyz, double localz=0.0) const
Fills the world coordinate of a point on the wire.
Definition: WireGeo.cxx:73
const char * cs
geo::GeometryCore const * geom
pointer to geometry service provider
QTextStream & endl(QTextStream &s)
void geo::GeometryTestAlg::printWiresInTPC ( const TPCGeo tpc,
std::string  indent = "" 
) const
private

Definition at line 525 of file GeometryTestAlg.cxx.

526  {
527  const unsigned int nPlanes = tpc.Nplanes();
528 
529  tpc.PrintTPCInfo(
530  mf::LogVerbatim("GeometryTest") << indent,
532  );
533 
534  for(unsigned int p = 0; p < nPlanes; ++p) {
535  const geo::PlaneGeo& plane = tpc.Plane(p);
536  const unsigned int nWires = plane.Nwires();
537 
538  plane.PrintPlaneInfo(
539  mf::LogVerbatim("GeometryTest") << indent << " ", indent + " ",
540  8 /* large verbosity */
541  );
542 
543  for(unsigned int w = 0; w < nWires; ++w) {
544  const geo::WireGeo& wire = plane.Wire(w);
545  // this additional check is preserved to test alternative transformation
546  // code paths; center is expected to match wire.GetCenter()
547  // BUG the double brace syntax is required to work around clang bug 21629
548  // (https://bugs.llvm.org/show_bug.cgi?id=21629)
549  std::array<double, 3U> const local = {{ 0.0, 0.0, 0.0 }};
550  std::array<double, 3U> center;
551  wire.LocalToWorld(local.data(), center.data());
552 
553  // the wire should be aligned on z axis, half on each side of 0,
554  // in its local frame
555  mf::LogVerbatim("GeometryTest") << indent
556  << " wire #" << w
557  << " at " << lar::dump::array<3>(center)
558  << "\n" << indent << " start at " << lar::dump::vector3D(wire.GetStart())
559  << "\n" << indent << " middle at " << lar::dump::vector3D(wire.GetCenter())
560  << "\n" << indent << " end at " << lar::dump::vector3D(wire.GetEnd())
561  ;
562  } // for wire
563  } // for plane
564  } // GeometryTestAlg::printWiresInTPC()
void GetStart(double *xyz) const
Definition: WireGeo.h:157
Geometry description of a TPC wireThe wire is a single straight segment on a wire plane...
Definition: WireGeo.h:65
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
auto vector3D(Vector3D const &v)
Returns a manipulator which will print the specified vector.
Definition: DumpUtils.h:301
WireGeo const & Wire(unsigned int iwire) const
Definition: PlaneGeo.cxx:506
void LocalToWorld(const double *wire, double *world) const
Transform point from local wire frame to world frame.
Definition: WireGeo.h:355
Geometry information for a single wire plane.The plane is represented in the geometry by a solid whic...
Definition: PlaneGeo.h:82
p
Definition: test.py:223
void PrintPlaneInfo(Stream &&out, std::string indent="", unsigned int verbosity=1) const
Prints information about this plane.
Definition: PlaneGeo.h:1539
void GetEnd(double *xyz) const
Definition: WireGeo.h:163
unsigned int Nwires() const
Number of wires in this plane.
Definition: PlaneGeo.h:269
def center(depos, point)
Definition: depos.py:117
void GetCenter(double *xyz, double localz=0.0) const
Fills the world coordinate of a point on the wire.
Definition: WireGeo.cxx:73
static constexpr unsigned int MaxVerbosity
Maximum verbosity supported by PrintTPCInfo().
Definition: TPCGeo.h:677
unsigned int geo::GeometryTestAlg::Run ( )
virtual

Runs the test, returns a number of errors (very unlikely!)

Definition at line 151 of file GeometryTestAlg.cxx.

152  {
153 
154  if (!geom) {
155  throw cet::exception("GeometryTestAlg")
156  << "GeometryTestAlg not configured: no valid geometry provided.\n";
157  }
158 
159  unsigned int nErrors = 0; // currently unused
160 
161  // change the printed version number when changing the "GeometryTest" output
162  //
163  // Version 1.1:
164  // more TPC information when printing all geometry
165  //
166  mf::LogVerbatim("GeometryTest") << "GeometryTest version 1.1";
167 
168  mf::LogInfo("GeometryTestInfo")
169  << "Running on detector: '" << geom->DetectorName() << "'";
170 
171  mf::LogVerbatim("GeometryTest")
172  << " Running on detector: '" << geom->DetectorName() << "'"
173  << "\nGeometry file: " << geom->ROOTFile();
174 
175  try{
176  if (shouldRunTests("DetectorIntro")) {
177  MF_LOG_INFO("GeometryTest") << "detector introduction:";
179  MF_LOG_INFO("GeometryTest") << "complete.";
180  }
181 
182  if (shouldRunTests("CheckOverlaps")) {
183  MF_LOG_INFO("GeometryTest") << "test for overlaps ...";
184  gGeoManager->CheckOverlaps(1e-5);
185  gGeoManager->PrintOverlaps();
186  if (!gGeoManager->GetListOfOverlaps()->IsEmpty()) {
187  mf::LogError("GeometryTest")
188  << gGeoManager->GetListOfOverlaps()->GetSize()
189  << " overlaps found in geometry during overlap test!";
190  ++nErrors;
191  }
192  MF_LOG_INFO("GeometryTest") << "complete.";
193  }
194 
195  if (shouldRunTests("ThoroughCheck")) {
196  MF_LOG_INFO("GeometryTest") << "thorough geometry test ...";
197  gGeoManager->CheckGeometryFull();
198  if (!gGeoManager->GetListOfOverlaps()->IsEmpty()) {
199  mf::LogError("GeometryTest")
200  << gGeoManager->GetListOfOverlaps()->GetSize()
201  << " overlaps found in geometry during thorough test!";
202  ++nErrors;
203  }
204  MF_LOG_INFO("GeometryTest") << "complete.";
205  }
206 
207  if (shouldRunTests("Cryostat")) {
208  MF_LOG_INFO("GeometryTest") << "test Cryostat methods ...";
209  testCryostat();
210  MF_LOG_INFO("GeometryTest") << "complete.";
211  }
212 
213  if (shouldRunTests("FindVolumes")) {
214  MF_LOG_INFO("GeometryTest") << "test FindAllVolumes method ...";
215  testFindVolumes();
216  MF_LOG_INFO("GeometryTest") << "complete.";
217  }
218 
219  if (shouldRunTests("PlaneDirections")) {
220  MF_LOG_INFO("GeometryTest") << "test plane directions...";
222  MF_LOG_INFO("GeometryTest") << "complete.";
223  }
224 
225  if (shouldRunTests("WireOrientations")) {
226  MF_LOG_INFO("GeometryTest") << "test wire orientations...";
228  MF_LOG_INFO("GeometryTest") << "complete.";
229  }
230 
231  if (shouldRunTests("ChannelToROP")) {
232  MF_LOG_INFO("GeometryTest") << "test channel to ROP and back ...";
234  MF_LOG_INFO("GeometryTest") << "complete.";
235  }
236 
237  if (shouldRunTests("ChannelToWire")) {
238  MF_LOG_INFO("GeometryTest") << "test channel to plane wire and back ...";
240  MF_LOG_INFO("GeometryTest") << "complete.";
241  }
242 
243  if (shouldRunTests("FindPlaneCenters")) {
244  MF_LOG_INFO("GeometryTest") << "test find plane centers...";
246  MF_LOG_INFO("GeometryTest") << "complete.";
247  }
248 
249  if (shouldRunTests("WireCoordFromPlane")) {
250  MF_LOG_INFO("GeometryTest") << "test PlaneGeo::WireCoordinate...";
252  MF_LOG_INFO("GeometryTest") << "complete.";
253  }
254 
255  if (shouldRunTests("ParallelWires")) {
256  MF_LOG_INFO("GeometryTest") << "test wire parallelism...";
258  MF_LOG_INFO("GeometryTest") << "complete.";
259  }
260 
261  if (shouldRunTests("PlanePointDecomposition")) {
262  MF_LOG_INFO("GeometryTest") << "test plane point decomposition...";
264  MF_LOG_INFO("GeometryTest") << "complete.";
265  }
266 
267  if (shouldRunTests("PlaneProjections")) {
268  MF_LOG_INFO("GeometryTest") << "test PlaneGeo::PointProjection...";
270  MF_LOG_INFO("GeometryTest") << "complete.";
271  }
272 
273  if (shouldRunTests("WireCoordAngle")) {
274  MF_LOG_INFO("GeometryTest") << "testWireCoordAngle...";
276  MF_LOG_INFO("GeometryTest") << "complete.";
277  }
278 
279  if (shouldRunTests("Projection")) {
280  MF_LOG_INFO("GeometryTest") << "testProject...";
281  testProject();
282  MF_LOG_INFO("GeometryTest") << "complete.";
283  }
284 
285  if (shouldRunTests("WirePos")) {
286  MF_LOG_INFO("GeometryTest") << "testWirePos...";
287  // There is a contradiction here, and these must be tested differently
288  // Testing based on detector ID should NOT become common practice
289  MF_LOG_INFO("GeometryTest") << "disabled.";
290  }
291 
292  if (shouldRunTests("NearestWire")) {
293  MF_LOG_INFO("GeometryTest") << "testNearestWire...";
294  testNearestWire();
295  MF_LOG_INFO("GeometryTest") << "complete.";
296  }
297 
298  if (shouldRunTests("WireIntersection")) {
299  MF_LOG_INFO("GeometryTest") << "testWireIntersection...";
301  MF_LOG_INFO("GeometryTest") << "testWireIntersection complete";
302  }
303 
304  if (shouldRunTests("ThirdPlane")) {
305  MF_LOG_INFO("GeometryTest") << "testThirdPlane...";
306  testThirdPlane();
307  MF_LOG_INFO("GeometryTest") << "complete.";
308  }
309 
310  if (shouldRunTests("ThirdPlaneSlope")) {
311  MF_LOG_INFO("GeometryTest") << "testThirdPlaneSlope...";
313  MF_LOG_INFO("GeometryTest") << "complete.";
314  }
315 
316  if (shouldRunTests("WirePitch")) {
317  MF_LOG_INFO("GeometryTest") << "testWirePitch...";
318  testWirePitch();
319  MF_LOG_INFO("GeometryTest") << "complete.";
320  }
321 
322  if (shouldRunTests("InterWireProjectedDistance")) {
323  MF_LOG_INFO("GeometryTest") << "testInterWireProjectedDistance...";
325  MF_LOG_INFO("GeometryTest") << "complete.";
326  }
327 
328  if (shouldRunTests("PlanePitch")) {
329  MF_LOG_INFO("GeometryTest") << "testPlanePitch...";
330  testPlanePitch();
331  MF_LOG_INFO("GeometryTest") << "complete.";
332  }
333 
334  if (shouldRunTests("Stepping")) {
335  MF_LOG_INFO("GeometryTest") << "testStepping...";
336  testStepping();
337  MF_LOG_INFO("GeometryTest") << "complete.";
338  }
339 
340  if (shouldRunTests("FindAuxDet")) {
341  MF_LOG_INFO("GeometryTest") << "testFindAuxDet...";
342  testFindAuxDet();
343  MF_LOG_INFO("GeometryTest") << "complete.";
344  }
345 
346  if (shouldRunTests("PrintWires")) {
347  MF_LOG_INFO("GeometryTest") << "printAllGeometry...";
349  MF_LOG_INFO("GeometryTest") << "complete.";
350  }
351  }
352  catch (cet::exception &e) {
353  mf::LogWarning("GeometryTest") << "exception caught: \n" << e;
354  if (fNonFatalExceptions.count(e.category()) == 0) throw;
355  }
356 
357  std::ostringstream out;
358  if (!fRunTests.CheckQueryRegistry(out)) {
359  throw cet::exception("GeometryTest")
360  << "(postumous) configuration error detected!\n"
361  << out.str();
362  }
363 
364  mf::LogInfo log("GeometryTest");
365  log << "Tests completed:";
366  auto tests_run = fRunTests.AcceptedNames();
367  if (tests_run.empty()) {
368  log << "\n no test run";
369  }
370  else {
371  log << "\n " << tests_run.size() << " tests run:\t ";
372  for (std::string const& test_name: tests_run) log << " " << test_name;
373  }
374  auto tests_skipped = fRunTests.RejectedNames();
375  if (!tests_skipped.empty()) {
376  log << "\n " << tests_skipped.size() << " tests skipped:\t ";
377  for (std::string const& test_name: tests_skipped) log << " " << test_name;
378  }
379 
380  return nErrors;
381  } // GeometryTestAlg::Run()
void testPlaneProjection() const
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
void testPlaneDirections() const
std::set< std::string > fNonFatalExceptions
std::string string
Definition: nybbler.cc:12
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
void testWireCoordAngle() const
Names_t RejectedNames() const
Returns a list of names that were rejected.
Definition: NameSelector.h:153
void printDetectorIntro() const
void testWireOrientations() const
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
bool shouldRunTests(std::string test_name) const
void testWireIntersection() const
std::string ROOTFile() const
Returns the full directory path to the geometry file source.
Names_t AcceptedNames() const
Returns a list of names that were accepted.
Definition: NameSelector.h:150
void testChannelToROP() const
const double e
testing::NameSelector fRunTests
test filter
void testThirdPlane_dTdW() const
std::string DetectorName() const
Returns a string with the name of the detector, as configured.
void testInterWireProjectedDistance() const
bool CheckQueryRegistry() const
Checks that no known element with valid response was left unqueried.
Definition: NameSelector.h:160
void testChannelToWire() const
#define MF_LOG_INFO(category)
void testWireCoordFromPlane() const
void printAllGeometry() const
void testParallelWires() const
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
void testFindAuxDet() const
void testThirdPlane() const
void testPlanePointDecomposition() const
geo::GeometryCore const * geom
pointer to geometry service provider
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
virtual void geo::GeometryTestAlg::Setup ( geo::GeometryCore const &  new_geo)
inlinevirtual

Runs the test.

Definition at line 117 of file GeometryTestAlg.h.

117 { geom = &new_geo; }
geo::GeometryCore const * geom
pointer to geometry service provider
bool geo::GeometryTestAlg::shouldRunTests ( std::string  test_name) const
inlineprivate

Definition at line 3597 of file GeometryTestAlg.cxx.

3597  {
3598  return fRunTests(test_name);
3599  } // GeometryTestAlg::shouldRunTests()
testing::NameSelector fRunTests
test filter
void geo::GeometryTestAlg::testAPAWirePos ( )
private

Definition at line 2044 of file GeometryTestAlg.cxx.

2045  {
2046  double origin[3] = {0.};
2047  double tpcworld[3] = {0.};
2048  double xyz[3] = {0.};
2049  double xyzprev[3] = {0.};
2050  for(size_t cs = 0; cs < geom->Ncryostats(); ++cs){
2051  for(size_t t = 0; t < geom->Cryostat(cs).NTPC(); ++t){
2052  const geo::TPCGeo* tpc = &geom->Cryostat(cs).TPC(t);
2053  tpc->LocalToWorld(origin, tpcworld);
2054 
2055  for (size_t i=0; i < tpc->Nplanes(); ++i) {
2056  const geo::PlaneGeo* plane = &tpc->Plane(i);
2057 
2058  for (size_t j = 1; j < plane->Nwires(); ++j) {
2059  geo::WireGeo const& wire = plane->Wire(j);
2060  geo::WireGeo const& wireprev = plane->Wire(j-1);
2061 
2062  wire.GetCenter(xyz);
2063  wireprev.GetCenter(xyzprev);
2064 
2065  // top TPC wires increase in -y
2066  if(tpcworld[1] > 0 && xyz[1] > xyzprev[1])
2067  throw cet::exception("WireOrderProblem") << "\n\ttop TPC wires do not increase in -y order in"
2068  << "Cryostat " << cs
2069  << ", TPC " << t
2070  << ", Plane " << i
2071  << "; at wire " << j << "\n";
2072  // bottom TPC wires increase in +y
2073  else if(tpcworld[1] < 0 && xyz[1] < xyzprev[1])
2074  throw cet::exception("WireOrderProblem") << "\n\tbottom TPC wires do not increase in +y order in"
2075  << "Cryostat " << cs
2076  << ", TPC " << t
2077  << ", Plane " << i
2078  << "; at wire " << j << "\n";
2079  }// end loop over wires
2080  }// end loop over planes
2081  }// end loop over tpcs
2082  }// end loop over cryostats
2083 
2084  }
Geometry description of a TPC wireThe wire is a single straight segment on a wire plane...
Definition: WireGeo.h:65
WireGeo const & Wire(unsigned int iwire) const
Definition: PlaneGeo.cxx:506
unsigned int Nplanes() const
Number of planes in this tpc.
Definition: TPCGeo.h:165
Geometry information for a single TPC.
Definition: TPCGeo.h:38
unsigned int Ncryostats() const
Returns the number of cryostats in the detector.
Geometry information for a single wire plane.The plane is represented in the geometry by a solid whic...
Definition: PlaneGeo.h:82
unsigned int NTPC() const
Number of TPCs in this cryostat.
Definition: CryostatGeo.h:181
CryostatGeo const & Cryostat(geo::CryostatID const &cryoid) const
Returns the specified cryostat.
const TPCGeo & TPC(unsigned int itpc) const
Return the itpc&#39;th TPC in the cryostat.
Definition: CryostatGeo.cxx:93
unsigned int Nwires() const
Number of wires in this plane.
Definition: PlaneGeo.h:269
PlaneGeo const & Plane(geo::View_t view) const
Return the plane in the tpc with View_t view.
Definition: TPCGeo.cxx:263
void GetCenter(double *xyz, double localz=0.0) const
Fills the world coordinate of a point on the wire.
Definition: WireGeo.cxx:73
const char * cs
void LocalToWorld(const double *tpc, double *world) const
Transform point from local TPC frame to world frame.
Definition: TPCGeo.h:563
geo::GeometryCore const * geom
pointer to geometry service provider
constexpr Point origin()
Returns a origin position with a point of the specified type.
Definition: geo_vectors.h:227
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
void geo::GeometryTestAlg::testChannelToROP ( ) const
private

Definition at line 1456 of file GeometryTestAlg.cxx.

1456  {
1457 
1458  // test that an invalid channel yields an invalid ROP
1459  try {
1461  if (invalidROP.isValid) {
1462  throw cet::exception("testChannelToROP")
1463  << "ROP from an invalid channel ("
1464  << raw::InvalidChannelID << ") is " << std::string(invalidROP)
1465  << " !?\n";
1466  } // if invalid rop is not invalid
1467  }
1468  catch (cet::exception const& e) {
1469  mf::LogWarning("testChannelToROP")
1470  << "Non-compilant ChannelToROP() throws on invalid channel.";
1471  }
1472 
1473  // for each channel, test that its ROP contains it;
1474  // we assume each ROP contains contiguous channel IDs
1475  for (raw::ChannelID_t channel = 0; channel < geom->Nchannels(); ++channel) {
1476 
1477  readout::ROPID const ropid = geom->ChannelToROP(channel);
1478 
1479  auto const firstChannel = geom->FirstChannelInROP(ropid);
1480  auto const lastChannel = firstChannel + geom->Nchannels(ropid);
1481 
1482  if ((channel < firstChannel) || (channel >= lastChannel)) {
1483  throw cet::exception("testChannelToROP")
1484  << "Channel " << channel << " comes from ROP " << std::string(ropid)
1485  << ", which contains only channels from " << firstChannel
1486  << " to " << lastChannel << " (excluded)\n";
1487  } // if
1488 
1489  } // for channel
1490 
1491 
1492  } // GeometryTestAlg::testChannelToROP()
std::string string
Definition: nybbler.cc:12
bool isValid
Whether this ID points to a valid element.
Definition: geo_types.h:211
uint8_t channel
Definition: CRTFragment.hh:201
unsigned int Nchannels() const
Returns the number of TPC readout channels in the detector.
constexpr ChannelID_t InvalidChannelID
ID of an invalid channel.
Definition: RawTypes.h:32
const double e
raw::ChannelID_t FirstChannelInROP(readout::ROPID const &ropid) const
Returns the ID of the first channel in the specified readout plane.
Class identifying a set of planes sharing readout channels.
readout::ROPID ChannelToROP(raw::ChannelID_t channel) const
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
unsigned int ChannelID_t
Type representing the ID of a readout channel.
Definition: RawTypes.h:28
geo::GeometryCore const * geom
pointer to geometry service provider
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
void geo::GeometryTestAlg::testChannelToWire ( ) const
private

Definition at line 1495 of file GeometryTestAlg.cxx.

1496  {
1497  using std::begin;
1498  using std::end;
1499 
1500  geo::PlaneID lastPlane; // starts invalid
1501  geo::View_t planeView = geo::kUnknown;
1502  geo::SigType_t planeSigType = geo::kMysteryType;
1503 
1504  for(geo::WireID testWireID: geom->IterateWireIDs()){
1505 
1507 
1508  if (!raw::isValidChannelID(channel)) {
1509  throw cet::exception("BadChannelLookup")
1510  << "Invalid channel returned for wire " << std::string(testWireID)
1511  << "\n";
1512  }
1513 
1514  auto const wireIDs = geom->ChannelToWire(channel);
1515 
1516  if (wireIDs.empty()) {
1517  throw cet::exception("BadChannelLookup")
1518  << "List of wires for channel #" << channel << " is empty;"
1519  << " it should have contained at least " << std::string(testWireID)
1520  << "\n";
1521  }
1522 
1523  if (std::count(begin(wireIDs), end(wireIDs), testWireID) != 1) {
1524  mf::LogError log("GeometryTest");
1525  log << wireIDs.size() << " wire IDs associated with channel #"
1526  << channel << ":";
1527  for (auto const& wid: wireIDs)
1528  log << "\n " << wid;
1529  throw cet::exception("BadChannelLookup")
1530  << "Wire " << std::string(testWireID) << " is on channel #" << channel
1531  << " but ChannelToWire() does not map the channel to that wire\n";
1532  }
1533 
1534  // currently (LArSoft 6.12) signal type from channel and from plane use
1535  // the same underlying code, so the following test is not very valuable
1536  auto const channelSigType = geom->SignalType(channel);
1537  if (channelSigType != geom->SignalType(testWireID.planeID())) {
1538  throw cet::exception("BadChannelLookup")
1539  << "Geometry service claims channel #" << channel
1540  << " to be of type " << channelSigType
1541  << " but that the plane of " << std::string(testWireID)
1542  << " is of type " << geom->SignalType(testWireID.planeID()) << "\n";
1543  }
1544 
1545  auto const channelView = geom->View(channel);
1546  if (channelView != geom->Plane(testWireID).View()) {
1547  throw cet::exception("BadChannelLookup")
1548  << "Geometry service claims channel #" << channel
1549  << " should be on view " << geom->View(channel)
1550  << " but the plane of " << std::string(testWireID)
1551  << " claims to be on view " << geom->Plane(testWireID).View() << "\n";
1552  }
1553 
1554  // check that all the channels on the same plane are consistent
1555  // (sort of: it's more like "all contiguous wires in the same plane")
1556  if (testWireID.planeID() == lastPlane) {
1557  if (planeView != channelView) {
1558  throw cet::exception("BadChannelLookup")
1559  << "Geometry service claims channel #" << channel
1560  << " is on view " << channelView
1561  << " but its plane " << std::string(lastPlane)
1562  << " claims to be on view " << planeView << "\n";
1563  }
1564  if (planeSigType != channelSigType) {
1565  throw cet::exception("BadChannelLookup")
1566  << "Geometry service claims channel #" << channel
1567  << " is if type " << channelSigType
1568  << " but its plane " << std::string(lastPlane)
1569  << " to be of type " << planeSigType << "\n";
1570  }
1571  }
1572  else {
1573  lastPlane = testWireID.planeID();
1574  planeView = channelView;
1575  planeSigType = channelSigType;
1576  }
1577 
1578  } // for wires in detector
1579 
1580  } // GeometryTestAlg::testChannelToWire()
end
while True: pbar.update(maxval-len(onlies[E][S])) #print iS, "/", len(onlies[E][S]) found = False for...
PlaneGeo const & Plane(unsigned int const p, unsigned int const tpc=0, unsigned int const cstat=0) const
Returns the specified wire.
Who knows?
Definition: geo_types.h:147
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
std::string string
Definition: nybbler.cc:12
Unknown view.
Definition: geo_types.h:136
The data type to uniquely identify a Plane.
Definition: geo_types.h:472
std::vector< geo::WireID > ChannelToWire(raw::ChannelID_t const channel) const
Returns a list of wires connected to the specified TPC channel.
uint8_t channel
Definition: CRTFragment.hh:201
SigType_t SignalType(geo::PlaneID const &pid) const
Returns the type of signal on the channels of specified TPC plane.
View_t View() const
Which coordinate does this plane measure.
Definition: PlaneGeo.h:184
IteratorBox< wire_id_iterator,&GeometryCore::begin_wire_id,&GeometryCore::end_wire_id > IterateWireIDs() const
Enables ranged-for loops on all wire IDs of the detector.
enum geo::_plane_sigtype SigType_t
constexpr bool isValidChannelID(raw::ChannelID_t channel)
Definition: RawTypes.h:37
View_t View(geo::PlaneID const &pid) const
Returns the view (wire orientation) on the channels of specified TPC plane.
raw::ChannelID_t PlaneWireToChannel(WireID const &wireid) const
Returns the ID of the TPC channel connected to the specified wire.
decltype(auto) constexpr begin(T &&obj)
ADL-aware version of std::begin.
Definition: StdUtils.h:72
unsigned int ChannelID_t
Type representing the ID of a readout channel.
Definition: RawTypes.h:28
geo::GeometryCore const * geom
pointer to geometry service provider
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
void geo::GeometryTestAlg::testCryostat ( )
private

Definition at line 674 of file GeometryTestAlg.cxx.

675  {
676  mf::LogVerbatim("GeometryTest") << "There are " << geom->Ncryostats() << " cryostats in the detector";
677 
678  for(geo::CryostatGeo const& cryo: geom->IterateCryostats()) {
679 
680  {
681  mf::LogVerbatim log("GeometryTest");
682 
683  log
684  << "\n\tCryostat " << cryo.ID()
685  << " " << cryo.Volume()->GetName()
686  << " Dimensions [cm]: " << cryo.Width()
687  << " x " << cryo.Height()
688  << " x " << cryo.Length()
689  ;
690  if (fComputeMass) {
691  log
692  << "\n\t\tmass [kg]: " << cryo.Mass();
693  }
694  log
695  << "\n\t\tCryostat boundaries:"
696  << " -x:" << cryo.MinX() << " +x:" << cryo.MaxX()
697  << " -y:" << cryo.MinY() << " +y:" << cryo.MaxY()
698  << " -z:" << cryo.MinZ() << " +z:" << cryo.MaxZ()
699  ;
700  }
701 
702  // pick a position in the middle of the cryostat in the world coordinates
703  double const worldLoc[3]
704  = { cryo.CenterX(), cryo.CenterY(), cryo.CenterZ() };
705 
706  MF_LOG_DEBUG("GeometryTest") << "\t testing GeometryCore::PoitionToCryostat....";
707  geo::CryostatID cid;
708  try{
709  geom->PositionToCryostat(worldLoc, cid);
710  }
711  catch(cet::exception &e){
712  mf::LogWarning("FailedToLocateCryostat") << "\n exception caught:" << e;
713  if (fNonFatalExceptions.count(e.category()) == 0) throw;
714  }
715  if (cid != cryo.ID()) {
716  throw cet::exception("CryostatTest")
717  << "Geometry points the middle of cryostat " << std::string(cryo.ID())
718  << " to a different one (" << std::string(cid) << ")\n";
719  }
720  MF_LOG_DEBUG("GeometryTest") << "done";
721 
722  MF_LOG_DEBUG("GeometryTest") << "\t Now test the TPCs associated with this cryostat";
723  testTPC(cryo.ID());
724  }
725 
726  return;
727  }
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
std::set< std::string > fNonFatalExceptions
CryostatGeo const & PositionToCryostat(geo::Point_t const &point) const
Returns the cryostat at specified location.
std::string string
Definition: nybbler.cc:12
Geometry information for a single cryostat.
Definition: CryostatGeo.h:43
unsigned int Ncryostats() const
Returns the number of cryostats in the detector.
const double e
void testTPC(geo::CryostatID const &cid)
IteratorBox< cryostat_iterator,&GeometryCore::begin_cryostat,&GeometryCore::end_cryostat > IterateCryostats() const
Enables ranged-for loops on all cryostats of the detector.
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
#define MF_LOG_DEBUG(id)
geo::GeometryCore const * geom
pointer to geometry service provider
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
The data type to uniquely identify a cryostat.
Definition: geo_types.h:190
bool fComputeMass
Whether to print the detector mass.
void geo::GeometryTestAlg::testFindAuxDet ( ) const
private

Definition at line 3543 of file GeometryTestAlg.cxx.

3543  {
3544 
3545  /*
3546  *
3547  * Picks the center of each sensitive detector and verifies that the
3548  * correct sensitive detector and auxiliary detector are found.
3549  *
3550  */
3551 
3552  unsigned int nErrors = 0;
3553 
3554  unsigned int const nAuxDets = geom->NAuxDets();
3555 
3556  for (unsigned int iDet = 0; iDet < nAuxDets; ++iDet) {
3557 
3558  geo::AuxDetGeo const& auxDet = geom->AuxDet(iDet);
3559  unsigned int const nSensitive = auxDet.NSensitiveVolume();
3560 
3561  if (nSensitive == 0) {
3562  std::array<double, 3U> center;
3563  auxDet.GetCenter(center.data());
3564 
3565  if (!CheckAuxDetAtPosition(center.data(), iDet)) ++nErrors;
3566 
3567  }
3568  else { // if one or more sensitive detectors
3569 
3570  for (unsigned int iDetSens = 0; iDetSens < nSensitive; ++iDetSens) {
3571 
3572  geo::AuxDetSensitiveGeo const& auxDetSens
3573  = auxDet.SensitiveVolume(iDetSens);
3574  std::array<double, 3U> center;
3575  auxDetSens.GetCenter(center.data());
3576 
3577  if (!CheckAuxDetAtPosition(center.data(), iDet)) ++nErrors;
3578  if (!CheckAuxDetSensitiveAtPosition(center.data(), iDet, iDetSens))
3579  ++nErrors;
3580 
3581  } // for sensitive detectors
3582 
3583  } // if ... else
3584 
3585  } // for all auxiliary detectors
3586 
3587  if (nErrors != 0) {
3588  throw cet::exception("FindAuxDet")
3589  << "Collected " << nErrors << " errors during testFindAuxDet() test!\n";
3590  }
3591 
3592 
3593  } // GeometryTestAlg::testFindAuxDet()
AuxDetSensitiveGeo const & SensitiveVolume(size_t sv) const
Definition: AuxDetGeo.h:171
size_t NSensitiveVolume() const
Definition: AuxDetGeo.h:172
bool CheckAuxDetAtPosition(double const pos[3], unsigned int expected) const
Returns whether the auxiliary detector at pos is the expected one.
AuxDetGeo const & AuxDet(unsigned int const ad=0) const
Returns the specified auxiliary detector.
def center(depos, point)
Definition: depos.py:117
void GetCenter(double *xyz, double localz=0.0) const
Return the center position of an AuxDet.
geo::GeometryCore const * geom
pointer to geometry service provider
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
unsigned int NAuxDets() const
Returns the number of auxiliary detectors.
void GetCenter(double *xyz, double localz=0.0) const
Return the center position of an AuxDet.
Definition: AuxDetGeo.cxx:62
bool CheckAuxDetSensitiveAtPosition(double const pos[3], unsigned int expectedDet, unsigned int expectedSens) const
Returns whether the auxiliary sensitive detector at pos is expected.
unsigned int geo::GeometryTestAlg::testFindCryostatVolumes ( )
private

Definition at line 761 of file GeometryTestAlg.cxx.

761  {
762 
763  unsigned int nErrors = 0;
764 
765  std::set<std::string> volume_names;
766  volume_names.insert(geom->GetWorldVolumeName());
767  volume_names.insert("volCryostat");
768 
769  std::vector<TGeoNode const*> nodes = geom->FindAllVolumes(volume_names);
770 
771  mf::LogVerbatim log("GeometryTest");
772  log << "Found " << nodes.size() << " world and cryostat volumes:";
773  for (TGeoNode const* node: nodes) {
774  TGeoVolume const* pVolume = node->GetVolume();
775  log << "\n - '" << pVolume->GetName() << "' (a "
776  << pVolume->GetShape()->GetName() << ")";
777  } // for
778  if (nodes.size() != (1 + geom->Ncryostats())) {
779  ++nErrors;
780  mf::LogError("GeometryTest")
781  << "Found " << nodes.size() << " world and cryostat volumes! "
782  "[expecting: 1 world and " << geom->Ncryostats() << " cryostats]";
783  } // if nodes
784 
785  return nErrors;
786  } // GeometryTestAlg::testFindCryostatVolumes()
const std::string GetWorldVolumeName() const
Return the name of the world volume (needed by Geant4 simulation)
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
unsigned int Ncryostats() const
Returns the number of cryostats in the detector.
#define nodes
std::vector< TGeoNode const * > FindAllVolumes(std::set< std::string > const &vol_names) const
Returns all the nodes with volumes with any of the specified names.
geo::GeometryCore const * geom
pointer to geometry service provider
void geo::GeometryTestAlg::testFindPlaneCenters ( )
private

Definition at line 1583 of file GeometryTestAlg.cxx.

1584  {
1585  double xyz[3] = {0.}, xyzW[3] = {0.};
1586  for(size_t i = 0; i < geom->Nplanes(); ++i){
1587  geom->Plane(i).LocalToWorld(xyz,xyzW);
1588  mf::LogVerbatim("GeometryTest") << "\n\tplane "
1589  << i << " is centered at (x,y,z) = ("
1590  << xyzW[0] << "," << xyzW[1]
1591  << "," << xyzW[2] << ")";
1592  }
1593  }
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
PlaneGeo const & Plane(unsigned int const p, unsigned int const tpc=0, unsigned int const cstat=0) const
Returns the specified wire.
unsigned int Nplanes(unsigned int tpc=0, unsigned int cstat=0) const
Returns the total number of wire planes in the specified TPC.
void LocalToWorld(const double *plane, double *world) const
Transform point from local plane frame to world frame.
Definition: PlaneGeo.h:1343
geo::GeometryCore const * geom
pointer to geometry service provider
unsigned int geo::GeometryTestAlg::testFindTPCvolumePaths ( )
private

Definition at line 789 of file GeometryTestAlg.cxx.

789  {
790 
791  unsigned int nErrors = 0;
792 
793  // search the full path of all TPCs
794  std::set<std::string> volume_names;
795  for (geo::TPCGeo const& TPC: geom->IterateTPCs())
796  volume_names.insert(TPC.TotalVolume()->GetName());
797 
798  // get the right answer: how many TPCs?
799  const unsigned int NTPCs = geom->TotalNTPC();
800 
801  std::vector<std::vector<TGeoNode const*>> node_paths
802  = geom->FindAllVolumePaths(volume_names);
803 
804  mf::LogVerbatim log("GeometryTest");
805  log << "Found " << node_paths.size() << " TPC volumes:";
806  for (auto const& path: node_paths) {
807  TGeoNode const* node = path.back();
808  TGeoVolume const* pVolume = node->GetVolume();
809  log << "\n - '" << pVolume->GetName() << "' (a "
810  << pVolume->GetShape()->GetName() << ") with " << (path.size()-1)
811  << " ancestors";
812  for (TGeoNode const* pNode: path) {
813  TGeoVolume const* pVolume = pNode->GetVolume();
814  log << "\n * '" << pVolume->GetName() << "' (a "
815  << pVolume->GetShape()->GetName() << ") with a "
816  << pNode->GetMatrix()->IsA()->GetName() << " that "
817  << (pNode->GetMatrix()->IsTranslation()? "is": "is not")
818  << " a simple translation";
819  } // for node
820  } // for path
821  if (node_paths.size() != NTPCs) {
822  ++nErrors;
823  mf::LogError("GeometryTest")
824  << "Found " << node_paths.size() << " TPC volumes! "
825  "[expecting: " << NTPCs << " TPCs]";
826  } // if missed some
827 
828  return nErrors;
829  } // GeometryTestAlg::testFindTPCvolumePaths()
unsigned int TotalNTPC() const
Returns the total number of TPCs in the detector.
Geometry information for a single TPC.
Definition: TPCGeo.h:38
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
IteratorBox< TPC_iterator,&GeometryCore::begin_TPC,&GeometryCore::end_TPC > IterateTPCs() const
Enables ranged-for loops on all TPCs of the detector.
std::vector< std::vector< TGeoNode const * > > FindAllVolumePaths(std::set< std::string > const &vol_names) const
Returns paths of all nodes with volumes with the specified names.
geo::GeometryCore const * geom
pointer to geometry service provider
void geo::GeometryTestAlg::testFindVolumes ( )
private

Definition at line 832 of file GeometryTestAlg.cxx.

832  {
833  /*
834  * Finds and checks a selected number of volumes by name:
835  * - world volume
836  * - cryostat volumes
837  * - TPCs (returns full paths)
838  */
839 
840  unsigned int nErrors = 0;
841 
842  nErrors += testFindWorldVolumes();
843  nErrors += testFindCryostatVolumes();
844  nErrors += testFindTPCvolumePaths();
845 
846  if (nErrors != 0) {
847  throw cet::exception("FindVolumes")
848  << "Collected " << nErrors << " errors during FindAllVolumes() test!\n";
849  }
850 
851  } // GeometryTestAlg::testFindVolumes()
unsigned int testFindCryostatVolumes()
unsigned int testFindTPCvolumePaths()
unsigned int testFindWorldVolumes()
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
unsigned int geo::GeometryTestAlg::testFindWorldVolumes ( )
private

Definition at line 730 of file GeometryTestAlg.cxx.

730  {
731 
732  unsigned int nErrors = 0;
733 
734  std::set<std::string> volume_names;
735  std::vector<TGeoNode const*> nodes;
736 
737  // world
738  volume_names.insert(geom->GetWorldVolumeName());
739  nodes = geom->FindAllVolumes(volume_names);
740  {
741  mf::LogVerbatim log("GeometryTest");
742  log << "Found " << nodes.size() << " world volumes '"
743  << geom->GetWorldVolumeName() << "':";
744  for (TGeoNode const* node: nodes) {
745  TGeoVolume const* pVolume = node->GetVolume();
746  log << "\n - '" << pVolume->GetName() << "' (a "
747  << pVolume->GetShape()->GetName() << ")";
748  } // for
749  } // anonymous block
750  if (nodes.size() != 1) {
751  ++nErrors;
752  mf::LogError("GeometryTest")
753  << "Found " << nodes.size() << " world volumes '"
754  << geom->GetWorldVolumeName() << "! [expecting: one!!]";
755  } // if nodes
756 
757  return nErrors;
758  } // GeometryTestAlg::testFindWorldVolumes()
const std::string GetWorldVolumeName() const
Return the name of the world volume (needed by Geant4 simulation)
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
#define nodes
std::vector< TGeoNode const * > FindAllVolumes(std::set< std::string > const &vol_names) const
Returns all the nodes with volumes with any of the specified names.
geo::GeometryCore const * geom
pointer to geometry service provider
void geo::GeometryTestAlg::testInterWireProjectedDistance ( ) const
private

Definition at line 3183 of file GeometryTestAlg.cxx.

3183  {
3184 
3185  /*
3186  * For each wire plane:
3187  * * we pick some projected directions; for each one:
3188  * * we manipulate it in some way that does not affect the result
3189  * * check that the projected distance is as expected
3190  * * add some arbitrary component on the drift direction; for each one:
3191  * * check that the projected distance is as expected
3192  * * check that the 3D distance is as expected
3193  *
3194  * We do not test directions parallel to the wires because they get
3195  * numerically unstable and the expectation may potentially differ a lot
3196  * being calculated with a different procedure.
3197  */
3198 
3199  constexpr lar::util::RealComparisons cmp { 1e-4 };
3200 
3201  static double const V3 = std::sqrt(3.0);
3202 
3203  // BUG the deduction guide for std::array seems not to be implemented yet
3204  // in Clang 5.0.0
3205  // BUG the double brace syntax is required to work around clang bug 21629
3206  // (https://bugs.llvm.org/show_bug.cgi?id=21629)
3207 // std::array const testProjections = {
3208  std::array<geo::PlaneGeo::WireCoordProjection_t, 5U> const testProjections {{
3214 // };
3215  }};
3216  // BUG the deduction guide for std::array seems not to be implemented yet
3217  // in Clang 5.0.0
3218  // BUG the double brace syntax is required to work around clang bug 21629
3219  // (https://bugs.llvm.org/show_bug.cgi?id=21629)
3220 // constexpr std::array testDriftOffsets { -20.0, -10.0, 0.0, +10.0, +20.0 };
3221  constexpr std::array<double, 5U> testDriftOffsets
3222  {{ -20.0, -10.0, 0.0, +10.0, +20.0 }};
3223 
3224  unsigned int nErrors = 0; // error count for the final report
3225 
3226  for (geo::PlaneGeo const& plane: geom->IteratePlanes()) {
3227 
3228  double const pitch = plane.WirePitch();
3229  auto const normalDir = plane.GetNormalDirection<geo::Vector_t>();
3230 
3231  for (auto const& testProjBase: testProjections) {
3232 
3233  //
3234  // expected result is kind of encoded in the chosen projections
3235  //
3236  double const expected = testProjBase.R() * pitch;
3237  double const expectedSqr = cet::square(expected);
3238 
3239  // we flip the projection around: result should not change
3240  for (double dirL: { -1.0, 1.0 }) for (double dirW: { -1.0, 1.0 }) for (double scale: { 0.5, 1.0, 3.0 }) {
3241 
3242  //
3243  // test the projection directly
3244  //
3246  { scale * dirL * testProjBase.X(), scale * dirW * testProjBase.Y() };
3247 
3248  double const interWireFromProj
3249  = plane.InterWireProjectedDistance(testProj);
3250  if (cmp.nonEqual(interWireFromProj, expected)) {
3251  mf::LogProblem("") << "ERROR: on plane " << plane.ID()
3252  << " distance between wires on projected direction " << testProj
3253  << " is " << interWireFromProj << " cm (expected: "
3254  << expected << ")"
3255  ;
3256  ++nErrors;
3257  } // if unexpected result
3258 
3259  // this is how much we needed to expand the test direction vector
3260  // (happens to work for the special case expected = 0 too)
3261  double const dScale = expected / testProj.R();
3262 
3263  //
3264  // test a 3D direction
3265  //
3266 
3267  auto const baseDir
3268  = plane.ComposeVector<geo::Vector_t>(0.0, testProj);
3269 
3270  for (double const driftOffset: testDriftOffsets) {
3271  auto const testDir = baseDir + driftOffset * normalDir;
3272 
3273  double const interWire = plane.InterWireProjectedDistance(testDir);
3274  if (cmp.nonEqual(interWire, expected)) {
3275  mf::LogProblem("") << "ERROR: on plane " << plane.ID()
3276  << " projected distance between wires on direction " << testDir
3277  << " (from projection " << testProj << " and offset "
3278  << driftOffset << ") is " << interWire << " cm (expected: "
3279  << expected << ")"
3280  ;
3281  ++nErrors;
3282  } // if unexpected result
3283 
3284  // build the expectation adding the drift component to the projected
3285  double const expected3D
3286  = std::sqrt(expectedSqr + cet::square(driftOffset * dScale));
3287 
3288  double const interWire3D = plane.InterWireDistance(testDir);
3289  if (cmp.nonEqual(interWire3D, expected3D)) {
3290  mf::LogProblem("") << "ERROR: on plane " << plane.ID()
3291  << " distance between wires on direction " << testDir
3292  << " (from projection " << testProj << " and offset "
3293  << driftOffset << ") is " << interWire3D << " cm (expected: "
3294  << expected3D << ")"
3295  ;
3296  ++nErrors;
3297  } // if unexpected result
3298 
3299 
3300 
3301  } // for drifts
3302 
3303  } // for dirL, dirW and scale
3304 
3305  } // for all test projection directions
3306 
3307  } // for all planes
3308 
3309 
3310  if (nErrors > 0U) {
3311  throw cet::exception("InterWireProjectedDistance")
3312  << "unexpected distances in " << nErrors << " tests!";
3313  } // end loop over planes
3314 
3315  } // GeometryTestAlg::testInterWireProjectedDistance()
IteratorBox< plane_iterator,&GeometryCore::begin_plane,&GeometryCore::end_plane > IteratePlanes() const
Enables ranged-for loops on all planes of the detector.
const char expected[]
Definition: Exception_t.cc:22
Provides simple real number checks.
ROOT::Math::DisplacementVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Vector_t
Type for representation of momenta in 3D space.
Definition: geo_vectors.h:164
constexpr T square(T x)
Definition: pow.h:21
MaybeLogger_< ELseverityLevel::ELsev_error, true > LogProblem
const double e
ROOT::Math::DisplacementVector2D< ROOT::Math::Cartesian2D< double >, WireCoordinateReferenceTag > WireCoordProjection_t
Type for projections in the wire base representation.
Definition: PlaneGeo.h:130
Geometry information for a single wire plane.The plane is represented in the geometry by a solid whic...
Definition: PlaneGeo.h:82
geo::GeometryCore const * geom
pointer to geometry service provider
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
void geo::GeometryTestAlg::testNearestWire ( )
private

Definition at line 2100 of file GeometryTestAlg.cxx.

2101  {
2102  // Even if you comment it out, please leave the TStopWatch code
2103  // in this code for additional testing. The NearestChannel routine
2104  // is the most frequently called in the simulation, so its execution time
2105  // is an important component of LArSoft's speed.
2106  TStopwatch stopWatch;
2107  stopWatch.Start();
2108 
2109  bool bTestWireCoordinate = true;
2110 
2111  // get a wire and find its center
2112  for (geo::PlaneGeo const& plane: geom->IteratePlanes()) {
2113 
2114  geo::PlaneID const& planeID = plane.ID();
2115  const unsigned int NWires = plane.Nwires();
2116 
2117  decltype(auto) IncreasingWireDir = plane.GetIncreasingWireDirection();
2118 
2119  MF_LOG_DEBUG("GeoTestWireCoordinate")
2120  << "The direction of increasing wires for plane " << planeID
2121  << " (theta=" << plane.Wire(0).ThetaZ() << " pitch="
2122  << plane.WirePitch() << " orientation="
2123  << (plane.Orientation() == geo::kHorizontal? "H": "V")
2124  << (plane.WireIDincreasesWithZ()? "+": "-")
2125  << ") is " << IncreasingWireDir;
2126 
2127  geo::WireID::WireID_t w = 0;
2128  for (geo::WireGeo const& wire: plane.IterateWires()) {
2129 
2130  geo::WireID wireID(planeID, w++);
2131 
2132  decltype(auto) wire_center = wire.GetCenter();
2133 
2134  uint32_t nearest = 0;
2135  std::vector< geo::WireID > wireIDs;
2136 
2137  try{
2138  // The double[] version tested here falls back on the
2139  // TVector3 version, so this test both.
2140  nearest = geom->NearestChannel(wire_center, planeID);
2141 
2142  // We also want to test the std::vector<double> version
2143  std::vector<double> posWorldV(3);
2144  for (int i=0; i<3; ++i) {
2145  posWorldV[i] = wire_center[i] + 0.001;
2146  }
2147  nearest = geom->NearestChannel(posWorldV, planeID);
2148  }
2149  catch(cet::exception &e){
2150  mf::LogWarning("GeoTestCaughtException") << e;
2151  if (fNonFatalExceptions.count(e.category()) == 0) throw;
2152  }
2153 
2154  try{
2155  wireIDs = geom->ChannelToWire(nearest);
2156 
2157  if ( wireIDs.empty() ) {
2158  throw cet::exception("BadPositionToChannel")
2159  << "test point is at " << wire_center
2160  << " nearest channel is " << nearest
2161  << " for " << std::string(wireID)
2162  << "\n";
2163  }
2164  }
2165  catch(cet::exception &e){
2166  mf::LogWarning("GeoTestCaughtException") << e;
2167  if (fNonFatalExceptions.count(e.category()) == 0) throw;
2168  }
2169 
2170  if(std::find(wireIDs.begin(), wireIDs.end(), wireID) == wireIDs.end()) {
2171  throw cet::exception("BadPositionToChannel")
2172  << "Current wire " << std::string(wireID)
2173  << " has a world position at " << wire_center
2174  << "\nNearestWire for this position is "
2175  << geom->NearestWire(wire_center, planeID)
2176  << "\nNearestChannel is " << nearest
2177  << "\nShould be channel " << geom->PlaneWireToChannel(wireID);
2178  } // if good lookup fails
2179 
2180 
2181  // nearest wire, integral and floating point
2182  try {
2183  // The test consists in sampling NStep (=5) points between the current
2184  // wire and the previous/next, following the normal to the wire.
2185  // We expect WireCoordinate() to reflect the same shift.
2186 
2187  // using absolute value just in case (what happens if w1 > w2?)
2188 
2189  const double pitch = std::abs(geom->WirePitch(planeID));
2190 
2191  TVector3 wire_shifted;
2192  TVector3 const step = pitch * IncreasingWireDir;
2193 
2194  constexpr int NSteps = 5; // odd value avoids testing half-way
2195  for (int i = -NSteps; i <= +NSteps; ++i) {
2196  // we move away by this fraction of wire:
2197  const double f = NSteps? (double(i) / NSteps): 0.0;
2198 
2199  // these are the actual shifts on the positive directions y and z
2200  TVector3 const delta = f * step;
2201  TVector3 const wire_shifted = wire_center + delta;
2202 
2203  // we expect this wire number
2204  const double expected_wire = wireID.Wire + f;
2205 
2206  double wire_from_wc = 0;
2207  if (bTestWireCoordinate) {
2208  try {
2209  wire_from_wc = geom->WireCoordinate(wire_shifted, planeID);
2210  }
2211  catch (cet::exception& e) {
2212  if (hasCategory(e, "NotImplemented")) {
2213  MF_LOG_ERROR("GeoTestWireCoordinate")
2214  << "WireCoordinate() is not implemented for your ChannelMap;"
2215  " skipping the test";
2216  bTestWireCoordinate = false;
2217  }
2218  else if (bTestWireCoordinate) throw;
2219  }
2220  }
2221  if (bTestWireCoordinate) {
2222  if (std::abs(wire_from_wc - expected_wire) > 1e-3) {
2223  // throw cet::exception("GeoTestErrorWireCoordinate")
2224  mf::LogError("GeoTestWireCoordinate")
2225  << "wire " << wireID
2226  << " [center: " << wire_center << "] on step of "
2227  << i << "/" << NSteps
2228  << " x" << step << " = " << delta
2229  << " cm shows " << wire_from_wc << ", " << expected_wire
2230  << " expected.\n";
2231  } // if mismatch
2232 
2233  } // if testing WireCoordinate
2234 
2235  if ((expected_wire > -0.5) && (expected_wire < NWires - 0.5)) {
2236  const unsigned int expected_wire_number = std::round(expected_wire);
2237  unsigned int wire_number_from_wc;
2238  try {
2239  wire_number_from_wc = geom->NearestWire(wire_shifted, planeID);
2240  }
2241  catch (cet::exception& e) {
2242  throw cet::exception("GeoTestErrorWireCoordinate", "", e)
2243  // MF_LOG_ERROR("GeoTestWireCoordinate")
2244  << "wire " << std::string(wireID)
2245  << " [center: " << wire_center << "] on step of "
2246  << i << "/" << NSteps
2247  << " x" << step << " = " << delta
2248  << " cm failed NearestWire(), " << expected_wire_number
2249  << " expected (more precisely, " << expected_wire << ").\n";
2250  }
2251 
2252  if (mf::isDebugEnabled()) {
2253  // In debug mode, we print a lot and we don't (fatally) complain
2254  std::stringstream e;
2255  e << "wire " << wireID
2256  << " [center: " << wire_center << "] on step of "
2257  << i << "/" << NSteps
2258  << " x" << step << " = " << delta << " cm near to "
2259  << wire_number_from_wc;
2260  if (wire_number_from_wc != expected_wire_number) {
2261  e << ", " << expected_wire_number
2262  << " expected (more precisely, " << expected_wire << ").";
2263  // throw e;
2264  MF_LOG_ERROR("GeoTestWireCoordinate") << e.str();
2265  }
2266  else {
2267  mf::LogVerbatim("GeoTestWireCoordinate") << e.str();
2268  }
2269  }
2270  else if (wire_number_from_wc != expected_wire_number) {
2271  // In production mode, we don't print anything and throw on error
2272  throw cet::exception("GeoTestErrorWireCoordinate")
2273  << "wire " << std::string(wireID)
2274  << " [center: " << wire_center << "] on step of "
2275  << i << "/" << NSteps
2276  << " x" << step << " = " << delta
2277  << " cm near to " << wire_number_from_wc
2278  << ", " << expected_wire_number
2279  << " expected (more precisely, " << expected_wire << ").";
2280  }
2281  } // if shifted wire not outside boundaries
2282 
2283  } // for i
2284 
2285  } // try
2286  catch(cet::exception &e) {
2287  mf::LogWarning("GeoTestCaughtException") << e;
2288  if (fNonFatalExceptions.count(e.category()) == 0) throw;
2289  }
2290 
2291  } // for all wires in the plane
2292  } // end loop over planes
2293 
2294  stopWatch.Stop();
2295  MF_LOG_DEBUG("GeoTestWireCoordinate") << "\tdone testing closest channel";
2296  stopWatch.Print();
2297 
2298  // trigger an exception with NearestChannel
2299  mf::LogVerbatim("GeoTestWireCoordinate") << "\tattempt to cause an exception to be caught "
2300  << "when looking for a nearest channel";
2301 
2302  // pick a position out of the world
2303  double posWorld[3];
2304  geom->WorldBox(nullptr, posWorld + 0,
2305  nullptr, posWorld + 1, nullptr, posWorld + 2);
2306  for (int i = 0; i < 3; ++i) posWorld[i] *= 2.;
2307 
2308  bool hasThrown = false;
2309  unsigned int nearest_to_what = 0;
2310  try{
2311  nearest_to_what = geom->NearestChannel(posWorld, 0, 0, 0);
2312  }
2313  catch(const geo::InvalidWireIDError& e){
2314  mf::LogWarning("GeoTestCaughtException") << e
2315  << "\nReturned wire would be: " << e.wire_number
2316  << ", suggested: " << e.better_wire_number;
2317  hasThrown = true;
2318  }
2319  catch(cet::exception& e){
2320  mf::LogWarning("GeoTestCaughtException") << e;
2321  hasThrown = true;
2322  }
2323  if (!hasThrown) {
2325  // ok, then why do we disable it?
2326  // an implementation might prefer to cap the wire number and go on
2327  // instead of throwing.
2328  MF_LOG_WARNING("GeoTestWireCoordinate")
2329  << "GeometryCore::NearestChannel() did not raise an exception"
2330  " on out-of-world position (" << posWorld[0] << "; "
2331  << posWorld[1] << "; " << posWorld[2] << "), and returned "
2332  << nearest_to_what << " instead.\n"
2333  "This is normally considered a failure.";
2334  }
2335  else {
2336  throw cet::exception("GeoTestErrorNearestChannel")
2337  << "GeometryCore::NearestChannel() did not raise an exception"
2338  " on out-of-world position (" << posWorld[0] << "; "
2339  << posWorld[1] << "; " << posWorld[2] << "), and returned "
2340  << nearest_to_what << " instead\n";
2341  }
2342  }
2343 
2344  }
geo::Length_t WireCoordinate(double YPos, double ZPos, geo::PlaneID const &planeid) const
Returns the index of the nearest wire to the specified position.
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
std::set< std::string > fNonFatalExceptions
IteratorBox< plane_iterator,&GeometryCore::begin_plane,&GeometryCore::end_plane > IteratePlanes() const
Enables ranged-for loops on all planes of the detector.
std::string string
Definition: nybbler.cc:12
void WorldBox(double *xlo, double *xhi, double *ylo, double *yhi, double *zlo, double *zhi) const
Fills the arguments with the boundaries of the world.
int wire_number
the invalid wire number
Definition: Exceptions.h:167
The data type to uniquely identify a Plane.
Definition: geo_types.h:472
#define MF_LOG_ERROR(category)
std::vector< geo::WireID > ChannelToWire(raw::ChannelID_t const channel) const
Returns a list of wires connected to the specified TPC channel.
STL namespace.
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
geo::Length_t WirePitch(geo::PlaneID const &planeid) const
Returns the distance between two consecutive wires.
T abs(T value)
Planes that are in the horizontal plane.
Definition: geo_types.h:140
int better_wire_number
a suggestion for a good wire number
Definition: Exceptions.h:168
const double e
IDparameter< geo::WireID > WireID
Member type of validated geo::WireID parameter.
geo::WireID::WireID_t NearestWire(geo::Point_t const &point, geo::PlaneID const &planeid) const
Returns the index of wire closest to position in the specified TPC.
Geometry information for a single wire plane.The plane is represented in the geometry by a solid whic...
Definition: PlaneGeo.h:82
bool fDisableValidWireIDcheck
disable test on out-of-world NearestWire()
bool isDebugEnabled()
raw::ChannelID_t PlaneWireToChannel(WireID const &wireid) const
Returns the ID of the TPC channel connected to the specified wire.
for(std::string line;std::getline(inFile, line);)
Definition: regex_t.cc:37
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
#define MF_LOG_DEBUG(id)
Exception thrown on invalid wire number (e.g. NearestWireID())
Definition: Exceptions.h:158
#define MF_LOG_WARNING(category)
raw::ChannelID_t NearestChannel(geo::Point_t const &worldLoc, geo::PlaneID const &planeid) const
Returns the ID of the channel nearest to the specified position.
LArSoft geometry interface.
Definition: ChannelGeo.h:16
geo::GeometryCore const * geom
pointer to geometry service provider
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
static std::array< double, 3 > GetIncreasingWireDirection(const geo::PlaneGeo &plane)
Returns the direction on plane orthogonal to wires where wire number increases.
void geo::GeometryTestAlg::testParallelWires ( ) const
private

Definition at line 1139 of file GeometryTestAlg.cxx.

1139  {
1140 
1141  //
1142  // checks that all the wires in the same plane are parallel
1143  //
1144  auto const vectorIs = lar::util::makeVector3DComparison(geom->coordIs);
1145 
1146  for (geo::PlaneGeo const& plane: geom->IteratePlanes()) {
1147 
1148  decltype(auto) genDir = plane.GetWireDirection();
1149 
1150  geo::WireID::WireID_t wireNo = 0;
1151  for (geo::WireGeo const& wire: plane.IterateWires()) {
1152 
1153  geo::WireID const wireID(plane.ID(), wireNo++);
1154 
1155  decltype(auto) wireDir = wire.Direction();
1156 
1157  if (vectorIs.nonEqual(wireDir, genDir)) {
1158  throw cet::exception("ParallelWires")
1159  << "Wire " << std::string(wireID) << " has direction " << wireDir
1160  << ", not parallel to wire direction " << genDir
1161  << " according to the plane " << std::string(plane.ID()) << "\n";
1162  }
1163 
1164  } // for wires in plane
1165 
1166  } // for planes
1167 
1168  } // GeometryTestAlg::testParallelWires()
IteratorBox< plane_iterator,&GeometryCore::begin_plane,&GeometryCore::end_plane > IteratePlanes() const
Enables ranged-for loops on all planes of the detector.
std::string string
Definition: nybbler.cc:12
static lar::util::RealComparisons< geo::Length_t > coordIs
Value of tolerance for equality comparisons.
auto makeVector3DComparison(RealType threshold)
Creates a Vector3DComparison from a RealComparisons object.
IDparameter< geo::WireID > WireID
Member type of validated geo::WireID parameter.
Geometry information for a single wire plane.The plane is represented in the geometry by a solid whic...
Definition: PlaneGeo.h:82
for(std::string line;std::getline(inFile, line);)
Definition: regex_t.cc:37
Direction
Definition: AssnsIter.h:13
if(!yymsg) yymsg
LArSoft geometry interface.
Definition: ChannelGeo.h:16
geo::GeometryCore const * geom
pointer to geometry service provider
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
void geo::GeometryTestAlg::testPlaneDirections ( ) const
private

Definition at line 941 of file GeometryTestAlg.cxx.

941  {
942  /*
943  * The test verifies that all the planes in the geometry respect the
944  * orientation requirements:
945  *
946  * { (wire direction) , (wire coordinate increase), (plane normal) }
947  *
948  * { (width direction) , (length direction), (plane normal) }
949  *
950  * both be a positively defined base.
951  */
952 
954 
955  unsigned int nErrors = 0;
956  for (geo::PlaneGeo const& plane: geom->IteratePlanes()) {
957 
958  //
959  // check the ( wire ; wire coordinate ; normal) base
960  //
961 
962  // this funny way declares a reference or not, depending on return type
963  decltype(auto) planeNormal = plane.GetNormalDirection();
964  decltype(auto) wireCoordDir = plane.GetIncreasingWireDirection();
965  decltype(auto) wireDir = plane.GetWireDirection();
966 
967  double const wireFrame = wireDir.Cross(wireCoordDir).Dot(planeNormal);
968 
969  if (coordIs.nonEqual(wireFrame, 1.0)) {
970  ++nErrors;
971 
972  mf::LogProblem("GeometryTestAlg")
973  << "Plane " << plane.ID()
974  << " has wire direction " << wireDir
975  << " wire coordinate direction " << wireCoordDir
976  << " and normal " << planeNormal
977  << " , yielding to a non-positive plane-coordinate definition"
978  << " (l x w . n = " << wireFrame << ")";
979  } // if error
980 
981  //
982  // check the ( width ; depth ; normal) base
983  //
984 
985  decltype(auto) widthDir = plane.WidthDir();
986  decltype(auto) depthDir = plane.DepthDir();
987  double const geoFrame = widthDir.Cross(depthDir).Dot(planeNormal);
988 
989  if (coordIs.nonEqual(geoFrame, 1.0)) {
990  ++nErrors;
991 
992  mf::LogProblem("GeometryTestAlg")
993  << "Plane " << plane.ID()
994  << " has width direction " << widthDir
995  << " depth direction " << depthDir
996  << " and normal " << planeNormal
997  << " , yielding to a non-positive plane-coordinate definition"
998  << " (w x d . n = " << geoFrame << ")";
999  } // if error
1000 
1001  } // for plane
1002 
1003  if (nErrors > 0) {
1004  throw cet::exception("GeometryTestAlg")
1005  << "testPlaneDirections() accumulated " << nErrors
1006  << " errors (see messages above)\n";
1007  }
1008 
1009  } // GeometryTestAlg::testPlaneDirections()
IteratorBox< plane_iterator,&GeometryCore::begin_plane,&GeometryCore::end_plane > IteratePlanes() const
Enables ranged-for loops on all planes of the detector.
Provides simple real number checks.
MaybeLogger_< ELseverityLevel::ELsev_error, true > LogProblem
const double e
Geometry information for a single wire plane.The plane is represented in the geometry by a solid whic...
Definition: PlaneGeo.h:82
if(!yymsg) yymsg
geo::GeometryCore const * geom
pointer to geometry service provider
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
static std::array< double, 3 > GetIncreasingWireDirection(const geo::PlaneGeo &plane)
Returns the direction on plane orthogonal to wires where wire number increases.
void geo::GeometryTestAlg::testPlanePitch ( )
private

Definition at line 3319 of file GeometryTestAlg.cxx.

3320  {
3321  // loop over all planes to be sure the pitch is consistent
3322 
3323  if (fExpectedPlanePitches.empty()) {
3324  // hard code the value we think it should be for each detector;
3325  // this is legacy and you should not add anything:
3326  // add the expectation to the FHiCL configuration of the test instead
3327  if(geom->DetectorName() == "bo") {
3328  fExpectedPlanePitches = { 0.65 };
3329  }
3330  if (!fExpectedPlanePitches.empty()) {
3331  mf::LogInfo("PlanePitch")
3332  << "Using legacy plane pitch parameters hard-coded for the detector '"
3333  << geom->DetectorName() << "'";
3334  }
3335  }
3336  if (fExpectedPlanePitches.empty()) {
3337  mf::LogWarning("PlanePitch")
3338  << "no expected plane pitch;"
3339  " I'll just check that they are all the same";
3340  }
3341  else {
3342  mf::LogInfo log("PlanePitch");
3343  log << "Expected plane pitch per plane pair, in centimetres:";
3344  for (double pitch: fExpectedPlanePitches) log << " " << pitch;
3345  log << " [...]";
3346  }
3347 
3348  unsigned int nPitchErrors = 0;
3349  for (geo::TPCID const& tpcid: geom->IterateTPCIDs()) {
3350 
3351  geo::TPCGeo const& TPC = geom->TPC(tpcid);
3352  const unsigned int nPlanes = TPC.Nplanes();
3353  if (nPlanes < 2) continue;
3354 
3355  double expectedPitch = 0.;
3356  if (fExpectedPlanePitches.empty()) {
3357  expectedPitch = TPC.PlanePitch(0, 1);
3358  MF_LOG_DEBUG("PlanePitch")
3359  << "Plane pitch between the first two planes of " << tpcid << ": "
3360  << expectedPitch << " cm";
3361  }
3362 
3363  geo::PlaneID::PlaneID_t p = 0; // plane number
3364  while (++p < nPlanes) {
3365  // which pitch to expect:
3366  // - if they did not tell us anything:
3367  // use the one from the first two planes (already in expectedPitch)
3368  // - if they did tell something, but not for this plane:
3369  // get the last pitch they told us
3370  // - if they told us about this plane: well, then use it!
3371  if (!fExpectedPlanePitches.empty()) {
3372  if (p - 1 < fExpectedPlanePitches.size())
3373  expectedPitch = fExpectedPlanePitches[p - 1];
3374  else
3375  expectedPitch = fExpectedPlanePitches.back();
3376  } // if we have directions about plane pitch
3377 
3378  const double thisPitch = std::abs(TPC.PlanePitch(p - 1, p));
3379  if (std::abs(thisPitch - expectedPitch) > 1e-5) {
3380  mf::LogProblem("PlanePitch") << "ERROR: pitch of planes P:" << (p - 1)
3381  << " and P: " << p << " in " << tpcid
3382  << " is " << thisPitch << " cm, not " << expectedPitch
3383  << " as expected!";
3384  ++nPitchErrors;
3385  } // if unexpected pitch
3386  } // while planes
3387 
3388  } // for TPCs
3389 
3390  if (nPitchErrors > 0) {
3391  throw cet::exception("UnexpectedPlanePitch")
3392  << "unexpected pitches between " << nPitchErrors << " planes!";
3393  } // end loop over planes
3394 
3395  } // GeometryTestAlg::testPlanePitch()
double PlanePitch(unsigned int p1=0, unsigned int p2=1) const
Definition: TPCGeo.cxx:388
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
unsigned int Nplanes() const
Number of planes in this tpc.
Definition: TPCGeo.h:165
unsigned int PlaneID_t
Type for the ID number.
Definition: geo_types.h:473
Geometry information for a single TPC.
Definition: TPCGeo.h:38
IteratorBox< TPC_id_iterator,&GeometryCore::begin_TPC_id,&GeometryCore::end_TPC_id > IterateTPCIDs() const
Enables ranged-for loops on all TPC IDs of the detector.
MaybeLogger_< ELseverityLevel::ELsev_error, true > LogProblem
std::vector< double > fExpectedPlanePitches
plane pitch on each plane
T abs(T value)
const double e
std::string DetectorName() const
Returns a string with the name of the detector, as configured.
p
Definition: test.py:223
The data type to uniquely identify a TPC.
Definition: geo_types.h:386
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
#define MF_LOG_DEBUG(id)
TPCGeo const & TPC(unsigned int const tpc=0, unsigned int const cstat=0) const
Returns the specified TPC.
geo::GeometryCore const * geom
pointer to geometry service provider
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
void geo::GeometryTestAlg::testPlanePointDecomposition ( ) const
private

Definition at line 1172 of file GeometryTestAlg.cxx.

1172  {
1173 
1174  //
1175  // For each plane:
1176  //
1177  // 1) create a plane point with arbitrary distance from the plane,
1178  // wire coordinate multiple (N) of wire pitch, and wire direction
1179  // coordinate 0 or half a wire length in either directions
1180  //
1181  // 2) compose into a 3D vector, and verify that the nearest wire is the one
1182  // expected (N)
1183  //
1184  // 3) decompose back the 3D vector, and verify that the result matches the
1185  // starting decomposition
1186  //
1187  // 4) also verify singly PointProjection() and DistanceFromPlane()
1188  //
1189  // 5) verify DriftPoint()
1190  //
1191  //
1192 
1194  auto vectorIs = lar::util::makeVector3DComparison(coordIs);
1195 
1196  unsigned int nErrors = 0;
1197  for (geo::PlaneGeo const& plane: geom->IteratePlanes()) {
1198 
1199  auto const& planeNorm = plane.GetNormalDirection();
1200  auto const& wirePitch = plane.WirePitch();
1201  auto const& refPoint = plane.ProjectionReferencePoint();
1202 
1203  unsigned int const lastWire = plane.Nwires() - 1;
1204  geo::WireID::WireID_t wireNo = 0;
1205  for (geo::WireGeo const& wire: plane.IterateWires()) {
1206 
1207  geo::WireID const wireID(plane.ID(), wireNo++);
1208 
1209  constexpr double distance = 5.0; // 5 cm from the plane
1210 
1211  auto const wireDirStep = wire.HalfL() / 2.0; // quarter of the length
1212 
1213  for (int iWDStep = -1; iWDStep <= 1; ++iWDStep) {
1214 
1215  //
1216  // prepare expectation
1217  //
1218  auto const wireDirOffset = iWDStep * wireDirStep;
1219 
1220  auto const expectedPoint = wire.GetCenter()
1221  + wireDirOffset * wire.Direction()
1222  + distance * planeNorm;
1223 
1224  auto const expectedWireCoord = wirePitch * wireID.Wire;
1225  auto const expectedWireDirCoord = wireDirOffset
1226  + wire.Direction().Dot(wire.GetCenter() - refPoint);
1227 
1228  geo::PlaneGeo::WireCoordProjection_t const expectedProj
1229  (expectedWireDirCoord, expectedWireCoord);
1230 
1231  //
1232  // composition
1233  //
1234  auto point = plane.ComposePoint(distance, expectedProj);
1235 
1236  if (vectorIs.nonEqual(point, expectedPoint)) {
1237  ++nErrors;
1238  mf::LogProblem("GeometryTestAlg")
1239  << "[testPlanePointDecomposition] ComposePoint(): "
1240  << "Point with projection " << expectedProj
1241  << " and distance from plane " << distance
1242  << " was reported as " << point
1243  << " while it is expected to be at " << expectedPoint
1244  << " (on wire " << std::string(wireID) << ")";
1245  } // if wrong point
1246 
1247  //
1248  // decomposition
1249  //
1250  auto const decomp = plane.DecomposePoint(point);
1251  if (coordIs.nonEqual(decomp.distance, distance)) {
1252  ++nErrors;
1253  mf::LogProblem("GeometryTestAlg")
1254  << "[testPlanePointDecomposition] DecomposePoint(): "
1255  << "Point " << point << " (on wire " << std::string(wireID) << ")"
1256  << " is reported to have distance from plane " << decomp.distance
1257  << " cm, while " << distance << " is expected";
1258  } // if wrong distance
1259  if (coordIs.nonEqual(decomp.projection.X(), expectedWireDirCoord)) {
1260  ++nErrors;
1261  mf::LogProblem("GeometryTestAlg")
1262  << "[testPlanePointDecomposition] DecomposePoint(): "
1263  << "Point " << point << " (on wire " << std::string(wireID) << ")"
1264  << " is reported to have wire direction coordinate "
1265  << decomp.projection.X()
1266  << " cm, while " << expectedWireDirCoord << " is expected";
1267  } // if wrong coordinate along the wire
1268  if (coordIs.nonEqual(decomp.projection.Y(), expectedWireCoord)) {
1269  ++nErrors;
1270  mf::LogProblem("GeometryTestAlg")
1271  << "[testPlanePointDecomposition] DecomposePoint(): "
1272  << "Point " << point << " (on wire " << std::string(wireID) << ")"
1273  << " is reported to have wire coordinate "
1274  << decomp.projection.Y()
1275  << " cm, while " << expectedWireCoord << " is expected";
1276  } // if wrong wire coordinate
1277 
1278  //
1279  // projection
1280  //
1281  auto const proj = plane.PointProjection(point);
1282  if (coordIs.nonEqual(proj.X(), expectedWireDirCoord)) {
1283  ++nErrors;
1284  mf::LogProblem("GeometryTestAlg")
1285  << "[testPlanePointDecomposition] PointProjection(): "
1286  << "Point " << point << " (on wire " << std::string(wireID) << ")"
1287  << " is reported to have wire direction coordinate " << proj.X()
1288  << " cm, while " << expectedWireDirCoord << " is expected";
1289  } // if wrong wire coordinate
1290  if (coordIs.nonEqual(proj.Y(), expectedWireCoord)) {
1291  ++nErrors;
1292  mf::LogProblem("GeometryTestAlg")
1293  << "[testPlanePointDecomposition] PointProjection(): "
1294  << "Point " << point << " (on wire " << std::string(wireID) << ")"
1295  << " is reported to have wire coordinate " << proj.Y()
1296  << " cm, while " << expectedWireCoord << " is expected";
1297  } // if wrong wire coordinate
1298 
1299  //
1300  // distance
1301  //
1302  auto const dist = plane.DistanceFromPlane(point);
1303  if (coordIs.nonEqual(dist, distance)) {
1304  ++nErrors;
1305  mf::LogProblem("GeometryTestAlg")
1306  << "[testPlanePointDecomposition] DistanceFromPlane(): "
1307  << "Point " << point << " (on wire " << std::string(wireID) << ")"
1308  << " is reported to have distance " << dist << " cm from plane "
1309  << plane.ID() << ", while " << distance << " is expected";
1310  } // if wrong distance
1311 
1312  //
1313  // drift
1314  //
1315  // BUG the double brace syntax is required to work around clang bug 21629
1316  // (https://bugs.llvm.org/show_bug.cgi?id=21629)
1317  std::array<double, 3> drifts {{ -distance, distance, 2.*distance }};
1318  for (double drift: drifts) {
1319 
1320  // DriftPoint() moves the point in the drift direction,
1321  // which is opposite to the plane normal:
1322  auto const expectedDistance = distance - drift;
1323 
1324  //
1325  // drift it by a known value
1326  //
1327  auto point = expectedPoint;
1328  plane.DriftPoint(point, drift);
1329 
1330  //
1331  // check the new distance
1332  //
1333  auto dist = plane.DistanceFromPlane(point);
1334  if (coordIs.nonEqual(dist, expectedDistance)) {
1335  ++nErrors;
1336  mf::LogProblem("GeometryTestAlg")
1337  << "[testPlanePointDecomposition] DriftPoint(): "
1338  << "Point " << expectedPoint
1339  << " (distant " << distance << " cm from the plane)"
1340  << " (on wire " << std::string(wireID) << ")"
1341  << " drifted by " << drift << " became " << point
1342  << " which has distance from plane " << dist
1343  << " cm, while " << expectedDistance << " was expected";
1344  } // if wrong distance
1345 
1346  } // for drifts
1347 
1348  //
1349  // containment
1350  //
1351  // skip this test for the first and last wire, since the containment
1352  // area might well exclude the whole wire if it's short enough
1353  bool const isExtremeWire
1354  = (wireID.Wire == 0) || (wireID.Wire == lastWire);
1355  if (!isExtremeWire && !plane.isProjectionOnPlane(expectedPoint)) {
1356  ++nErrors;
1357  mf::LogProblem("GeometryTestAlg")
1358  << "[testPlanePointDecomposition] isProjectionOnPlane(): "
1359  << "Point " << expectedPoint << " (center of " << wireID
1360  << ") is not believed to be on the plane, but it should.";
1361  }
1362 
1363  } // for different wire direction steps
1364 
1365  } // for wires in plane
1366 
1367  } // for planes
1368 
1369  if (nErrors > 0) {
1370  throw cet::exception("GeometryTestAlg")
1371  << "testPlanePointDecomposition() accumulated " << nErrors
1372  << " errors (see messages above)\n";
1373  }
1374 
1375  } // GeometryTestAlg::testPlanePointDecomposition()
Geometry description of a TPC wireThe wire is a single straight segment on a wire plane...
Definition: WireGeo.h:65
IteratorBox< plane_iterator,&GeometryCore::begin_plane,&GeometryCore::end_plane > IteratePlanes() const
Enables ranged-for loops on all planes of the detector.
std::string string
Definition: nybbler.cc:12
Provides simple real number checks.
MaybeLogger_< ELseverityLevel::ELsev_error, true > LogProblem
const double e
auto makeVector3DComparison(RealType threshold)
Creates a Vector3DComparison from a RealComparisons object.
ROOT::Math::DisplacementVector2D< ROOT::Math::Cartesian2D< double >, WireCoordinateReferenceTag > WireCoordProjection_t
Type for projections in the wire base representation.
Definition: PlaneGeo.h:130
Geometry information for a single wire plane.The plane is represented in the geometry by a solid whic...
Definition: PlaneGeo.h:82
double distance(double x1, double y1, double z1, double x2, double y2, double z2)
constexpr double dist(const TReal *x, const TReal *y, const unsigned int dimension)
unsigned int WireID_t
Type for the ID number.
Definition: geo_types.h:561
geo::GeometryCore const * geom
pointer to geometry service provider
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
void geo::GeometryTestAlg::testPlanePointDecompositionFrame ( ) const
private

Definition at line 1635 of file GeometryTestAlg.cxx.

1635  {
1636 
1637  //
1638  // For each plane:
1639  //
1640  // 1) create a plane point with arbitrary distance from the plane,
1641  // width and depth coordinates all across the plane
1642  //
1643  // 2) compose into a 3D vector
1644  //
1645  // 3) decompose back the 3D vector, and verify that the result matches the
1646  // starting decomposition
1647  //
1648  // 4) also verify singly PointProjection() and DistanceFromPlane()
1649  //
1650  // 5) verify DriftPoint()
1651  //
1652  //
1653 
1655  auto vectorIs = lar::util::makeVector3DComparison(coordIs);
1656 
1657  // steps on each side of the center, within the plane:
1658  constexpr int steps = 5;
1659  static_assert(steps > 0, "Steps must be a positive number.");
1660  // how many width units we go far from the center (1 means stay inside)
1661  constexpr int nOutsides = 1;
1662 
1663  unsigned int nErrors = 0;
1664  for (geo::PlaneGeo const& plane: geom->IteratePlanes()) {
1665 
1666  auto const& planeNorm = plane.GetNormalDirection();
1667  auto const& widthDir = plane.WidthDir();
1668  auto const& depthDir = plane.DepthDir();
1669  auto const& refPoint = plane.GetCenter();
1670 
1671  double const halfWidth = plane.Width() / 2;
1672  double const halfDepth = plane.Depth() / 2;
1673  double const dW_2 = halfWidth / (steps * 2); // half width step
1674  double const dD_2 = halfDepth / (steps * 2); // half depth step
1675 
1676  for (int iW = -nOutsides * steps; iW <= +nOutsides * steps; ++iW) {
1677 
1678  double const expected_w = dW_2 * (iW * 2 + 1); // width coordinate
1679  bool const inWidth = std::abs(expected_w) <= halfWidth;
1680 
1681  for (int iD = -nOutsides * steps; iD <= +nOutsides * steps; ++iD) {
1682 
1683  double const expected_d = dD_2 * (iD * 2 + 1); // depth coordinate
1684  bool const inDepth = std::abs(expected_d) <= halfDepth;
1685 
1686  constexpr double distance = 5.0; // we might test this too...
1687 
1688  //
1689  // prepare expectation
1690  //
1691  auto const expectedPoint = refPoint
1692  + expected_w * widthDir
1693  + expected_d * depthDir
1694  + distance * planeNorm;
1695 
1696  geo::PlaneGeo::WidthDepthProjection_t const expectedProj
1697  (expected_w, expected_d);
1698 
1699  //
1700  // composition
1701  //
1702  auto point = plane.ComposePoint(distance, expectedProj);
1703 
1704  if (vectorIs.nonEqual(point, expectedPoint)) {
1705  ++nErrors;
1706  mf::LogProblem("GeometryTestAlg")
1707  << "[testPlanePointDecompositionFrame] ComposePoint(): "
1708  << "Point with projection " << expectedProj
1709  << " (width: " << expected_w << ", depth: " << expected_d
1710  << ") and distance from plane " << distance
1711  << " was reported as " << point
1712  << " while it is expected to be at " << expectedPoint;
1713  } // if wrong point
1714 
1715  //
1716  // decomposition
1717  //
1718  auto const decomp = plane.DecomposePointWidthDepth(point);
1719  if (coordIs.nonEqual(decomp.distance, distance)) {
1720  ++nErrors;
1721  mf::LogProblem("GeometryTestAlg")
1722  << "[testPlanePointDecomposition] DecomposePointWidthDepth(): "
1723  << "Point " << point
1724  << " (width: " << expected_w << ", depth: " << expected_d
1725  << ") is reported to have distance from plane " << decomp.distance
1726  << " cm, while " << distance << " is expected";
1727  } // if wrong distance
1728  if (coordIs.nonEqual(decomp.projection.X(), expected_w)) {
1729  ++nErrors;
1730  mf::LogProblem("GeometryTestAlg")
1731  << "[testPlanePointDecomposition] DecomposePointWidthDepth(): "
1732  << "Point " << point
1733  << " (width: " << expected_w << ", depth: " << expected_d
1734  << ") is reported to have width direction coordinate "
1735  << decomp.projection.X()
1736  << " cm, while " << expected_w << " is expected";
1737  } // if wrong coordinate along the wire
1738  if (coordIs.nonEqual(decomp.projection.Y(), expected_d)) {
1739  ++nErrors;
1740  mf::LogProblem("GeometryTestAlg")
1741  << "[testPlanePointDecomposition] DecomposePointWidthDepth(): "
1742  << "Point " << point
1743  << " (width: " << expected_w << ", depth: " << expected_d
1744  << ") is reported to have depth direction coordinate "
1745  << decomp.projection.Y()
1746  << " cm, while " << expected_d << " is expected";
1747  } // if wrong wire coordinate
1748 
1749  //
1750  // projection
1751  //
1752  auto const proj = plane.PointWidthDepthProjection(point);
1753  if (coordIs.nonEqual(proj.X(), expected_w)) {
1754  ++nErrors;
1755  mf::LogProblem("GeometryTestAlg")
1756  << "[testPlanePointDecomposition] PointProjection(): "
1757  << "Point " << point
1758  << " (width: " << expected_w << ", depth: " << expected_d
1759  << ") is reported to have width direction coordinate " << proj.X()
1760  << " cm, while " << expected_w << " is expected";
1761  } // if wrong wire coordinate
1762  if (coordIs.nonEqual(proj.Y(), expected_d)) {
1763  ++nErrors;
1764  mf::LogProblem("GeometryTestAlg")
1765  << "[testPlanePointDecomposition] PointProjectionWidthDepth(): "
1766  << "Point " << point
1767  << " (width: " << expected_w << ", depth: " << expected_d
1768  << ") is reported to have wire coordinate " << proj.Y()
1769  << " cm, while " << expected_d << " is expected";
1770  } // if wrong wire coordinate
1771 
1772  //
1773  // distance
1774  //
1775  auto const dist = plane.DistanceFromPlane(point);
1776  if (coordIs.nonEqual(dist, distance)) {
1777  ++nErrors;
1778  mf::LogProblem("GeometryTestAlg")
1779  << "[testPlanePointDecomposition] DistanceFromPlane(): "
1780  << "Point " << point
1781  << " (width: " << expected_w << ", depth: " << expected_d
1782  << ") is reported to have distance from plane " << dist
1783  << " cm, while " << distance << " is expected";
1784  } // if wrong distance
1785 
1786  //
1787  // drift
1788  //
1789  // BUG the double brace syntax is required to work around clang bug 21629
1790  // (https://bugs.llvm.org/show_bug.cgi?id=21629)
1791  std::array<double, 3> drifts {{ -distance, distance, 2.*distance }};
1792  for (double drift: drifts) {
1793 
1794  // DriftPoint() moves the point in the drift direction,
1795  // which is opposite to the plane normal:
1796  auto const expectedDistance = distance - drift;
1797 
1798  //
1799  // drift it by a known value
1800  //
1801  auto point = expectedPoint;
1802  plane.DriftPoint(point, drift);
1803 
1804  //
1805  // check the new distance
1806  //
1807  auto dist = plane.DistanceFromPlane(point);
1808  if (coordIs.nonEqual(dist, expectedDistance)) {
1809  ++nErrors;
1810  mf::LogProblem("GeometryTestAlg")
1811  << "[testPlanePointDecomposition] DriftPoint(): "
1812  << "Point " << expectedPoint
1813  << " (distant " << distance << " cm from the plane)"
1814  << " (width: " << expected_w << ", depth: " << expected_d
1815  << ") drifted by " << drift << " became " << point
1816  << " which has distance from plane " << dist
1817  << " cm, while " << expectedDistance << " was expected";
1818  } // if wrong distance
1819 
1820  } // for drifts
1821 
1822  //
1823  // containment
1824  //
1825  const bool expected_onPlane = inWidth && inDepth;
1826  const bool onPlane = plane.isProjectionOnPlane(expectedPoint);
1827  if (expected_onPlane != onPlane) {
1828  // always
1829  ++nErrors;
1830  mf::LogProblem("GeometryTestAlg")
1831  << "[testPlanePointDecompositionFrame] isProjectionOnPlane(): "
1832  << "Point " << expectedPoint
1833  << " (width: " << expected_w << ", depth: " << expected_d
1834  << ") is" << (onPlane? "": " not")
1835  << " believed to be on plane " << plane.ID()
1836  << ", but it should" << (expected_onPlane? "": " not be") << ".";
1837  }
1838 
1839  } // for different wire direction steps
1840 
1841  } // for wires in plane
1842 
1843  } // for planes
1844 
1845  if (nErrors > 0) {
1846  throw cet::exception("GeometryTestAlg")
1847  << "testPlanePointDecomposition() accumulated " << nErrors
1848  << " errors (see messages above)\n";
1849  }
1850 
1851  } // GeometryTestAlg::testPlanePointDecompositionFrame()
IteratorBox< plane_iterator,&GeometryCore::begin_plane,&GeometryCore::end_plane > IteratePlanes() const
Enables ranged-for loops on all planes of the detector.
Provides simple real number checks.
MaybeLogger_< ELseverityLevel::ELsev_error, true > LogProblem
T abs(T value)
ROOT::Math::DisplacementVector2D< ROOT::Math::Cartesian2D< double >, WidthDepthReferenceTag > WidthDepthProjection_t
Definition: PlaneGeo.h:151
const double e
auto makeVector3DComparison(RealType threshold)
Creates a Vector3DComparison from a RealComparisons object.
Geometry information for a single wire plane.The plane is represented in the geometry by a solid whic...
Definition: PlaneGeo.h:82
double distance(double x1, double y1, double z1, double x2, double y2, double z2)
constexpr double dist(const TReal *x, const TReal *y, const unsigned int dimension)
geo::GeometryCore const * geom
pointer to geometry service provider
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
void geo::GeometryTestAlg::testPlaneProjection ( ) const
private

Definition at line 1986 of file GeometryTestAlg.cxx.

1986  {
1987 
1988  //
1989  // Check the definition of the reference
1990  //
1991 
1993 
1994  //
1995  // Check the projections and point composition in the plane frame reference
1996  //
1998 
1999  //
2000  // Check containment
2001  //
2003 
2004 
2005  } // GeometryTestAlg::testPlaneProjection()
void testPlaneProjectionReference() const
void testPlaneProjectionOnFrame() const
void testPlanePointDecompositionFrame() const
void geo::GeometryTestAlg::testPlaneProjectionOnFrame ( ) const
private

Definition at line 1855 of file GeometryTestAlg.cxx.

1855  {
1856 
1857  //
1858  // Tests:
1859  //
1860  // 1. containment (isProjectionOnPlane())
1861  //
1862  //
1863  // 2. capping by closest point
1864  //
1865 
1867  auto vectorIs = lar::util::makeVector3DComparison(coordIs);
1868  auto const& vector2Dis = vectorIs.comp2D();
1869 
1870  // steps on each side of the center, within the plane:
1871  constexpr int steps = 5;
1872  static_assert(steps > 0, "Steps must be a positive number.");
1873  // how many width units we go far from the center (1 means stay inside)
1874  constexpr int nOutsides = 2;
1875 
1876  unsigned int nErrors = 0;
1877  for (geo::PlaneGeo const& plane: geom->IteratePlanes()) {
1878 
1879  double const halfWidth = plane.Width() / 2;
1880  double const halfDepth = plane.Depth() / 2;
1881  double const dW_2 = halfWidth / (steps * 2); // half width step
1882  double const dD_2 = halfDepth / (steps * 2); // half depth step
1883 
1884  for (int iW = -nOutsides * steps; iW <= +nOutsides * steps; ++iW) {
1885 
1886  double const expected_w = dW_2 * (iW * 2 + 1); // width coordinate
1887  bool const inWidth = std::abs(expected_w) <= halfWidth;
1888 
1889  for (int iD = -nOutsides * steps; iD <= +nOutsides * steps; ++iD) {
1890 
1891  double const expected_d = dD_2 * (iD * 2 + 1); // depth coordinate
1892  bool const inDepth = std::abs(expected_d) <= halfDepth;
1893 
1894  for (double distance: { -30., 0.0, +30.0 }) {
1895 
1896  //
1897  // prepare expectation
1898  //
1899  // definition of the test point
1900  geo::PlaneGeo::WidthDepthProjection_t const expected_proj
1901  (expected_w, expected_d);
1902 
1903  auto const expected_point
1904  = plane.ComposePoint(distance, expected_proj);
1905 
1906  //
1907  // 1. Containment test
1908  //
1909  const bool expected_onPlane = inWidth && inDepth;
1910  const bool onPlane
1911  = plane.isProjectionOnPlane(expected_point);
1912  if (expected_onPlane != onPlane) {
1913  ++nErrors;
1914  mf::LogProblem("GeometryTestAlg")
1915  << "[testPlaneProjectionOnFrame] isProjectionOnPlane(): "
1916  << "Point " << expected_point
1917  << " (width: " << expected_w << ", depth: " << expected_d
1918  << ") is" << (onPlane? "": " not")
1919  << " believed to be on plane " << plane.ID()
1920  << ", but it should" << (expected_onPlane? "": " not be")
1921  << ".";
1922  }
1923 
1924  //
1925  // 2. capping by closest point
1926  //
1927  // 2.1. capping projection
1928  //
1929  geo::PlaneGeo::WidthDepthProjection_t expected_movedProjection(
1930  (inWidth?
1931  expected_w: (expected_w < 0)? -halfWidth: +halfWidth),
1932  (inDepth?
1933  expected_d: (expected_d < 0)? -halfDepth: +halfDepth)
1934  );
1935  auto movedProjection = plane.MoveProjectionToPlane(expected_proj);
1936  if (vector2Dis.nonEqual(movedProjection, expected_movedProjection))
1937  {
1938  ++nErrors;
1939  mf::LogProblem("GeometryTestAlg")
1940  << "[testPlaneProjectionOnFrame] moveProjectionOnPlane():"
1941  << "Projection of ooint " << expected_point
1942  << " (width: " << expected_w << ", depth: " << expected_d
1943  << ") (" << (onPlane? "on": "off")
1944  << "-plane " << plane.ID() << ") was moved to "
1945  << movedProjection << " while it should have moved to "
1946  << expected_movedProjection
1947  << ".";
1948  }
1949 
1950  //
1951  // 2.2. capping point
1952  //
1953  auto expected_movedPoint
1954  = plane.ComposePoint(distance, expected_movedProjection);
1955 
1956  auto movedPoint = plane.MovePointOverPlane(expected_point);
1957  if (vectorIs.nonEqual(movedPoint, expected_movedPoint)) {
1958  ++nErrors;
1959  mf::LogProblem("GeometryTestAlg")
1960  << "[testPlaneProjectionOnFrame] movePointOnPlane(): "
1961  << "Point " << expected_point
1962  << " (width: " << expected_w << ", depth: " << expected_d
1963  << ") (" << (onPlane? "on": "off")
1964  << "-plane " << plane.ID() << ") was moved to "
1965  << movedPoint << " while it should have moved to "
1966  << expected_movedPoint
1967  << ".";
1968  }
1969 
1970  } // for distance
1971 
1972  } // for depth step
1973 
1974  } // for width step
1975 
1976  } // for planes
1977 
1978  if (nErrors > 0) {
1979  throw cet::exception("GeoTestPlaneProjection")
1980  << "Accumulated " << nErrors << " errors (see messages above)\n";
1981  }
1982 
1983  } // testPlaneProjectionOnFrame()
IteratorBox< plane_iterator,&GeometryCore::begin_plane,&GeometryCore::end_plane > IteratePlanes() const
Enables ranged-for loops on all planes of the detector.
Provides simple real number checks.
MaybeLogger_< ELseverityLevel::ELsev_error, true > LogProblem
T abs(T value)
ROOT::Math::DisplacementVector2D< ROOT::Math::Cartesian2D< double >, WidthDepthReferenceTag > WidthDepthProjection_t
Definition: PlaneGeo.h:151
const double e
auto makeVector3DComparison(RealType threshold)
Creates a Vector3DComparison from a RealComparisons object.
Geometry information for a single wire plane.The plane is represented in the geometry by a solid whic...
Definition: PlaneGeo.h:82
double distance(double x1, double y1, double z1, double x2, double y2, double z2)
geo::GeometryCore const * geom
pointer to geometry service provider
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
void geo::GeometryTestAlg::testPlaneProjectionReference ( ) const
private

Definition at line 1596 of file GeometryTestAlg.cxx.

1596  {
1597 
1598  //
1599  // Check the definition of the projection reference
1600  //
1601  unsigned int nErrors = 0;
1602  for (auto const& plane: geom->IteratePlanes()) {
1603 
1604  TVector3 reference = plane.ProjectionReferencePoint();
1605 
1606  auto decomp = plane.DecomposePoint(reference);
1607 
1608  if (geom->coordIs.nonZero(decomp.distance)) {
1609  MF_LOG_ERROR("GeometryTest")
1610  << "Plane " << plane.ID() << " reference point " << reference
1611  << " has distance " << decomp.distance << " cm (should be 0)";
1612  ++nErrors;
1613  }
1614 
1615  if (geom->coordIs.nonZero(decomp.projection.X())
1616  || geom->coordIs.nonZero(decomp.projection.Y())
1617  )
1618  {
1619  MF_LOG_ERROR("GeometryTest")
1620  << "Plane " << plane.ID() << " reference point " << reference
1621  << " has projection ( " << decomp.projection.X()
1622  << " ; " << decomp.projection.Y() << " ) cm (should be (0;0) )";
1623  ++nErrors;
1624  }
1625 
1626  } // for planes
1627  if (nErrors > 0) {
1628  throw cet::exception("GeoTestPlaneProjectionReference")
1629  << "Accumulated " << nErrors << " errors (see messages above)\n";
1630  }
1631  } // GeometryTestAlg::testPlaneProjectionReference()
IteratorBox< plane_iterator,&GeometryCore::begin_plane,&GeometryCore::end_plane > IteratePlanes() const
Enables ranged-for loops on all planes of the detector.
#define MF_LOG_ERROR(category)
static lar::util::RealComparisons< geo::Length_t > coordIs
Value of tolerance for equality comparisons.
constexpr bool nonZero(Value_t value) const
Returns whether the value is farther from 0 than the threshold.
geo::GeometryCore const * geom
pointer to geometry service provider
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
void geo::GeometryTestAlg::testProject ( )
private

Definition at line 3452 of file GeometryTestAlg.cxx.

3453  {
3454  double xlo, xhi;
3455  double ylo, yhi;
3456  double zlo, zhi;
3457  geom->WorldBox(&xlo, &xhi, &ylo, &yhi, &zlo, &zhi);
3458 
3459  double xyz[3] = { 0.0, 0.0, 0.0};
3460  double dxyz1[3] = { 1.0, 0.0, 0.0};
3461  double dxyz2[3] = {-1.0, 0.0, 0.0};
3462  double dxyz3[3] = { 0.0, 1.0, 0.0};
3463  double dxyz4[3] = { 0.0,-1.0, 0.0};
3464  double dxyz5[3] = { 0.0, 0.0, 1.0};
3465  double dxyz6[3] = { 0.0, 0.0,-1.0};
3466 
3467  double xyzo[3];
3468  geo::ProjectToBoxEdge(xyz, dxyz1, xlo, xhi, ylo, yhi, zlo, zhi, xyzo);
3469  if (std::abs(xyzo[0]-xhi)>1.E-6) abort();
3470 
3471  geo::ProjectToBoxEdge(xyz, dxyz2, xlo, xhi, ylo, yhi, zlo, zhi, xyzo);
3472  if (std::abs(xyzo[0]-xlo)>1.E-6) abort();
3473 
3474  geo::ProjectToBoxEdge(xyz, dxyz3, xlo, xhi, ylo, yhi, zlo, zhi, xyzo);
3475  if (std::abs(xyzo[1]-yhi)>1.E-6) abort();
3476 
3477  geo::ProjectToBoxEdge(xyz, dxyz4, xlo, xhi, ylo, yhi, zlo, zhi, xyzo);
3478  if (std::abs(xyzo[1]-ylo)>1.E-6) abort();
3479 
3480  geo::ProjectToBoxEdge(xyz, dxyz5, xlo, xhi, ylo, yhi, zlo, zhi, xyzo);
3481  if (std::abs(xyzo[2]-zhi)>1.E-6) abort();
3482 
3483  geo::ProjectToBoxEdge(xyz, dxyz6, xlo, xhi, ylo, yhi, zlo, zhi, xyzo);
3484  if (std::abs(xyzo[2]-zlo)>1.E-6) abort();
3485  }
void WorldBox(double *xlo, double *xhi, double *ylo, double *yhi, double *zlo, double *zhi) const
Fills the arguments with the boundaries of the world.
T abs(T value)
void ProjectToBoxEdge(const double xyz[], const double dxyz[], double xlo, double xhi, double ylo, double yhi, double zlo, double zhi, double xyzout[])
Definition: geo.h:49
geo::GeometryCore const * geom
pointer to geometry service provider
void geo::GeometryTestAlg::testStandardWirePos ( )
private

Definition at line 2009 of file GeometryTestAlg.cxx.

2010  {
2011  double xyz[3] = {0.};
2012  double xyzprev[3] = {0.};
2013  for(size_t cs = 0; cs < geom->Ncryostats(); ++cs){
2014  for(size_t t = 0; t < geom->Cryostat(cs).NTPC(); ++t){
2015  const geo::TPCGeo* tpc = &geom->Cryostat(cs).TPC(t);
2016 
2017  for (size_t i=0; i < tpc->Nplanes(); ++i) {
2018  const geo::PlaneGeo* plane = &tpc->Plane(i);
2019 
2020  for (size_t j = 1; j < plane->Nwires(); ++j) {
2021 
2022  geo::WireGeo const& wire = plane->Wire(j);
2023  geo::WireGeo const& wireprev = plane->Wire(j-1);
2024 
2025  wire.GetCenter(xyz);
2026  wireprev.GetCenter(xyzprev);
2027 
2028  // wires increase in +z order
2029  if(xyz[2] < xyzprev[2])
2030  throw cet::exception("WireOrderProblem") << "\n\twires do not increase in +z order in"
2031  << "Cryostat " << cs
2032  << ", TPC " << t
2033  << ", Plane " << i
2034  << "; at wire " << j << "\n";
2035 
2036  }// end loop over wires
2037  }// end loop over planes
2038  }// end loop over tpcs
2039  }// end loop over cryostats
2040 
2041 }
Geometry description of a TPC wireThe wire is a single straight segment on a wire plane...
Definition: WireGeo.h:65
WireGeo const & Wire(unsigned int iwire) const
Definition: PlaneGeo.cxx:506
unsigned int Nplanes() const
Number of planes in this tpc.
Definition: TPCGeo.h:165
Geometry information for a single TPC.
Definition: TPCGeo.h:38
unsigned int Ncryostats() const
Returns the number of cryostats in the detector.
Geometry information for a single wire plane.The plane is represented in the geometry by a solid whic...
Definition: PlaneGeo.h:82
unsigned int NTPC() const
Number of TPCs in this cryostat.
Definition: CryostatGeo.h:181
CryostatGeo const & Cryostat(geo::CryostatID const &cryoid) const
Returns the specified cryostat.
const TPCGeo & TPC(unsigned int itpc) const
Return the itpc&#39;th TPC in the cryostat.
Definition: CryostatGeo.cxx:93
unsigned int Nwires() const
Number of wires in this plane.
Definition: PlaneGeo.h:269
PlaneGeo const & Plane(geo::View_t view) const
Return the plane in the tpc with View_t view.
Definition: TPCGeo.cxx:263
void GetCenter(double *xyz, double localz=0.0) const
Fills the world coordinate of a point on the wire.
Definition: WireGeo.cxx:73
const char * cs
geo::GeometryCore const * geom
pointer to geometry service provider
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
void geo::GeometryTestAlg::testStepping ( )
private

Definition at line 3399 of file GeometryTestAlg.cxx.

3400  {
3401  //
3402  // Test stepping. Example is similar to what one would do for photon
3403  // transport. Rattles photons around inside the scintillator
3404  // bouncing them off walls.
3405  //
3406  double xyz[3] = {0.};
3407  double xyzWire[3] = {0.};
3408  double dxyz[3] = {0.};
3409  double dxyzWire[3] = {0, sin(0.1), cos(0.1)};
3410 
3411  geom->Plane(1).Wire(0).LocalToWorld(xyzWire,xyz);
3412  geom->Plane(1).Wire(0).LocalToWorldVect(dxyzWire,dxyz);
3413 
3414  mf::LogVerbatim("GeometryTest") << "\n\t" << xyz[0] << "\t" << xyz[1] << "\t" << xyz[2] ;
3415  mf::LogVerbatim("GeometryTest") << "\t" << dxyz[0] << "\t" << dxyz[1] << "\t" << dxyz[2];
3416 
3417  gGeoManager->InitTrack(xyz, dxyz);
3418  for (int i=0; i<10; ++i) {
3419  const double* pos = gGeoManager->GetCurrentPoint();
3420  const double* dir = gGeoManager->GetCurrentDirection();
3421  mf::LogVerbatim("GeometryTest") << "\tnode = "
3422  << gGeoManager->GetCurrentNode()->GetName()
3423  << "\n\t\tpos=" << "\t"
3424  << pos[0] << "\t"
3425  << pos[1] << "\t"
3426  << pos[2]
3427  << "\n\t\tdir=" << "\t"
3428  << dir[0] << "\t"
3429  << dir[1] << "\t"
3430  << dir[2]
3431  << "\n\t\tmat = "
3432  << gGeoManager->GetCurrentNode()->GetVolume()->GetMaterial()->GetName();
3433 
3434  gGeoManager->FindNextBoundary();
3435  gGeoManager->FindNormal();
3436  gGeoManager->Step(kTRUE,kTRUE);
3437  }
3438 
3439  xyz[0] = 306.108; xyz[1] = -7.23775; xyz[2] = 856.757;
3440  gGeoManager->InitTrack(xyz, dxyz);
3441  mf::LogVerbatim("GeometryTest") << "\tnode = "
3442  << gGeoManager->GetCurrentNode()->GetName()
3443  << "\n\tmat = "
3444  << gGeoManager->GetCurrentNode()->GetVolume()->GetMaterial()->GetName();
3445 
3446  gGeoManager->GetCurrentNode()->GetVolume()->GetMaterial()->Print();
3447 
3448  }
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
PlaneGeo const & Plane(unsigned int const p, unsigned int const tpc=0, unsigned int const cstat=0) const
Returns the specified wire.
WireGeo const & Wire(unsigned int iwire) const
Definition: PlaneGeo.cxx:506
string dir
void LocalToWorldVect(const double *wire, double *world) const
Transform direction vector from local to world.
Definition: WireGeo.h:363
void LocalToWorld(const double *wire, double *world) const
Transform point from local wire frame to world frame.
Definition: WireGeo.h:355
geo::GeometryCore const * geom
pointer to geometry service provider
void geo::GeometryTestAlg::testThirdPlane ( ) const
private

Definition at line 2729 of file GeometryTestAlg.cxx.

2729  {
2730  /*
2731  * This is a test for ThirdPlane() function, that returns the plane that is
2732  * not specified in the input.
2733  * Currently, the only implemented signature is designed for TPCs with 3
2734  * planes.
2735  *
2736  * The test strategy is to check all the TPC one by one:
2737  * - for all combinations of two planes, if the result is the expected one
2738  *
2739  * All tests are performed; at the end, the test is considered a failure
2740  * if any of the single tests failed.
2741  */
2742 
2743  unsigned int nErrors = 0;
2744  for (geo::GeometryCore::TPC_id_iterator iTPC(&*geom); iTPC; ++iTPC) {
2745  const geo::TPCID tpcid = *iTPC;
2746  const geo::TPCGeo& TPC = geom->TPC(tpcid);
2747 
2748  const unsigned int nPlanes = TPC.Nplanes();
2749  MF_LOG_DEBUG("GeometryTest") << tpcid << " (" << nPlanes << " planes)";
2750 
2751 
2752  for (geo::PlaneID::PlaneID_t iPlane1 = 0; iPlane1 < nPlanes; ++iPlane1) {
2753  geo::PlaneID pid1(tpcid, iPlane1);
2754 
2755  for (geo::PlaneID::PlaneID_t iPlane2 = 0; iPlane2 < nPlanes; ++iPlane2)
2756  {
2757  geo::PlaneID pid2(tpcid, iPlane2);
2758 
2759  const bool isValidInput = (nPlanes == 3) && (iPlane1 != iPlane2);
2760  bool bError = false;
2761  geo::PlaneID third_plane;
2762  try {
2763  third_plane = geom->ThirdPlane(pid1, pid2);
2764  }
2765  catch (cet::exception const& e) {
2766  if (isValidInput) throw;
2767  // we have gotten the error we were looking for
2768  // if "GeometryCore" is included in the categories of the exception
2769  bError = hasCategory(e, "GeometryCore");
2770  } // try...catch
2771 
2772  MF_LOG_TRACE("GeometryTest")
2773  << " [" << pid1 << "], [" << pid2 << "] => "
2774  << (bError? "error": std::string(third_plane));
2775  if (bError) continue; // we got what we wanted
2776 
2777  if (!bError && !isValidInput) {
2778  MF_LOG_ERROR("GeometryTest") << "ThirdPlane() on " << pid1
2779  << " and " << pid2 << " returned " << third_plane
2780  << ", while should have thrown an exception";
2781  ++nErrors;
2782  continue;
2783  } // if no error
2784 
2785  if (third_plane != tpcid) {
2786  MF_LOG_ERROR("GeometryTest") << "ThirdPlane() on " << pid1
2787  << " and " << pid2 << " returned " << third_plane
2788  << ", on a different TPC!!!";
2789  ++nErrors;
2790  }
2791  else if (!third_plane.isValid) {
2792  MF_LOG_ERROR("GeometryTest") << "ThirdPlane() on " << pid1
2793  << " and " << pid2 << " returned an invalid " << third_plane;
2794  ++nErrors;
2795  }
2796  else if (third_plane.Plane >= nPlanes) {
2797  MF_LOG_ERROR("GeometryTest") << "ThirdPlane() on " << pid1
2798  << " and " << pid2 << " returned " << third_plane
2799  << " with plane out of range";
2800  ++nErrors;
2801  }
2802  else if (third_plane == pid1) {
2803  MF_LOG_ERROR("GeometryTest") << "ThirdPlane() on " << pid1
2804  << " and " << pid2 << " returned " << third_plane
2805  << ", same as the first input";
2806  ++nErrors;
2807  }
2808  else if (third_plane == pid2) {
2809  MF_LOG_ERROR("GeometryTest") << "ThirdPlane() on " << pid1
2810  << " and " << pid2 << " returned " << third_plane
2811  << ", same as the second input";
2812  ++nErrors;
2813  }
2814 
2815  } // for plane 2
2816 
2817  } // for plane 1
2818 
2819  } // for TPC
2820 
2821  if (nErrors > 0) {
2822  throw cet::exception("GeoTestThirdPlane")
2823  << "Accumulated " << nErrors << " errors (see messages above)\n";
2824  }
2825 
2826  } // GeometryTestAlg::testThirdPlane()
std::string string
Definition: nybbler.cc:12
Base forward iterator browsing all TPC IDs in the detector.
Definition: GeometryCore.h:292
unsigned int Nplanes() const
Number of planes in this tpc.
Definition: TPCGeo.h:165
unsigned int PlaneID_t
Type for the ID number.
Definition: geo_types.h:473
The data type to uniquely identify a Plane.
Definition: geo_types.h:472
bool isValid
Whether this ID points to a valid element.
Definition: geo_types.h:211
Geometry information for a single TPC.
Definition: TPCGeo.h:38
#define MF_LOG_ERROR(category)
geo::PlaneID ThirdPlane(geo::PlaneID const &pid1, geo::PlaneID const &pid2) const
Returns the plane that is not in the specified arguments.
const double e
#define MF_LOG_TRACE(id)
The data type to uniquely identify a TPC.
Definition: geo_types.h:386
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:493
#define MF_LOG_DEBUG(id)
TPCGeo const & TPC(unsigned int const tpc=0, unsigned int const cstat=0) const
Returns the specified TPC.
geo::GeometryCore const * geom
pointer to geometry service provider
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
void geo::GeometryTestAlg::testThirdPlane_dTdW ( ) const
private

Definition at line 2830 of file GeometryTestAlg.cxx.

2830  {
2831  /*
2832  * This is a test for ThirdPlane_dTdW() function, that returns the apparent
2833  * slope on a wire plane, given the ones observed on other two planes.
2834  *
2835  * The test strategy is to check all the TPC one by one:
2836  * - if a query for planes on different cryostats fails
2837  * - if a query for planes on different TPCs fails
2838  * - for selected 3D points, compute the three dT/dW and verify them by
2839  * test these slopes by testThirdPlane__dTdW_at() function (see)
2840  *
2841  * All tests are performed; at the end, the test is considered a failure
2842  * if any of the single tests failed.
2843  */
2844 
2845  unsigned int nErrors = 0;
2846  for (geo::GeometryCore::TPC_id_iterator iTPC(&*geom); iTPC; ++iTPC) {
2847  const geo::TPCID tpcid = *iTPC;
2848  const geo::TPCGeo& TPC = geom->TPC(tpcid);
2849 
2850  const double driftVelocity = 0.1
2851  * ((TPC.DriftDirection() == geo::kNegX)? -1.: +1);
2852 
2853  const unsigned int NPlanes = TPC.Nplanes();
2854  MF_LOG_DEBUG("GeometryTest") << tpcid << " (" << NPlanes << " planes)";
2855 
2856  // sanity: planes on different cryostats
2857  if (tpcid.Cryostat < geom->Ncryostats() - 1) {
2858  geo::PlaneID p1 { tpcid, 0 }, p2 { tpcid.Cryostat + 1, tpcid.TPC, 1 };
2859  bool bError = false;
2860  double slope;
2861  try {
2862  slope = geom->ThirdPlane_dTdW(p1, +1.0, p2, -1.0);
2863  }
2864  catch (cet::exception const& e) {
2865  // we have gotten the error we were looking for
2866  // if "GeometryCore" is included in the categories of the exception
2867  bError = hasCategory(e, "GeometryCore");
2868  } // try...catch
2869  if (!bError) {
2870  MF_LOG_ERROR("GeometryTest") << "ThirdPlane_dTdW() on " << p1
2871  << " and " << p2 << " returned " << slope
2872  << ", while should have thrown an exception";
2873  ++nErrors;
2874  } // if no error
2875  } // if not the last cryostat
2876 
2877  // sanity: planes on different TPC
2878  if (tpcid.TPC < geom->NTPC(tpcid.Cryostat) - 1) {
2879  geo::PlaneID p1 { tpcid, 0 }, p2 { tpcid.Cryostat, tpcid.TPC + 1, 1 };
2880  bool bError = false;
2881  double slope;
2882  try {
2883  slope = geom->ThirdPlane_dTdW(p1, +1.0, p2, -1.0);
2884  }
2885  catch (cet::exception const& e) {
2886  // we have gotten the error we were looking for
2887  // if "GeometryCore" is included in the categories of the exception
2888  bError = hasCategory(e, "GeometryCore");
2889  } // try...catch
2890  if (!bError) {
2891  MF_LOG_ERROR("GeometryTest") << "ThirdPlane_dTdW() on " << p1
2892  << " and " << p2 << " returned " << slope
2893  << ", while should have thrown an exception";
2894  ++nErrors;
2895  } // if no error
2896  } // if not the last TPC in its cryostat
2897 
2898 
2899  // pick a point in the very middle of the TPC:
2900  // BUG the double brace syntax is required to work around clang bug 21629
2901  // (https://bugs.llvm.org/show_bug.cgi?id=21629)
2902  const std::array<double, 3> A
2903  = {{ TPC.CenterX(), TPC.CenterY(), TPC.CenterZ() }};
2904  // pick a radius half the way to the closest border
2905  const double radius
2906  = std::min({ TPC.HalfWidth(), TPC.HalfHeight(), TPC.Length()/2. }) / 2.;
2907 
2908  // I arbitrary decide that the second point will have dX equal size as
2909  // the radius, and on the positive x direction (may be negative dT)
2910  const double dX = radius;
2911  const double dT = driftVelocity * dX;
2912 
2913  // sample a circle of SplitAngles directions around A
2914  constexpr unsigned int NAngles = 19;
2915  const double start_angle = 0.05; // radians
2916  const double step_angle = 2. * util::pi<double>() / NAngles; // radians
2917 
2918  for (unsigned int iAngle = 0; iAngle < NAngles; ++iAngle) {
2919  const double angle = start_angle + iAngle * step_angle;
2920 
2921  // define B as a point "radius" far from A in the angle direction,
2922  // with some arbitrary and fixed dx offset
2923  // BUG the double brace syntax is required to work around clang bug 21629
2924  // (https://bugs.llvm.org/show_bug.cgi?id=21629)
2925  std::array<double, 3> B = {{
2926  A[0] + dX,
2927  A[1] + radius * std::sin(angle),
2928  A[2] + radius * std::cos(angle)
2929  }};
2930 
2931  // get the expectation; this function assumes a drift velocity of
2932  // 1 mm per tick by default; for the test, it does not matter
2933  std::vector<std::pair<geo::PlaneID, double>> dTdWs
2934  = ExpectedPlane_dTdW(A, B, driftVelocity);
2935 
2936  if (mf::isDebugEnabled()) {
2937  mf::LogTrace log("GeometryTest");
2938  log << "Expected dT/dW for a segment with " << radius
2939  << " cm long projection at "
2940  << angle << " rad, and dT " << dT << " cm:";
2941  for (auto const& dTdW_info: dTdWs)
2942  log << " " << dTdW_info.first << " slope:" << dTdW_info.second;
2943  } // if debug
2944 
2945  // run the test
2946  nErrors += testThirdPlane_dTdW_at(dTdWs);
2947 
2948  } // for angle
2949 
2950  } // for TPC
2951 
2952  if (nErrors > 0) {
2953  throw cet::exception("GeoTestThirdPlane_dTdW")
2954  << "Accumulated " << nErrors << " errors (see messages above)\n";
2955  }
2956 
2957  } // GeometryTestAlg::testThirdPlane_dTdW()
Base forward iterator browsing all TPC IDs in the detector.
Definition: GeometryCore.h:292
unsigned int Nplanes() const
Number of planes in this tpc.
Definition: TPCGeo.h:165
The data type to uniquely identify a Plane.
Definition: geo_types.h:472
Geometry information for a single TPC.
Definition: TPCGeo.h:38
double CenterX() const
Returns the world x coordinate of the center of the box.
Definition: BoxBoundedGeo.h:94
#define MF_LOG_ERROR(category)
double CenterZ() const
Returns the world z coordinate of the center of the box.
CryostatID_t Cryostat
Index of cryostat.
Definition: geo_types.h:212
std::vector< std::pair< geo::PlaneID, double > > ExpectedPlane_dTdW(std::array< double, 3 > const &A, std::array< double, 3 > const &B, const double driftVelocity=-0.1) const
Returns dT/dW expected from the specified segment A-to-B.
Drift towards negative X values.
Definition: geo_types.h:162
unsigned int Ncryostats() const
Returns the number of cryostats in the detector.
unsigned int testThirdPlane_dTdW_at(std::vector< std::pair< geo::PlaneID, double >> const &plane_dTdW) const
Performs the third plane slope test with a single configuration.
double Length() const
Length is associated with z coordinate [cm].
Definition: TPCGeo.h:115
const double e
bool isDebugEnabled()
The data type to uniquely identify a TPC.
Definition: geo_types.h:386
DriftDirection_t DriftDirection() const
Returns an enumerator value describing the drift direction.
Definition: TPCGeo.h:144
unsigned int NTPC(unsigned int cstat=0) const
Returns the total number of TPCs in the specified cryostat.
double HalfHeight() const
Height is associated with y coordinate [cm].
Definition: TPCGeo.h:111
T min(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:55
double CenterY() const
Returns the world y coordinate of the center of the box.
#define MF_LOG_DEBUG(id)
TPCGeo const & TPC(unsigned int const tpc=0, unsigned int const cstat=0) const
Returns the specified TPC.
TPCID_t TPC
Index of the TPC within its cryostat.
Definition: geo_types.h:406
double ThirdPlane_dTdW(geo::PlaneID const &pid1, double slope1, geo::PlaneID const &pid2, double slope2, geo::PlaneID const &output_plane) const
Returns dT/dW on the third plane, given it in the other two.
double HalfWidth() const
Width is associated with x coordinate [cm].
Definition: TPCGeo.h:107
geo::GeometryCore const * geom
pointer to geometry service provider
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
unsigned int geo::GeometryTestAlg::testThirdPlane_dTdW_at ( std::vector< std::pair< geo::PlaneID, double >> const &  plane_dTdW) const
private

Performs the third plane slope test with a single configuration.

Definition at line 3028 of file GeometryTestAlg.cxx.

3029  {
3030  /*
3031  * This function tests that for every combination of planes, the expected
3032  * slope is returned within some tolerance.
3033  * It returns the number of mistakes found.
3034  *
3035  * The parameter is a list if pair of expected slope on the paired plane.
3036  */
3037 
3038  unsigned int nErrors = 0;
3039  for (std::pair<geo::PlaneID, double> const& input1: plane_dTdW) {
3040  for (std::pair<geo::PlaneID, double> const& input2: plane_dTdW) {
3041 
3042  const bool bValidInput = input1.first != input2.first;
3043 
3044  for (std::pair<geo::PlaneID, double> const& output: plane_dTdW) {
3045  bool bError = false;
3046  double output_slope = 0.;
3047  try {
3048  output_slope = geom->ThirdPlane_dTdW(
3049  input1.first, input1.second,
3050  input2.first, input2.second,
3051  output.first
3052  );
3053  }
3054  catch (cet::exception const& e) {
3055  // catch only "GeometryCore" category, and only if input is faulty;
3056  // otherwise, rethrow
3057  if (bValidInput) throw;
3058  bError = hasCategory(e, "GeometryCore");
3059  if (!bError) throw;
3060  MF_LOG_TRACE("GeometryTest")
3061  << input1.first << " slope:" << input1.second
3062  << " " << input2.first << " slope:" << input2.second
3063  << " => exception";
3064  continue;
3065  }
3066 
3067  if (!bValidInput && !bError) {
3068  MF_LOG_ERROR("GeometryTest")
3069  << "GeometryCore::ThirdPlane_dTdW() on "
3070  << input1.first << " and " << input2.first
3071  << " should have thrown an exception, it returned "
3072  << output_slope << " instead";
3073  ++nErrors;
3074  continue;
3075  } // if faulty input and no error
3076 
3077  MF_LOG_TRACE("GeometryTest")
3078  << input1.first << " slope:" << input1.second
3079  << " " << input2.first << " slope:" << input2.second
3080  << " => " << output.first << " slope:" << output_slope;
3081  if (((output.second == 0.) && (output_slope > 1e-3))
3082  || std::abs(output_slope/output.second - 1.) > 1e-3)
3083  {
3084  MF_LOG_ERROR("testThirdPlane_dTdW_at")
3085  << "GeometryCore::ThirdPlane_dTdW(): "
3086  << input1.first << " slope:" << input1.second
3087  << " " << input2.first << " slope:" << input2.second
3088  << " => " << output.first << " slope:" << output_slope
3089  << " (expected: " << output.second << ")";
3090  ++nErrors;
3091  } // if too far
3092 
3093  // now test the automatic detection of the other plane
3094 
3095  } // for output
3096  } // for second input
3097  } // for first input
3098  return nErrors;
3099  } // GeometryTestAlg::testThirdPlane_dTdW_at()
#define MF_LOG_ERROR(category)
T abs(T value)
const double e
#define MF_LOG_TRACE(id)
double ThirdPlane_dTdW(geo::PlaneID const &pid1, double slope1, geo::PlaneID const &pid2, double slope2, geo::PlaneID const &output_plane) const
Returns dT/dW on the third plane, given it in the other two.
geo::GeometryCore const * geom
pointer to geometry service provider
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
void geo::GeometryTestAlg::testTPC ( geo::CryostatID const &  cid)
private

Definition at line 855 of file GeometryTestAlg.cxx.

856  {
857  geo::CryostatGeo const& cryo = geom->Cryostat(cid);
858 
859  mf::LogVerbatim("GeometryTest") << "\tThere are " << cryo.NTPC()
860  << " TPCs in the detector";
861 
862  for(size_t t = 0; t < cryo.NTPC(); ++t){
863  geo::TPCID const tpcid(cid, t);
864  geo::TPCGeo const& tpc = cryo.TPC(tpcid);
865  decltype(auto) activeCenter = tpc.GetActiveVolumeCenter();
866 
867  {
868  mf::LogVerbatim log { "GeometryTest" };
869  log
870  << "\n\t\tTPC " << tpcid
871  << " " << geom->GetLArTPCVolumeName(tpcid)
872  << " has " << tpc.Nplanes() << " planes."
873  << "\n\t\tTPC location: ( "
874  << tpc.MinX() << " ; " << tpc.MinY() << " ; "<< tpc.MinZ()
875  << " ) => ( "
876  << tpc.MaxX() << " ; " << tpc.MaxY() << " ; "<< tpc.MaxZ()
877  << " ) [cm]"
878  << "\n\t\tTPC Dimensions (W x H x L, cm): "
879  << tpc.Width() << " (" << directionName(tpc.WidthDir()) << ")"
880  << " x " << tpc.Height() << " (" << directionName(tpc.HeightDir()) << ")"
881  << " x " << tpc.Length() << " (" << directionName(tpc.LengthDir()) << ")"
882  << "\n\t\tTPC Active Dimensions: "
883  << 2.*tpc.ActiveHalfWidth() << " x " << 2.*tpc.ActiveHalfHeight() << " x " << tpc.ActiveLength()
884  << " around ( " << activeCenter.X() << " ; " << activeCenter.Y()
885  << " ; "<< activeCenter.Z() << " ) cm"
886  ;
887  if (fComputeMass)
888  log << "\n\t\tTPC mass: " << tpc.ActiveMass();
889  log
890  << "\n\t\tTPC drift distance: " << tpc.DriftDistance()
891  << ", direction: " << tpc.DriftDir();
892  }
893 
894  for(size_t p = 0; p < tpc.Nplanes(); ++p) {
895  geo::PlaneGeo const& plane = tpc.Plane(p);
896 
897  // first line indented with two tabs, the others with two more spaces;
898  // very verbose (8)
899  plane.PrintPlaneInfo
900  (mf::LogVerbatim("GeometryTest") << "\t\t", "\t\t ", 8);
901 
902  mf::LogVerbatim("GeometryTest")
903  << "\t\t pitch from plane 0 is " << tpc.Plane0Pitch(p);
904 
905  } // for plane
907  if (dir == geo::kNegX) {
908  mf::LogVerbatim("GeometryTest")
909  << "\t\tdrift direction is towards negative values: "
910  << tpc.DriftDir();
911  }
912  else if(dir == geo::kPosX) {
913  mf::LogVerbatim("GeometryTest")
914  << "\t\tdrift direction is towards positive values: "
915  << tpc.DriftDir();
916  }
917  else{
918  throw cet::exception("UnknownDriftDirection") << "\t\tdrift direction is unknown\n";
919  }
920 
921  MF_LOG_DEBUG("GeometryTest") << "\t testing PositionToTPC...";
922  // pick a position in the middle of the TPC in the world coordinates
923  double worldLoc[3] = {0.};
924  double localLoc[3] = {0.};
925  tpc.LocalToWorld(localLoc, worldLoc);
926 
927  const unsigned int tpcNo = cryo.FindTPCAtPosition(worldLoc, 1+1.e-4);
928 
929  if(tpcNo != t)
930  throw cet::exception("BadTPCLookupFromPosition") << "TPC look up returned tpc = "
931  << tpcNo << " should be " << t << "\n";
932 
933  MF_LOG_DEBUG("GeometryTest") << "done.";
934  } // for TPC
935 
936  return;
937  }
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
Vector DriftDir() const
Returns the direction of the drift (vector pointing toward the planes).
Definition: TPCGeo.h:773
decltype(auto) LengthDir() const
Returns the direction Length() is measured on.
Definition: TPCGeo.h:138
double ActiveHalfHeight() const
Half height (associated with y coordinate) of active TPC volume [cm].
Definition: TPCGeo.h:99
unsigned int Nplanes() const
Number of planes in this tpc.
Definition: TPCGeo.h:165
double MinX() const
Returns the world x coordinate of the start of the box.
Definition: BoxBoundedGeo.h:88
Geometry information for a single TPC.
Definition: TPCGeo.h:38
double MaxX() const
Returns the world x coordinate of the end of the box.
Definition: BoxBoundedGeo.h:91
double ActiveMass() const
Definition: TPCGeo.h:118
string dir
Drift towards negative X values.
Definition: geo_types.h:162
Geometry information for a single cryostat.
Definition: CryostatGeo.h:43
double Width() const
Width is associated with x coordinate [cm].
Definition: TPCGeo.h:109
double Height() const
Height is associated with y coordinate [cm].
Definition: TPCGeo.h:113
double Length() const
Length is associated with z coordinate [cm].
Definition: TPCGeo.h:115
const double e
enum geo::driftdir DriftDirection_t
Drift direction: positive or negative.
double ActiveHalfWidth() const
Half width (associated with x coordinate) of active TPC volume [cm].
Definition: TPCGeo.h:95
Geometry information for a single wire plane.The plane is represented in the geometry by a solid whic...
Definition: PlaneGeo.h:82
double MinZ() const
Returns the world z coordinate of the start of the box.
p
Definition: test.py:223
unsigned int NTPC() const
Number of TPCs in this cryostat.
Definition: CryostatGeo.h:181
CryostatGeo const & Cryostat(geo::CryostatID const &cryoid) const
Returns the specified cryostat.
The data type to uniquely identify a TPC.
Definition: geo_types.h:386
double Plane0Pitch(unsigned int p) const
Definition: TPCGeo.cxx:324
DriftDirection_t DriftDirection() const
Returns an enumerator value describing the drift direction.
Definition: TPCGeo.h:144
double ActiveLength() const
Length (associated with z coordinate) of active TPC volume [cm].
Definition: TPCGeo.h:103
double MaxY() const
Returns the world y coordinate of the end of the box.
geo::TPCID::TPCID_t FindTPCAtPosition(double const worldLoc[3], double const wiggle) const
Returns the index of the TPC at specified location.
double DriftDistance() const
Definition: TPCGeo.h:155
decltype(auto) HeightDir() const
Returns the direction Height() is measured on.
Definition: TPCGeo.h:131
Drift towards positive X values.
Definition: geo_types.h:161
void PrintPlaneInfo(Stream &&out, std::string indent="", unsigned int verbosity=1) const
Prints information about this plane.
Definition: PlaneGeo.h:1539
const TPCGeo & TPC(unsigned int itpc) const
Return the itpc&#39;th TPC in the cryostat.
Definition: CryostatGeo.cxx:93
double MaxZ() const
Returns the world z coordinate of the end of the box.
#define MF_LOG_DEBUG(id)
PlaneGeo const & Plane(geo::View_t view) const
Return the plane in the tpc with View_t view.
Definition: TPCGeo.cxx:263
double MinY() const
Returns the world y coordinate of the start of the box.
void LocalToWorld(const double *tpc, double *world) const
Transform point from local TPC frame to world frame.
Definition: TPCGeo.h:563
geo::GeometryCore const * geom
pointer to geometry service provider
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
std::string GetLArTPCVolumeName(geo::TPCID const &tpcid) const
Return the name of specified LAr TPC volume.
bool fComputeMass
Whether to print the detector mass.
decltype(auto) WidthDir() const
Returns the direction Width() is measured on.
Definition: TPCGeo.h:124
void geo::GeometryTestAlg::testWireCoordAngle ( ) const
private

Definition at line 1379 of file GeometryTestAlg.cxx.

1379  {
1380  /*
1381  * Tests that the angle PhiZ() actually points to the next wire.
1382  *
1383  * The test, for each plane, performs the following:
1384  * - pick the middle wire, verify that we can get the expected wire
1385  * coordinate for its centre
1386  * - move one wire pitch away from the centre in the direction determined
1387  * by PhiZ(), verify that the coordinate increases by 1
1388  */
1389 
1390  for (geo::PlaneID const& planeid: geom->IteratePlaneIDs()) {
1391 
1392  geo::PlaneGeo const& plane = geom->Plane(planeid);
1393 
1394  // define the wires to work with
1395  const unsigned int nWires = plane.Nwires();
1396 
1397  geo::WireID middle_wire_id(planeid, nWires / 2);
1398  geo::WireID next_wire_id(planeid, nWires / 2 + 1);
1399 
1400  if (next_wire_id.Wire >= nWires) {
1401  throw cet::exception("WeirdGeometry")
1402  << "Plane " << std::string(planeid) << " has only " << nWires
1403  << " wires?!?\n";
1404  }
1405 
1406 
1407  geo::WireGeo const& middle_wire = geom->Wire(middle_wire_id);
1408  decltype(auto) middle_wire_center = middle_wire.GetCenter();
1409  MF_LOG_TRACE("GeometryTest")
1410  << "Center of " << middle_wire_id << " at " << middle_wire_center;
1411 
1412  // cross check: we can find the middle wire
1413  const double middle_coord = geom->WireCoordinate
1414  (middle_wire_center, planeid);
1415 
1416  if (std::abs(middle_coord - double(middle_wire_id.Wire)) > 2e-3) {
1417  throw cet::exception("WireCoordAngle")
1418  << "Center of " << std::string(middle_wire_id) << " at ("
1419  << middle_wire_center[0]
1420  << "; " << middle_wire_center[1] << "; " << middle_wire_center[2]
1421  << ") has wire coordinate " << middle_coord
1422  << " (" << middle_wire_id.Wire << " expected)\n";
1423  } // if
1424 
1425  // the check: this coordinate should lie on the next wire
1426  const double pitch = plane.WirePitch();
1427  decltype(auto) wireCoordDir = plane.GetIncreasingWireDirection();
1428 
1429  MF_LOG_TRACE("GeometryTest")
1430  << " pitch: " << pitch << " wire coord dir: cos(phi_z): "
1431  << wireCoordDir;
1432 
1433  auto on_next_wire = middle_wire_center + pitch * wireCoordDir;
1434 
1435  const double next_coord = geom->WireCoordinate(on_next_wire, planeid);
1436 
1437  if (std::abs(next_coord - double(next_wire_id.Wire)) > 2e-3) {
1438  std::cerr
1439  << " pitch: " << pitch << " wire coord dir: " << wireCoordDir
1440  << "\n start wire: " << middle_wire_id
1441  << " centered at " << middle_wire_center
1442  << "\n querying point " << on_next_wire
1443  << std::endl;
1444  throw cet::exception("WireCoordAngle")
1445  << "Position " << on_next_wire
1446  << " is expected to be on wire " << std::string(next_wire_id)
1447  << " but it has wire coordinate " << next_coord << "\n";
1448  } // if
1449 
1450  } // for
1451 
1452  } // GeometryTestAlg::testWireCoordAngle()
Geometry description of a TPC wireThe wire is a single straight segment on a wire plane...
Definition: WireGeo.h:65
PlaneGeo const & Plane(unsigned int const p, unsigned int const tpc=0, unsigned int const cstat=0) const
Returns the specified wire.
std::string string
Definition: nybbler.cc:12
auto coord(Vector &v, unsigned int n) noexcept
Returns an object to manage the coordinate n of a vector.
The data type to uniquely identify a Plane.
Definition: geo_types.h:472
STL namespace.
string dir
IteratorBox< plane_id_iterator,&GeometryCore::begin_plane_id,&GeometryCore::end_plane_id > IteratePlaneIDs() const
Enables ranged-for loops on all plane IDs of the detector.
T abs(T value)
const double e
Geometry information for a single wire plane.The plane is represented in the geometry by a solid whic...
Definition: PlaneGeo.h:82
#define MF_LOG_TRACE(id)
WireGeo const & Wire(geo::WireID const &wireid) const
Returns the specified wire.
unsigned int Nwires() const
Number of wires in this plane.
Definition: PlaneGeo.h:269
if(!yymsg) yymsg
geo::GeometryCore const * geom
pointer to geometry service provider
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
double WirePitch() const
Return the wire pitch (in centimeters). It is assumed constant.
Definition: PlaneGeo.h:411
QTextStream & endl(QTextStream &s)
static std::array< double, 3 > GetIncreasingWireDirection(const geo::PlaneGeo &plane)
Returns the direction on plane orthogonal to wires where wire number increases.
void geo::GeometryTestAlg::testWireCoordFromPlane ( ) const
private

Definition at line 1065 of file GeometryTestAlg.cxx.

1065  {
1066 
1067  //
1068  // For each wire:
1069  //
1070  // * picks points lying on the planes including a wire and the normal to the
1071  // wire plane (which have all the same wire coordinate)
1072  //
1073  // * tests that the coordinates are as expected (wire number times pitch)
1074  //
1075 
1076  unsigned int nErrors = 0;
1077  for (geo::PlaneGeo const& plane: geom->IteratePlanes()) {
1078 
1079  auto const nWires = plane.Nwires();
1080  auto const wirePitch = plane.WirePitch();
1081 
1082  double const driftDistance = geom->TPC(plane.ID()).DriftDistance();
1083 
1084  decltype(auto) planeNormal = plane.GetNormalDirection();
1085 
1086  for (geo::WireID::WireID_t wireNo = 0; wireNo < nWires; ++wireNo) {
1087 
1088  geo::WireGeo const& wire = plane.Wire(wireNo);
1089 
1090  double const expected = wireNo * wirePitch;
1091 
1092  // sample 7 points on wire
1093  constexpr int shifts = 3;
1094  double const step = wire.HalfL() / (std::abs(shifts) + 1);
1095  for (int iOfs = -shifts; iOfs <= shifts; ++iOfs) {
1096 
1097  double const offset = iOfs * step;
1098 
1099  auto const basePoint = wire.GetPositionFromCenter(offset);
1100 
1101  // at 4 different distances from the plane
1102  constexpr int quotas = 4;
1103  double const jump = driftDistance / (std::abs(quotas) + 1);
1104  for (int iQuota = 0; iQuota < quotas; ++iQuota) {
1105 
1106  // translate the point along the normal to the plane;
1107  // this should not change the result
1108  auto const point = basePoint + iQuota * jump * planeNormal;
1109 
1110  double const distance = plane.PlaneCoordinate(point);
1111 
1112  if (std::abs(distance - expected) > 1e-4) {
1113  mf::LogProblem("GeometryTestAlg") << "Point " << point
1114  << " (offset: " << iOfs << "x" << step << ", at " << iQuota
1115  << "x" << jump << " from plane) is reported to be " << distance
1116  << " cm far from wire " << plane.ID() << " W: " << wireNo
1117  << " (" << expected << " expected)";
1118  ++nErrors;
1119  } // if unexpected
1120 
1121  } // for quotas
1122 
1123  } // for iOfs
1124 
1125  } // for wires
1126 
1127  } // for planes
1128 
1129  if (nErrors > 0) {
1130  throw cet::exception("GeometryTestAlg")
1131  << "testWireCoordFromPlane() accumulated " << nErrors
1132  << " errors (see messages above)\n";
1133  }
1134 
1135  } // GeometryTestAlg::testWireCoordFromPlane()
Geometry description of a TPC wireThe wire is a single straight segment on a wire plane...
Definition: WireGeo.h:65
IteratorBox< plane_iterator,&GeometryCore::begin_plane,&GeometryCore::end_plane > IteratePlanes() const
Enables ranged-for loops on all planes of the detector.
const char expected[]
Definition: Exception_t.cc:22
Point GetPositionFromCenter(double localz) const
Returns the position (world coordinate) of a point on the wire.
Definition: WireGeo.h:182
MaybeLogger_< ELseverityLevel::ELsev_error, true > LogProblem
T abs(T value)
const double e
IDparameter< geo::WireID > WireID
Member type of validated geo::WireID parameter.
Geometry information for a single wire plane.The plane is represented in the geometry by a solid whic...
Definition: PlaneGeo.h:82
double distance(double x1, double y1, double z1, double x2, double y2, double z2)
double HalfL() const
Returns half the length of the wire [cm].
Definition: WireGeo.h:128
for(std::string line;std::getline(inFile, line);)
Definition: regex_t.cc:37
TPCGeo const & TPC(unsigned int const tpc=0, unsigned int const cstat=0) const
Returns the specified TPC.
LArSoft geometry interface.
Definition: ChannelGeo.h:16
geo::GeometryCore const * geom
pointer to geometry service provider
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
void geo::GeometryTestAlg::testWireIntersection ( ) const
private

Definition at line 2371 of file GeometryTestAlg.cxx.

2371  {
2372  /*
2373  * This is a test for geo::GeometryCore::WireIDsIntersect() and
2374  * geo::WireGeo::IntersectionWith() methods, that return whether two wires
2375  * intersect, and where.
2376  *
2377  * The test strategy is to check all the TPC one by one:
2378  * - if a query for wires on different cryostats fails
2379  * - if a query for wires on different TPCs fails
2380  * - if a query for wires on the same plane fails
2381  * - for points at the centre of a grid SplitY x SplitZ on the wire planes,
2382  * test these point by testWireIntersectionAt() function (see)
2383  * All tests are performed; at the end, the test is considered a failure
2384  * if any of the single tests failed.
2385  */
2386 
2387  unsigned int nErrors = 0;
2388  for (geo::TPCGeo const& TPC: geom->IterateTPCs()) {
2389 
2390  MF_LOG_DEBUG("GeometryTest") << "Wire intersection test on " << TPC.ID();
2391 
2392  // sanity: wires on different cryostats
2393  if (TPC.ID().Cryostat < geom->Ncryostats() - 1) {
2394 
2395  geo::WireID const w1 { geo::PlaneID{ TPC.ID(), 0 }, 0 };
2396  geo::WireGeo const& wire1 = geom->Wire(w1);
2397  geo::Vector_t const& wireDir = wire1.Direction<geo::Vector_t>();
2398 
2399  geo::CryostatGeo const& otherCryo
2400  = geom->Cryostat(TPC.ID().Cryostat + 1);
2401  geo::PlaneGeo const* otherPlane = nullptr;
2402  for (geo::PlaneGeo const& plane: geom->IteratePlanes(otherCryo.ID())) {
2403  if (isWireAlignedToPlaneDirections(plane, wireDir)) continue;
2404  otherPlane = &plane;
2405  break;
2406  }
2407  if (otherPlane) {
2408  geo::WireID const w2 { otherPlane->ID(), 1 };
2409  MF_LOG_TRACE("GeometryTest")
2410  << "Off cryostat test (" << w1 << "): chosen wire " << w2;
2411  geo::Point_t xingPoint;
2412  if (geom->WireIDsIntersect(w1, w2, xingPoint)) {
2413  MF_LOG_ERROR("GeometryTest") << "WireIDsIntersect() on " << w1
2414  << " and " << w2 << " returned " << xingPoint
2415  << " cm, while should have reported no intersection at all";
2416  ++nErrors;
2417  } // if intersect
2418  try {
2419  // the value of result is not checked here
2420  wire1.IntersectionWith(geom->Wire(w2));
2421  }
2422  catch(...) {
2423  MF_LOG_ERROR("GeometryTest") << "WiresIntersect() on " << w1
2424  << " and " << w2 << " threw an exception, which should not have";
2425  ++nErrors;
2426  }
2427  }
2428  else {
2429  MF_LOG_WARNING("GeometryTest")
2430  << "No wire plane found in " << otherCryo.ID()
2431  << " with wires not aligned with " << w1.asPlaneID()
2432  << ", " << wireDir << "; off-cryostat sanity check skipped.";
2433  }
2434 
2435  } // if not the last cryostat
2436 
2437  // sanity: wires on different TPC
2438  if (TPC.ID().TPC < geom->NTPC(TPC.ID().asCryostatID()) - 1) {
2439 
2440  geo::PlaneID const refPlaneID { TPC.ID(), 0 };
2441  geo::WireID const w1 { refPlaneID, 0 };
2442  geo::WireGeo const& wire1 = geom->Wire(w1);
2443  geo::Vector_t const& wireDir = wire1.Direction<geo::Vector_t>();
2444 
2445  geo::CryostatGeo const& cryo = geom->Cryostat(TPC.ID());
2446  geo::PlaneGeo const* otherPlane = nullptr;
2447  for (geo::PlaneGeo const& plane: geom->IteratePlanes(cryo.ID())) {
2448  if (plane.ID().asTPCID() == TPC.ID()) continue; // on the same TPC
2449  if (isWireAlignedToPlaneDirections(plane, wireDir)) continue;
2450  otherPlane = &plane;
2451  break;
2452  }
2453  if (otherPlane) {
2454  geo::WireID const w2 { otherPlane->ID(), 1 };
2455  MF_LOG_TRACE("GeometryTest")
2456  << "Off TPC test (" << w1 << "): chosen wire " << w2;
2457  geo::Point_t xingPoint;
2458  if (geom->WireIDsIntersect(w1, w2, xingPoint)) {
2459  MF_LOG_ERROR("GeometryTest") << "WireIDsIntersect() on " << w1
2460  << " and " << w2 << " returned " << xingPoint
2461  << ", while should have reported no intersection at all";
2462  ++nErrors;
2463  } // if intersect
2464  try {
2465  // the value of result is not checked here
2466  wire1.IntersectionWith(geom->Wire(w2));
2467  }
2468  catch(...) {
2469  MF_LOG_ERROR("GeometryTest") << "WiresIntersect() on " << w1
2470  << " and " << w2 << " threw an exception, which should not have";
2471  ++nErrors;
2472  }
2473  }
2474  else {
2475  MF_LOG_WARNING("GeometryTest")
2476  << "No wire plane found in any TPC of " << cryo.ID()
2477  << " with wires not aligned with " << refPlaneID
2478  << ", " << wireDir << "; off-TPC sanity check skipped.";
2479  }
2480  } // if not the last TPC
2481 
2482  // sanity: wires on same plane
2483  for (geo::PlaneGeo const& plane: TPC.IteratePlanes()) {
2484  geo::WireID const w1 { plane.ID(), 0 }, w2 { plane.ID(), 1 };
2485  geo::Point_t xingPoint;
2486  if (geom->WireIDsIntersect(w1, w2, xingPoint)) {
2487  MF_LOG_ERROR("GeometryTest") << "WireIDsIntersect() on " << w1
2488  << " and " << w2 << " returned " << xingPoint
2489  << ", while should have reported no intersection at all";
2490  ++nErrors;
2491  } // if intersect
2492  // prerequisites of WiresIntersect() are not met here, no test possible
2493  } // for all planes
2494 
2495 
2496  // sample the area covered by all planes, split into SplitW and SplitD
2497  // rectangles; drift coordinate is chosen roughly in the middle of the TPC
2498  geo::PlaneGeo const& refPlane = TPC.SmallestPlane();
2499  constexpr unsigned int SplitW = 19, SplitD = 17;
2500 
2501  auto const driftOffset = -TPC.DriftDistance() / 2.0 * TPC.DriftDir();
2502  auto const refPoint = refPlane.GetCenter() + driftOffset;
2503 
2504  decltype(auto) coverage = refPlane.ActiveArea();
2505  const double stepW = coverage.width.length() / SplitW;
2506  const double stepD = coverage.depth.length() / SplitD;
2507  const int stepsW = SplitW / 2;
2508  const int stepsD = SplitD / 2;
2509 
2510  // let's pick a point:
2511  for (int iW = -stepsW; iW <= +stepsW; ++iW) {
2512 
2513  auto const widthOffset = (iW * stepW) * refPlane.WidthDir();
2514 
2515  for (int iD = -stepsD; iD < +stepsD; ++iD) {
2516 
2517  auto const depthOffset = (iD * stepD) * refPlane.DepthDir();
2518 
2519  auto const point = refPoint + widthOffset + depthOffset;
2520 
2521  // finally, let's test this point...
2522  nErrors += testWireIntersectionAt(TPC, point);
2523  } // for y
2524  } // for z
2525  } // for TPC
2526 
2527  if (nErrors > 0) {
2528  throw cet::exception("GeoTestWireIntersection")
2529  << "Accumulated " << nErrors << " errors (see messages above)\n";
2530  }
2531 
2532  } // GeometryTestAlg::testWireIntersection()
Geometry description of a TPC wireThe wire is a single straight segment on a wire plane...
Definition: WireGeo.h:65
IteratorBox< plane_iterator,&GeometryCore::begin_plane,&GeometryCore::end_plane > IteratePlanes() const
Enables ranged-for loops on all planes of the detector.
The data type to uniquely identify a Plane.
Definition: geo_types.h:472
Geometry information for a single TPC.
Definition: TPCGeo.h:38
#define MF_LOG_ERROR(category)
ROOT::Math::DisplacementVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Vector_t
Type for representation of momenta in 3D space.
Definition: geo_vectors.h:164
Geometry information for a single cryostat.
Definition: CryostatGeo.h:43
unsigned int Ncryostats() const
Returns the number of cryostats in the detector.
bool isWireAlignedToPlaneDirections(geo::PlaneGeo const &plane, geo::Vector_t const &wireDir) const
Helper function for testWireIntersection().
Point IntersectionWith(geo::WireGeo const &other) const
Returns the point of this wire that is closest to other wire.
Definition: WireGeo.h:637
IteratorBox< TPC_iterator,&GeometryCore::begin_TPC,&GeometryCore::end_TPC > IterateTPCs() const
Enables ranged-for loops on all TPCs of the detector.
Point GetCenter() const
Returns the centre of the wire plane in world coordinates [cm].
Definition: PlaneGeo.h:479
Geometry information for a single wire plane.The plane is represented in the geometry by a solid whic...
Definition: PlaneGeo.h:82
#define MF_LOG_TRACE(id)
CryostatGeo const & Cryostat(geo::CryostatID const &cryoid) const
Returns the specified cryostat.
Vector Direction() const
Definition: WireGeo.h:587
unsigned int NTPC(unsigned int cstat=0) const
Returns the total number of TPCs in the specified cryostat.
Vector DepthDir() const
Return the direction of plane depth.
Definition: PlaneGeo.h:236
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Point_t
Type for representation of position in physical 3D space.
Definition: geo_vectors.h:184
bool WireIDsIntersect(WireID const &wid1, WireID const &wid2, geo::Point_t &intersection) const
Computes the intersection between two wires.
WireGeo const & Wire(geo::WireID const &wireid) const
Returns the specified wire.
for(std::string line;std::getline(inFile, line);)
Definition: regex_t.cc:37
unsigned int testWireIntersectionAt(geo::TPCGeo const &TPC, TVector3 const &point) const
Performs the wire intersection test at a single point.
geo::PlaneID const & ID() const
Returns the identifier of this plane.
Definition: PlaneGeo.h:203
#define MF_LOG_DEBUG(id)
Vector WidthDir() const
Return the direction of plane width.
Definition: PlaneGeo.h:221
#define MF_LOG_WARNING(category)
geo::GeometryCore const * geom
pointer to geometry service provider
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
unsigned int geo::GeometryTestAlg::testWireIntersectionAt ( geo::TPCGeo const &  TPC,
TVector3 const &  point 
) const
private

Performs the wire intersection test at a single point.

Definition at line 2536 of file GeometryTestAlg.cxx.

2537  {
2538  /* Tests WireIDsIntersect() and WiresIntersect() on the specified point on
2539  * the wire planes of a given TPC.
2540  *
2541  * The test follows this strategy:
2542  * - find the ID of the wires closest to the point on each plane
2543  * - for all wire plane pairing, ask for the intersection between the wires
2544  * - fail if the returned point is farther than half a pitch from the
2545  * original point
2546  *
2547  * This function returns the number of accumulated failures.
2548  */
2549 
2550  using geo::vect::dot;
2551 
2552  unsigned int nErrors = 0;
2553 
2554  const unsigned int NPlanes = TPC.Nplanes();
2555 
2556  bool const bDriftOnX = (TPC.DriftDir<geo::Vector_t>() == geo::Xaxis())
2557  || (TPC.DriftDir<geo::Vector_t>() == -geo::Xaxis());
2558 
2559  // collect information per plane:
2560  std::vector<double> WirePitch(NPlanes); // for convenience
2561  std::vector
2562  <decltype(std::declval<geo::PlaneGeo>().GetIncreasingWireDirection())>
2563  WireCoordDirs(NPlanes);
2564  std::vector<geo::WireID> WireIDs; // ID of the closest wire
2565  WireIDs.reserve(NPlanes);
2566  std::vector<double> WireDistances(NPlanes); // distance from the closest wire
2567  for (unsigned int iPlane = 0; iPlane < NPlanes; ++iPlane) {
2568  const geo::PlaneGeo& plane = TPC.Plane(iPlane);
2569  WireCoordDirs[iPlane] = plane.GetIncreasingWireDirection();
2570  WirePitch[iPlane] = plane.WirePitch();
2571 
2572  const double WireDistance = geom->WireCoordinate(point, plane.ID());
2573  WireIDs.emplace_back(plane.ID(), (unsigned int) std::round(WireDistance));
2574  WireDistances[iPlane]
2575  = (WireDistance - std::round(WireDistance)) * WirePitch[iPlane];
2576 
2577  MF_LOG_DEBUG("GeometryTest") << "Nearest wire to " << point
2578  << " on plane " << std::string(plane.ID())
2579  << " (pitch: " << WirePitch[iPlane]
2580  << ", coord.dir.=" << WireCoordDirs[iPlane]
2581  << ") is " << WireIDs[iPlane] << " (position: " << WireDistance << ")";
2582  } // for planes
2583 
2584  // test all the combinations
2585 
2587  auto vectorIs = lar::util::makeVector3DComparison(coordIs);
2588  for (unsigned int iPlane1 = 0; iPlane1 < NPlanes; ++iPlane1) {
2589 
2590  const geo::WireID& w1 = WireIDs[iPlane1];
2591  geo::PlaneGeo const& plane1 = TPC.Plane(w1);
2592  geo::WireGeo const& w1obj = plane1.Wire(w1);
2593 
2594  for (unsigned int iPlane2 = iPlane1 + 1; iPlane2 < NPlanes; ++iPlane2) {
2595  const geo::WireID& w2 = WireIDs[iPlane2];
2596  geo::PlaneGeo const& plane2 = TPC.Plane(w2);
2597  geo::WireGeo const& w2obj = plane2.Wire(w2);
2598 
2599  geo::Point_t xingPoint;
2600  if (!geom->WireIDsIntersect(w1, w2, xingPoint)) {
2601  MF_LOG_ERROR("GeometryTest") << "Wires " << w1 << " and " << w2
2602  << " should intersect around " << point << " of TPC " << TPC.ID()
2603  << ", but they seem not to intersect at all!";
2604  ++nErrors;
2605  continue;
2606  }
2607  geo::Point_t xingPoint2 = xingPoint; // matching point on plane 2
2608  plane2.DriftPoint(xingPoint2);
2609 
2610  if (bDriftOnX) { // legacy code check
2611 
2612  geo::WireIDIntersection widIntersect;
2613  if (!geom->WireIDsIntersect(w1, w2, widIntersect)) {
2614  MF_LOG_ERROR("GeometryTest") << "Legacy check: wires " << w1 << " and "
2615  << w2 << " should intersect around " << point << " of TPC "
2616  << TPC.ID() << ", but they seem not to intersect at all!";
2617  ++nErrors;
2618  }
2619 
2620  if (coordIs.nonEqual(widIntersect.y, xingPoint.Y())
2621  || coordIs.nonEqual(widIntersect.z, xingPoint.Z()))
2622  {
2623  MF_LOG_ERROR("GeometryTest") << "Legacy check: wires " << w1 << " and "
2624  << w2 << " should intersect around " << point << " of TPC "
2625  << TPC.ID() << ", but legacy code says (?, "
2626  << widIntersect.y << ", " << widIntersect.z << ")!";
2627  ++nErrors;
2628  }
2629 
2630  } // bDriftOnX
2631 
2632 
2633  geo::Point_t xingPointInv;
2634  if (!geom->WireIDsIntersect(w2, w1, xingPointInv)) {
2635  MF_LOG_ERROR("GeometryTest") << "Wires " << w2 << " and " << w1
2636  << " (reversed test) should intersect around " << point
2637  << " of TPC " << TPC.ID()
2638  << ", but they seem not to intersect at all!";
2639  ++nErrors;
2640  continue;
2641  }
2642  if (vectorIs.nonEqual(xingPointInv, xingPoint2)) {
2643  MF_LOG_ERROR("GeometryTest")
2644  << "WireIDsIntersect() gives different intersections for "
2645  << w1 << " and " << w2
2646  << ": " << xingPoint2 << " (direct+shift) and " << xingPointInv
2647  << " (reversed)";
2648  ++nErrors;
2649  continue;
2650  }
2651 
2652  // the expected distance between the probe point and the
2653  // intersection point is geometrically determined, given the distances
2654  // of the probe point from the two wires and the angle between the wires
2655  // the formula is a mix between the Carnot theorem and sine definition;
2656  // the definition of the angles is tricky though, and we rely on the
2657  // strong vector definition enforced in the geometry to get the proper
2658  // "cosine" and corresponding sine
2659  const double d1 = WireDistances[iPlane1], d2 = WireDistances[iPlane2],
2660  cosAlpha = dot(WireCoordDirs[iPlane1], WireCoordDirs[iPlane2]);
2661  const double expected_d = std::sqrt(
2662  (cet::square(d1) + cet::square(d2) - 2.0 * d1 * d2 * cosAlpha) / (1 - cet::square(cosAlpha))
2663  );
2664  // the actual distance we have found:
2665  double const d = plane1.VectorProjection(xingPoint - geo::vect::toPoint(point)).R();
2666  MF_LOG_DEBUG("GeometryTest")
2667  << " - wires " << w1 << " and " << w2 << " intersect at " << xingPoint
2668  << ", " << d << " cm far from starting point (expected: "
2669  << expected_d << ")";
2670 
2671  // precision of the test is an issue; the 10^-3 x pitch threshold
2672  // is roughly tuned so that we don't get errors
2674  (std::max(WirePitch[iPlane1], WirePitch[iPlane2]) * 1e-3); // cm
2675  if (wireCoordIs.nonEqual(d, expected_d)) {
2676  MF_LOG_ERROR("GeometryTest")
2677  << "wires " << w1 << " and " << w2 << " intersect at " << xingPoint
2678  << ", "
2679  << d << " cm far from starting point: too far from the expected "
2680  << expected_d << " cm!";
2681  ++nErrors;
2682  continue;
2683  } // if too far
2684 
2685  geo::Point_t objXingPoint;
2686 
2687  // test that geo::WiresIntersection() gives the same result as
2688  // the already validated one from geom->WireIDsIntersect()
2689  objXingPoint = geo::WiresIntersection(w1obj, w2obj);
2690  if (vectorIs.nonEqual(objXingPoint, xingPoint)) {
2691  MF_LOG_ERROR("GeometryTest")
2692  << "geo::WiresIntersection() gives wrong intersection for "
2693  << w1 << " and " << w2
2694  << ": " << objXingPoint << " vs. " << xingPoint << " (expected)";
2695  ++nErrors;
2696  continue;
2697  }
2698 
2699  // test that geo::WireGeo::IntersectionWith() gives the same result as
2700  // the already validated one from geom->WireIDsIntersect()
2701  objXingPoint = w1obj.IntersectionWith(w2obj);
2702  if (vectorIs.nonEqual(objXingPoint, xingPoint)) {
2703  MF_LOG_ERROR("GeometryTest")
2704  << "geo::WireGeo[" << w1
2705  << "]::IntersectionWith() gives wrong intersection with " << w2
2706  << ": " << objXingPoint << " vs. " << xingPoint << " (expected)";
2707  ++nErrors;
2708  continue;
2709  }
2710 
2711  objXingPoint = w2obj.IntersectionWith(w1obj);
2712  if (vectorIs.nonEqual(objXingPoint, xingPointInv)) {
2713  MF_LOG_ERROR("GeometryTest")
2714  << "geo::WireGeo[" << w2
2715  << "]::IntersectionWith() gives wrong intersection with " << w1
2716  << ": " << objXingPoint << " vs. " << xingPoint << " (expected)";
2717  ++nErrors;
2718  continue;
2719  }
2720 
2721  } // for iPlane2
2722  } // for iPlane1
2723 
2724  return nErrors;
2725  } // GeometryTestAlg::testWireIntersectionAt()
geo::Length_t WireCoordinate(double YPos, double ZPos, geo::PlaneID const &planeid) const
Returns the index of the nearest wire to the specified position.
Geometry description of a TPC wireThe wire is a single straight segment on a wire plane...
Definition: WireGeo.h:65
double z
z position of intersection
Definition: geo_types.h:805
constexpr auto dot(Vector const &a, Vector const &b)
Return cross product of two vectors.
WireGeo const & Wire(unsigned int iwire) const
Definition: PlaneGeo.cxx:506
std::string string
Definition: nybbler.cc:12
Provides simple real number checks.
struct vector vector
::geo::Point_t toPoint(Point const &p)
Convert the specified point into a geo::Point_t.
#define MF_LOG_ERROR(category)
ROOT::Math::DisplacementVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Vector_t
Type for representation of momenta in 3D space.
Definition: geo_vectors.h:164
constexpr T square(T x)
Definition: pow.h:21
Vector GetIncreasingWireDirection() const
Returns the direction of increasing wires.
Definition: PlaneGeo.h:457
void DriftPoint(geo::Point_t &position, double distance) const
Shifts the position of an electron drifted by a distance.
Definition: PlaneGeo.h:643
geo::Point_t WiresIntersection(geo::WireGeo const &wireA, geo::WireGeo const &wireB)
Returns the point of wireA that is closest to wireB.
Definition: WireGeo.h:654
Point IntersectionWith(geo::WireGeo const &other) const
Returns the point of this wire that is closest to other wire.
Definition: WireGeo.h:637
const double e
auto makeVector3DComparison(RealType threshold)
Creates a Vector3DComparison from a RealComparisons object.
Geometry information for a single wire plane.The plane is represented in the geometry by a solid whic...
Definition: PlaneGeo.h:82
constexpr Vector Xaxis()
Returns a x axis vector of the specified type.
Definition: geo_vectors.h:215
static int max(int a, int b)
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Point_t
Type for representation of position in physical 3D space.
Definition: geo_vectors.h:184
bool WireIDsIntersect(WireID const &wid1, WireID const &wid2, geo::Point_t &intersection) const
Computes the intersection between two wires.
WireCoordProjection_t VectorProjection(geo::Vector_t const &v) const
Definition: PlaneGeo.h:963
double y
y position of intersection
Definition: geo_types.h:804
geo::PlaneID const & ID() const
Returns the identifier of this plane.
Definition: PlaneGeo.h:203
#define MF_LOG_DEBUG(id)
geo::GeometryCore const * geom
pointer to geometry service provider
double WirePitch() const
Return the wire pitch (in centimeters). It is assumed constant.
Definition: PlaneGeo.h:411
static std::array< double, 3 > GetIncreasingWireDirection(const geo::PlaneGeo &plane)
Returns the direction on plane orthogonal to wires where wire number increases.
void geo::GeometryTestAlg::testWireOrientations ( ) const
private

Definition at line 1013 of file GeometryTestAlg.cxx.

1013  {
1014  /*
1015  * The test verifies that all the wires in the geometry respect the
1016  * orientation requirement described in geo::WireGeo documentation:
1017  *
1018  * { (wire direction) , (wire coordinate increase), (plane normal) }
1019  *
1020  * be a positively defined base.
1021  *
1022  */
1023 
1024  unsigned int nErrors = 0;
1025  for (geo::PlaneGeo const& plane: geom->IteratePlanes()) {
1026 
1027  // this funny way declares a reference or not, depending on return type
1028  decltype(auto) planeNormal = plane.GetNormalDirection();
1029  decltype(auto) wireCoordDir = plane.GetIncreasingWireDirection();
1030 
1031  unsigned int nWires = plane.Nwires();
1032  for (unsigned int wireNo = 0; wireNo < nWires; ++wireNo) {
1033 
1034  geo::WireGeo const& wire = plane.Wire(wireNo);
1035 
1036  double positive = wire.Direction().Cross(wireCoordDir).Dot(planeNormal);
1037 
1038  if (positive < 0.5) {
1039  ++nErrors;
1040 
1041  // detail the problem; details of the plane should be read in the
1042  // output from other tests
1043  decltype(auto) wireDir = wire.Direction();
1045  << "Wire " << plane.ID() << " W: " << wireNo
1046  << " has direction ( " << wireDir
1047  << " , yielding to a non-positive plane-coordinate definition"
1048  << " (l x w . n = " << positive << ")";
1049  } // if error
1050 
1051  } // for wire
1052 
1053  } // for plane
1054 
1055  if (nErrors > 0) {
1056  throw cet::exception("GeometryTestAlg")
1057  << "testWireOrientations() accumulated " << nErrors
1058  << " errors (see messages above)\n";
1059  }
1060 
1061  } // GeometryTestAlg::testWireOrientations()
Geometry description of a TPC wireThe wire is a single straight segment on a wire plane...
Definition: WireGeo.h:65
GeometryTestAlg(fhicl::ParameterSet const &pset)
IteratorBox< plane_iterator,&GeometryCore::begin_plane,&GeometryCore::end_plane > IteratePlanes() const
Enables ranged-for loops on all planes of the detector.
unsigned int ID
static QStrList * l
Definition: config.cpp:1044
MaybeLogger_< ELseverityLevel::ELsev_error, true > LogProblem
std::void_t< T > n
const double a
Geometry information for a single wire plane.The plane is represented in the geometry by a solid whic...
Definition: PlaneGeo.h:82
Vector Direction() const
Definition: WireGeo.h:587
for(std::string line;std::getline(inFile, line);)
Definition: regex_t.cc:37
list x
Definition: train.py:276
Direction
Definition: AssnsIter.h:13
if(!yymsg) yymsg
geo::GeometryCore const * geom
pointer to geometry service provider
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
static std::array< double, 3 > GetIncreasingWireDirection(const geo::PlaneGeo &plane)
Returns the direction on plane orthogonal to wires where wire number increases.
void geo::GeometryTestAlg::testWirePitch ( )
private

Definition at line 3103 of file GeometryTestAlg.cxx.

3104  {
3105  // loop over all planes and wires to be sure the pitch is consistent
3106  unsigned int nPitchErrors = 0;
3107 
3108  if (fExpectedWirePitches.empty()) {
3109  // hard code the value we think it should be for each detector;
3110  // this is legacy and you should not add anything:
3111  // add the expectation to the FHiCL configuration of the test instead
3112  if(geom->DetectorName() == "bo") {
3113  fExpectedWirePitches = { 0.46977, 0.46977, 0.46977 };
3114  }
3115  if (!fExpectedWirePitches.empty()) {
3116  mf::LogInfo("WirePitch")
3117  << "Using legacy wire pitch parameters hard-coded for the detector '"
3118  << geom->DetectorName() << "'";
3119  }
3120  }
3121  if (fExpectedWirePitches.empty()) {
3122  mf::LogWarning("WirePitch")
3123  << "no expected wire pitch; I'll just check that they are all the same";
3124  }
3125  else {
3126  mf::LogInfo log("WirePitch");
3127  log << "Expected wire pitch per plane, in centimetres:";
3128  for (double pitch: fExpectedWirePitches) log << " " << pitch;
3129  log << " [...]";
3130  }
3131 
3132  for (geo::PlaneID const& planeid: geom->IteratePlaneIDs()) {
3133 
3134  geo::PlaneGeo const& plane = geom->Plane(planeid);
3135  const unsigned int nWires = plane.Nwires();
3136  if (nWires < 2) continue;
3137 
3138  geo::WireGeo const* pWire = &(plane.Wire(0));
3139 
3140  // which pitch to expect:
3141  // - if they did not tell us anything:
3142  // get the one from the first two wires
3143  // - if they did tell something, but not for this plane:
3144  // get the last pitch they told us
3145  // - if they told us about this plane: well, then use it!
3146  double expectedPitch = 0.;
3147  if (fExpectedWirePitches.empty()) {
3148  geo::WireGeo const& wire1 = plane.Wire(1); // pWire now points to wire0
3149  expectedPitch = geo::WireGeo::WirePitch(*pWire, wire1);
3150  MF_LOG_DEBUG("WirePitch")
3151  << "Wire pitch on " << planeid << ": " << expectedPitch << " cm";
3152  }
3153  else if (planeid.Plane < fExpectedWirePitches.size())
3154  expectedPitch = fExpectedWirePitches[planeid.Plane];
3155  else
3156  expectedPitch = fExpectedWirePitches.back();
3157 
3158  geo::WireID::WireID_t w = 0; // wire number
3159  while (++w < nWires) {
3160  geo::WireGeo const* pPrevWire = pWire;
3161  pWire = &(plane.Wire(w));
3162 
3163  const double thisPitch = std::abs(pWire->DistanceFrom(*pPrevWire));
3164  if (std::abs(thisPitch - expectedPitch) > 1e-5) {
3165  mf::LogProblem("WirePitch") << "ERROR: on plane " << planeid
3166  << " pitch between wires W:" << (w-1) << " and W:" << w
3167  << " is " << thisPitch << " cm, not " << expectedPitch
3168  << " as expected!";
3169  ++nPitchErrors;
3170  } // if unexpected pitch
3171  } // while
3172 
3173  } // for
3174 
3175  if (nPitchErrors > 0) {
3176  throw cet::exception("UnexpectedWirePitch")
3177  << "unexpected pitches between " << nPitchErrors << " wires!";
3178  } // end loop over planes
3179 
3180  } // GeometryTestAlg::testWirePitch()
Geometry description of a TPC wireThe wire is a single straight segment on a wire plane...
Definition: WireGeo.h:65
PlaneGeo const & Plane(unsigned int const p, unsigned int const tpc=0, unsigned int const cstat=0) const
Returns the specified wire.
std::vector< double > fExpectedWirePitches
wire pitch on each plane
WireGeo const & Wire(unsigned int iwire) const
Definition: PlaneGeo.cxx:506
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
double DistanceFrom(geo::WireGeo const &wire) const
Returns 3D distance from the specified wire.
Definition: WireGeo.cxx:113
The data type to uniquely identify a Plane.
Definition: geo_types.h:472
static double WirePitch(geo::WireGeo const &w1, geo::WireGeo const &w2)
Returns the pitch (distance on y/z plane) between two wires, in cm.
Definition: WireGeo.h:476
IteratorBox< plane_id_iterator,&GeometryCore::begin_plane_id,&GeometryCore::end_plane_id > IteratePlaneIDs() const
Enables ranged-for loops on all plane IDs of the detector.
MaybeLogger_< ELseverityLevel::ELsev_error, true > LogProblem
T abs(T value)
const double e
std::string DetectorName() const
Returns a string with the name of the detector, as configured.
Geometry information for a single wire plane.The plane is represented in the geometry by a solid whic...
Definition: PlaneGeo.h:82
unsigned int Nwires() const
Number of wires in this plane.
Definition: PlaneGeo.h:269
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
#define MF_LOG_DEBUG(id)
unsigned int WireID_t
Type for the ID number.
Definition: geo_types.h:561
geo::GeometryCore const * geom
pointer to geometry service provider
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33

Member Data Documentation

bool geo::GeometryTestAlg::fComputeMass = true
private

Whether to print the detector mass.

Definition at line 134 of file GeometryTestAlg.h.

bool geo::GeometryTestAlg::fDisableValidWireIDcheck
private

disable test on out-of-world NearestWire()

Definition at line 129 of file GeometryTestAlg.h.

std::vector<double> geo::GeometryTestAlg::fExpectedPlanePitches
private

plane pitch on each plane

Definition at line 132 of file GeometryTestAlg.h.

std::vector<double> geo::GeometryTestAlg::fExpectedWirePitches
private

wire pitch on each plane

Definition at line 131 of file GeometryTestAlg.h.

std::set<std::string> geo::GeometryTestAlg::fNonFatalExceptions
private

Definition at line 130 of file GeometryTestAlg.h.

testing::NameSelector geo::GeometryTestAlg::fRunTests
private

test filter

Definition at line 137 of file GeometryTestAlg.h.

geo::GeometryCore const* geo::GeometryTestAlg::geom
private

pointer to geometry service provider

Definition at line 127 of file GeometryTestAlg.h.


The documentation for this class was generated from the following files: