GeometryTest_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 /// \file GeometryTest_module.cc
3 /// \brief unit tests for the Geometry service
4 /// \author brebel@fnal.gov
5 ////////////////////////////////////////////////////////////////////////
6 #include <cmath>
7 #include <string>
8 #include <vector>
9 
10 // ROOT includes
11 #include "TH1D.h"
12 #include "TH2D.h"
13 #include "TVector3.h"
14 #include "TNtuple.h"
15 #include "TGeoManager.h"
16 
17 // NOvA includes
18 #include "Geometry/Geometry.h"
19 #include "GeometryObjects/Geo.h"
23 
24 // Framework includes
32 #include "fhiclcpp/ParameterSet.h"
34 
35 namespace geo{
36  class GeometryTest : public art::EDAnalyzer {
37  public:
38  explicit GeometryTest(fhicl::ParameterSet const& pset);
39  virtual ~GeometryTest();
40 
41  void analyze(art::Event const& evt);
42 
43  private:
44 
45  void testCellId();
46  void testProject();
47  int testUniqueId();
48  int testFindCell();
49  void testStepping();
50  void testCellPos();
51  void testFindPlanes();
52  void testCellIdFromPos();
53  };
54 
55  //......................................................................
57  : EDAnalyzer(pset)
58  {
59  }
60 
61  //......................................................................
63  {
64  }
65 
66  //......................................................................
68  {
70 
71  // output some basic information
72  std::cerr << "Looking at detector: " << geom->DetId()
73  << " with base file name: " << geom->FileBaseName()
74  << "\nGDML file: " << geom->GDMLFile()
75  << std::endl;
76 
77  // Generate an exception for testing
78  try { geom->Plane(999999); }
79  catch (...) {
80  std::cerr <<
81  "Exception generated for illegal plane and caught." << std::endl;
82  }
83 
84  double wpos[3] = {0.};
85  for(unsigned int i = 0; i < geom->Plane(0)->Ncells(); ++i){
86  geom->Plane(0)->Cell(i)->GetCenter(wpos,0.);
87  std::cerr << "horizontal cell " << i << " "
88  << wpos[0] << " " << wpos[1] << " " << wpos[2]
89  << " " << geom->DetHalfWidth()
90  << " " << geom->Plane(0)->Cell(i)->HalfW()
91  << " " << geom->Plane(0)->Cell(i)->HalfL()
92  << " " << geom->Plane(0)->Cell(i)->HalfD()
93  << std::endl;
94  }
95  for(unsigned int i = 0; i < geom->Plane(1)->Ncells(); ++i){
96  geom->Plane(1)->Cell(i)->GetCenter(wpos,0.);
97  std::cerr << " vertical cell " << i << " "
98  << wpos[0] << " " << wpos[1] << " " << wpos[2]
99  << " " << geom->DetHalfHeight()
100  << " " << geom->Plane(1)->Cell(0)->HalfW()
101  << " " << geom->Plane(1)->Cell(0)->HalfL()
102  << " " << geom->Plane(1)->Cell(0)->HalfD()
103  << std::endl;
104  }
105 
106  try {
107  std::cerr << "Cell width " << geom->Plane(0)->Cell(0)->HalfW() << std::endl;
108  std::cerr << "Cell length " << geom->Plane(0)->Cell(0)->HalfL() << std::endl;
109  std::cerr << "Cell depth " << geom->Plane(0)->Cell(0)->HalfD() << std::endl;
110 
111  std::cerr << "Total mass " << geom->TotalMass() << std::endl;
112 
113  std::cerr << "testFindPlanes...";
114  testFindPlanes();
115  std::cerr << "complete." << std::endl;
116 
117  std::cerr << "testProject...";
118  testProject();
119  std::cerr << "complete." << std::endl;
120 
121  std::cerr << "testCellPos...";
122  testCellPos();
123  std::cerr << "complete." << std::endl;
124 
125  std::cerr << "testFindCell...";
126  testFindCell();
127  std::cerr << "complete." << std::endl;
128 
129  std::cerr << "testUniqueId...";
130  testUniqueId();
131  std::cerr << "complete." << std::endl;
132 
133  std::cerr << "testStepping...";
134  testStepping();
135  std::cerr << "complete." << std::endl;
136 
137  std::cerr << "testCellId...";
138  testCellId();
139  std::cerr << "complete." << std::endl;
140 
141  std::cerr << "testCellIdFromPos...";
143  std::cerr << "complete." << std::endl;
144  }
145  catch (...) {
146  abort();
147  }
148 
149  return;
150  }
151 
152  //......................................................................
153 
155  {
156  // Print a table of the planes and their orientations
158 
159  const std::set<unsigned int>& allpl = geom->GetPlanesByView(geo::kXorY);
160  const std::set<unsigned int>& xpl = geom->GetPlanesByView(geo::kX);
161  const std::set<unsigned int>& ypl = geom->GetPlanesByView(geo::kY);
162 
163  std::set<unsigned int>::const_iterator itrA;
164  std::set<unsigned int>::const_iterator itrX;
165  std::set<unsigned int>::const_iterator itrY;
166  std::set<unsigned int>::const_iterator itrEndA(allpl.end());
167  std::set<unsigned int>::const_iterator itrEndX(xpl.end());
168  std::set<unsigned int>::const_iterator itrEndY(ypl.end());
169 
170  // Test that x planes appear once in the all planes and x planes
171  // lists and zero times in the y planes list. Then do the same for
172  // the y planes.
173  for (itrY=ypl.begin(); itrY!=itrEndY; ++itrY) {
174  if (allpl.count(*itrY)!=1) abort();
175  if (xpl. count(*itrY)!=0) abort();
176  if (ypl. count(*itrY)!=1) abort();
177  assert(geom->Plane(*itrY)->View() == geo::kY);
178  }
179  for (itrX=xpl.begin(); itrX!=itrEndX; ++itrX) {
180  if (allpl.count(*itrX)!=1) abort();
181  if (xpl. count(*itrX)!=1) abort();
182  if (ypl. count(*itrX)!=0) abort();
183  assert(geom->Plane(*itrX)->View() == geo::kX);
184  }
185 
186  std::cout << std::endl;
187  itrA = allpl.begin();
188  for (; itrA!=itrEndA; ++itrA) {
189  std::cout << (*itrA) << " : ";
190  if (geom->Plane(*itrA)->View()==geo::kX) std::cout << " X ";
191  if (geom->Plane(*itrA)->View()==geo::kY) std::cout << " Y ";
192 
193  std::cout << " DS,same v=" << geom->NextPlaneInView(*itrA,+1);
194  std::cout << " US,same v=" << geom->NextPlaneInView(*itrA,-1);
195  std::cout << " DS,opos v=" << geom->NextPlaneOtherView(*itrA,+1);
196  std::cout << " US,opos v=" << geom->NextPlaneOtherView(*itrA,-1);
197  std::cout << std::endl;
198  }
199 
200  // Test the find routines
201  unsigned int p0 = 10;
202  unsigned int p1 = geom->NextPlaneInView(p0,+1);
203  unsigned int p2 = geom->NextPlaneInView(p0,-1);
204  unsigned int p3 = geom->NextPlaneOtherView(p0,+1);
205  unsigned int p4 = geom->NextPlaneOtherView(p0,-1);
206 
207  assert(p1>p0);
208  assert(p2<p0);
209  // in our geometry should never have to go futher than three steps
210  assert(p1-p0<=3);
211  assert(p0-p2<=3);
212  assert(geom->Plane(p0)->View() == geom->Plane(p1)->View());
213  assert(geom->Plane(p0)->View() == geom->Plane(p2)->View());
214 
215  assert(p3>p0);
216  assert(p4<p0);
217  // in our geometry should never have to go futher than three steps
218  assert(p3-p0<=3);
219  assert(p0-p4<=3);
220  assert(geom->Plane(p0)->View() != geom->Plane(p3)->View());
221  assert(geom->Plane(p0)->View() != geom->Plane(p4)->View());
222  }
223 
224  //......................................................................
225 
227  {
228  double xyz[3] = {0.};
229  int iview = 0;
230 
232  TNtuple *fnt = tfs->make<TNtuple>("cellnt","cellpos","ipl:icell:iview:x:y:z");
233 
235 
236  for (int i=0; ; ++i) {
237  const geo::PlaneGeo* plane;
238  try { plane = geom->Plane(i); }
239  catch (...) {
240  // An exception is generated when we run off the end of the
241  // detector. Use it as a flag to say "stop" when we reach the
242  // end of which ever detector geometry is loaded
243  break;
244  }
245  iview = plane->View();
246  for (unsigned int j=0; j<plane->Ncells(); ++j) {
247  const geo::CellGeo* cell = plane->Cell(j);
248  cell->GetCenter(xyz);
249  fnt->Fill(i,j,iview,xyz[0],xyz[1],xyz[2]);
250  }
251  }
252 
253  fnt->Write();
254  }
255 
256  //......................................................................
257 
259  {
260  //
261  // Test stepping. Example is similar to what one would do for photon
262  // transport. Rattles photons around inside the scintillator
263  // bouncing them off walls.
264  //
266 
267  double xyz[3], xyzCell[3] = {0,0,0};
268  double dxyz[3], dxyzCell[3] = {0,sin(0.1),cos(0.1)};
269 
270  geom->Plane(10)->Cell(0)->LocalToWorld(xyzCell,xyz);
271  geom->Plane(10)->Cell(0)->LocalToWorldVect(dxyzCell,dxyz);
272 
273  std::cerr << xyz[0] << "\t" << xyz[1] << "\t" << xyz[2] << std::endl;
274  std::cerr << dxyz[0] << "\t" << dxyz[1] << "\t" << dxyz[2] << std::endl;
275 
276  gGeoManager->InitTrack(xyz, dxyz);
277  for (int i=0; i<10; ++i) {
278  const double* pos = gGeoManager->GetCurrentPoint();
279  const double* dir = gGeoManager->GetCurrentDirection();
280  std::cerr << "pos=" << "\t"
281  << pos[0] << "\t"
282  << pos[1] << "\t"
283  << pos[2] << std::endl;
284  std::cerr << "dir=" << "\t"
285  << dir[0] << "\t"
286  << dir[1] << "\t"
287  << dir[2] << std::endl;
288  std::cerr << "mat = " <<
289  gGeoManager->GetCurrentNode()->
290  GetVolume()->GetMaterial()->GetName() << std::endl;
291 
292  gGeoManager->FindNextBoundary();
293  /* const double* norm = */ gGeoManager->FindNormal();
294  gGeoManager->Step(kTRUE,kTRUE);
295 
296  // Reflect -- bogus but OK for testing...
297  const TGeoNode* node = gGeoManager->GetCurrentNode();
298  if (strncmp("PVC",node->GetVolume()->GetMaterial()->GetName(),3)==0) {
299  gGeoManager->SetCurrentDirection(-dir[0],-dir[1],-dir[2]);
300 
301  // And step back into the scintillator
302  gGeoManager->FindNextBoundary();
303  gGeoManager->Step(kTRUE,kTRUE);
304  }
305  }
306  }
307  //......................................................................
308 
310  {
311  //
312  // Test that I can associate planes and cells with track steps
313  //
315 
316  double xyz[3], xyzCell[3] = {0,0,0};
317  double dxyz[3], dxyzCell[3] = {1/sqrt(3),
318  1/sqrt(3),
319  1/sqrt(3)};
320 
321  // Start off in plane 10, cell 10
322  geom->Plane(10)->Cell(10)->LocalToWorld(xyzCell,xyz);
323  geom->Plane(10)->Cell(10)->LocalToWorldVect(dxyzCell,dxyz);
324 
325  int aok, ip, ic;
326  gGeoManager->InitTrack(xyz, dxyz);
327  aok = geom->CurrentCell(&ip, &ic);
328 
329  clock_t start = clock();
330  for (int i=0; i<10000; ++i) {
331  const TGeoNode* node = gGeoManager->GetCurrentNode();
332  if (strncmp(node->GetVolume()->GetMaterial()->GetName(),"Sci",3)==0) {
333  aok = geom->CurrentCell(&ip, &ic);
334  if (aok==1) {
335  const double* pos = gGeoManager->GetCurrentPoint();
336  std::cout << "Stepping through (plane,cell)=\t"
337  << ip << "\t" << ic << "\t"
338  << pos[0] << "\t" << pos[1] << "\t" << pos[2] << "\t"
339  << gGeoManager->GetCurrentNodeId()
340  << std::endl;
341  }
342  else {
343  std::cout << "Failed to find plane/cell address for node"
344  << std::endl;
345  return 0;
346  }
347  }
348  gGeoManager->FindNextBoundary();
349  gGeoManager->Step(kTRUE,kTRUE);
350  }
351  clock_t end = clock();
352  std::cout << "Completed steps in "
353  << (end-start)/(float)CLOCKS_PER_SEC
354  << " seconds" << std::endl;
355  return 1;
356  }
357 
358  //......................................................................
359 
361  {
363 
364  // After the root6 upgrade most of these functions don't exist anymore
365  /*
366  // Test the basic functions used for setting UID's
367  const char* path1 = "vCellV_1112";
368  std::vector<int> id1;
369  geo::GetNodeNumbers(path1, id1);
370  if (id1[0]!=1112) abort();
371 
372  const char* path2 =
373  "vWorld_1/vMother_2/vDetEnclosure_3/"
374  "vSuperBlock_44/vBlock_55/vPlane_66/vModule_777/vCellV_888";
375  std::vector<int> id2;
376  geo::GetNodeNumbers(path2, id2);
377  if (id2[0]!=1) abort();
378  if (id2[1]!=2) abort();
379  if (id2[2]!=3) abort();
380  if (id2[3]!=44) abort();
381  if (id2[4]!=55) abort();
382  if (id2[5]!=66) abort();
383  if (id2[6]!=777) abort();
384  if (id2[7]!=888) abort();
385 
386  std::cout << "\tGetNodeNumbers : OK" << std::endl;
387 
388  // need 6 levels to get the unique ID
389  geo::CellUniqueId id3 =
390  888ULL+1000ULL*(777ULL+1000ULL*(66ULL+1000ULL*(55ULL+1000ULL*(44ULL+(1000ULL*3ULL)))));
391  if (id3 != geo::PathToUniqueId(path2)) abort();
392  std::cout << "\tPathToUniqueId : OK" << std::endl;
393  */
394 
395  // Make a test going from cell to id and id to cell to make sure the
396  // loop completes. If it does not, we have big problems.
397  int ipin = 11;
398  int icin = 12;
399  geo::CellUniqueId idin = geom->Plane(ipin)->Cell(icin)->Id();
400 
401  int ipout, icout;
402  const geo::CellGeo* cell = geom->IdToCell(idin,&ipout,&icout);
403 
404  if (!((cell->Id()==idin) && (ipin==ipout) && (icin==icout))) {
405  std::cout << "Id mismatch" << std::endl;
406  abort();
407  }
408 
409  std::cout << "\tCell->ID, ID->Cell : OK" << std::endl;
410 
411  return 1;
412  }
413 
414  //......................................................................
415 
417  {
419 
420  double xlo, xhi;
421  double ylo, yhi;
422  double zlo, zhi;
423  geom->WorldBox(&xlo, &xhi, &ylo, &yhi, &zlo, &zhi);
424 
425  double xyz[3] = {0.0, 0.0, 0.0};
426  double dxyz1[3] = { 1.0, 0.0, 0.0};
427  double dxyz2[3] = {-1.0, 0.0, 0.0};
428  double dxyz3[3] = { 0.0, 1.0, 0.0};
429  double dxyz4[3] = { 0.0,-1.0, 0.0};
430  double dxyz5[3] = { 0.0, 0.0, 1.0};
431  double dxyz6[3] = { 0.0, 0.0,-1.0};
432 
433  double xyzo[3];
434  geo::ProjectToBoxEdge(xyz, dxyz1, xlo, xhi, ylo, yhi, zlo, zhi, xyzo);
435  if (fabs(xyzo[0]-xhi)>1.E-6) abort();
436 
437  geo::ProjectToBoxEdge(xyz, dxyz2, xlo, xhi, ylo, yhi, zlo, zhi, xyzo);
438  if (fabs(xyzo[0]-xlo)>1.E-6) abort();
439 
440  geo::ProjectToBoxEdge(xyz, dxyz3, xlo, xhi, ylo, yhi, zlo, zhi, xyzo);
441  if (fabs(xyzo[1]-yhi)>1.E-6) abort();
442 
443  geo::ProjectToBoxEdge(xyz, dxyz4, xlo, xhi, ylo, yhi, zlo, zhi, xyzo);
444  if (fabs(xyzo[1]-ylo)>1.E-6) abort();
445 
446  geo::ProjectToBoxEdge(xyz, dxyz5, xlo, xhi, ylo, yhi, zlo, zhi, xyzo);
447  if (fabs(xyzo[2]-zhi)>1.E-6) abort();
448 
449  geo::ProjectToBoxEdge(xyz, dxyz6, xlo, xhi, ylo, yhi, zlo, zhi, xyzo);
450  if (fabs(xyzo[2]-zlo)>1.E-6) abort();
451  }
452 
453  //......................................................................
454 
456  {
458 
459  double xyz[100][3];
460  unsigned int iplane[100];
461  unsigned int icell[100];
462  geo::CellUniqueId id[100];
463 
464  // Get a list of 100 'random' cells, their positions and their cell
465  // id's
466  unsigned int ip = 0;
467  unsigned int ic = 0;
468  for (int i=0; i<100; ++i) {
469  ip += 2437;
470  ic += 6577;
471 
472  iplane[i] = ip%geom->NPlanes();
473  icell[i] = ic%geom->Plane(iplane[i])->Ncells();
474 
475  geom->Plane(iplane[i])->Cell(icell[i])->GetCenter(xyz[i]);
476  id[i] = geom->Plane(iplane[i])->Cell(icell[i])->Id();
477  }
478 
479  // Now check that we get the same results using the "CellId"
480  // function
481  for (int i=0; i<100; ++i) {
482  geo::CellUniqueId idtest;
483 
484  idtest = geom->CellId(xyz[i][0], xyz[i][1], xyz[i][2]);
485  if (idtest != id[i]) {
486  std::cerr << "testCellId failed. Id's do not match for position "
487  << id[i] << " != "
488  << idtest << " "
489  << "plane,cell = "
490  << iplane[i] << " "
491  << icell[i] << " "
492  << std::endl;
493 
494  double pos[3];
495  geom->Plane(iplane[i])->Cell(icell[i])->GetCenter(pos);
496  std::cerr << "Positions: " << std::endl
497  << xyz[i][0] << "\t"
498  << xyz[i][1] << "\t"
499  << xyz[i][2] << std::endl
500  << pos[0] << "\t"
501  << pos[1] << "\t"
502  << pos[2] << std::endl;
503 
504  abort();
505  }
506  }
507  }
508 
509  //......................................................................
510 
512  {
514  unsigned int p = 0, c = 0;
515 
516 // const geo::CellUniqueId testid = geom->CellId(25.0808, 57.9953, 301.284);
517 // geom->IdToCell(testid, &p, &c);
518 
519 // std::cerr << "cell id from test point is "
520 // << testid
521 // << " corresponds to " << p << " " << c
522 // << std::endl;
523 
524 
525  // get the center of a cell
526  for(p = 0; p < geom->NPlanes(); ++p){
527  for(c = 0; c < geom->Plane(p)->Ncells(); ++c){
528 
529  const geo::CellGeo* cell = geom->Plane(p)->Cell(c);
530  const geo::CellUniqueId id = cell->Id();
531 
532  double cellCenter[3] = {0.};
533  cell->GetCenter(cellCenter);
534 
535  // use the geometry method for getting cell id based on position and direction
536  const geo::CellUniqueId posid = geom->CellId(cellCenter[0], cellCenter[1], cellCenter[2]);
537 
538  if(posid != id){
539  std::cerr << "cell id from position does not match cell id from CellGeo "
540  << posid << " vs " << id
541  << std::endl;
542  abort();
543  } // end if ids are the same
544  } // end loop over cells
545  }// end loop over planes
546  }
547 
549 }//end namespace
double HalfL() const
Definition: CellGeo.cxx:198
double HalfD() const
Definition: CellGeo.cxx:205
void GetCenter(double *xyz, double localz=0.0) const
Definition: CellGeo.cxx:159
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
X or Y views.
Definition: PlaneGeo.h:30
void LocalToWorld(const double *local, double *world) const
Definition: CellGeo.cxx:80
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
const char * p
Definition: xmltok.h:285
T sqrt(T number)
Definition: d0nt_math.hpp:156
Vertical planes which measure X.
Definition: PlaneGeo.h:28
unsigned int Ncells() const
Number of cells in this plane.
Definition: PlaneGeo.h:43
double HalfW() const
Definition: CellGeo.cxx:191
OStream cerr
Definition: OStream.cxx:7
TString ip
Definition: loadincs.C:5
const PlaneGeo * Plane(unsigned int i) const
DEFINE_ART_MODULE(TestTMapFile)
GeometryTest(fhicl::ParameterSet const &pset)
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
Give every cell in the geometry a unique ID number based on the TGeo path to the node.
int CurrentCell(int *ip, int *ic) const
const CellUniqueId & Id() const
Definition: CellGeo.h:55
void LocalToWorldVect(const double *local, double *world) const
Definition: CellGeo.cxx:90
View_t View() const
Which coordinate does this plane measure.
Definition: PlaneGeo.h:53
Float_t E
Definition: plot.C:20
void WorldBox(double *xlo_cm, double *xhi_cm, double *ylo_cm, double *yhi_cm, double *zlo_cm, double *zhi_cm) const
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
Definition: GeometryBase.h:243
Geometry information for a single readout plane.
Definition: PlaneGeo.h:36
int evt
const CellUniqueId CellId(const double &x, const double &y, const double &z, double dxds=0., double dyds=0., double dzds=1., double step=0.01) const
Collect Geo headers and supply basic geometry functions.
const double j
Definition: BetheBloch.cxx:29
void analyze(art::Event const &evt)
double DetHalfHeight() const
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
std::string GDMLFile() const
Definition: GeometryBase.h:235
OStream cout
Definition: OStream.cxx:6
unsigned long long int CellUniqueId
Definition: CellUniqueId.h:15
double TotalMass(const char *volume="vDetEnclosure") const
T * make(ARGS...args) const
T sin(T number)
Definition: d0nt_math.hpp:132
double DetHalfWidth() const
TDirectory * dir
Definition: macro.C:5
const unsigned int NextPlaneInView(unsigned int p, int d=+1) const
const std::set< unsigned int > & GetPlanesByView(View_t v=kXorY) const
void geom(int which=0)
Definition: geom.C:163
T cos(T number)
Definition: d0nt_math.hpp:78
assert(nhit_max >=nhit_nbins)
void ProjectToBoxEdge(const double xyz[], const double dxyz[], double xlo, double xhi, double ylo, double yhi, double zlo, double zhi, double xyzout[])
Project along a direction from a particular starting point to the edge of a box.
Definition: Geo.cxx:38
unsigned int NPlanes() const
const CellGeo * IdToCell(const CellUniqueId &id, int *iplane, int *icell) const
Encapsulate the cell geometry.
Definition: CellGeo.h:25
Helper for AttenCurve.
Definition: Path.h:10
Encapsulate the geometry of one entire detector (near, far, ndos)
const unsigned int NextPlaneOtherView(unsigned int p, int d=+1) const
std::string FileBaseName() const