GeometryBase.cxx
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 /// \file GeometryBase.cxx
3 /// \brief
4 /// \version $Id: GeometryBase.cxx,v 1.12 2012-12-19 05:34:46 gsdavies Exp $
5 /// \author messier@indiana.edu
6 ////////////////////////////////////////////////////////////////////////
7 
9 
10 #include <sys/stat.h>
11 #include <sys/types.h>
12 
13 #include <cassert>
14 #include <iostream>
15 #include <stdio.h>
16 
17 //ROOT includes
18 #include "TGeoManager.h"
19 #include "TGeoNode.h"
20 #include "TStopwatch.h"
21 #include "TVirtualGeoPainter.h"
22 
23 // Framework includes
24 #include "cetlib_except/exception.h"
25 #include "cetlib/search_path.h"
27 #include "CLHEP/Random/RandFlat.h"
28 
29 // NOvA includes
31 #include "GeometryObjects/Geo.h"
33 #include "Utilities/func/MathUtil.h"
34 
35 namespace geo
36 {
37  //....................................................................
38  static bool plane_sort(const PlaneGeo* p1, const PlaneGeo* p2)
39  {
40  double xyz1[3], xyz2[3];
41  p1->Cell(0)->GetCenter(xyz1);
42  p2->Cell(0)->GetCenter(xyz2);
43  return xyz1[2]<xyz2[2];
44  }
45 
46  //......................................................................
48  : fParams(params)
49  , fDetId(novadaq::cnv::kUNKNOWN_DET)
50  , fRunNumber(0)
51  {
53 
54  // If they want to force FCL usage only, they better have supplied a file
55  assert(!(fParams.ForceUseFCLOnly() && fParams.GDMLFile().empty()));
56  if(fParams.GDMLFile().empty()){
57  // That's fine. Normal operating mode. We'll get GDML information in
58  // beginRun().
59  return;
60  }
61 
62  //NOTE: Setting fROOTFile to the same thing as GDMLFile, in the future
63  //this should be removed but for now EventGeneratorBase has a print
64  //statement requiring the root file name
66 
68 
69  // constructor decides if initialized value is a path or an environment variable
70  cet::search_path sp("FW_SEARCH_PATH");
71 
72  if ( !sp.find_file(fParams.GDMLFile(), fGDMLFile) ) {
73  throw cet::exception("GeometryFileLoad")
74  << "cannot find GDML file " << fGDMLFile << " bail ungracefully\n"
75  << __FILE__ << ":" << __LINE__ << "\n";
76  }
77  mf::LogWarning("LoadNewGeometry") << "loading new geometry files\n"
78  << fGDMLFile << "\n";
79 
81  }
82 
83  //......................................................................
85  {
86  for (unsigned int i=0; i<fPlanes.size(); ++i) {
87  if (fPlanes[i]) { delete fPlanes[i]; fPlanes[i] = 0; }
88  }
89 
90  //clean up an existing temporary file
91  if (fGDMLFile.find("TEMP_GEO_GDML") != std::string::npos){
92  this->RemoveTmpFile(fGDMLFile);
93  }
94 
95 
96  }
97 
98  //....................................................................
99  ///
100  /// Method to make a temporary file holding Detector geometry (gdml) file
101  ///
103  {
104  int fStat = -1; //file status variable, -1 indicates error
105  FILE *tmpFile; //file variable
106 
107  //hold temporary gdml file name
108  std::string tmpFileName;
109 
110  //Use mkstemp to make a temporary file, mkstemp takes in a template
111  //character array as an argument. The last 6 characters must be XXXXXX
112  //mkstemp will then attempt to open a file with that template and the
113  //last 6 characters will be changed to be unique. If the file cannot
114  //be successfully opened an error code will be returned.
115  //If a file is properly opened it will be unique. File is set to be
116  //opened in /temp.
117  char tmp[256];
118  sprintf(tmp,"%s/TEMP_GEO_GDML_XXXXXX", fParams.StoreTempGeo().c_str());
119  //check to see if file opens properly
120  if ( ((fStat = mkstemp(tmp)) == -1) ||
121  ((tmpFile = fopen(tmp, "w+")) == NULL) ) {
122 
123  if (fStat != -1){
124  fclose(tmpFile);
125  remove(tmp);
126  std::cerr<<"Error making temporary file "
127  <<tmp<<std::endl;
128  }
129  //if file could not be created properly exit and return -1
130  std::cerr<<"Failed to produce temporary geometry file"<<std::endl;
131  return -1;
132  }
133 
134  //this means a unique file has been properly created and opened
135  fclose(tmpFile); //properly close file
136 
137  //add temp file name to string with proper .gdml extension
138  tmpFileName += tmp;
139  tmpFileName += ".gdml";
140 
141  //now attempt to open new file with .gdml extension
142  FILE *gdmlFile;
143  gdmlFile = fopen(tmpFileName.c_str(), "w+");
144  if (gdmlFile == NULL){
145  //file was not opened properly, remove temp files
146  //and exit with a status -1
147  remove(tmp);
148  remove(tmpFileName.c_str());
149  std::cerr<<"Error making temporary file "
150  <<tmpFileName<<std::endl;
151  return -1;
152  }
153 
154  //temporary gdml file was properly opened, now put our gdml info in file
155  fputs (gdmlInfo.c_str(), gdmlFile);
156 
157  //now remove initial temp file without the .gdml extension
158  remove(tmp);
159 
160  //close temporary gdml file
161  fclose(gdmlFile);
162 
163  //now store name of temporary gdml file in GDMLFile
164  //this way the temp file can also be accessed by GEANT4,
165  //which loads the geometry directly
166  fGDMLFile = tmpFileName;
167 
168  //return 0 indicating successful temporary file creation
169  return 0;
170  }
171 
172  //......................................................................
173  ///
174  /// Method to remove temporary file that holds detector geometry file
176  remove(fileName.c_str());
177  }
178 
179  //......................................................................
181  {
182  return ExtractGDML(fGDMLFile, true);
183  }
184 
185  //......................................................................
187  {
188  if(!fullpath) fname = FindGDMLFile(fname);
189 
190  std::string gdmlInfo;
191 
192  std::ifstream geoFile(fname.c_str());
194 
195  if (geoFile.is_open()){
196  while (geoFile.good()){
197  getline(geoFile, line);
198  gdmlInfo += line;
199  gdmlInfo += '\n';
200  }
201  geoFile.close();
202  }
203 
204  return gdmlInfo;
205  }
206 
207  //......................................................................
209  {
210  cet::search_path sp("FW_SEARCH_PATH");
211  size_t lastSlash = fname.find_last_of("/");
212  std::string basename = fname.substr(lastSlash+1);
213 
214  fname = "Geometry/gdml/"+basename;
215 
217  if(!sp.find_file(fname, ret)){
218  throw cet::exception("FindGDMLFile")
219  << "cannot find GDML file " << fname << " bail ungracefully\n"
220  << __FILE__ << ":" << __LINE__ << "\n";
221  }
222  return ret;
223  }
224 
225  //......................................................................
226  ///
227  /// Method to load detector geometry file
228  ///
231 
232  int old_verbosity = gGeoManager->GetVerboseLevel();
233 
234  // TGeoManager is too verbose when loading geometry.
235  // Make it quiet.
236  gGeoManager->SetVerboseLevel(0);
237 
238  // Open the GDML file
239  struct stat sb;
240  if (gdmlfile.empty() || stat(gdmlfile.c_str(), &sb)!=0) {
241  // failed to resolve the file name
242  throw cet::exception("GeometryFileLoad")
243  << "empty GDML string " << gdmlfile << " bail ungracefully\n"
244  << __FILE__ << ":" << __LINE__ << "\n";
245  }
246  // clear the plane array and id map
247  for (unsigned int i=0; i<fPlanes.size(); ++i) {
248  if (fPlanes[i]) { delete fPlanes[i]; fPlanes[i] = 0; }
249  }
250  fPlanes.clear();
251  fPlanesInMuonCatcher.clear();
252  fIdMap.clear();
253  CellGeo::clear();
254 
255  TGeoManager::Import(gdmlfile.c_str());
256 
257  this->SetDrawOptions();
258 
259  std::vector<const TGeoNode*> n(16);
260  n[0] = gGeoManager->GetTopNode();
261  this->FindPlanes(n, 0);
262  sort(fPlanes.begin(), fPlanes.end(), plane_sort);
264 
265  this->SetDetectorSize();
266  this->BuildMaps();
267 
268  if(detId == novadaq::cnv::kUNKNOWN_DET){
269  // have to set the fDetId based on the input file then
270  if( gdmlfile.find("ndos") != std::string::npos) detId = novadaq::cnv::kNDOS;
271  else if(gdmlfile.find("near") != std::string::npos) detId = novadaq::cnv::kNEARDET;
272  else if(gdmlfile.find("far") != std::string::npos) detId = novadaq::cnv::kFARDET;
273  else if(gdmlfile.find("testbeam-1x1-2block") != std::string::npos) detId = novadaq::cnv::kTESTBEAM;
274  else if(gdmlfile.find("testbeam-1x1-3block") != std::string::npos) detId = novadaq::cnv::kTESTBEAM;
275  else if(gdmlfile.find("testbeam-2x2-2block") != std::string::npos) detId = novadaq::cnv::kTESTBEAM;
276  else if(gdmlfile.find("testbeam-2x2-3block") != std::string::npos) detId = novadaq::cnv::kTESTBEAM;
277  }
278  setDetectorID(detId);
279 
282  else{
288  fDetectorBigBoxZLo = 0.;
289  }
290 
291  mf::LogInfo("LoadGeometryFile")
292  << "Detector geometry loaded from " << gdmlfile;
293 
294 
295  // Restore verbosity level before we leave because
296  // we have to be polite and listen
297  gGeoManager->SetVerboseLevel(old_verbosity);
298 
299  return true;
300  }// end of LoadGeometryFile
301 
302 
303  //......................................................................
304  ///
305  /// Method to extract gdml file basename and store as an std::string
306  ///
308  assert(fGDMLFile.find(".gdml") != std::string::npos);
309 
310  // to find the base name of the file
311  // use string operations to find the last
312  // instance of "/" in the fROOTFile string
313  // and the last instance of ".", the base
314  // string should be between those
315  size_t lastSlash = fGDMLFile.find_last_of("/");
316  std::string almost = fGDMLFile.substr(lastSlash+1);
317 
318  // resize almost to drop the file ".gdml"
319  almost.resize(almost.length() - 5);
320 
321  return almost;
322  }// end of FileBaseName
323 
324  //......................................................................
325  ///
326  /// Return the plane and cell number for the current tracking position
327  /// held by the gGeoManager
328  ///
329  /// \param planeid : pointer to plane number, >=0 on success, -1 on failure
330  /// \param cellid : pointer to cell number, >=0 on success, -1 on failure
331  /// \returns >0 if cell is found, <0 if not found.
332  ///
333  int GeometryBase::CurrentCell(int* planeid, int* cellid) const {
334  if (! this->IdToCell(this->CurrentCellId(), planeid, cellid)) {
335  *planeid = -1;
336  *cellid = -1;
337  throw cet::exception("UNKNOWN_ID")
338  << "unknown id " << this->CurrentCellId() << "\n"
339  << __FILE__ << ":" << __LINE__ << "\n";
340  }
341  return 1;
342  }// end of
343 
344  //......................................................................
345  ///
346  /// Return the unique cell id for the current cell -- zero if not in any cell
347  ///
349  const std::string path = gGeoManager->GetPath();
350  //std::cout << "path to cell is " << path << std::endl;
351 
352  // check that the path actually contains a cell
353  if(path.find("Cell") == std::string::npos){
354  LOG_DEBUG("Geometry") << "no Cell in the path to the current volume:" << path << std::endl;
355  return 0;
356  }
357  if(path.find("HAirBubble") != std::string::npos){
358  LOG_DEBUG("Geometry") << "This an air bubble, not a cell:" << path << std::endl;
359  return 0;
360  }
361 
362  // We used to return PathToUniqueId(path.c_str()) at this point, but as of
363  // the root6 upgrade, the paths are all like
364  // vPlaneV_0/vModuleV_0/vExtruV_0/vCellV_0 rather than
365  // vPlaneV_3/vModuleV_1/vExtruV_4/vCellV_1, so use the list of parent
366  // pointers to determine the unique ID instead. The code here is derived
367  // from what TGeoManager::GetPath() does.
368  TGeoNodeCache* cache = gGeoManager->GetCurrentNavigator()->GetCache();
369  std::vector<const TGeoNode*> ns;
370  for(int i = 0; i <= cache->GetLevel(); ++i){
371  ns.push_back(((TGeoNode**)cache->GetBranch())[i]);
372  }
373  return NodesToUniqueId(ns, ns.size()-1);
374  }// end of CurrentCellId
375 
376  //......................................................................
377  ///
378  /// Return the unique cel id for a specified location and direction in the detector
379  ///
380  /// \param x - x position coordinate (cm, world coordinates)
381  /// \param y - y position coordinate (cm, world coordinates)
382  /// \param z - z position coordinate (cm, world coordinates)
383  /// \param dxds - dxds direction in which to explore a step
384  /// \param dyds - dyds direction in which to explore a step
385  /// \param dzds - dzds direction in which to explore a step
386  /// \param step - step size to take
387  ///
388  /// For the sake of efficiency, this method does no checking to ensure
389  /// that the location is in fact inside a cell. The user is expected
390  /// to perform that check before using the position (x,y,z).
391  ///
392  /// Note further, that this uses the gGeoManager->InitTrack method and
393  /// will disrupt TGeo-based tracking and should be called once any
394  /// TGeo-based tracking is completed.
395  ///
396  const CellUniqueId GeometryBase::CellId(const double& x, const double& y, const double& z,
397  double dxds, double dyds, double dzds,
398  double step) const {
399  // The strategy for finding the cell ID is to try the exact point
400  // we've been handed. If that does not return a valid cell ID, then
401  // try shifting the point around slightly to account for tracking
402  // and edge finding uncertainties. If after trying that we still
403  // don't have a valid ID, return 0.
404 
405  gGeoManager->InitTrack(x, y, z, dxds, dyds, dzds);
406  gGeoManager->FindNextBoundary(1.E-6);
407  CellUniqueId ID = this->CurrentCellId();
408  if (ID != 0) return ID;
409 
410  double eps = 1.0E-1;
411  const double s = std::max(step, 1.e-3);
412 
413  // walk the step points backwards and forwards to see if we can find the
414  // cell
415  while(eps < 1.){
416 
417  const double eps_s = eps*s;
418 
419  // std::cout << "step back pos " << x-dxds*eps*s
420  // << "," << y-dyds*eps*s << "," << z-dzds*eps*s
421  // << " " << eps << " " << step << std::endl;
422  gGeoManager->InitTrack(x-dxds*eps_s, y-dyds*eps_s, z-dzds*eps_s, dxds, dyds, dzds);
423  gGeoManager->FindNextBoundary(1.E-6);
424  ID = this->CurrentCellId();
425  if (ID != 0) return ID;
426 
427  // std::cout << "step forward pos " << x+dxds*eps*s
428  // << "," << y+dyds*eps*s << "," << z+dzds*eps*s << std::endl;
429  gGeoManager->InitTrack(x+dxds*eps_s, y+dyds*eps_s, z+dzds*eps_s, dxds, dyds, dzds);
430  ID = this->CurrentCellId();
431  if (ID != 0) return ID;
432 
433  eps += 1.e-1;
434  }
435 
436  // No valid ID found
437  return 0;
438  }// end of CellId
439 
440  //......................................................................
441  ///
442  /// Method to obtain current detector material at given (x,y,z) coordinates
443  ///
444  TGeoMaterial* GeometryBase::Material(double x, double y, double z) const {
445  const TGeoNode* node = gGeoManager->FindNode(x,y,z);
446  return node->GetMedium()->GetMaterial();
447  }// end of Material
448 
449  //......................................................................
450  ///
451  /// Is the material m an active or passive material?
452  ///
453  bool GeometryBase::IsActive(const TGeoMaterial* m) const {
454  const std::string scint = "Scintillator";
455  if (scint == m->GetName()) return true;
456  return false;
457  }
458 
459  //......................................................................
460  unsigned int GeometryBase::NPlanes() const
461  {
463  return fPlanes.size();
464  }
465 
466  //......................................................................
468  {
470  return fDetHalfWidth;
471  }
472 
473  //......................................................................
475  {
477  return fDetHalfHeight;
478  }
479 
480  //......................................................................
481  double GeometryBase::DetLength() const
482  {
484  return fDetLength;
485  }
486 
487  //......................................................................
488 
489  ///
490  /// A typical y-position value at the surface (where earth meets air)
491  /// for this detector site
492  ///
493  /// \returns typical y position at surface in units of cm
494  ///
495  /// \throws geo::Exception if detector ID is not set properly
496  ///
497  /// \todo: Only far det number is reasonable accurate. Near and IPND
498  /// need some thought
499  ///
500  double GeometryBase::SurfaceY() const {
501 
502  switch(fDetId){
504  return 9.94E2; // far : surface ~10m above det. center
505 
506  case novadaq::cnv::kNDOS:
507  return -1.966E2; // NDOS: surface ~2m below det. center
508 
510  return 130.0E2; // near: surface ~130m above det. center
511 
513  return -1.966E2; // NDOS: surface ~2m below det. center
514 
515  default:
516  throw cet::exception("BadGeoConfig")
517  << "Bad geometry configuration: " << (int)fDetId << "\n"
518  << __FILE__ << ":" << __LINE__ << "\n";
519  }
520  }// end of SurfaceY
521 
522  //---------------------------------------------------------------------------
524  {
525  // These are obtained by plotting rec.mc.nu.p.fP.fZ/rec.mc.nu.p.fE etc from
526  // CAF files and fitting the resulting distribution to a gaussian.
527  // A true vertex inside of the detector is require as well as the NuE group FA reconstruction
528  // and containment cuts.
529  // CAF Files used:
530  // NearDet:
531  // prod_caf_S15-05-22a_nd_genie_fhc_nonswap_ndnewpos
532  // FarDet:
533  // prod_caf_S15-05-22a_fd_genie_fhc_nonswap_fdfluxv08
534  // Unknown for NDOS
535 
536  switch(fDetId){
537  case novadaq::cnv::kFARDET: return TVector3(-6.833e-05,
538  +6.388e-02,
539  +9.980e-1).Unit();
540  case novadaq::cnv::kNDOS: return TVector3(-3.845e-4,
541  +6.517e-2,
542  +9.967e-1).Unit();
543  case novadaq::cnv::kNEARDET: return TVector3(-8.424e-04,
544  -6.174e-02,
545  +9.981e-1).Unit();
546  case novadaq::cnv::kTESTBEAM: return TVector3(0,
547  0,
548  1).Unit();
549  default:
550  throw cet::exception("BadGeoConfig")
551  << "Bad geometry configuration: " << fDetId << "\n"
552  << __FILE__ << ":" << __LINE__ << "\n";
553  }
554  }
555 
556  //---------------------------------------------------------------------------
558  {
559  // No simulation available for NDOS Booster. These values copied from the Cana
560  // package, where the comment read:
561  //
562  // NDOS BNB direction from Zelimir Djurcic <zdjurcic@hep.anl.gov> 2/26/2011,
563  // private communication.
564  //
565  // ND BNB direction from Ryan Murphy <rwmurphy@indiana.edu> 10/14/2016
566  // Used same method as described in NuMIBeamDirection()
567  // CAF File used:
568  // NearDet:
569  // /nova/ana/users/rwmurphy/bnb_mc_all_hadd_caf.caf.root
570 
571  switch(fDetId){
572  case novadaq::cnv::kNDOS: return TVector3(-0.367, 0.012, 0.930).Unit();
573 
574  case novadaq::cnv::kNEARDET: return TVector3(-2.961e-01, -1.206e-01,+9.475e-01).Unit();
575 
576  default:
577  throw cet::exception("BadGeoConfig")
578  << "Bad geometry configuration: " << fDetId << "\n"
579  << __FILE__ << ":" << __LINE__ << "\n";
580  }
581  }
582 
583  ///
584  /// Return the ranges of x,y and z for the "world volume" that the
585  /// entire geometry lives in. If any pointers are 0, then those
586  /// coordinates are ignored.
587  ///
588  /// \param xlo : On return, lower bound on x positions
589  /// \param xhi : On return, upper bound on x positions
590  /// \param ylo : On return, lower bound on y positions
591  /// \param yhi : On return, upper bound on y positions
592  /// \param zlo : On return, lower bound on z positions
593  /// \param zhi : On return, upper bound on z positions
594  ///
595  void GeometryBase::WorldBox(double* xlo, double* xhi,
596  double* ylo, double* yhi,
597  double* zlo, double* zhi) const {
598  const TGeoShape* s = gGeoManager->GetVolume("vWorld")->GetShape();
599  assert(s);
600 
601  if (xlo || xhi) {
602  double x1, x2;
603  s->GetAxisRange(1,x1,x2); if (xlo) *xlo = x1; if (xhi) *xhi = x2;
604  }
605  if (ylo || yhi) {
606  double y1, y2;
607  s->GetAxisRange(2,y1,y2); if (ylo) *ylo = y1; if (yhi) *yhi = y2;
608  }
609  if (zlo || zhi) {
610  double z1, z2;
611  s->GetAxisRange(3,z1,z2); if (zlo) *zlo = z1; if (zhi) *zhi = z2;
612  }
613  }// end of WorldBox
614 
615  ///
616  /// Return the ranges of x,y and z for the "Detector volume" that the
617  /// entire geometry lives in. If any pointers are 0, then those
618  /// coordinates are ignored.
619  ///
620  /// \param xlo : On return, lower bound on x positions
621  /// \param xhi : On return, upper bound on x positions
622  /// \param ylo : On return, lower bound on y positions
623  /// \param yhi : On return, upper bound on y positions
624  /// \param zlo : On return, lower bound on z positions
625  /// \param zhi : On return, upper bound on z positions
626  ///
628  double* ylo, double* yhi,
629  double* zlo, double* zhi) const {
630  const TGeoShape* s = gGeoManager->GetVolume("vDetEnclosure")->GetShape();
631  assert(s);
632 
633  if (xlo || xhi) {
634  double x1, x2;
635  s->GetAxisRange(1,x1,x2); if (xlo) *xlo = x1; if (xhi) *xhi = x2;
636  }
637  if (ylo || yhi) {
638  double y1, y2;
639  s->GetAxisRange(2,y1,y2); if (ylo) *ylo = y1; if (yhi) *yhi = y2;
640  }
641  if (zlo || zhi) {
642  double z1, z2;
643  s->GetAxisRange(3,z1,z2); if (zlo) *zlo = z1; if (zhi) *zhi = z2;
644  }
645  }// end ofDetectorEnclosureBox
646 
647  ///---------------------------------------------------------------------------
648  /// Return the ranges of x,y and z for the "Detector volume" that the
649  /// entire geometry lives in. If any pointers are 0, then those
650  /// coordinates are ignored.
651  /// \param xlo : On return, lower bound on x positions
652  /// \param xhi : On return, upper bound on x positions
653  /// \param ylo : On return, lower bound on y positions
654  /// \param yhi : On return, upper bound on y positions
655  /// \param zlo : On return, lower bound on z positions
656  /// \param zhi : On return, upper bound on z positions
657  ///
658  void GeometryBase::DetectorBigBox(double* xlo, double* xhi,
659  double* ylo, double* yhi,
660  double* zlo, double* zhi) const {
661  //we don't want to return without passing along the big box coordinates
662  //note the big box is set by default to just be the detector size
663  //if 'isDetectorBigBox' is false
664  //if(!this->isDetectorBigBoxUsed())return;
665  *xlo = fDetectorBigBoxXLo;
666  *xhi = fDetectorBigBoxXHi;
667  *ylo = fDetectorBigBoxYLo;
668  *yhi = fDetectorBigBoxYHi;
669  *zlo = fDetectorBigBoxZLo;
670  *zhi = fDetectorBigBoxZHi;
671  }// end of DetectorBigBox
672 
673  //......................................................................
674  bool GeometryBase::IntersectsDetector(double *xyz, double* dxyz) const {
675  return geo::IntersectsBox(xyz, dxyz,
678  0., fDetLength);
679  }// end of IntersectsDetector
680 
681  //......................................................................
682  bool GeometryBase::IntersectsBigBox(double *xyz, double* dxyz) const {
683  return geo::IntersectsBox(xyz, dxyz,
687  }// end of IntersectsBigBox
688 
689  //......................................................................
690 
692 
693  //......................................................................
694 
695  ///
696  /// Recursively search for planes in the geometry
697  ///
698  /// \param n - Array holding the path of nodes traversed
699  /// \param depth - Depth of the search to this point
700  /// \param inMuonCatcher - internal flag, recursed into the muon catcher?
701  /// \throws geo::Exception is maximum depth is exceeded.
702  ///
703  void GeometryBase::FindPlanes(std::vector<const TGeoNode*>& n,
704  unsigned int depth,
705  bool inMuonCatcher)
706  {
707  const char* nm = n[depth]->GetName();
708  if (strncmp(nm,"vPlane", 6)==0) {
709  this->MakePlane(n, depth, inMuonCatcher);
710  return;
711  }
712 
713  if(strncmp(nm, "vBlockMuon", 10) == 0) inMuonCatcher = true;
714 
715  // Explore the next layer down
716  unsigned int deeper = depth+1;
717  if (deeper>=n.size()) {
718  throw cet::exception("BAD_NODE")
719  << "Exceeded maximum depth (" << n.size() << ") " << deeper << "\n"
720  << __FILE__ << ":" << __LINE__ << "\n";
721  }
722  const TGeoVolume* v = n[depth]->GetVolume();
723  int nd = v->GetNdaughters();
724  for (int i=0; i<nd; ++i) {
725  n[deeper] = v->GetNode(i);
726  this->FindPlanes(n, deeper, inMuonCatcher);
727  }
728  }// end of FindPlanes
729 
730  //........................................................................
731  void GeometryBase::MakePlane(std::vector<const TGeoNode*>& n,
732  unsigned int depth,
733  bool inMuonCatcher)
734  {
735  fPlanes.push_back(new PlaneGeo(n, depth));
736 
737  if(inMuonCatcher) fPlanesInMuonCatcher.push_back(fPlanes.back());
738  }
739 
740  //......................................................................
741  ///
742  /// Return the geometry description of the ith plane in the detector.
743  ///
744  /// \param i : input plane number, starting from 0
745  /// \returns plane geometry for ith plane
746  ///
747  /// \throws geo::Exception if iplane is outside allowed range
748  ///
749  const PlaneGeo* GeometryBase::Plane(unsigned int i) const {
750  if (i>=fPlanes.size()) {
751  throw cet::exception("BAD_PLANE_NUMBER")
752  << "Plane index " << i << " out of range " << fPlanes.size() << "\n"
753  << __FILE__ << ":" << __LINE__ << "\n";
754  return 0;
755  }
756  return fPlanes[i];
757  }//end of Plane
758 
759  //........................................................................
760  ///
761  /// Return the transverse position of the cell in which ever view its
762  /// in
763  ///
764  /// \param ip - plane number
765  /// \param ic - cell number
766  /// \param w - optional position along the local cell z-axis
767  ///
768  /// \returns - transverse position in units of cm
769  ///
770  /// \throws - geo::Exception is plane, cell, or view are out of range
771  ///
772  double GeometryBase::CellTpos(unsigned int ip, unsigned int ic, double w) const {
773  if (ip>=this->NPlanes()) {
774  throw cet::exception("BAD_PLANE_NUMBER")
775  << " iplane " << ip << " [0," << this->NPlanes() << ")\n"
776  << __FILE__ << ":" << __LINE__ << "\n";
777  }
778  const PlaneGeo* p = this->Plane(ip);
779  if (p) {
780  if (ic>=p->Ncells()) {
781  throw cet::exception("BAD_CELL_NUMBER")
782  << " iplane " << ic << " [0," << p->Ncells() << ")\n"
783  << __FILE__ << ":" << __LINE__ << "\n";
784  }
785  const CellGeo* c = p->Cell(ic);
786  if (c) {
787  double pos[3];
788  c->GetCenter(pos,w);
789  if (p->View() == kX) { return pos[0]; }
790  if (p->View() == kY) { return pos[1]; }
791  throw cet::exception("BAD_PLANE_VIEW")
792  << "view " << p->View() << "\n"
793  << __FILE__ << ":" << __LINE__ << "\n";
794  return 0.0;
795  }
796  }
797  return 0.0;
798  }// end of CellTpos
799 
800  //......................................................................
801  ///
802  /// Given a plane and cell number, return which view the cell
803  /// measures, its nominal center position, and nominal size about that
804  /// center position.
805  ///
806  /// \param ip - plane number
807  /// \param ic - cell number
808  /// \param view - on return, the view this cell measures in
809  /// \param pos - the nominal (xyz) center position of this cell (cm)
810  /// \param dpos - half-dimensions of the cell in xyz (cm)
811  /// \throw geo::Exception
812  ///
813  void GeometryBase::CellInfo(unsigned int ip, unsigned int ic,
814  View_t* view, double* pos, double* dpos) const {
815  if (ip>=this->NPlanes()) {
816  throw cet::exception("BAD_PLANE_NUMBER")
817  << " iplane " << ip << " [0," << this->NPlanes() << ")\n"
818  << __FILE__ << ":" << __LINE__ << "\n";
819  }
820  const PlaneGeo* p = this->Plane(ip);
821  if (p) {
822  if (view) *view = p->View();
823 
824  if (ic>=p->Ncells()) {
825  throw cet::exception("BAD_CELL_NUMBER")
826  << " iplane " << ic << " [0," << p->Ncells() << ")\n"
827  << __FILE__ << ":" << __LINE__ << "\n";
828  }
829  const CellGeo* c = p->Cell(ic);
830  if (pos) c->GetCenter(pos);
831 
832  if (dpos) {
833  if (p->View() == kX) {
834  dpos[0] = c->HalfW();
835  dpos[1] = c->HalfL();
836  dpos[2] = c->HalfD();
837  }
838  else if (p->View() == kY) {
839  dpos[0] = c->HalfL();
840  dpos[1] = c->HalfW();
841  dpos[2] = c->HalfD();
842  }
843  }
844  }
845  }// end of CellInfo
846 
847  //......................................................................
848  ///
849  /// Given a unique cell identifier, look up the plane and cell
850  /// numbers.
851  ///
852  /// \param id : Cell identifier
853  /// \param iplane : pointer to return plane number
854  /// \param icell : pointer to return cell number
855  ///
856  /// \returns cell geometry description
857  /// \throws geo::Exception if id is not found
858  ///
860  int* iplane,
861  int* icell) const {
862  const OfflineChan offline_chan = getPlaneCellMap(id);
863  *iplane = (int)offline_chan.Plane();
864  *icell = (int)offline_chan.Cell();
865  return this->Plane(*iplane)->Cell(*icell);
866  }// end of IdToCell
867 
868  //......................................................................
869  ///
870  /// Given a unique cell identifier, look up the plane number.
871  ///
872  /// \param id : Cell identifier
873  /// \param iplane : optional, pointer to return plane number
874  ///
875  /// \returns Plane geometry description
876  /// \throws geo::Exception if id is not found
877  ///
878 
879  const PlaneGeo* GeometryBase::IdToPlane(const CellUniqueId& id, int* iplane) const {
880  *iplane = getPlaneCellMap(id).Plane();
881  return this->Plane(*iplane);
882  }// end of IdToCell
883  //......................................................................
885  return getPlaneCellMap(id).Plane();
886  }
887  //......................................................................
888  int GeometryBase::getCellID (const CellUniqueId& id) const{
889  return getPlaneCellMap(id).Cell();
890  }
891  //......................................................................
892  bool GeometryBase::getPlaneAndCellID(const CellUniqueId& id, int& plane, int& cell) const{
893  const OfflineChan offline_chan = getPlaneCellMap(id);
894  plane = (int)offline_chan.Plane();
895  cell = (int)offline_chan.Cell();
896  return true;
897  }
898  //......................................................................
900 
901  UIDMap::const_iterator itr = fIdMap.find(id);
902  if (itr == fIdMap.end()) {
903  throw cet::exception("BAD_CELL_UNIQUE_ID")
904  << "CellUniqueId " << id << "\n"
905  << __FILE__ << ":" << __LINE__ << "\n";
906  }
907 
908  return itr->second;
909  }
910 
911  //......................................................................
912  /// Return list of planes which measure the requested projection
913  ///
914  /// \param v : X/Y or All (see View_t typedef in PlaneGeo.h)
915  /// \returns : a set of plane indices satisfying the request
916  /// \throws : geo::Exception if initialization fails
917  ///
918  const std::set<unsigned int>& GeometryBase::GetPlanesByView(View_t v) const {
919  // Return the user's choice, fall back on all planes
920  if (v==kX) return fXplanes;
921  else if (v==kY) return fYplanes;
922  return fAllPlanes;
923  }// end of GetPlanesByView
924 
925  //......................................................................
926  /// Return the index of the next plane in a particular view
927  ///
928  /// \param p1 : Index of current plane
929  /// \param d : Direction to step in (>0 next down stream, <0 upstream)
930  ///
931  /// \returns Index of next plane in same view. If no plane can be
932  /// found (eg. reach end of detector) return kPLANE_NOT_FOUND.
933  ///
934  /// \throws : Asserts that p1 is in valid range
935  ///
936  const unsigned int GeometryBase::NextPlaneInView(unsigned int p1, int d) const {
937  assert(p1<fPlanes.size());
938 
939  // March along from p1 in the direction d until we find a match
940  int s = (d>0 ? 1 : -1);
941  if (s<0 && p1==0) return kPLANE_NOT_FOUND;
942  if (s>0 && p1==fPlanes.size()-1) return kPLANE_NOT_FOUND;
943  for (unsigned int p2=p1+s; 1; p2+=s) {
944  if (fPlanes[p1]->View() == fPlanes[p2]->View()) return p2;
945  if (p2==fPlanes.size()-1) return kPLANE_NOT_FOUND;
946  if (p2==0) return kPLANE_NOT_FOUND;
947  }
948  }// end of NextPlaneInView
949 
950  //......................................................................
951  /// Return the index of the next plane in a particular view
952  ///
953  /// \param p1 : Index of current plane
954  /// \param d : Direction to step in (>0 next down stream, <0 upstream)
955  ///
956  /// \returns index of next plane in same view.
957  ///
958  /// \throws : Asserts that p1 is in valid range
959  ///
960  const unsigned int GeometryBase::NextPlaneOtherView(unsigned int p1, int d) const
961  {
962  assert(p1<fPlanes.size());
963 
964  // March along from p1 in the direction d until we find a match
965  int s = (d>0 ? 1 : -1);
966  if (s<0 && p1==0) return kPLANE_NOT_FOUND;
967  if (s>0 && p1==fPlanes.size()-1) return kPLANE_NOT_FOUND;
968  for (unsigned int p2=p1+s; 1; p2+=s) {
969  if (fPlanes[p1]->View() != fPlanes[p2]->View()) return p2;
970  if (p2==fPlanes.size()) return kPLANE_NOT_FOUND;
971  if (p2==0) return kPLANE_NOT_FOUND;
972  }
973  }// end of NextPlaneOtherView
974 
975  //......................................................................
976  /// \brief Returns the index of the first plane contained in the muon catcher
977  ///
978  /// Or kPLANE_NOT_FOUND if there is no muon catcher
979  const unsigned int GeometryBase::FirstPlaneInMuonCatcher() const
980  {
981  if(fPlanesInMuonCatcher.empty()){
982  mf::LogWarning("GeometryBase") << "Requesting muon catcher info for a detector without one";
983  return kPLANE_NOT_FOUND;
984  }
985 
986  auto it = std::find(fPlanes.begin(), fPlanes.end(), fPlanesInMuonCatcher[0]);
987 
988  assert(it != fPlanes.end());
989  return it-fPlanes.begin();
990  }
991 
992  //......................................................................
993  ///
994  /// Build maps used for quick access to the geometry
995  ///
997  {
998  fAllPlanes.clear();
999  fXplanes.clear();
1000  fYplanes.clear();
1001 
1002  // Sets which list planes by view
1003  for (unsigned int i=0; i<fPlanes.size(); ++i) {
1004  fAllPlanes.insert(i);
1005  if (fPlanes[i]->View()==kX) fXplanes.insert(i);
1006  else if (fPlanes[i]->View()==kY) fYplanes.insert(i);
1007  else {
1008  throw cet::exception("BAD_GEO_CONFIG")
1009  << "Bad plane view" << fPlanes[i]->View() << ", plane " << i << "\n"
1010  << __FILE__ << ":" << __LINE__ << "\n";
1011  }
1012  }
1013 
1014  // Unique Id -> cell map
1015  for (unsigned int i=0; i<fPlanes.size(); ++i) {
1016  const PlaneGeo* p = fPlanes[i];
1017  for (unsigned int j=0; j<p->Ncells(); ++j) {
1018  // Assert that the cell ids are in fact unique
1019  if (fIdMap.find(this->Plane(i)->Cell(j)->Id())!=fIdMap.end()) {
1020  throw cet::exception("BAD_CELL_UNIQUE_ID")
1021  << "Duplicate Cell UIDs" << this->Plane(i)->Cell(j)->Id()
1022  << ", plane " << i << ", cell " << j << "\n"
1023  << __FILE__ << ":" << __LINE__ << "\n";
1024  }
1025  OfflineChan c(i,j);
1026  fIdMap[this->Plane(i)->Cell(j)->Id()] = c;
1027  }
1028  }
1029  }// end of BuildMaps
1030 
1031  //......................................................................
1032  ///
1033  /// Return the total mass of the detector
1034  ///
1035  ///
1036  double GeometryBase::TotalMass(const char *volume) const
1037  {
1038  double mass = 0.;
1039 
1040  // the detector resides in the vDetEnclosure volume so start there
1041  const TGeoVolume *enclosure = gGeoManager->FindVolumeFast(volume);
1042 
1043  int old_verbosity = gGeoManager->GetVerboseLevel();
1044  gGeoManager->SetVerboseLevel(0);
1045  // loop over the daughters in the volume and let ROOT calculate the
1046  // mass in kg for you for the daughter nodes that have vBlock in
1047  // their names this only goes one level down, so no worries of
1048  // double counting non-superblock blocks
1049  //
1050  for(int d = 0; d < enclosure->GetNdaughters(); ++d){
1051  const TGeoNode* node = enclosure->GetNode(d);
1052  const char* nm = node->GetName();
1053  if(strncmp(nm,"vSuperBlock",11) == 0){
1054  std::string name = nm;
1055  name.erase(name.find('_'));
1056  mass += gGeoManager->FindVolumeFast(name.c_str())->Weight();
1057  }
1058  }
1059  gGeoManager->SetVerboseLevel(old_verbosity);
1060  return mass;
1061  }// end of TotalMass
1062 
1063 
1064  //---------------------------------------------------------------------------
1065  ///
1066  /// Method to calculate masses of the detector
1067  ///
1068  bool GeometryBase::calculateMassesLong(const unsigned int number_of_points,
1069  CLHEP::HepRandomEngine& engine) const {
1070 
1071  std::map <std::string, double> density; // Map of densities of all materials in the detector. Dimensionless. Mapped by material names.
1072  std::map <std::string, double> volume; // Map of volumes of all materials in the detector. Dimensionless. Mapped by material names.
1073  std::map <std::string, double> mass; // Map of total masses of all materials the detector. Dimensionless. Mapped by material names.
1074  std::map <std::string, double> volume_fiducial; // Map of volumes of all materials in the fiducial region of the detector. Dimensionless. Mapped by material names.
1075  std::map <std::string, double> mass_fiducial; // Map of total masses of all materials in the fiducial region of the detector. Dimensionless. Mapped by material names.
1076 
1077  CLHEP::RandFlat flat(engine);
1078 
1079  TStopwatch timer;
1080  timer.Start();
1081 
1082  // Effective cell volume for a point
1083  const double cell_volume = (2.0*fDetHalfWidth*CLHEP::cm)*(2.0*fDetHalfHeight*CLHEP::cm)*(fDetLength*CLHEP::cm) / double(number_of_points);
1084 
1085  // Loop over points
1086  for(unsigned int ipoint = 0; ipoint < number_of_points; ++ipoint){
1087 
1088  // Set the progress bar
1089  ROOTGeoManager()->GetGeomPainter()->OpProgress("Calculating masses", ipoint, number_of_points-1, &timer, kFALSE);
1090 
1091  // X, Y, Z of the point in cm
1092  const double x = flat.fire(-fDetHalfWidth , fDetHalfWidth ) * CLHEP::cm;
1093  const double y = flat.fire(-fDetHalfHeight, fDetHalfHeight) * CLHEP::cm;
1094  const double z = flat.fire(0.0 , fDetLength ) * CLHEP::cm;
1095  const double current_position[3] = {x,y,z};
1096 
1097  // Get the material of the point
1098  const TGeoMaterial* current_material = Material(x/CLHEP::cm,
1099  y/CLHEP::cm,
1100  z/CLHEP::cm
1101  );
1102  // Get the material name
1103  const std::string current_material_name = current_material->GetName();
1104 
1105  // See if we have not defined the material
1106  if(density.find(current_material_name) == density.end()){
1107  density [current_material_name] = current_material->GetDensity() * (CLHEP::g/CLHEP::cm3); //http://indra.in2p3.fr/KaliVedaDoc/1.8.1/src/KVedaLossMaterial.h.html
1108  // Initialize these to zero
1109  volume [current_material_name] = 0.0;
1110  volume_fiducial[current_material_name] = 0.0;
1111  mass [current_material_name] = 0.0;
1112  mass_fiducial [current_material_name] = 0.0;
1113  }// end of searching current_material_name
1114 
1115  volume[current_material_name] += cell_volume;
1116  mass [current_material_name] += cell_volume * density[current_material_name];
1117 
1118  // See if the point is inside fiducial volume
1119  if(isInsideFiducialVolume(current_position)){
1120  volume_fiducial[current_material_name] += cell_volume;
1121  mass_fiducial [current_material_name] += cell_volume * density[current_material_name];
1122  }// end of checking if the point is inside fiducial volume
1123 
1124  }// end of looping over points
1125 
1126  std::cout<<"Calculated masses, volumes etc.:\n";
1127 
1128  // Printing the mass table for each detector material
1129  for(std::map<std::string, double>::iterator iter = density.begin(); iter != density.end(); ++iter) {
1130 
1131  const std::string key_str = iter->first;
1132 
1133  std::cout<<std::setw(9)<<"\t\tMaterial="<<std::setw(20)<<key_str
1134  <<std::setw(10)<<" density="<<std::setw(10)<<density[key_str]/CLHEP::kg*CLHEP::m3<<std::setw(6)<<"kg/m3"
1135  <<std::setw(9)<<" volume="<<std::setw(10)<<volume[key_str]/CLHEP::m3<<std::setw(2)<<"m3"
1136  <<std::setw(7)<<" mass="<<std::setw(10)<<mass[key_str] / CLHEP::kg<<std::setw(2)<<"kg"
1137  <<std::setw(17)<<" volume_fiducial="<<std::setw(10)<<volume_fiducial[key_str] / CLHEP::m3<<std::setw(1)<<"m3"
1138  <<std::setw(15)<<" mass_fiducial="<<std::setw(10)<<mass_fiducial[key_str] / CLHEP::kg<<std::setw(3)<<"kg\n";
1139 
1140  }// end of printing the mass table for each detector material
1141 
1142  return true;
1143  }
1144 
1145  //......................................................................
1146  ///
1147  /// Return the column density between 2 points in gm/cm^2
1148  ///
1149  /// \param p1 : pointer to array holding xyz of first point in world
1150  /// coordinates
1151  ///
1152  /// \param p2 : pointer to array holding xyz of second point in world
1153  /// coordinates
1154  ///
1155  double GeometryBase::MassBetweenPoints(double *p1, double *p2) const
1156  {
1157  // The purpose of this method is to determine the column density
1158  // between the two points given. Do that by starting at p1 and
1159  // stepping until you get to the node of p2. calculate the distance
1160  // between the point just inside that node and p2 to get the last
1161  // bit of column density
1162  double columnD = 0.;
1163 
1164  // first initialize a track - get the direction cosines
1165  double length = util::pythag(p2[0]-p1[0],
1166  p2[1]-p1[1],
1167  p2[2]-p1[2]);
1168  double dxyz[3] = {(p2[0]-p1[0])/length,
1169  (p2[1]-p1[1])/length,
1170  (p2[2]-p1[2])/length};
1171 
1172  gGeoManager->InitTrack(p1,dxyz);
1173  gGeoManager->FindNextBoundary(1.E-6);
1174 
1175  // might be helpful to have a point to a TGeoNode
1176  TGeoNode *node = gGeoManager->GetCurrentNode();
1177 
1178  // check that the points are not in the same volume already.
1179  // if they are in different volumes, keep stepping until you
1180  // are in the same volume as the second point
1181  while(!gGeoManager->IsSameLocation(p2[0], p2[1], p2[2])){
1182  gGeoManager->FindNextBoundary();
1183 
1184  columnD +=
1185  gGeoManager->GetStep() *
1186  node->GetMedium()->GetMaterial()->GetDensity();
1187 
1188  // the act of stepping puts you in the next node and returns that
1189  // node
1190  node = gGeoManager->Step();
1191  } //end loop to get to volume of second point
1192 
1193  // now you are in the same volume as the last point, but not at that
1194  // point. get the distance between the current point and the last
1195  // one
1196  const double *current = gGeoManager->GetCurrentPoint();
1197  length = util::pythag(p2[0]-current[0],
1198  p2[1]-current[1],
1199  p2[2]-current[2]);
1200  columnD += length*node->GetMedium()->GetMaterial()->GetDensity();
1201 
1202  return columnD;
1203  }// end of MassBetweenPoints
1204 
1205  //......................................................................
1206 
1208  (const double *p1,
1209  const double *p2,
1210  std::vector<double>& ds,
1211  std::vector<const TGeoMaterial*>& m) const
1212  {
1213  double d[3], dmag;
1214  dmag = util::pythag(p2[0]-p1[0],p2[1]-p1[1],p2[2]-p1[2]);
1215  d[0] = (p2[0]-p1[0])/dmag;
1216  d[1] = (p2[1]-p1[1])/dmag;
1217  d[2] = (p2[2]-p1[2])/dmag;
1218 
1219  gGeoManager->InitTrack(p1,d);
1220  gGeoManager->FindNextBoundary(1.E-6);
1221  TGeoNode *node = gGeoManager->GetCurrentNode();
1222 
1223  while (!gGeoManager->IsSameLocation(p2[0], p2[1], p2[2])) {
1224  gGeoManager->FindNextBoundary();
1225 
1226  // If we take a step and haven't gone anywhere, then we will be
1227  // stuck in a very long loop. So we will quit if the next step
1228  // takes us nowhere.
1229  double step = gGeoManager->GetStep();
1230  if(step < 1.0e-12) break;
1231 
1232  ds.push_back(step);
1233  m.push_back(node->GetMedium()->GetMaterial());
1234 
1235  node = gGeoManager->Step();
1236  if (node==0) {
1237  std::cerr << __FILE__ << ":" << __LINE__
1238  << " Cannot continue tracking MaterialsBetweenPoints(). "
1239  << "Null node pointer." << std::endl;
1240  return;
1241  }
1242  }
1243  //
1244  // Handle final step by hand
1245  //
1246  const double* p3 = gGeoManager->GetCurrentPoint();
1247  ds.push_back(util::pythag(p3[0]-p2[0],p3[1]-p2[1],p3[2]-p2[2]));
1248  m.push_back(node->GetMedium()->GetMaterial());
1249  }
1250 
1251  ///
1252  ///Overloaded MaterialsBetweenPoints to for using TVector3 instead of
1253  ///arrays because rb::Track trajectory points are TVector3 and there
1254  ///is a desire to track a 'track' through the geometry. This will make
1255  ///that easier. Function is really just a wrapper around the array version
1256  ///
1258  (TVector3 v1,
1259  TVector3 v2,
1260  std::vector<double>& ds,
1261  std::vector<const TGeoMaterial*>& m) const
1262  {
1263  double p1[3], p2[3];
1264  v1.GetXYZ(p1);
1265  v2.GetXYZ(p2);
1266 
1267  MaterialsBetweenPoints(p1, p2, ds, m);
1268  }
1269 
1270  ///
1271  /// Calculate the detector size based on the geometry. Requires a
1272  /// little bit of work, so do this once when the geometry is updated
1273  ///
1275  {
1276  double dummy;
1277 
1278  // Compute the height and width based on the sizes of the vertical
1279  // and horizontal planes. Remember, in the local plane frame z goes
1280  // along the length of the modules. Hence the use of axis=3 to get
1281  // the detector height (from the vertical planes) and width (from
1282  // the horizontal planes).
1283  const TGeoVolume* vvp = gGeoManager->GetVolume("vPlaneV");
1284  const TGeoVolume* hvp = gGeoManager->GetVolume("vPlaneH");
1285 
1286  // Handle the case of the IPND where the volumes are names slightly
1287  // different accounting for the different widths of the modules.
1288  if (vvp==0) vvp = gGeoManager->GetVolume("vPlaneVNDOS");
1289  if (hvp==0) hvp = gGeoManager->GetVolume("vPlaneHNDOS");
1290  if (vvp==0 && hvp==0) {
1291  vvp = gGeoManager->GetVolume("vPlaneVND");
1292  hvp = gGeoManager->GetVolume("vPlaneHND");
1293  }
1294 
1295  /// If there no Vertical or horizontal planes, then we go get the size from the Detector Enclosure
1296  /// It used to be that once this condition is true, then it we would just throw an exception.
1297  /// However, now Giuseppe Ferone and Denis Perevalov are studying interactions in simple detector,
1298  /// which don't have planes in it. Thus, they don't want the code breaking.
1299  if (vvp==0 || hvp==0) {
1300  mf::LogWarning("SetDetectorSize")
1301  << "Warning! BAD_GEO_CONFIG. Unable to find shapes to set detector size\n"
1302  << "Now getting the detector size from vDetEnclosure\n";
1303 
1304  /// Getting Detector enclosure volume
1305  const TGeoVolume* vdet_enclosure = gGeoManager->GetVolume("vDetEnclosure");
1306 
1307  /// Lo and Hi ranges of each of the axes. To be filled later
1308  double x_lo, x_hi;
1309  double y_lo, y_hi;
1310  double z_lo, z_hi;
1311 
1312  /// Fill them
1313  vdet_enclosure->GetShape()->GetAxisRange(1, x_lo, x_hi);
1314  vdet_enclosure->GetShape()->GetAxisRange(2, y_lo, y_hi);
1315  vdet_enclosure->GetShape()->GetAxisRange(3, z_lo, z_hi);
1316 
1317  /// Set detector size
1318  fDetHalfWidth = x_hi;
1319  fDetHalfHeight = y_hi;
1320  fDetLength = z_hi - z_lo;
1321 
1322  return;
1323  }
1324 
1325  const TGeoShape* vp = vvp->GetShape();
1326  const TGeoShape* hp = hvp->GetShape();
1327  vp->GetAxisRange(3,dummy,fDetHalfHeight);
1328  hp->GetAxisRange(3,dummy,fDetHalfWidth);
1329 
1330  // Remember -- in the local plane coordinate frame, the "depth" is
1331  // along y, hence axis=2. In the world coordinate this is along z.
1332  double vplanez1, vplanez2;
1333  vp->GetAxisRange(2,vplanez1,vplanez2);
1334 
1335  // Compute the detector length based on the center positions of
1336  // cells in the first and last plane. Adjust to account for the
1337  // plane widths. This ajustment works if the vertical and horizontal
1338  // planes have the same depth, or if the detector begins and ends
1339  // with vertical planes. As of the TDR both of these conditions was
1340  // true for all the NOvA detectors.
1341  double xyz1[3] = {0,0,0};
1342  double xyz2[3] = {0,0,0};
1343  this->Plane(0)-> Cell(0)->GetCenter(xyz1);
1344  this->Plane(fPlanes.size()-1)->Cell(0)->GetCenter(xyz2);
1345  fDetLength = (xyz2[2]-xyz1[2])+(vplanez2-vplanez1);
1346  }// end of SetDetectorSize
1347 
1348  //......................................................................
1349 
1350  void GeometryBase::setDetectorBigBox(double detector_bigbox_range) {
1351  fIsDetectorBigBoxUsed = true;
1352 
1353  double x1, x2;
1354  double y1, y2;
1355  double z1, z2;
1356 
1357  /// Get the Detector size
1358  this->DetectorEnclosureBox(&x1, &x2, &y1, &y2, &z1, &z2);
1359 
1360  /// Dimensions of the DetectorBigBox in cm
1361  fDetectorBigBoxXHi = x2 + detector_bigbox_range;
1362  fDetectorBigBoxXLo = x1 - detector_bigbox_range;
1363  fDetectorBigBoxYHi = y2 + detector_bigbox_range;
1364  fDetectorBigBoxYLo = y1 - detector_bigbox_range;
1365  fDetectorBigBoxZHi = z2 + detector_bigbox_range;
1366  fDetectorBigBoxZLo = z1 - detector_bigbox_range;
1367 
1368  }// end of setDetectorBigBox
1369 
1370  //......................................................................
1371  bool GeometryBase::isInsideDetectorBigBox(const double x, const double y, const double z) const {
1372  if( x < fDetectorBigBoxXHi && x > fDetectorBigBoxXLo &&
1373  y < fDetectorBigBoxYHi && y > fDetectorBigBoxYLo &&
1374  z < fDetectorBigBoxZHi && z > fDetectorBigBoxZLo)
1375  return true;
1376  return false;
1377  }// end of GeometryBase::isInsideDetectorBigBox
1378 
1379  //......................................................................
1380  ///
1381  /// Is the particle inside the fiducial volume?
1382  ///
1383  /// \param x_cm : x coordinate in cm
1384  ///
1385  /// \param y_cm : y coordinate in cm
1386  ///
1387  /// \param z_cm : z coordinate in cm
1388  ///
1389  bool GeometryBase::isInsideFiducialVolume(const double x_cm, const double y_cm, const double z_cm) const {
1390 
1391  return (x_cm < fFiducialVolumeXHi && x_cm > fFiducialVolumeXLo &&
1392  y_cm < fFiducialVolumeYHi && y_cm > fFiducialVolumeYLo &&
1393  z_cm < fFiducialVolumeZHi && z_cm > fFiducialVolumeZLo);
1394  }// end of GeometryBase::isInsideFiducialVolume
1395 
1396  //......................................................................
1397  ///
1398  /// Is the particle inside the fiducial volume?
1399  ///
1400  /// \param xyz : pointer to array holding xyz in world coordinates
1401  ///
1402  bool GeometryBase::isInsideFiducialVolume(const double* xyz) const {
1403 
1404  return isInsideFiducialVolume(xyz[0]/CLHEP::cm,
1405  xyz[1]/CLHEP::cm,
1406  xyz[2]/CLHEP::cm);
1407  }// end of GeometryBase::isInsideFiducialVolume
1408 
1409  //......................................................................
1410  void GeometryBase::FiducialBox(TVector3& r0, TVector3& r1) const
1411  {
1414  }
1415 
1416  //----------------------------------------------------------------------
1417  ///
1418  /// Compute the distance from a point to the detector wall along a
1419  /// particular direction.
1420  ///
1421  /// \param x0 - Input, [x,y,z] start point (cm)
1422  /// \param dx - Input, direction vector
1423  /// \param exitp - Ouput, Calculated [x,y,z] of exit point (cm)
1424  ///
1425  /// \returns distance to exit point in cm
1426  ///
1427  double GeometryBase::DEdge(const double* x0,
1428  const double* dx,
1429  double* exitp) const
1430  {
1431  double xyzout[3];
1432  geo::ProjectToBoxEdge(x0, dx,
1435  0, fDetLength,
1436  xyzout);
1437  if (exitp!=0) {
1438  exitp[0] = xyzout[0];
1439  exitp[1] = xyzout[1];
1440  exitp[2] = xyzout[2];
1441  }
1442  return util::pythag(x0[0]-xyzout[0],
1443  x0[1]-xyzout[1],
1444  x0[2]-xyzout[2]);
1445  }// end of GeometryBase::DEdge
1446 
1447  //--------------------------------------------------------------------------
1448  /// \brief Find the distance of closest approach between a cell and a track
1449  ///
1450  /// @param ip - plane number
1451  /// @param ic - cell number
1452  /// @param vtx - track start point
1453  /// @param dir - track direction
1454  /// @param dx - offset within the cell (in cm, dx=0 is center of cell)
1455  /// @param dy - offset within the cell (in cm, dy=0 is center of cell)
1456  /// @param w - Distance from read-out end where the closest approach occurs (cm)
1457  /// @param s - Distance down the track where the closest approach occurs (cm)
1458  /// @param point_in_cell - Point inside the cell where the closest approach occurs
1459  /// @param point_on_track - Point on track where the closest approach occurs
1460  /// @returns The distance^2 of closest approach
1461  ///
1462  double GeometryBase::ClosestApproach(unsigned int ip,
1463  unsigned int ic,
1464  const TVector3& vtx,
1465  const TVector3& dir,
1466  double dx,
1467  double dy,
1468  double* w,
1469  double* s,
1470  TVector3* point_in_cell,
1471  TVector3* point_on_track) const
1472  {
1473  // ==
1474  // Find the line that runs down the cell at the request position
1475  //
1476  const PlaneGeo* plane = this->Plane(ip);
1477  View_t view = plane->View();;
1478 
1479  const CellGeo* cell = plane->Cell(ic);
1480 
1481  double p0[3], p1[3];
1482  double dz = cell->HalfL();
1483  cell->GetCenter(p0, +dz);
1484  cell->GetCenter(p1, -dz);
1485 
1486  // Shift to the ray requested
1487  p0[2] += dy;
1488  p1[2] += dy;
1489  if (view==kX) { p0[0] += dx; p1[0] += dx; }
1490  else if (view==kY) { p0[1] += dx; p1[1] += dx; }
1491 
1492  TVector3 P0(p0[0], p0[1], p0[2]);
1493  TVector3 P1(p1[0], p1[1], p1[2]);
1494  // ==
1495 
1496  // Get two points on the track, the vertex and a point at vertex+dir
1497  TVector3 Q1(vtx); Q1 += dir;
1498 
1499  // Use the closest approach function to compute the rest
1500  return geo::ClosestApproach(P0, P1,
1501  vtx, Q1,
1502  w, s,
1503  point_in_cell,
1504  point_on_track);
1505  }
1506 
1507  //--------------------------------------------------------------------------
1508  TGeoManager* GeometryBase::ROOTGeoManager() const { return gGeoManager;}
1509 
1510 
1511  //--------------------------------------------------------------------------
1512  ///
1513  /// Method to set DetectorID
1514  ///
1516  {
1517  LOG_DEBUG("Geometry")
1518  <<"GeometryBase::setDetectorID\n";
1519 
1520  fDetId = det_id;
1521 
1522  // Change coordinate transformation class as well
1524  //fCoordinateTransformation.print();
1525 
1526  // Work out which config to use
1528  switch(fDetId){
1529  case novadaq::cnv::kNDOS: config = fParams.ndos(); break;
1530  case novadaq::cnv::kNEARDET: config = fParams.nd(); break;
1531  case novadaq::cnv::kFARDET: config = fParams.fd(); break;
1532  case novadaq::cnv::kTESTBEAM: config = fParams.tb(); break;
1533  default: return; // Probably unknown or testbeam. Nothing sane to set the
1534  // below variables to.
1535  }
1536 
1537  // Set fiducial volume
1544 
1545  fIsDetectorBigBoxUsed = config.BigBoxUsed();
1546  fBigBoxRange = config.BigBoxRange();
1547 
1548  }// end of GeometryBase::setDetectorID
1549 
1550 
1551  //--------------------------------------------------------------------------
1552  ///
1553  /// Method to count the number of cells passed through by a line segment
1554  /// between two points in the detector.
1555  ///
1556  void GeometryBase::CountCellsOnLine(double x1, double y1, double z1, double x2, double y2, double z2,
1557  std::vector<OfflineChan>& xhitsonline,
1558  std::vector<OfflineChan>& yhitsonline)
1559  {
1560  std::vector<CellOnLine> xhs, yhs;
1561  CountCellsOnLine(x1, y1, z1, x2, y2, z2, xhs, yhs);
1562 
1563  for(const CellOnLine& c: xhs) xhitsonline.push_back(c.chan);
1564  for(const CellOnLine& c: yhs) yhitsonline.push_back(c.chan);
1565  } // end CountCellsOnLine
1566 
1567  //---------------------------------------------------------------------------
1568  ///
1569  /// Method to count the number of cells passed through by a line segment
1570  /// between two points in the detector.
1571  ///
1572  void GeometryBase::CountCellsOnLine(TVector3 r1, TVector3 r2,
1573  std::vector<OfflineChan>& xhitsonline,
1574  std::vector<OfflineChan>& yhitsonline)
1575  {
1576  CountCellsOnLine(r1.X(), r1.Y(), r1.Z(), r2.X(), r2.Y(), r2.Z(),
1577  xhitsonline, yhitsonline);
1578  }
1579 
1580  //---------------------------------------------------------------------------
1581  ///
1582  /// Method to count the number of cells passed through by a line segment
1583  /// between two points in the detector and keep track of the entry and exit
1584  /// points of the line on the cells
1585  ///
1586  void GeometryBase::CountCellsOnLine(double x1, double y1, double z1, double x2, double y2, double z2,
1587  std::vector<CellOnLine>& xhitsonline, std::vector<CellOnLine>& yhitsonline)
1588  {
1589 
1590  // This function returns a vector of cells crossed by the line and the
1591  // entry and exit points of the line through the cells for each view. If
1592  // the end point is within a cell, the entry or exit point is set to be
1593  // the end of the line.
1594  //
1595  // WARNING: This function makes use of gGeoManager and is therefore
1596  // accessing global variables! Take caution to properly
1597  // initialize these global variables when accessing them after
1598  // using this function!
1599 
1600  // define initial position and normalized direction vector
1601  double xyz[3] = {x1, y1, z1};
1602  double d = util::pythag(x2-x1, y2-y1, z2-z1);
1603  double dxyz[3] = {(x2-x1)/d, (y2-y1)/d, (z2-z1)/d};
1604 
1605  gGeoManager->InitTrack(xyz, dxyz);
1606  const double* pos = gGeoManager->GetCurrentPoint();
1607 
1608  double dtemp = util::pythag(pos[0]-x1, pos[1]-y1, pos[2]-z1);
1609 
1610  bool enter = false;
1611 
1612  // We update the properties of this cell intersection as we go, and store
1613  // it one every exit.
1614  CellOnLine col;
1615 
1616  // create a list of positions in scintillator
1617  while(dtemp <= d) {
1618  // if the current material is scintillator, count that as a cell passed through
1619  if(strncmp("Scintillator", gGeoManager->GetCurrentNode()->GetVolume()->GetMaterial()->GetName(), 12) == 0) {
1620  enter = true;
1621 
1622  col.entry = TVector3(pos[0],pos[1],pos[2]);
1623 
1624  int plane, cell;
1625  IdToCell(CurrentCellId(), &plane, &cell);
1626 
1627  col.chan = OfflineChan(plane, cell);
1628  }
1629 
1630  // step forward to the next material boundary
1631  gGeoManager->FindNextBoundary();
1632  gGeoManager->Step(kTRUE,kTRUE);
1633 
1634  dtemp = util::pythag(pos[0]-x1, pos[1]-y1, pos[2]-z1);
1635 
1636  if(enter){
1637  enter = false;
1638 
1639  if(dtemp <= d){
1640  col.exit = TVector3(pos[0], pos[1], pos[2]);
1641  }
1642  else{
1643  col.exit = TVector3(x2, y2, z2);
1644  }
1645 
1646  const View_t view = Plane(col.chan.Plane())->View();
1647  if(view == geo::kX){ xhitsonline.push_back(col); }
1648  if(view == geo::kY){ yhitsonline.push_back(col); }
1649 
1650  }
1651  }
1652 
1653  } // end CountCellsOnLine
1654 
1655  //---------------------------------------------------------------------------
1656  ///
1657  /// Method to count the number of cells passed through by a line segment
1658  /// between two points in the detector and keep track of the entry and exit points
1659  /// of the line on the cells
1660  ///
1661  void GeometryBase::CountCellsOnLine(TVector3 r1, TVector3 r2,
1662  std::vector<CellOnLine>& Xhitsonline,
1663  std::vector<CellOnLine>& Yhitsonline)
1664  {
1665  CountCellsOnLine(r1.X(), r1.Y(), r1.Z(), r2.X(), r2.Y(), r2.Z(),
1666  Xhitsonline, Yhitsonline);
1667  }
1668 
1669  //---------------------------------------------------------------------------
1670  ///
1671  /// Fast method to count the number of cells passed through by a line
1672  /// segment between two points in the detector.
1673  ///
1675  double v1, double z1,
1676  double v2, double z2,
1677  std::vector<OfflineChan>& hitsonline)
1678  {
1679  assert(view == geo::kX || view == geo::kY);
1680 
1681  std::vector<OfflineChan> junk;
1682  if(view == geo::kX)
1683  CountCellsOnLineFast(v1, 0, z1, v2, 0, z2, hitsonline, junk);
1684  else
1685  CountCellsOnLineFast(0, v1, z1, 0, v2, z2, junk, hitsonline);
1686  }
1687 
1688  //---------------------------------------------------------------------------
1689  ///
1690  /// Fast method to count the number of cells passed through by a line
1691  /// segment between two points in the detector.
1692  ///
1693  void GeometryBase::CountCellsOnLineFast(double x1, double y1, double z1,
1694  double x2, double y2, double z2,
1695  std::vector<OfflineChan>& xhitsonline,
1696  std::vector<OfflineChan>& yhitsonline)
1697  {
1698  CountCellsOnLineFast(TVector3(x1, y1, z1), TVector3(x2, y2, z2),
1699  xhitsonline, yhitsonline);
1700  }
1701 
1702  //---------------------------------------------------------------------------
1703  ///
1704  /// Fast method to count the number of cells passed through by a line
1705  /// segment between two points in the detector.
1706  ///
1707  void GeometryBase::CountCellsOnLineFast(TVector3 r1, TVector3 r2,
1708  std::vector<OfflineChan>& xhitsonline,
1709  std::vector<OfflineChan>& yhitsonline)
1710  {
1711  // Otherwise we would count when the ray was outside the detector in the
1712  // other view.
1713  ClampRayToDetectorVolume(&r1, &r2, this);
1714 
1715  // We assume downstream-going below
1716  if(r1.Z() > r2.Z()) std::swap(r1, r2);
1717 
1718  const double z1 = r1.Z();
1719  const double z2 = r2.Z();
1720 
1721  const unsigned int planeMax = fPlanes.size();
1722  for(unsigned int plane = 0; plane < planeMax; ++plane){
1723  const geo::PlaneGeo* geoplane = fPlanes[plane];
1724  const geo::CellGeo* geocell0 = geoplane->Cell(0);
1725  // Assuming here perfect alignment and identical cells across plane
1726  double xyzp[3];
1727  geocell0->GetCenter(xyzp);
1728  const double zc = xyzp[2];
1729  const double dz = geocell0->HalfD();
1730 
1731  if(z1 > zc+dz) continue;
1732  // Plane positions increase monotonically
1733  if(z2 < zc-dz) break;
1734 
1735  const geo::View_t view = geoplane->View();
1736 
1737  const double v1 = r1[view];
1738  const double v2 = r2[view];
1739 
1740  const double denom = z2-z1;
1741  const double dvdz = denom ? (v2-v1)/denom : 0;
1742 
1743  // Position entering the cell or starting
1744  const double zin = std::max(zc-dz, z1);
1745  double vlo = v1+(zin-z1)*dvdz;
1746  // Position exiting or ending
1747  const double zout = std::min(zc+dz, z2);
1748  double vhi = v1+(zout-z1)*dvdz;
1749  // Ensure they're from low to high instead
1750  if(vlo > vhi) std::swap(vlo, vhi);
1751 
1752  for(unsigned int cell = 0; cell < geoplane->Ncells(); ++cell){
1753  double xyz[3];
1754  geoplane->Cell(cell)->GetCenter(xyz);
1755  const double vc = xyz[view];
1756  const double dv = geoplane->Cell(cell)->HalfW();
1757 
1758  if(vlo > vc+dv) continue;
1759  // Cell positions increase monotonically
1760  if(vhi < vc-dv) break;
1761 
1762  if(view == geo::kX) xhitsonline.push_back(geo::OfflineChan(plane, cell));
1763  if(view == geo::kY) yhitsonline.push_back(geo::OfflineChan(plane, cell));
1764  } // end for cell
1765  } // end for plane
1766  }
1767 
1768  //......................................................................
1769  /// Return distance from a point in the cell to readout electronics, including pigtail
1770  /// \param localz : Z position in cm in local cell frame
1771  double GeometryBase::DistToElectronics(double localz, const CellGeo& cell) const
1772  {
1773  return (cell.HalfL() - localz + this->GetPigtail(cell.Id()));
1774  }
1775 
1776  //......................................................................
1777  double GeometryBase::GetPigtail(const CellUniqueId& id) const
1778  {
1779  int planeid,cellid;
1780  this->IdToCell(id, &planeid, &cellid);
1781  View_t view = (this->IdToPlane(id,&planeid))->View();
1782  //NOTE: Do we think 32 cells per module can ever change
1783  //This number can be found in DAQChannelMap but did not want to add that
1784  //dependency to geometry unless needed.
1785  const int kCellsPerModule = 32;
1786  cellid = cellid % kCellsPerModule;
1787  // In vertical planes, high cell numbers have longer fibres.
1788  // In horizontal planes it's the opposite, so flip it round.
1789  // For Test Beam, vertical planes actually follow the same pattern as
1790  // the horizontals, so flip both.
1791  if( view == geo::kY || this->DetId() == novadaq::cnv::kTESTBEAM ) cellid = kCellsPerModule-cellid-1;
1792  // This should really never happen, but just to be safe...
1793  if(cellid < 0 || cellid >= kCellsPerModule) return 100;
1794  // Email from Tom Chase 2011-04-29
1795  // NB: this isn't just a ~3.8cm change per cell. Every 8 cells something
1796  // different happens.
1797  const double kPigtails[kCellsPerModule] = {
1798  34.5738, 38.4379, 42.3020, 46.1660, 50.0301, 53.8942, 57.7582, 61.6223,
1799  64.7504, 68.6144, 72.4785, 76.3426, 80.2067, 84.0707, 87.9348, 91.0790,
1800  95.3301, 99.1941, 103.058, 106.922, 110.786, 114.650, 118.514, 122.379,
1801  125.507, 129.371, 133.235, 137.099, 140.963, 144.827, 148.691, 150.751
1802  };
1803  return kPigtails[cellid];
1804  }
1805 
1806 } // end namespace geo
1807 ////////////////////////////////////////////////////////////////////////
novadaq::cnv::DetId fDetId
id for the detector being used
Definition: GeometryBase.h:285
double fFiducialVolumeZLo
Definition: GeometryBase.h:328
virtual ~GeometryBase()
T max(const caf::Proxy< T > &a, T b)
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
const XML_Char * name
Definition: expat.h:151
double HalfL() const
Definition: CellGeo.cxx:198
double CellTpos(unsigned int ip, unsigned int ic, double w=0.0) const
fileName
Definition: plotROC.py:78
double HalfD() const
Definition: CellGeo.cxx:205
void GetCenter(double *xyz, double localz=0.0) const
Definition: CellGeo.cxx:159
TVector3 BoosterBeamDirection() const
Direction of neutrinos from the Booster beam (unit vector)
Atom< bool > ForceUseFCLOnly
Definition: GeometryBase.h:71
bool getPlaneAndCellID(const CellUniqueId &id, int &plane, int &cell) const
bool IntersectsBox(const double xyz[], const double dxyz[], double xlo, double xhi, double ylo, double yhi, double zlo, double zhi)
Determine if a particle starting at xyz with direction dxyz will intersect a box defined by xlo...
Definition: Geo.cxx:97
set< int >::iterator it
double fDetectorBigBoxXHi
Dimensions of the DetectorBigBox in cm.
Definition: GeometryBase.h:314
static constexpr double cm3
TGeoMaterial * Material(double x, double y, double z) const
GeometryBase(const Params &params)
double fDetLength
Detector length (cm)
Definition: GeometryBase.h:299
bool setDetectorAndBeam(uint32_t detector_id, uint32_t beam_type)
Try to set the detector and beam.
Float_t y1[n_points_granero]
Definition: compare.C:5
std::string fROOTFile
root file holding the geometry
Definition: GeometryBase.h:282
static constexpr double g
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
void MakePlane(std::vector< const TGeoNode * > &n, unsigned int depth, bool inMuonCatcher)
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
const PlaneGeo * IdToPlane(const CellUniqueId &id, int *iplane) const
Float_t x1[n_points_granero]
Definition: compare.C:5
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
const char * p
Definition: xmltok.h:285
bool calculateMassesLong(const unsigned int number_of_points, CLHEP::HepRandomEngine &engine) const
const CellUniqueId CurrentCellId() const
Vertical planes which measure X.
Definition: PlaneGeo.h:28
unsigned int Ncells() const
Number of cells in this plane.
Definition: PlaneGeo.h:43
std::string fGDMLFile
gdml file holding the geometry
Definition: GeometryBase.h:281
Definition: scint.py:1
static constexpr Double_t nm
Definition: Munits.h:133
Table< DetectorParams > tb
Definition: GeometryBase.h:77
Table< DetectorParams > nd
Definition: GeometryBase.h:75
double HalfW() const
Definition: CellGeo.cxx:191
bool fIsDetectorBigBoxUsed
Do we need to use the BigBox cut?
Definition: GeometryBase.h:310
double DetLength() const
TGeoManager * ROOTGeoManager() const
std::set< unsigned int > fAllPlanes
List of all planes.
Definition: GeometryBase.h:295
double fFiducialVolumeXLo
Definition: GeometryBase.h:324
OStream cerr
Definition: OStream.cxx:7
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
Table< DetectorParams > ndos
Definition: GeometryBase.h:74
double MassBetweenPoints(double *p1, double *p2) const
TString ip
Definition: loadincs.C:5
Float_t tmp
Definition: plot.C:36
bool ClampRayToDetectorVolume(TVector3 *p0, TVector3 *p1, const GeometryBase *geom)
If either endpoint is outside the detector move it to the edge.
Definition: Geo.cxx:640
const PlaneGeo * Plane(unsigned int i) const
double fDetHalfHeight
Detector 1/2 height (cm)
Definition: GeometryBase.h:300
std::vector< PlaneGeo * > fPlanes
The detector planes.
Definition: GeometryBase.h:291
void MaterialsBetweenPoints(const double *p1, const double *p2, std::vector< double > &ds, std::vector< const TGeoMaterial * > &mat) const
double fDetHalfWidth
Detector 1/2 width (cm)
Definition: GeometryBase.h:301
std::string find_file(std::string const &filename) const
TVector3 exit
Exit point from cell.
Definition: GeometryBase.h:45
static constexpr double kg
CellUniqueId NodesToUniqueId(const std::vector< const TGeoNode * > &n, unsigned int depth)
double fFiducialVolumeXHi
Definition: GeometryBase.h:323
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
void CellInfo(unsigned int ip, unsigned int ic, View_t *view=0, double *pos=0, double *dpos=0) const
void CountCellsOnLineFast(double x1, double y1, double z1, double x2, double y2, double z2, std::vector< OfflineChan > &Xhitsonline, std::vector< OfflineChan > &Yhitsonline)
Make list of cells in each view traversed by a line segment.
int CurrentCell(int *ip, int *ic) const
const CellUniqueId & Id() const
Definition: CellGeo.h:55
void swap(art::HLTGlobalStatus &lhs, art::HLTGlobalStatus &rhs)
fclose(fg1)
const XML_Char * s
Definition: expat.h:262
static constexpr double cm
Definition: SystemOfUnits.h:99
void CountCellsOnLine(double X1, double Y1, double Z1, double X2, double Y2, double Z2, std::vector< OfflineChan > &Xhitsonline, std::vector< OfflineChan > &Yhitsonline)
double DEdge(const double *x0, const double *dx, double *exitp=0) const
View_t View() const
Which coordinate does this plane measure.
Definition: PlaneGeo.h:53
Far Detector at Ash River, MN.
double dy[NP][NC]
length
Definition: demo0.py:21
double dx[NP][NC]
double dz[NP][NC]
Float_t E
Definition: plot.C:20
Encapsulate the geometry of one entire detector (near, far, ndos)
void WorldBox(double *xlo_cm, double *xhi_cm, double *ylo_cm, double *yhi_cm, double *zlo_cm, double *zhi_cm) const
Int_t col[ntarg]
Definition: Style.C:29
Prototype Near Detector on the surface at FNAL.
TVector3 entry
Entry point into cell.
Definition: GeometryBase.h:44
bool LoadGeometryFile(std::string gdmlfile, novadaq::cnv::DetId det_id=novadaq::cnv::kUNKNOWN_DET)
In case of unknown ID, guesses from the filename.
std::set< unsigned int > fYplanes
List of Y measuring planes.
Definition: GeometryBase.h:297
UIDMap fIdMap
Unique ID -> Plane,Cell.
Definition: GeometryBase.h:293
void setDetectorBigBox(double detector_bigbox_range)
int getPlaneID(const CellUniqueId &id) 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
double fDetectorBigBoxZHi
Definition: GeometryBase.h:318
static std::string FindGDMLFile(std::string fname)
Search for Geometry/gdml/fname in FW_SEARCH_PATH, return full path.
Near Detector in the NuMI cavern.
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
bool isInsideDetectorBigBox(const double x_cm, const double y_cm, const double z_cm) const
Is the particle inside the detector Big Box?
Collect Geo headers and supply basic geometry functions.
Float_t d
Definition: plot.C:236
CoordinateTransformation fCoordinateTransformation
Coordinate Transformation class.
Definition: GeometryBase.h:289
double fDetectorBigBoxXLo
Definition: GeometryBase.h:315
double fDetectorBigBoxYLo
Definition: GeometryBase.h:317
const double j
Definition: BetheBloch.cxx:29
double fFiducialVolumeZHi
Definition: GeometryBase.h:327
double DetHalfHeight() const
double fBigBoxRange
Range of big box.
Definition: GeometryBase.h:311
static constexpr double m3
Atom< std::string > StoreTempGeo
Definition: GeometryBase.h:72
static const double ns
Module that plots metrics from reconstructed cosmic ray data.
static bool plane_sort(const PlaneGeo *p1, const PlaneGeo *p2)
Definition: View.py:1
z
Definition: test.py:28
A very simple service to remember what detector we&#39;re working in.
void FindPlanes(std::vector< const TGeoNode * > &n, unsigned int depth, bool inMuonCatcher=false)
OStream cout
Definition: OStream.cxx:6
unsigned short Plane() const
Definition: OfflineChan.h:31
std::vector< PlaneGeo * > fPlanesInMuonCatcher
Same pointers as fPlanes.
Definition: GeometryBase.h:292
double GetPigtail(const CellUniqueId &id) const
Return length of fiber in cm from end of cell to apd.
const std::string path
Definition: plot_BEN.C:43
TVector3 NuMIBeamDirection() const
Direction of neutrinos from the NuMI beam (unit vector)
unsigned long long int CellUniqueId
Definition: CellUniqueId.h:15
std::set< unsigned int > fXplanes
List of X measuring planes.
Definition: GeometryBase.h:296
double fDetectorBigBoxZLo
Definition: GeometryBase.h:319
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
OfflineChan chan
Channel passed through by line.
Definition: GeometryBase.h:43
double TotalMass(const char *volume="vDetEnclosure") const
static void clear()
Definition: CellGeo.h:57
unsigned short Cell() const
Definition: OfflineChan.h:32
int MakeTmpFile(std::string gdmlInfo)
double DetHalfWidth() const
Table< DetectorParams > fd
Definition: GeometryBase.h:76
double ClosestApproach(const double point[], const double intercept[], const double slopes[], double closest[])
Find the distance of closest approach between point and line.
Definition: Geo.cxx:222
double DistToElectronics(double localz, const CellGeo &cell) const
get distance from local z position in cell to apd in cm, including pigtail
TDirectory * dir
Definition: macro.C:5
double SurfaceY() const
double pythag(double x, double y)
2D Euclidean distance
Definition: MathUtil.h:29
void FiducialBox(TVector3 &r0, TVector3 &r1) const
const unsigned int NextPlaneInView(unsigned int p, int d=+1) const
A (plane, cell) pair.
Definition: OfflineChan.h:17
double ClosestApproach(unsigned int ip, unsigned int ic, const TVector3 &vtx, const TVector3 &dir, double dx=0.0, double dy=0.0, double *w=0, double *s=0, TVector3 *point_in_cell=0, TVector3 *point_on_track=0) const
Find the distance of closest approach between a cell and a track.
const std::set< unsigned int > & GetPlanesByView(View_t v=kXorY) const
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
assert(nhit_max >=nhit_nbins)
int getCellID(const CellUniqueId &id) const
std::string ExtractGDML() const
Extract contents from fGDMLFile and return as a string.
bool IntersectsDetector(double *xyz_cm, double *dxyz) const
void DetectorEnclosureBox(double *xlo_cm, double *xhi_cm, double *ylo_cm, double *yhi_cm, double *zlo_cm, double *zhi_cm) const
double fDetectorBigBoxYHi
Definition: GeometryBase.h:316
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
Atom< std::string > GDMLFile
Definition: GeometryBase.h:70
tuple density
Definition: barite.py:33
unsigned int NPlanes() const
T min(const caf::Proxy< T > &a, T b)
double fFiducialVolumeYHi
Definition: GeometryBase.h:325
std::string fGDMLFromFCL
keep track of original fcl parameter
Definition: GeometryBase.h:283
Float_t e
Definition: plot.C:35
bool IntersectsBigBox(double *xyz_cm, double *dxyz) 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
double fFiducialVolumeYLo
Definition: GeometryBase.h:326
Float_t w
Definition: plot.C:20
geo::OfflineChan getPlaneCellMap(const CellUniqueId &id) const
bool IsActive(const TGeoMaterial *m) const
bool isInsideFiducialVolume(const double x_cm, const double y_cm, const double z_cm) const
Is the particle inside the detector Fiducial Volume?
virtual void setDetectorID(novadaq::cnv::DetId)
Method to set DetectorID.
const unsigned int NextPlaneOtherView(unsigned int p, int d=+1) const
void DetectorBigBox(double *xlo_cm, double *xhi_cm, double *ylo_cm, double *yhi_cm, double *zlo_cm, double *zhi_cm) const
std::string FileBaseName() const
void RemoveTmpFile(std::string fileName)
Method to remove temporary file that holds detector geometry file.
const unsigned int FirstPlaneInMuonCatcher() const
Returns the index of the first plane contained in the muon catcher.