GeometryTestAlg.cxx
Go to the documentation of this file.
1 /**
2  * @file GeometryTestAlg.cxx
3  * @brief Unit test for geometry functionalities: implementation file
4  * @date 2011/02/17
5  * @author brebel@fnal.gov
6  * @see GeometryTestAlg.h
7  */
8 
9 // our header
10 #include "test/Geometry/GeometryTestAlg.h"
11 
12 // LArSoft includes
13 #include "larcorealg/Geometry/Decomposer.h"
14 #include "larcorealg/Geometry/Exceptions.h"
15 #include "larcorealg/Geometry/SimpleGeo.h"
16 #include "larcorealg/Geometry/GeometryCore.h"
17 #include "larcorealg/Geometry/CryostatGeo.h"
18 #include "larcorealg/Geometry/TPCGeo.h"
19 #include "larcorealg/Geometry/PlaneGeo.h"
20 #include "larcorealg/Geometry/WireGeo.h"
21 #include "larcorealg/Geometry/AuxDetGeo.h"
22 #include "larcorealg/Geometry/AuxDetSensitiveGeo.h"
23 #include "larcorealg/Geometry/geo.h"
24 #include "larcorealg/Geometry/geo_vectors_utils.h"
25 #include "larcorealg/CoreUtils/RealComparisons.h"
26 #include "larcorealg/CoreUtils/DumpUtils.h" // lar::dump::vector3D(), ...
27 #include "larcoreobj/SimpleTypesAndConstants/geo_vectors.h"
28 #include "larcoreobj/SimpleTypesAndConstants/readout_types.h" // readout::ROPID
29 #include "larcoreobj/SimpleTypesAndConstants/RawTypes.h" // raw::ChannelID_t
30 #include "larcoreobj/SimpleTypesAndConstants/geo_types.h"
31 #include "larcoreobj/SimpleTypesAndConstants/PhysicalConstants.h" // util::pi<>
32 
33 // Framework includes
34 #include "cetlib/pow.h"
35 #include "cetlib_except/exception.h"
36 #include "fhiclcpp/ParameterSet.h"
37 #include "messagefacility/MessageLogger/MessageLogger.h"
38 
39 // ROOT includes
40 #include "TGeoManager.h"
41 #include "TStopwatch.h"
42 #include "TGeoNode.h"
43 #include "TGeoVolume.h"
44 #include "TClass.h"
45 #include "Math/GenVector/DisplacementVector2D.h"
46 #include "Math/GenVector/DisplacementVector3D.h"
47 #include "Math/GenVector/PositionVector3D.h"
48 #include "RtypesCore.h" // for kTRUE
49 #include "TGeoMaterial.h"
50 #include "TGeoMatrix.h"
51 #include "TGeoShape.h"
52 #include "TObjArray.h"
53 #include "TVector2.h"
54 
55 // C/C++ standard libraries
56 #include <stdint.h>
57 #include <stdlib.h> // for abort
58 #include <memory>
59 #include <cmath>
60 #include <vector>
61 #include <iterator> // std::inserter()
62 #include <algorithm> // std::copy()
63 #include <set>
64 #include <array>
65 #include <string>
66 #include <sstream>
67 #include <iostream>
68 #include <limits> // std::numeric_limits<>
69 #include <initializer_list>
70 
71 
72 //------------------------------------------------------------------------------
73 // custom stream insertion operators; they need to be defined in the "right"
74 // namespace
75 std::ostream& operator<< (std::ostream& out, TVector3 const& v) {
76  out << "( " << v.X() << " ; " << v.Y() << " ; " << v.Z() << " )";
77  return out;
78 } // operator<< (TVector3)
79 
80 std::ostream& operator<< (std::ostream& out, TVector2 const& v) {
81  out << "( " << v.X() << " ; " << v.Y() << " )";
82  return out;
83 } // operator<< (TVector2)
84 
85 
86 //------------------------------------------------------------------------------
87 namespace {
88 
89  /// Returns whether the CET exception e contains the specified category cat
90  bool hasCategory(cet::exception const& e, std::string const& cat) {
91  for (auto const& e_category: e.history())
92  if (e_category == cat) return true;
93  return false;
94  } // hasCategory()
95 
96 
97  /// Returns a convenience string for the specified direction
98  template <typename Vector>
99  std::string directionName(Vector const& v) {
100  if (v == geo::Xaxis<Vector>()) return "x";
101  if (v == geo::Yaxis<Vector>()) return "y";
102  if (v == geo::Zaxis<Vector>()) return "z";
103  if (v == -geo::Xaxis<Vector>()) return "-x";
104  if (v == -geo::Yaxis<Vector>()) return "-y";
105  if (v == -geo::Zaxis<Vector>()) return "-z";
106  std::ostringstream sstr;
107  sstr << v;
108  return sstr.str();
109  } // directionName()
110 
111 } // local namespace
112 
113 
114 namespace geo{
115 
116 
117  //......................................................................
118  GeometryTestAlg::GeometryTestAlg(fhicl::ParameterSet const& pset)
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;
144  fRunTests.PrintConfiguration(sstr);
145 
146  mf::LogInfo("GeometryTestAlg") << "Test selection:" << sstr.str();
147 
148  } // GeometryTestAlg::GeometryTestAlg()
149 
150  //......................................................................
151  unsigned int GeometryTestAlg::Run()
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:";
178  printDetectorIntro();
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...";
221  testPlaneDirections();
222  MF_LOG_INFO("GeometryTest") << "complete.";
223  }
224 
225  if (shouldRunTests("WireOrientations")) {
226  MF_LOG_INFO("GeometryTest") << "test wire orientations...";
227  testWireOrientations();
228  MF_LOG_INFO("GeometryTest") << "complete.";
229  }
230 
231  if (shouldRunTests("ChannelToROP")) {
232  MF_LOG_INFO("GeometryTest") << "test channel to ROP and back ...";
233  testChannelToROP();
234  MF_LOG_INFO("GeometryTest") << "complete.";
235  }
236 
237  if (shouldRunTests("ChannelToWire")) {
238  MF_LOG_INFO("GeometryTest") << "test channel to plane wire and back ...";
239  testChannelToWire();
240  MF_LOG_INFO("GeometryTest") << "complete.";
241  }
242 
243  if (shouldRunTests("FindPlaneCenters")) {
244  MF_LOG_INFO("GeometryTest") << "test find plane centers...";
245  testFindPlaneCenters();
246  MF_LOG_INFO("GeometryTest") << "complete.";
247  }
248 
249  if (shouldRunTests("WireCoordFromPlane")) {
250  MF_LOG_INFO("GeometryTest") << "test PlaneGeo::WireCoordinate...";
251  testWireCoordFromPlane();
252  MF_LOG_INFO("GeometryTest") << "complete.";
253  }
254 
255  if (shouldRunTests("ParallelWires")) {
256  MF_LOG_INFO("GeometryTest") << "test wire parallelism...";
257  testParallelWires();
258  MF_LOG_INFO("GeometryTest") << "complete.";
259  }
260 
261  if (shouldRunTests("PlanePointDecomposition")) {
262  MF_LOG_INFO("GeometryTest") << "test plane point decomposition...";
263  testPlanePointDecomposition();
264  MF_LOG_INFO("GeometryTest") << "complete.";
265  }
266 
267  if (shouldRunTests("PlaneProjections")) {
268  MF_LOG_INFO("GeometryTest") << "test PlaneGeo::PointProjection...";
269  testPlaneProjection();
270  MF_LOG_INFO("GeometryTest") << "complete.";
271  }
272 
273  if (shouldRunTests("WireCoordAngle")) {
274  MF_LOG_INFO("GeometryTest") << "testWireCoordAngle...";
275  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...";
300  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...";
312  testThirdPlane_dTdW();
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...";
324  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...";
348  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()
382 
383 
384 
385  //......................................................................
386  void GeometryTestAlg::printDetectorIntro() const {
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()
413 
414 
415  //......................................................................
416  void GeometryTestAlg::printChannelSummary()
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  }
442 
443  //......................................................................
444  void GeometryTestAlg::printVolBounds()
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  }
473 
474 
475 
476  //......................................................................
477  // great sanity check for geometry, only call in analyze when debugging
478  void GeometryTestAlg::printDetDim()
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: "
486  << geom->CryostatHalfHeight(c);
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  }
499 
500  //......................................................................
501  // great sanity check for volume sorting, only call in analyze when debugging
502  void GeometryTestAlg::printWirePos()
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  }
521 
522  //......................................................................
523  // great insanity: print all wires in a TPC
524  void GeometryTestAlg::printWiresInTPC
525  (const geo::TPCGeo& tpc, std::string indent /* = "" */) const
526  {
527  const unsigned int nPlanes = tpc.Nplanes();
528 
529  tpc.PrintTPCInfo(
530  mf::LogVerbatim("GeometryTest") << indent,
531  indent, geo::TPCGeo::MaxVerbosity
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()
565 
566 
567  void GeometryTestAlg::printAllGeometry() const {
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
582  printAuxiliaryDetectors();
583  mf::LogVerbatim("GeometryTest") << "End of detector "
584  << geom->DetectorName() << " geometry.";
585  } // GeometryTestAlg::printAllGeometry()
586 
587 
588  //......................................................................
589  void GeometryTestAlg::printAuxiliaryDetectors() const {
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()
601 
602 
603  //......................................................................
604  template <typename Stream>
605  void GeometryTestAlg::printAuxDetGeo(
606  Stream&& out, geo::AuxDetGeo const& auxDet,
607  std::string indent, std::string firstIndent
608  ) const
609  {
610 
611  lar::util::RealComparisons<double> coordIs(1e-4);
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 ";
630  printAuxDetSensitiveGeo(
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()
647 
648 
649  //......................................................................
650  template <typename Stream>
651  void GeometryTestAlg::printAuxDetSensitiveGeo(
652  Stream&& out, geo::AuxDetSensitiveGeo const& auxDetSens,
653  std::string /* indent */, std::string firstIndent
654  ) const
655  {
656 
657  lar::util::RealComparisons<double> coordIs(1e-4);
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()
672 
673  //......................................................................
674  void GeometryTestAlg::testCryostat()
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  }
728 
729  //......................................................................
730  unsigned int GeometryTestAlg::testFindWorldVolumes() {
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()
759 
760 
761  unsigned int GeometryTestAlg::testFindCryostatVolumes() {
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()
787 
788 
789  unsigned int GeometryTestAlg::testFindTPCvolumePaths() {
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()
830 
831 
832  void GeometryTestAlg::testFindVolumes() {
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()
852 
853 
854  //......................................................................
855  void GeometryTestAlg::testTPC(geo::CryostatID const& cid)
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
906  geo::DriftDirection_t dir = tpc.DriftDirection();
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  }
938 
939 
940  //......................................................................
941  void GeometryTestAlg::testPlaneDirections() const {
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 
953  lar::util::RealComparisons<double> coordIs(1e-5);
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()
1010 
1011 
1012  //......................................................................
1013  void GeometryTestAlg::testWireOrientations() const {
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();
1044  mf::LogProblem("GeometryTestAlg")
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()
1062 
1063 
1064  //......................................................................
1065  void GeometryTestAlg::testWireCoordFromPlane() const {
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()
1136 
1137 
1138  //......................................................................
1139  void GeometryTestAlg::testParallelWires() const {
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()
1169 
1170 
1171  //......................................................................
1172  void GeometryTestAlg::testPlanePointDecomposition() const {
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 
1193  lar::util::RealComparisons<double> coordIs(1e-5);
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()
1376 
1377 
1378  //......................................................................
1379  void GeometryTestAlg::testWireCoordAngle() const {
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()
1453 
1454 
1455  //......................................................................
1456  void GeometryTestAlg::testChannelToROP() const {
1457 
1458  // test that an invalid channel yields an invalid ROP
1459  try {
1460  readout::ROPID invalidROP = geom->ChannelToROP(raw::InvalidChannelID);
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()
1493 
1494  //......................................................................
1495  void GeometryTestAlg::testChannelToWire() const
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 
1506  raw::ChannelID_t channel = geom->PlaneWireToChannel(testWireID);
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()
1581 
1582  //......................................................................
1583  void GeometryTestAlg::testFindPlaneCenters()
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  }
1594 
1595  //......................................................................
1596  void GeometryTestAlg::testPlaneProjectionReference() const {
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()
1632 
1633 
1634  //......................................................................
1635  void GeometryTestAlg::testPlanePointDecompositionFrame() const {
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 
1654  lar::util::RealComparisons<double> coordIs(1e-5);
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()
1852 
1853 
1854  //......................................................................
1855  void GeometryTestAlg::testPlaneProjectionOnFrame() const {
1856 
1857  //
1858  // Tests:
1859  //
1860  // 1. containment (isProjectionOnPlane())
1861  //
1862  //
1863  // 2. capping by closest point
1864  //
1865 
1866  lar::util::RealComparisons<double> coordIs(1e-5);
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()
1984 
1985  //......................................................................
1986  void GeometryTestAlg::testPlaneProjection() const {
1987 
1988  //
1989  // Check the definition of the reference
1990  //
1991 
1992  testPlaneProjectionReference();
1993 
1994  //
1995  // Check the projections and point composition in the plane frame reference
1996  //
1997  testPlanePointDecompositionFrame();
1998 
1999  //
2000  // Check containment
2001  //
2002  testPlaneProjectionOnFrame();
2003 
2004 
2005  } // GeometryTestAlg::testPlaneProjection()
2006 
2007 
2008  //......................................................................
2009  void GeometryTestAlg::testStandardWirePos()
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 }
2042 
2043  //......................................................................
2044  void GeometryTestAlg::testAPAWirePos()
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  }
2085 
2086 
2087  //......................................................................
2088  inline std::array<double, 3> GeometryTestAlg::GetIncreasingWireDirection
2089  (const geo::PlaneGeo& plane)
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()
2097 
2098 
2099  //......................................................................
2100  void GeometryTestAlg::testNearestWire()
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) {
2324  if (fDisableValidWireIDcheck) {
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  }
2345 
2346  //......................................................................
2347  bool GeometryTestAlg::isWireAlignedToPlaneDirections
2348  (geo::PlaneGeo const& plane, geo::Vector_t const& wireDir) const
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()
2369 
2370 
2371  void GeometryTestAlg::testWireIntersection() const {
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()
2533 
2534 
2535  unsigned int GeometryTestAlg::testWireIntersectionAt
2536  (const geo::TPCGeo& TPC, TVector3 const& point) const
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 
2586  lar::util::RealComparisons<double> coordIs(1e-3);
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
2673  lar::util::RealComparisons<double> wireCoordIs
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()
2726 
2727 
2728  //......................................................................
2729  void GeometryTestAlg::testThirdPlane() const {
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()
2827 
2828 
2829  //......................................................................
2830  void GeometryTestAlg::testThirdPlane_dTdW() const {
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()
2958 
2959 
2960  std::vector<std::pair<geo::PlaneID, double>>
2961  GeometryTestAlg::ExpectedPlane_dTdW(
2962  std::array<double, 3> const& A, std::array<double, 3> const& B,
2963  const double driftVelocity /* = 0.1 */
2964  ) const
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()
3025 
3026 
3027  unsigned int GeometryTestAlg::testThirdPlane_dTdW_at
3028  (std::vector<std::pair<geo::PlaneID, double>> const& plane_dTdW) const
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()
3100 
3101 
3102  //......................................................................
3103  void GeometryTestAlg::testWirePitch()
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()
3181 
3182  //......................................................................
3183  void GeometryTestAlg::testInterWireProjectedDistance() const {
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 {{
3209  geo::PlaneGeo::WireCoordProjection_t{ 0.00, 1.00 },
3210  geo::PlaneGeo::WireCoordProjection_t{ 0.75, 1.00 },
3211  geo::PlaneGeo::WireCoordProjection_t{ 1.00, 1.00 },
3212  geo::PlaneGeo::WireCoordProjection_t{ V3, 1.00 },
3213  geo::PlaneGeo::WireCoordProjection_t{ 2.40, 1.00 }
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  //
3245  geo::PlaneGeo::WireCoordProjection_t const testProj
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()
3316 
3317 
3318  //......................................................................
3319  void GeometryTestAlg::testPlanePitch()
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()
3396 
3397  //......................................................................
3398 
3399  void GeometryTestAlg::testStepping()
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  }
3449 
3450  //......................................................................
3451 
3452  void GeometryTestAlg::testProject()
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  }
3486 
3487 
3488  //......................................................................
3489  bool GeometryTestAlg::CheckAuxDetAtPosition
3490  (double const pos[3], unsigned int expected) const
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()
3512 
3513  //......................................................................
3514  bool GeometryTestAlg::CheckAuxDetSensitiveAtPosition
3515  (double const pos[3], unsigned int expectedDet, unsigned int expectedSens)
3516  const
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()
3541 
3542  //......................................................................
3543  void GeometryTestAlg::testFindAuxDet() const {
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()
3594 
3595  //......................................................................
3596 
3597  inline bool GeometryTestAlg::shouldRunTests(std::string test_name) const {
3598  return fRunTests(test_name);
3599  } // GeometryTestAlg::shouldRunTests()
3600 
3601 }//end namespace
geo::TPCGeo::ID
geo::TPCID const & ID() const
Returns the identifier of this TPC.
Definition: TPCGeo.h:333
geo::GeometryCore::WireCoordinate
geo::Length_t WireCoordinate(double YPos, double ZPos, geo::PlaneID const &planeid) const
Returns the index of the nearest wire to the specified position.
Definition: GeometryCore.cxx:1128
ValidateOpDetReco.end
end
while True: pbar.update(maxval-len(onlies[E][S])) #print iS, "/", len(onlies[E][S]) found = False for...
Definition: ValidateOpDetReco.py:548
geo::GeometryCore::FindAuxDetSensitiveAtPosition
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.
Definition: GeometryCore.cxx:525
geo::WireGeo::GetStart
void GetStart(double *xyz) const
Definition: WireGeo.h:157
geo::WireGeo
Geometry description of a TPC wireThe wire is a single straight segment on a wire plane...
Definition: WireGeo.h:65
geo::GeometryTestAlg::testPlaneProjection
void testPlaneProjection() const
Definition: GeometryTestAlg.cxx:1986
mf::LogVerbatim
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
Definition: MessageLogger.h:216
QString::operator<<
QDataStream & operator<<(QDataStream &s, const QString &str)
Definition: qstring.cpp:15066
geo_vectors_utils.h
Utilities to extend the interface of geometry vectors.
geo::CryostatGeo::GetCenter
geo::Point_t GetCenter() const
Returns the geometrical center of the cryostat.
Definition: CryostatGeo.h:124
geo::GeometryTestAlg::GeometryTestAlg
GeometryTestAlg(fhicl::ParameterSet const &pset)
Definition: GeometryTestAlg.cxx:118
geo::GeometryCore::CryostatHalfHeight
geo::Length_t CryostatHalfHeight(geo::CryostatID const &cid) const
Returns the height of the cryostat (y direction)
Definition: GeometryCore.cxx:778
geo::TPCGeo::DriftDir
Vector DriftDir() const
Returns the direction of the drift (vector pointing toward the planes).
Definition: TPCGeo.h:773
geo::WireIDIntersection::z
double z
z position of intersection
Definition: geo_types.h:805
geo::GeometryTestAlg::testPlaneDirections
void testPlaneDirections() const
Definition: GeometryTestAlg.cxx:941
geo::TPCGeo::PlanePitch
double PlanePitch(unsigned int p1=0, unsigned int p2=1) const
Definition: TPCGeo.cxx:388
geo::GeometryTestAlg::fNonFatalExceptions
std::set< std::string > fNonFatalExceptions
Definition: GeometryTestAlg.h:130
geo::TPCGeo::PrintTPCInfo
void PrintTPCInfo(Stream &&out, std::string indent="", unsigned int verbosity=1) const
Prints information about this TPC.
Definition: TPCGeo.h:789
lar::dump::vector3D
auto vector3D(Vector3D const &v)
Returns a manipulator which will print the specified vector.
Definition: DumpUtils.h:301
geo::TPCGeo::GetActiveVolumeCenter
Point GetActiveVolumeCenter() const
Returns the center of the TPC active volume in world coordinates [cm].
Definition: TPCGeo.h:288
geo::vect::dot
constexpr auto dot(Vector const &a, Vector const &b)
Return cross product of two vectors.
Definition: geo_vectors_utils.h:989
keras_to_tensorflow.f
f
Definition: keras_to_tensorflow.py:162
geo::GeometryCore::Plane
PlaneGeo const & Plane(unsigned int const p, unsigned int const tpc=0, unsigned int const cstat=0) const
Returns the specified wire.
Definition: GeometryCore.h:2994
make_norm_csv.output
list output
Definition: make_norm_csv.py:48
geo::GeometryCore::DetHalfWidth
geo::Length_t DetHalfWidth(geo::TPCID const &tpcid) const
Returns the half width of the active volume of the specified TPC.
Definition: GeometryCore.cxx:752
geo::GeometryCore::TotalNTPC
unsigned int TotalNTPC() const
Returns the total number of TPCs in the detector.
Definition: GeometryCore.cxx:868
geo::kMysteryType
Who knows?
Definition: geo_types.h:147
lar::util::RealComparisons::nonEqual
constexpr bool nonEqual(Value_t a, Value_t b) const
Returns whether a and b are farther than the threshold.
Definition: RealComparisons.h:85
CryostatGeo.h
Encapsulate the construction of a single cyostat.
geo::GeometryCore::FindAuxDetAtPosition
unsigned int FindAuxDetAtPosition(double const worldLoc[3], double tolerance=0) const
Returns the index of the auxiliary detector at specified location.
Definition: GeometryCore.cxx:505
geo::GeometryCore::PositionToCryostat
CryostatGeo const & PositionToCryostat(geo::Point_t const &point) const
Returns the cryostat at specified location.
Definition: GeometryCore.cxx:460
geo::GeometryCore::IteratePlanes
IteratorBox< plane_iterator,&GeometryCore::begin_plane,&GeometryCore::end_plane > IteratePlanes() const
Enables ranged-for loops on all planes of the detector.
Definition: GeometryCore.h:3209
geo::TPCGeo::LengthDir
decltype(auto) LengthDir() const
Returns the direction Length() is measured on.
Definition: TPCGeo.h:138
geo.h
Collect all the geometry header files together.
geo::GeometryTestAlg::fExpectedWirePitches
std::vector< double > fExpectedWirePitches
wire pitch on each plane
Definition: GeometryTestAlg.h:131
readout_types.h
Classes identifying readout-related concepts.
geo::AuxDetSensitiveGeo::GetNormalVector
geo::Vector_t GetNormalVector() const
Returns the unit normal vector to the detector.
Definition: AuxDetSensitiveGeo.cxx:57
expected
const char expected[]
Definition: Exception_t.cc:22
mf::MaybeLogger_
Definition: MessageLogger.h:113
geo::GeometryTestAlg::testPlaneProjectionReference
void testPlaneProjectionReference() const
Definition: GeometryTestAlg.cxx:1596
geo::PlaneGeo::Wire
WireGeo const & Wire(unsigned int iwire) const
Definition: PlaneGeo.cxx:506
geo::View_t
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
string
std::string string
Definition: nybbler.cc:12
geo::TPCGeo::ActiveHalfHeight
double ActiveHalfHeight() const
Half height (associated with y coordinate) of active TPC volume [cm].
Definition: TPCGeo.h:99
geo::details::TPC_id_iterator_base
Base forward iterator browsing all TPC IDs in the detector.
Definition: GeometryCore.h:292
geo::GeometryTestAlg::testPlaneProjectionOnFrame
void testPlaneProjectionOnFrame() const
Definition: GeometryTestAlg.cxx:1855
mf::LogInfo
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
Definition: MessageLogger.h:209
geo::AuxDetGeo::SensitiveVolume
AuxDetSensitiveGeo const & SensitiveVolume(size_t sv) const
Definition: AuxDetGeo.h:171
geo::WireGeo::DistanceFrom
double DistanceFrom(geo::WireGeo const &wire) const
Returns 3D distance from the specified wire.
Definition: WireGeo.cxx:113
geo::kUnknown
Unknown view.
Definition: geo_types.h:136
geo::TPCGeo::Nplanes
unsigned int Nplanes() const
Number of planes in this tpc.
Definition: TPCGeo.h:165
geo::GeometryCore::WorldBox
void WorldBox(double *xlo, double *xhi, double *ylo, double *yhi, double *zlo, double *zhi) const
Fills the arguments with the boundaries of the world.
Definition: GeometryCore.cxx:926
lar::util::RealComparisons
Provides simple real number checks.
Definition: RealComparisons.h:66
testing::NameSelector::AddToDefinition
void AddToDefinition(std::string set_name, NAMES...names)
Parses a list of names and add them to the specified definition.
Definition: NameSelector.h:127
geo::InvalidWireIDError::wire_number
int wire_number
the invalid wire number
Definition: Exceptions.h:167
geo::BoxBoundedGeo::MinX
double MinX() const
Returns the world x coordinate of the start of the box.
Definition: BoxBoundedGeo.h:88
geo::PlaneID::PlaneID_t
unsigned int PlaneID_t
Type for the ID number.
Definition: geo_types.h:473
geo::PlaneID
The data type to uniquely identify a Plane.
Definition: geo_types.h:472
geo::CryostatID::isValid
bool isValid
Whether this ID points to a valid element.
Definition: geo_types.h:211
geo::TPCGeo
Geometry information for a single TPC.
Definition: TPCGeo.h:38
geo::GeometryTestAlg::testWireCoordAngle
void testWireCoordAngle() const
Definition: GeometryTestAlg.cxx:1379
vector
struct vector vector
bump_copyright.d
string d
Definition: bump_copyright.py:70
geo::BoxBoundedGeo::CenterX
double CenterX() const
Returns the world x coordinate of the center of the box.
Definition: BoxBoundedGeo.h:94
geo::vect::toPoint
::geo::Point_t toPoint(Point const &p)
Convert the specified point into a geo::Point_t.
Definition: geo_vectors_utils.h:1280
geo::GeometryCore::GetWorldVolumeName
const std::string GetWorldVolumeName() const
Return the name of the world volume (needed by Geant4 simulation)
Definition: GeometryCore.cxx:619
MF_LOG_ERROR
#define MF_LOG_ERROR(category)
Definition: MessageLogger.h:229
testing::NameSelector::RejectedNames
Names_t RejectedNames() const
Returns a list of names that were rejected.
Definition: NameSelector.h:153
geo::BoxBoundedGeo::CenterZ
double CenterZ() const
Returns the world z coordinate of the center of the box.
Definition: BoxBoundedGeo.h:124
geo::GeometryCore::ChannelToWire
std::vector< geo::WireID > ChannelToWire(raw::ChannelID_t const channel) const
Returns a list of wires connected to the specified TPC channel.
Definition: GeometryCore.cxx:1107
geo::CryostatID::Cryostat
CryostatID_t Cryostat
Index of cryostat.
Definition: geo_types.h:212
std
STL namespace.
geo::GeometryTestAlg::printDetectorIntro
void printDetectorIntro() const
Definition: GeometryTestAlg.cxx:386
geo::GeometryTestAlg::ExpectedPlane_dTdW
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.
Definition: GeometryTestAlg.cxx:2961
geo::GeometryCore::MaxWires
unsigned int MaxWires() const
Returns the largest number of wires among all planes in this detector.
Definition: GeometryCore.cxx:888
test_nxdot.dot
dot
Definition: test_nxdot.py:12
geo::BoxBoundedGeo::MaxX
double MaxX() const
Returns the world x coordinate of the end of the box.
Definition: BoxBoundedGeo.h:91
geo::PlaneGeo::IterateWires
ElementIteratorBox IterateWires() const
Definition: PlaneGeo.h:401
geo::WireID::Wire
WireID_t Wire
Index of the wire within its plane.
Definition: geo_types.h:580
geo::TPCGeo::ActiveMass
double ActiveMass() const
Definition: TPCGeo.h:118
geo::GeometryTestAlg::testWireOrientations
void testWireOrientations() const
Definition: GeometryTestAlg.cxx:1013
testing::NameSelector::ParseNames
void ParseNames(NAMES...names)
Parses a list of names and adds them to the selector.
Definition: NameSelector.h:97
geo::Vector_t
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
testing::NameSelector::PrintConfiguration
void PrintConfiguration(std::ostream &) const
Prints the configuration into a stream.
Definition: NameSelector.cxx:49
channel
uint8_t channel
Definition: CRTFragment.hh:201
dir
string dir
Definition: INukeNucleonCorr.cxx:61
geo::GeometryCore::SignalType
SigType_t SignalType(geo::PlaneID const &pid) const
Returns the type of signal on the channels of specified TPC plane.
Definition: GeometryCore.cxx:578
mf::LogError
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
Definition: MessageLogger.h:211
geo::kNegX
Drift towards negative X values.
Definition: geo_types.h:162
geo::PlaneGeo::ActiveArea
Rect const & ActiveArea() const
Returns an area covered by the wires in the plane.
Definition: PlaneGeo.h:758
geo::CryostatGeo
Geometry information for a single cryostat.
Definition: CryostatGeo.h:43
geo::PlaneGeo::GetNormalDirection
Vector GetNormalDirection() const
Returns the direction normal to the plane.
Definition: PlaneGeo.h:442
geo::GeometryTestAlg::Run
virtual unsigned int Run()
Runs the test, returns a number of errors (very unlikely!)
Definition: GeometryTestAlg.cxx:151
geo::GeometryCore::Ncryostats
unsigned int Ncryostats() const
Returns the number of cryostats in the detector.
Definition: GeometryCore.h:1994
geo::WireGeo::WirePitch
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
cet::square
constexpr T square(T x)
Definition: pow.h:21
geo::GeometryCore::ThirdPlane
geo::PlaneID ThirdPlane(geo::PlaneID const &pid1, geo::PlaneID const &pid2) const
Returns the plane that is not in the specified arguments.
Definition: GeometryCore.cxx:1482
geo::PlaneGeo::GetIncreasingWireDirection
Vector GetIncreasingWireDirection() const
Returns the direction of increasing wires.
Definition: PlaneGeo.h:457
geo::WireGeo::ThetaZ
double ThetaZ() const
Returns angle of wire with respect to z axis in the Y-Z plane in radians.
Definition: WireGeo.h:250
geo::GeometryCore::WirePitch
geo::Length_t WirePitch(geo::PlaneID const &planeid) const
Returns the distance between two consecutive wires.
Definition: GeometryCore.cxx:824
geo::GeometryTestAlg::testThirdPlane_dTdW_at
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.
Definition: GeometryTestAlg.cxx:3028
RealComparisons.h
Class for approximate comparisons.
geo::GeometryCore::IterateTPCIDs
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.
Definition: GeometryCore.h:2597
geo::TPCGeo::Width
double Width() const
Width is associated with x coordinate [cm].
Definition: TPCGeo.h:109
geo::GeometryTestAlg::shouldRunTests
bool shouldRunTests(std::string test_name) const
Definition: GeometryTestAlg.cxx:3597
geo::TPCGeo::Height
double Height() const
Height is associated with y coordinate [cm].
Definition: TPCGeo.h:113
geo::GeometryTestAlg::testWireIntersection
void testWireIntersection() const
Definition: GeometryTestAlg.cxx:2371
geo::PlaneGeo::InterWireProjectedDistance
double InterWireProjectedDistance(WireCoordProjection_t const &projDir) const
Returns the distance between wires along the specified direction.
Definition: PlaneGeo.cxx:705
geo::GeometryCore::ROOTFile
std::string ROOTFile() const
Returns the full directory path to the geometry file source.
Definition: GeometryCore.h:1740
geo::PlaneGeo::DriftPoint
void DriftPoint(geo::Point_t &position, double distance) const
Shifts the position of an electron drifted by a distance.
Definition: PlaneGeo.h:643
geo::GeometryTestAlg::isWireAlignedToPlaneDirections
bool isWireAlignedToPlaneDirections(geo::PlaneGeo const &plane, geo::Vector_t const &wireDir) const
Helper function for testWireIntersection().
Definition: GeometryTestAlg.cxx:2348
geo::GeometryCore::Nchannels
unsigned int Nchannels() const
Returns the number of TPC readout channels in the detector.
Definition: GeometryCore.cxx:208
geo::TPCGeo::Length
double Length() const
Length is associated with z coordinate [cm].
Definition: TPCGeo.h:115
geo::GeometryCore::TotalMass
double TotalMass() const
Returns the total mass [kg] of the specified volume (default: world).
Definition: GeometryCore.h:1891
geo::WiresIntersection
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
geo::GeometryCore::coordIs
static lar::util::RealComparisons< geo::Length_t > coordIs
Value of tolerance for equality comparisons.
Definition: GeometryCore.h:1513
geo::GeometryCore::IteratePlaneIDs
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.
Definition: GeometryCore.h:3127
geo::GeometryCore::FindTPCAtPosition
geo::TPCID FindTPCAtPosition(double const worldLoc[3]) const
Returns the ID of the TPC at specified location.
Definition: GeometryCore.h:2454
geo::WireGeo::LocalToWorldVect
void LocalToWorldVect(const double *wire, double *world) const
Transform direction vector from local to world.
Definition: WireGeo.h:363
geo::AuxDetGeo::HalfWidth2
double HalfWidth2() const
Definition: AuxDetGeo.h:104
geo::WireGeo::GetPositionFromCenter
Point GetPositionFromCenter(double localz) const
Returns the position (world coordinate) of a point on the wire.
Definition: WireGeo.h:182
mf::LogProblem
MaybeLogger_< ELseverityLevel::ELsev_error, true > LogProblem
Definition: MessageLogger.h:218
geo::GeometryTestAlg::fExpectedPlanePitches
std::vector< double > fExpectedPlanePitches
plane pitch on each plane
Definition: GeometryTestAlg.h:132
abs
T abs(T value)
Definition: sparse_vector_test.cc:30
testing::NameSelector::AcceptedNames
Names_t AcceptedNames() const
Returns a list of names that were accepted.
Definition: NameSelector.h:150
geo::GeometryTestAlg::testChannelToROP
void testChannelToROP() const
Definition: GeometryTestAlg.cxx:1456
geo::PlaneGeo::View
View_t View() const
Which coordinate does this plane measure.
Definition: PlaneGeo.h:184
geo::GeometryCore::IterateWireIDs
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.
Definition: GeometryCore.h:3554
geo::PlaneGeo::WidthDepthProjection_t
ROOT::Math::DisplacementVector2D< ROOT::Math::Cartesian2D< double >, WidthDepthReferenceTag > WidthDepthProjection_t
Definition: PlaneGeo.h:151
geo::WireGeo::LocalToWorld
void LocalToWorld(const double *wire, double *world) const
Transform point from local wire frame to world frame.
Definition: WireGeo.h:355
geo::GeometryCore::Nplanes
unsigned int Nplanes(unsigned int tpc=0, unsigned int cstat=0) const
Returns the total number of wire planes in the specified TPC.
Definition: GeometryCore.h:2848
geo::kHorizontal
Planes that are in the horizontal plane.
Definition: geo_types.h:140
raw::InvalidChannelID
constexpr ChannelID_t InvalidChannelID
ID of an invalid channel.
Definition: RawTypes.h:32
geo::GeometryCore::MaxPlanes
unsigned int MaxPlanes() const
Returns the largest number of planes among all TPCs in this detector.
Definition: GeometryCore.cxx:878
geo::WireGeo::IntersectionWith
Point IntersectionWith(geo::WireGeo const &other) const
Returns the point of this wire that is closest to other wire.
Definition: WireGeo.h:637
geo::InvalidWireIDError::better_wire_number
int better_wire_number
a suggestion for a good wire number
Definition: Exceptions.h:168
e
const double e
Definition: gUpMuFluxGen.cxx:165
geo::GeometryCore::DetHalfHeight
geo::Length_t DetHalfHeight(geo::TPCID const &tpcid) const
Returns the half height of the active volume of the specified TPC.
Definition: GeometryCore.cxx:758
lar::util::makeVector3DComparison
auto makeVector3DComparison(RealType threshold)
Creates a Vector3DComparison from a RealComparisons object.
Definition: RealComparisons.h:252
AuxDetSensitiveGeo.h
Encapsulate the geometry of the sensitive portion of an auxiliary detector.
geo::GeometryTestAlg::fRunTests
testing::NameSelector fRunTests
test filter
Definition: GeometryTestAlg.h:137
nodes
#define nodes
Decomposer.h
Classes to project and compose a vector on a plane.
geo::GeometryTestAlg::testThirdPlane_dTdW
void testThirdPlane_dTdW() const
Definition: GeometryTestAlg.cxx:2830
geo::GeometryCore::DetectorName
std::string DetectorName() const
Returns a string with the name of the detector, as configured.
Definition: GeometryCore.h:1761
geo::DriftDirection_t
enum geo::driftdir DriftDirection_t
Drift direction: positive or negative.
geo::GeometryTestAlg::testInterWireProjectedDistance
void testInterWireProjectedDistance() const
Definition: GeometryTestAlg.cxx:3183
geo::AuxDetGeo::HalfHeight
double HalfHeight() const
Definition: AuxDetGeo.h:105
geo::WireGeo::RMax
double RMax() const
Returns the outer half-size of the wire [cm].
Definition: WireGeo.cxx:91
Exceptions.h
Collection of exceptions for Geometry system.
testing::NameSelector::Parse
void Parse(LIST const &items)
Parses the names in the list and adds them to the selector.
Definition: NameSelector.h:77
geo::GeometryCore::FirstChannelInROP
raw::ChannelID_t FirstChannelInROP(readout::ROPID const &ropid) const
Returns the ID of the first channel in the specified readout plane.
Definition: GeometryCore.cxx:1758
geo::PlaneGeo::WireCoordProjection_t
ROOT::Math::DisplacementVector2D< ROOT::Math::Cartesian2D< double >, WireCoordinateReferenceTag > WireCoordProjection_t
Type for projections in the wire base representation.
Definition: PlaneGeo.h:130
geo::SigType_t
enum geo::_plane_sigtype SigType_t
geo::GeometryCore::IterateTPCs
IteratorBox< TPC_iterator,&GeometryCore::begin_TPC,&GeometryCore::end_TPC > IterateTPCs() const
Enables ranged-for loops on all TPCs of the detector.
Definition: GeometryCore.h:2655
reco_momentum_tuples.t
t
Definition: reco_momentum_tuples.py:25
geo::TPCGeo::ActiveHalfWidth
double ActiveHalfWidth() const
Half width (associated with x coordinate) of active TPC volume [cm].
Definition: TPCGeo.h:95
geo::GeometryTestAlg::printAuxiliaryDetectors
void printAuxiliaryDetectors() const
Method to print the auxiliary detectors on screen.
Definition: GeometryTestAlg.cxx:589
geo::PlaneGeo::GetCenter
Point GetCenter() const
Returns the centre of the wire plane in world coordinates [cm].
Definition: PlaneGeo.h:479
geo::GeometryCore::NearestWire
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.
Definition: GeometryCore.cxx:1141
make_pion_angle.shifts
list shifts
Definition: make_pion_angle.py:17
geo::CryostatGeo::HalfWidth
double HalfWidth() const
Half width of the cryostat [cm].
Definition: CryostatGeo.cxx:190
testing::NameSelector::CheckQueryRegistry
bool CheckQueryRegistry() const
Checks that no known element with valid response was left unqueried.
Definition: NameSelector.h:160
geo::PlaneGeo
Geometry information for a single wire plane.The plane is represented in the geometry by a solid whic...
Definition: PlaneGeo.h:82
geo::PlaneGeo::InterWireDistance
double InterWireDistance(geo::Vector_t const &dir) const
Returns the distance between wires along the specified direction.
Definition: PlaneGeo.cxx:713
MF_LOG_TRACE
#define MF_LOG_TRACE(id)
Definition: MessageLogger.h:270
geo::BoxBoundedGeo::MinZ
double MinZ() const
Returns the world z coordinate of the start of the box.
Definition: BoxBoundedGeo.h:118
geo::GeometryTestAlg::fDisableValidWireIDcheck
bool fDisableValidWireIDcheck
disable test on out-of-world NearestWire()
Definition: GeometryTestAlg.h:129
test.p
p
Definition: test.py:223
raw::isValidChannelID
constexpr bool isValidChannelID(raw::ChannelID_t channel)
Definition: RawTypes.h:37
geo::CryostatGeo::NTPC
unsigned int NTPC() const
Number of TPCs in this cryostat.
Definition: CryostatGeo.h:181
geo::Xaxis
constexpr Vector Xaxis()
Returns a x axis vector of the specified type.
Definition: geo_vectors.h:215
AuxDetGeo.h
Encapsulate the geometry of an auxiliary detector.
geo::GeometryCore::DetLength
geo::Length_t DetLength(geo::TPCID const &tpcid) const
Returns the length of the active volume of the specified TPC.
Definition: GeometryCore.cxx:764
geo::GeometryCore::Cryostat
CryostatGeo const & Cryostat(geo::CryostatID const &cryoid) const
Returns the specified cryostat.
Definition: GeometryCore.cxx:321
geo::AuxDetGeo::NSensitiveVolume
size_t NSensitiveVolume() const
Definition: AuxDetGeo.h:172
mf::isDebugEnabled
bool isDebugEnabled()
Definition: MessageLogger.cc:611
geo::GeometryTestAlg::CheckAuxDetAtPosition
bool CheckAuxDetAtPosition(double const pos[3], unsigned int expected) const
Returns whether the auxiliary detector at pos is the expected one.
Definition: GeometryTestAlg.cxx:3490
geo::GeometryTestAlg::testChannelToWire
void testChannelToWire() const
Definition: GeometryTestAlg.cxx:1495
lar::util::simple_geo::Rectangle::width
Range_t width
Range along width direction.
Definition: SimpleGeo.h:393
protoana::distance
double distance(double x1, double y1, double z1, double x2, double y2, double z2)
Definition: truepionXsection_module.cc:217
geo::AuxDetGeo::GetNormalVector
geo::Vector_t GetNormalVector() const
Returns the unit normal vector to the detector.
Definition: AuxDetGeo.cxx:72
max
static int max(int a, int b)
Definition: fortranscanner.cpp:60366
geo::TPCID
The data type to uniquely identify a TPC.
Definition: geo_types.h:386
geo::PlaneID::Plane
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:493
test_nxjs.indent
indent
Definition: test_nxjs.py:18
MF_LOG_INFO
#define MF_LOG_INFO(category)
Definition: MessageLogger.h:225
geo::GeometryCore::View
View_t View(geo::PlaneID const &pid) const
Returns the view (wire orientation) on the channels of specified TPC plane.
Definition: GeometryCore.cxx:601
geo_types.h
Definition of data types for geometry description.
geo::TPCGeo::Plane0Pitch
double Plane0Pitch(unsigned int p) const
Definition: TPCGeo.cxx:324
geo::GeometryTestAlg::testWireCoordFromPlane
void testWireCoordFromPlane() const
Definition: GeometryTestAlg.cxx:1065
geo::TPCGeo::DriftDirection
DriftDirection_t DriftDirection() const
Returns an enumerator value describing the drift direction.
Definition: TPCGeo.h:144
geo::GeometryTestAlg::printWiresInTPC
void printWiresInTPC(const TPCGeo &tpc, std::string indent="") const
Definition: GeometryTestAlg.cxx:525
readout::ROPID
Class identifying a set of planes sharing readout channels.
Definition: readout_types.h:168
geo::GeometryCore::AuxDet
AuxDetGeo const & AuxDet(unsigned int const ad=0) const
Returns the specified auxiliary detector.
Definition: GeometryCore.cxx:340
geo::TPCGeo::ActiveLength
double ActiveLength() const
Length (associated with z coordinate) of active TPC volume [cm].
Definition: TPCGeo.h:103
geo::WireGeo::Direction
Vector Direction() const
Definition: WireGeo.h:587
geo::BoxBoundedGeo::MaxY
double MaxY() const
Returns the world y coordinate of the end of the box.
Definition: BoxBoundedGeo.h:106
geo::GeometryTestAlg::testTPC
void testTPC(geo::CryostatID const &cid)
Definition: GeometryTestAlg.cxx:855
keras_to_tensorflow.int
int
Definition: keras_to_tensorflow.py:69
geo::PlaneGeo::Orientation
Orient_t Orientation() const
What is the orientation of the plane.
Definition: PlaneGeo.h:187
geo::TPCID::asTPCID
constexpr TPCID const & asTPCID() const
Conversion to TPCID (for convenience of notation).
Definition: geo_types.h:446
WireGeo.h
Encapsulate the geometry of a wire.
geo::GeometryCore::FindAllVolumes
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.
Definition: GeometryCore.cxx:706
geo::GeometryCore::NTPC
unsigned int NTPC(unsigned int cstat=0) const
Returns the total number of TPCs in the specified cryostat.
Definition: GeometryCore.h:2303
geo::CryostatGeo::HalfHeight
double HalfHeight() const
Half height of the cryostat [cm].
Definition: CryostatGeo.cxx:196
geo::GeometryCore::ChannelToROP
readout::ROPID ChannelToROP(raw::ChannelID_t channel) const
Definition: GeometryCore.cxx:1113
geo::TPCGeo::HalfHeight
double HalfHeight() const
Height is associated with y coordinate [cm].
Definition: TPCGeo.h:111
geo::PlaneGeo::DepthDir
Vector DepthDir() const
Return the direction of plane depth.
Definition: PlaneGeo.h:236
geo::WireGeo::HalfL
double HalfL() const
Returns half the length of the wire [cm].
Definition: WireGeo.h:128
geo::Point_t
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
geo::GeometryCore::IterateCryostats
IteratorBox< cryostat_iterator,&GeometryCore::begin_cryostat,&GeometryCore::end_cryostat > IterateCryostats() const
Enables ranged-for loops on all cryostats of the detector.
Definition: GeometryCore.h:2198
geo::GeometryCore::WireIDsIntersect
bool WireIDsIntersect(WireID const &wid1, WireID const &wid2, geo::Point_t &intersection) const
Computes the intersection between two wires.
geo::CryostatGeo::FindTPCAtPosition
geo::TPCID::TPCID_t FindTPCAtPosition(double const worldLoc[3], double const wiggle) const
Returns the index of the TPC at specified location.
Definition: CryostatGeo.cxx:127
geo::GeometryCore::Wire
WireGeo const & Wire(geo::WireID const &wireid) const
Returns the specified wire.
Definition: GeometryCore.h:3419
geo::GeometryTestAlg::testFindCryostatVolumes
unsigned int testFindCryostatVolumes()
Definition: GeometryTestAlg.cxx:761
geo::CryostatGeo::LocalToWorld
void LocalToWorld(const double *cryo, double *world) const
Transform point from local cryostat frame to world frame.
Definition: CryostatGeo.h:387
geo::GeometryCore::PlaneWireToChannel
raw::ChannelID_t PlaneWireToChannel(WireID const &wireid) const
Returns the ID of the TPC channel connected to the specified wire.
Definition: GeometryCore.cxx:1226
geo::GeometryTestAlg::testPlanePointDecompositionFrame
void testPlanePointDecompositionFrame() const
Definition: GeometryTestAlg.cxx:1635
geo::TPCGeo::DriftDistance
double DriftDistance() const
Definition: TPCGeo.h:155
geo::TPCGeo::HeightDir
decltype(auto) HeightDir() const
Returns the direction Height() is measured on.
Definition: TPCGeo.h:131
geo::GeometryTestAlg::testFindTPCvolumePaths
unsigned int testFindTPCvolumePaths()
Definition: GeometryTestAlg.cxx:789
geo::PlaneGeo::VectorProjection
WireCoordProjection_t VectorProjection(geo::Vector_t const &v) const
Definition: PlaneGeo.h:963
cet::sqlite::min
T min(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:55
geo::kPosX
Drift towards positive X values.
Definition: geo_types.h:161
PlaneGeo.h
Encapsulate the construction of a single detector plane.
geo::AuxDetGeo::Length
double Length() const
Definition: AuxDetGeo.h:102
geo::AuxDetGeo::Name
std::string Name() const
Definition: AuxDetGeo.h:154
geo::PlaneGeo::PrintPlaneInfo
void PrintPlaneInfo(Stream &&out, std::string indent="", unsigned int verbosity=1) const
Prints information about this plane.
Definition: PlaneGeo.h:1539
larg4::dist
constexpr double dist(const TReal *x, const TReal *y, const unsigned int dimension)
Definition: OpFastScintillation.hh:563
geo::BoxBoundedGeo::CenterY
double CenterY() const
Returns the world y coordinate of the center of the box.
Definition: BoxBoundedGeo.h:109
geo::GeometryCore::CryostatHalfWidth
geo::Length_t CryostatHalfWidth(geo::CryostatID const &cid) const
Returns the half width of the cryostat (x direction)
Definition: GeometryCore.cxx:771
geo::CryostatGeo::TPC
const TPCGeo & TPC(unsigned int itpc) const
Return the itpc&#39;th TPC in the cryostat.
Definition: CryostatGeo.cxx:93
geo::WireGeo::GetEnd
void GetEnd(double *xyz) const
Definition: WireGeo.h:163
geo::GeometryTestAlg::printAllGeometry
void printAllGeometry() const
Definition: GeometryTestAlg.cxx:567
geo::WireIDIntersection::y
double y
y position of intersection
Definition: geo_types.h:804
geo::GeometryTestAlg::testWireIntersectionAt
unsigned int testWireIntersectionAt(geo::TPCGeo const &TPC, TVector3 const &point) const
Performs the wire intersection test at a single point.
Definition: GeometryTestAlg.cxx:2536
geo::BoxBoundedGeo::MaxZ
double MaxZ() const
Returns the world z coordinate of the end of the box.
Definition: BoxBoundedGeo.h:121
geo::PlaneGeo::ID
geo::PlaneID const & ID() const
Returns the identifier of this plane.
Definition: PlaneGeo.h:203
geo::GeometryTestAlg::testParallelWires
void testParallelWires() const
Definition: GeometryTestAlg.cxx:1139
lar::util::RealComparisons::nonZero
constexpr bool nonZero(Value_t value) const
Returns whether the value is farther from 0 than the threshold.
Definition: RealComparisons.h:77
geo::PlaneGeo::Nwires
unsigned int Nwires() const
Number of wires in this plane.
Definition: PlaneGeo.h:269
mf::LogWarning
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
Definition: MessageLogger.h:210
wirecell.gen.depos.center
def center(depos, point)
Definition: depos.py:117
geo::PlaneGeo::WireIDincreasesWithZ
bool WireIDincreasesWithZ() const
Returns whether the higher z wires have higher wire ID.
Definition: PlaneGeo.cxx:525
MF_LOG_DEBUG
#define MF_LOG_DEBUG(id)
Definition: MessageLogger.h:268
SimpleGeo.h
Some simple functions to represent geometry entities.
geo::GeometryTestAlg::testFindAuxDet
void testFindAuxDet() const
Definition: GeometryTestAlg.cxx:3543
geo::PlaneGeo::WidthDir
Vector WidthDir() const
Return the direction of plane width.
Definition: PlaneGeo.h:221
copy
T copy(T const &v)
Definition: RangeForWrapper_test.cc:98
geo::GeometryTestAlg::testThirdPlane
void testThirdPlane() const
Definition: GeometryTestAlg.cxx:2729
geo::GeometryCore::CryostatLength
geo::Length_t CryostatLength(geo::CryostatID const &cid) const
Returns the length of the cryostat (z direction)
Definition: GeometryCore.cxx:784
reco_momentum_tuples.count
int count
Definition: reco_momentum_tuples.py:202
geo::WireID::WireID_t
unsigned int WireID_t
Type for the ID number.
Definition: geo_types.h:561
geo_vectors.h
Definitions of geometry vector data types.
GeometryCore.h
Access the description of detector geometry.
geo::GeometryCore::FindAllVolumePaths
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.
Definition: GeometryCore.cxx:723
geo::GeometryCore::TPC
TPCGeo const & TPC(unsigned int const tpc=0, unsigned int const cstat=0) const
Returns the specified TPC.
Definition: GeometryCore.h:2423
util::begin
decltype(auto) constexpr begin(T &&obj)
ADL-aware version of std::begin.
Definition: StdUtils.h:72
geo::GeometryTestAlg::testPlanePointDecomposition
void testPlanePointDecomposition() const
Definition: GeometryTestAlg.cxx:1172
geo::InvalidWireIDError
Exception thrown on invalid wire number (e.g. NearestWireID())
Definition: Exceptions.h:158
geo::ProjectToBoxEdge
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::TPCGeo::Plane
PlaneGeo const & Plane(geo::View_t view) const
Return the plane in the tpc with View_t view.
Definition: TPCGeo.cxx:263
geo::AuxDetGeo::HalfWidth1
double HalfWidth1() const
Definition: AuxDetGeo.h:103
geo::WireGeo::GetCenter
void GetCenter(double *xyz, double localz=0.0) const
Fills the world coordinate of a point on the wire.
Definition: WireGeo.cxx:73
geo::GeometryCore::Nviews
unsigned int Nviews() const
Returns the number of views (different wire orientations)
Definition: GeometryCore.cxx:307
raw::ChannelID_t
unsigned int ChannelID_t
Type representing the ID of a readout channel.
Definition: RawTypes.h:28
cs
const char * cs
Definition: qtextcodec.cpp:1350
art::get
auto const & get(AssnsNode< L, R, D > const &r)
Definition: AssnsNode.h:115
MF_LOG_WARNING
#define MF_LOG_WARNING(category)
Definition: MessageLogger.h:227
geo::TPCID::TPC
TPCID_t TPC
Index of the TPC within its cryostat.
Definition: geo_types.h:406
PhysicalConstants.h
Collection of Physical constants used in LArSoft.
geo::GeometryCore::NearestChannel
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.
Definition: GeometryCore.cxx:1207
geo::GeometryTestAlg::testFindWorldVolumes
unsigned int testFindWorldVolumes()
Definition: GeometryTestAlg.cxx:730
geo::GeometryTestAlg::printAuxDetSensitiveGeo
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.
Definition: GeometryTestAlg.cxx:651
geo::TPCGeo::MaxVerbosity
static constexpr unsigned int MaxVerbosity
Maximum verbosity supported by PrintTPCInfo().
Definition: TPCGeo.h:677
geo
LArSoft geometry interface.
Definition: ChannelGeo.h:16
geo::BoxBoundedGeo::MinY
double MinY() const
Returns the world y coordinate of the start of the box.
Definition: BoxBoundedGeo.h:103
geo::GeometryCore::ThirdPlane_dTdW
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.
Definition: GeometryCore.cxx:1559
geo::GeometryTestAlg::printAuxDetGeo
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.
Definition: GeometryTestAlg.cxx:605
bool
int bool
Definition: qglobal.h:345
geo::PlaneGeo::ComposeVector
Vector ComposeVector(WireDecomposedVector_t const &decomp) const
Returns the 3D vector from composition of projection and distance.
Definition: PlaneGeo.h:982
geo::GeometryCore::MaxTPCs
unsigned int MaxTPCs() const
Returns the largest number of TPCs a cryostat in the detector has.
Definition: GeometryCore.cxx:858
geo::TPCGeo::LocalToWorld
void LocalToWorld(const double *tpc, double *world) const
Transform point from local TPC frame to world frame.
Definition: TPCGeo.h:563
geo::WireGeo::RMin
double RMin() const
Returns the inner radius of the wire (usually 0) [cm].
Definition: WireGeo.cxx:95
geo::AuxDetSensitiveGeo::GetCenter
void GetCenter(double *xyz, double localz=0.0) const
Return the center position of an AuxDet.
Definition: AuxDetSensitiveGeo.cxx:47
geo::PlaneGeo::LocalToWorld
void LocalToWorld(const double *plane, double *world) const
Transform point from local plane frame to world frame.
Definition: PlaneGeo.h:1343
geo::TPCGeo::HalfWidth
double HalfWidth() const
Width is associated with x coordinate [cm].
Definition: TPCGeo.h:107
geo::GeometryTestAlg::geom
geo::GeometryCore const * geom
pointer to geometry service provider
Definition: GeometryTestAlg.h:127
geo::origin
constexpr Point origin()
Returns a origin position with a point of the specified type.
Definition: geo_vectors.h:227
fhicl::exception
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
geo::PlaneGeo::WirePitch
double WirePitch() const
Return the wire pitch (in centimeters). It is assumed constant.
Definition: PlaneGeo.h:411
lar::util::simple_geo::Range::length
Data_t length() const
Returns the distance between upper and lower bounds.
Definition: SimpleGeo.h:341
geo::CryostatGeo::Length
double Length() const
Length of the cryostat [cm].
Definition: CryostatGeo.h:107
endl
QTextStream & endl(QTextStream &s)
Definition: qtextstream.cpp:2030
geo::GeometryCore::NAuxDets
unsigned int NAuxDets() const
Returns the number of auxiliary detectors.
Definition: GeometryCore.h:4549
geo::GeometryTestAlg::GetIncreasingWireDirection
static std::array< double, 3 > GetIncreasingWireDirection(const geo::PlaneGeo &plane)
Returns the direction on plane orthogonal to wires where wire number increases.
Definition: GeometryTestAlg.cxx:2089
TPCGeo.h
Encapsulate the construction of a single detector plane.
geo::CryostatID
The data type to uniquely identify a cryostat.
Definition: geo_types.h:190
geo::GeometryCore::GetLArTPCVolumeName
std::string GetLArTPCVolumeName(geo::TPCID const &tpcid) const
Return the name of specified LAr TPC volume.
Definition: GeometryCore.cxx:740
geo::AuxDetGeo::GetCenter
void GetCenter(double *xyz, double localz=0.0) const
Return the center position of an AuxDet.
Definition: AuxDetGeo.cxx:62
geo::GeometryTestAlg::fComputeMass
bool fComputeMass
Whether to print the detector mass.
Definition: GeometryTestAlg.h:134
geo::TPCGeo::WidthDir
decltype(auto) WidthDir() const
Returns the direction Width() is measured on.
Definition: TPCGeo.h:124
geo::GeometryTestAlg::CheckAuxDetSensitiveAtPosition
bool CheckAuxDetSensitiveAtPosition(double const pos[3], unsigned int expectedDet, unsigned int expectedSens) const
Returns whether the auxiliary sensitive detector at pos is expected.
Definition: GeometryTestAlg.cxx:3515
Generated by   doxygen 1.8.11