ROOTGeomAnalyzer.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*
3  Copyright (c) 2003-2019, The GENIE Collaboration
4  For the full text of the license visit http://copyright.genie-mc.org
5  or see $GENIE/LICENSE
6 
7  Author: Anselmo Meregaglia <anselmo.meregaglia \at cern.ch>, ETH Zurich
8  Costas Andreopoulos <costas.andreopoulos \at stfc.ac.uk>, STFC - Rutherford Lab
9  Robert Hatcher <rhatcher \at fnal.gov>, Fermilab
10 
11  May 24, 2005
12 
13  For the class documentation see the corresponding header file.
14 
15  Important revisions after version 2.0.0 :
16  @ Feb 28, 2008 - CA
17  Slight code restructuring to make it easier keeping track of conversions
18  between the SI and the current geometry system of units.
19  Fixed a bug in density unit conversions reported by Elaine Schulte.
20  @ Mar 28, 2008 - Pawel Guzowksi, Jim Dobson
21  Fixed bug preventing the code crossing more than a fixed number of volume
22  boundaries. The problem was not seen before as the code was tested with
23  relatively simpler geometries (small number of volumes).
24  Fixed problem with SetTopVolume() not actually setting it, so that the
25  option to specify sub-volumes of the master volume wasn't working properly.
26  @ Nov 12, 2008 - CA, Jim Dobson
27  Cleaned-up geometry navigation code. Improve the vertex placement method
28  GenerateVertex() to avoid geo volume overshots.
29  @ May 30, 2008 - CA, Jim Dobson
30  Fixed path-lengths units problem in Local2SI(). The problem went undetected
31  for long as the relevant quantities were used in ratios & the extra
32  multiplicative factor was cancelled out. It was spotted when we tried to
33  work out T2K event sample normalization in terms of POTs.
34  @ June 4, 2008 - CA
35  Modified code forming pdg codes in case of mixtures: Getting A from
36  TGeoMixture::GetAmixt()[ielement] rather than TGeoElement::GetA().
37  Updated code enables GENIE to see all isotopes in the fixed nd280 geometry.
38  @ August 29, 2008 - Pawel Guzowksi
39  Fixed a long-standing limitation: Now the top volume can be any geometry
40  volume and not just the master geometry volume. ROOT is placing the origin
41  of the coordinatee system at the centre of the top volume. On the other
42  hand GENIE assumed the the origin of the coordinate system is at the centre
43  of the master geometry volume. So, if the top volume is to be set to anything
44  other than the master geometry volume, an appropriate transformation is
45  required for the input flux neutrino positions/directions and the generated
46  neutrino interaction vertices.
47  Added Master2Top(v), Master2TopDir(v) and Top2Master(v) methods and the
48  fMasterToTop TGeoHMatrix to store the transformation matrix between the
49  master and top volume coordinates.
50  @ June 4, 2009 - RWH
51  Refactorized code so swimming the same start point and direction only
52  ever happens once for all the target PDG codes; do so generates a
53  PathSegmentList from whence the individual PDG lengths can be derived.
54  This also allows for future culling or trimming of path segments by an
55  externally supplied function as well as being more efficient.
56  @ July 17, 2009 - RWH
57  StepUntilEntering() wasn't accumulating total distance stepped if the
58  StepToNextBoudary() didn't leave the position IsEntering().
59  Also, small tweaks to LOG messages.
60  @ July 27, 2009 - RWH
61  The method StepToNextBoundary() doesn't actually move the current point
62  but just finds the boundary so the returned "step" size shouldn't be part
63  of the sum in StepUntilEntering().
64  @ August 3, 2009 - RWH
65  Handle case where initial step is from outside the geometry into the top
66  volume. It appears that the volume name is given as the top volume
67  (perhaps incorrectly), but the material (correctly) is null. Such steps
68  don't contribute to any material path length but document the total path.
69  Provision for keeping track of maximum ray distance and step size errors.
70  @ August 28, 2009 - RWH
71  Adjust to PathSegmentList move to genie::geometry namespace.
72  Add AdoptGeomVolSelector() function with takes ownership of GeomVolSelectorI*.
73  This selector can generate an trimmed PathSegmentList* which can be swapped
74  in for the original to limit what material is considered.
75  @ February 4, 2010 - RWH
76  Allow forcing fetch of geometry hierarchy path (or'ed with desire of the
77  GeomVolSelector, if any) but to get it by default as it is costly.
78  Fill the PathSegment medium and material upon entry along with volume name,
79  not at exit.
80  @ February 4, 2011 - JD
81  Change the way the topvolume is matched in SetTopVolName(string name).
82  Previously used TString::Contains("vol2match") which did not require the string
83  length to be the same and sometime lead to degeneracies and selection of
84  incorrect top volume. Bug and fix were found by Kevin Connolly.
85 
86 */
87 //____________________________________________________________________________
88 
89 #include <cassert>
90 #include <cstdlib>
91 #include <iomanip>
92 #include <set>
93 
94 #include <TGeoVolume.h>
95 #include <TGeoManager.h>
96 #include <TGeoShape.h>
97 #include <TGeoMedium.h>
98 #include <TGeoMaterial.h>
99 #include <TGeoMatrix.h>
100 #include <TGeoNode.h>
101 #include <TObjArray.h>
102 #include <TLorentzVector.h>
103 #include <TList.h>
104 #include <TSystem.h>
105 #include <TMath.h>
106 #include <TPolyMarker3D.h>
107 #include <TGeoBBox.h>
108 
122 
123 using namespace genie;
124 using namespace genie::geometry;
125 using namespace genie::controls;
126 
127 //#define RWH_DEBUG
128 //#define RWH_DEBUG_2
129 //#define RWH_COUNTVOLS
130 
131 #ifdef RWH_COUNTVOLS
132 // keep some statistics about how many volumes traversed for each box face
133 long int mxsegments = 0; //rwh
134 long int curface = 0; //rwh
135 long int nswims[6] = { 0, 0, 0, 0, 0, 0}; //rwh
136 long int nnever[6] = { 0, 0, 0, 0, 0, 0}; //rwh
137 double dnvols[6] = { 0, 0, 0, 0, 0, 0}; //rwh
138 double dnvols2[6] = { 0, 0, 0, 0, 0, 0}; //rwh
139 bool accum_vol_stat = false;
140 #endif
141 
142 //___________________________________________________________________________
143 ROOTGeomAnalyzer::ROOTGeomAnalyzer(string geometry_filename)
144  : GeomAnalyzerI()
145 {
146 ///
147 /// Constructor from a geometry file
148 ///
149  LOG("GROOTGeom", pDEBUG)
150  << "ROOTGeomAnalyzer ctor \"" << geometry_filename << "\"";
151  this->Initialize();
152  this->Load(geometry_filename);
153 }
154 
155 //___________________________________________________________________________
157  : GeomAnalyzerI()
158 {
159 ///
160 /// Constructor from a TGeoManager
161 ///
162  LOG("GROOTGeom", pDEBUG)
163  << "ROOTGeomAnalyzer ctor passed TGeoManager*";
164  this->Initialize();
165  this->Load(gm);
166 }
167 
168 //___________________________________________________________________________
170 {
171  this->CleanUp();
172 
173  if ( fmxddist > 0 || fmxdstep > 0 )
174  LOG("GROOTGeom",pNOTICE)
175  << "ROOTGeomAnalyzer "
176  << " mxddist " << fmxddist
177  << " mxdstep " << fmxdstep;
178 }
179 
180 //===========================================================================
181 // Geometry driver interface implementation:
182 
183 //___________________________________________________________________________
185 {
186  return *fCurrPDGCodeList;
187 }
188 
189 //___________________________________________________________________________
191 {
192 /// Computes the maximum path lengths for all materials in the input
193 /// geometry. The computed path lengths are in SI units (kgr/m^2, if
194 /// density weighting is enabled)
195 
196  LOG("GROOTGeom", pNOTICE)
197  << "Computing the maximum path lengths for all materials";
198 
199  if (!fGeometry) {
200  LOG("GROOTGeom", pFATAL) << "No ROOT geometry is loaded!!";
201  exit(1);
202  }
203 
204  //-- initialize max path lengths
206 
207  //-- select maximum path length calculation method
208  if ( fFlux ) {
209  this->MaxPathLengthsFluxMethod();
210  // clear any accumulated exposure accounted generated
211  // while exploring the geometry
212  fFlux->Clear("CycleHistory");
213  } else {
214  this->MaxPathLengthsBoxMethod();
215  }
216 
217  return *fCurrMaxPathLengthList;
218 }
219 
220 //___________________________________________________________________________
222  const TLorentzVector & x, const TLorentzVector & p)
223 {
224 /// Computes the path-length within each detector material for a
225 /// neutrino starting from point x (master coord) and travelling along
226 /// the direction of p (master coord).
227 /// The computed path lengths are in SI units (kgr/m^2, if density
228 /// weighting is enabled)
229 
230  //LOG("GROOTGeom", pDEBUG)
231  // << "Computing path-lengths for the input neutrino";
232 
233 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
234  LOG("GROOTGeom", pDEBUG)
235  << "\nInput nu: 4p (GeV) = " << utils::print::P4AsShortString(&p)
236  << ", 4x (m,s) = " << utils::print::X4AsString(&x);
237 #endif
238 
239  // if trimming configure with neutrino ray's info
240  if ( fGeomVolSelector ) {
243  }
244 
245  TVector3 udir = p.Vect().Unit(); // unit vector along direction
246  TVector3 pos = x.Vect(); // initial position
247  this->SI2Local(pos); // SI -> curr geom units
248 
249  if (!fMasterToTopIsIdentity) {
250  this->Master2Top(pos); // transform position (master -> top)
251  this->Master2TopDir(udir); // transform direction (master -> top)
252  }
253 
254  // reset current list of path-lengths
256 
257  //loop over materials & compute the path-length
258  vector<int>::iterator itr;
259  for (itr=fCurrPDGCodeList->begin();itr!=fCurrPDGCodeList->end();itr++) {
260 
261  int pdgc = *itr;
262 
263  Double_t pl = this->ComputePathLengthPDG(pos,udir,pdgc);
265 
266 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
267  LOG("GROOTGeom", pINFO)
268  <<"Calculated path length for material: " << pdgc << " = " << pl;
269 #endif
270 
271  } // loop over materials
272 
273  this->Local2SI(*fCurrPathLengthList); // curr geom units -> SI
274 
275  return *fCurrPathLengthList;
276 }
277 
278 //___________________________________________________________________________
280  const TLorentzVector & x, const TLorentzVector & p, int tgtpdg)
281 {
282 /// Generates a random vertex, within the detector material with the input
283 /// PDG code, for a neutrino starting from point x (master coord) and
284 /// travelling along the direction of p (master coord).
285 
286  LOG("GROOTGeom", pNOTICE)
287  << "Generating vtx in material: " << tgtpdg
288  << " along the input neutrino direction";
289 
290  int nretry = 0;
291  retry: // goto label in case of abject failure
292  nretry++;
293 
294  // reset current interaction vertex
295  fCurrVertex->SetXYZ(0.,0.,0.);
296 
297 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
298  LOG("GROOTGeom", pDEBUG)
299  << "\nInput nu: 4p (GeV) = " << utils::print::P4AsShortString(&p)
300  << ", 4x (m,s) = " << utils::print::X4AsString(&x);
301 #endif
302 
303  if (!fGeometry) {
304  LOG("GROOTGeom", pFATAL) << "No ROOT geometry is loaded!!";
305  exit(1);
306  }
307 
308  // calculate the max path length for the selected material starting from
309  // x and looking along the direction of p
310  TVector3 udir = p.Vect().Unit();
311  TVector3 pos = x.Vect();
312  this->SI2Local(pos); // SI -> curr geom units
313 
314  if (!fMasterToTopIsIdentity) {
315  this->Master2Top(pos); // transform position (master -> top)
316  this->Master2TopDir(udir); // transform direction (master -> top)
317  }
318 
319  double maxwgt_dist = this->ComputePathLengthPDG(pos,udir,tgtpdg);
320  if ( maxwgt_dist <= 0 ) {
321  LOG("GROOTGeom", pERROR)
322  << "The current trajectory does not cross the selected material!!";
323  return *fCurrVertex;
324  }
325 
326  // generate random number between 0 and max_dist
328  double genwgt_dist(maxwgt_dist * rnd->RndGeom().Rndm());
329 
330  LOG("GROOTGeom", pINFO)
331  << "Swim mass: Top Vol dir = " << utils::print::P3AsString(&udir)
332  << ", pos = " << utils::print::Vec3AsString(&pos);
333  LOG("GROOTGeom", pINFO)
334  << "Max {L x Density x Weight} given (init,dir) = " << maxwgt_dist;
335  LOG("GROOTGeom", pINFO)
336  << "Generated 'distance' in selected material = " << genwgt_dist;
337 #ifdef RWH_DEBUG
338  if ( ( fDebugFlags & 0x01 ) ) {
340  LOG("GROOTGeom", pINFO) << *fCurrPathSegmentList; //RWH
341  double mxddist = 0, mxdstep = 0;
342  fCurrPathSegmentList->CrossCheck(mxddist,mxdstep);
343  fmxddist = TMath::Max(fmxddist,mxddist);
344  fmxdstep = TMath::Max(fmxdstep,mxdstep);
345  }
346 #endif
347 
348  // compute the pdg weight for each material just once, then use a stl map
354  // loop over map to get tgt weight for each material (once)
355  // steps outside the geometry may have no assigned material
356  for ( ; mitr != mitr_end; ++mitr ) {
357  const TGeoMaterial* mat = mitr->first;
358  double wgt = ( mat ) ? this->GetWeight(mat,tgtpdg) : 0;
359  wgtmap[mat] = wgt;
360 #ifdef RWH_DEBUG
361  if ( ( fDebugFlags & 0x02 ) ) {
362  LOG("GROOTGeom", pINFO)
363  << " wgtmap[" << mat->GetName() << "] pdg " << tgtpdg << " wgt " << Form("%.6f",wgt);
364  }
365 #endif
366  }
367 
368  // walk down the path to pick the vertex
372  double walked = 0;
373  for ( sitr = segments.begin(); sitr != segments.end(); ++sitr) {
374  const genie::geometry::PathSegment& seg = *sitr;
375  const TGeoMaterial* mat = seg.fMaterial;
376  double trimmed_step = seg.GetSummedStepRange();
377  double wgtstep = trimmed_step * wgtmap[mat];
378  double beyond = walked + wgtstep;
379 #ifdef RWH_DEBUG
380  if ( ( fDebugFlags & 0x04 ) ) {
381  LOG("GROOTGeom", pINFO)
382  << " beyond " << beyond << " genwgt_dist " << genwgt_dist
383  << " trimmed_step " << trimmed_step << " wgtstep " << wgtstep;
384  }
385 #endif
386  if ( beyond > genwgt_dist ) {
387  // the end of this segment is beyond our generation point
388 #ifdef RWH_DEBUG
389  if ( ( fDebugFlags & 0x08 ) ) {
390  LOG("GROOTGeom", pINFO)
391  << "Choose vertex pos walked=" << walked
392  << " beyond=" << beyond
393  << " wgtstep " << wgtstep
394  << " ( " << trimmed_step << "*" << wgtmap[mat] << ")"
395  << " look for " << genwgt_dist
396  << " in " << seg.fVolume->GetName() << " "
397  << mat->GetName();
398  }
399 #endif
400  // choose a vertex in this segment (possibly multiple steps)
401  double frac = ( genwgt_dist - walked ) / wgtstep;
402  if ( frac > 1.0 ) {
403  LOG("GROOTGeom", pWARN)
404  << "Hey, frac = " << frac << " ( > 1.0 ) "
405  << genwgt_dist << " " << walked << " " << wgtstep;
406  }
407  pos = seg.GetPosition(frac);
408  fGeometry -> SetCurrentPoint (pos[0],pos[1],pos[2]);
409  fGeometry -> FindNode();
410  LOG("GROOTGeom", pINFO)
411  << "Choose vertex position in " << seg.fVolume->GetName() << " "
413  break;
414  }
415  walked = beyond;
416  }
417 
418  LOG("GROOTGeom", pNOTICE)
419  << "The vertex was placed in volume: "
420  << fGeometry->GetCurrentVolume()->GetName()
421  << ", path: " << fGeometry->GetPath();
422 
423  // warn for any volume overshoots
424  bool ok = this->FindMaterialInCurrentVol(tgtpdg);
425  if (!ok) {
426  LOG("GROOTGeom", pWARN)
427  << "Geometry volume was probably overshot";
428  LOG("GROOTGeom", pWARN)
429  << "No material with code = " << tgtpdg << " could be found at genwgt_dist="
430  << genwgt_dist << " (maxwgt_dist=" << maxwgt_dist << ")";
431  if ( nretry < 10 ) {
432  LOG("GROOTGeom", pWARN)
433  << "retry placing vertex";
434  goto retry; // yeah, I know! MyBad.
435  }
436  }
437 
438  if (!fMasterToTopIsIdentity) {
439  this->Top2Master(pos); // transform position (top -> master)
440  }
441 
442  this->Local2SI(pos); // curr geom units -> SI
443 
444  fCurrVertex->SetXYZ(pos[0],pos[1],pos[2]);
445 
446 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
447  LOG("GROOTGeom", pDEBUG)
448  << "Vtx (m) = " << utils::print::Vec3AsString(&pos);
449 #endif
450 
451  return *fCurrVertex;
452 }
453 
454 //===========================================================================
455 // Driver configuration methods:
456 
457 //___________________________________________________________________________
459 {
460 /// Use the units of the input geometry,
461 /// e.g. SetLengthUnits(genie::units::centimeter)
462 /// GENIE uses the physical system of units (hbar=c=1) almost throughtout
463 /// so everything is expressed in GeV but when analyzing detector geometries
464 /// use meters. Setting input geometry units will allow the code to compute
465 /// the conversion factor.
466 /// As input, use one of the constants in $GENIE/src/Conventions/Units.h
467 
469  LOG("GROOTGeom", pNOTICE)
470  << "Geometry length units scale factor (geom units -> m): "
471  << fLengthScale;
472 }
473 
474 //___________________________________________________________________________
476 {
477 /// Like SetLengthUnits, but for density (default units = kgr/m3)
478 
480  LOG("GROOTGeom", pNOTICE)
481  << "Geometry density units scale factor (geom units -> kgr/m3): "
482  << fDensityScale;
483 }
484 
485 //___________________________________________________________________________
487 {
488 /// Set a factor that can multiply the computed max path lengths.
489 /// The maximum path lengths are computed by performing an MC scanning of
490 /// the input geometry. If you configure the scanner with a low number of
491 /// points or rays you might understimate the path lengths, so you might
492 /// want to 'inflate' them a little bit using this method.
493 /// Do not set this number too high, because the max interaction probability
494 /// will be grossly overestimated and you would need lots of attempts before
495 /// getting a flux neutrino to interact...
496 
497  fMaxPlSafetyFactor = sf;
498 
499  LOG("GROOTGeom", pNOTICE)
500  << "Max path length safety factor: " << fMaxPlSafetyFactor;
501 }
502 
503 //___________________________________________________________________________
505 {
506 /// Set it to x, if the relative weight proportions of elements in a mixture
507 /// add up to x (eg x=1, 100, etc). Set it to a negative value to explicitly
508 /// compute the correct weight normalization.
509 
510  fMixtWghtSum = sum;
511 }
512 
513 //___________________________________________________________________________
515 {
516 /// Set the name of the top volume.
517 /// This driver would ask the TGeoManager::GetTopVolume() for the top volume.
518 /// Use this method for changing this if for example you want to set a smaller
519 /// volume as the top one so as to generate events only in a specific part of
520 /// your detector.
521 
522  if (name.size() == 0) return;
523 
525  LOG("GROOTGeom",pNOTICE) << "Geometry Top Volume name: " << fTopVolumeName;
526 
527  TGeoVolume * gvol = fGeometry->GetVolume(fTopVolumeName.c_str());
528  if (!gvol) {
529  LOG("GROOTGeom",pWARN) << "Could not find volume: " << name.c_str();
530  LOG("GROOTGeom",pWARN) << "Will not change the current top volume";
531  fTopVolumeName = "";
532  return;
533  }
534 
535  // Get a matrix connecting coordinates of master and top volumes.
536  // The matrix will be used for transforming the coordinates of incoming
537  // flux neutrinos & generated interaction vertices.
538  // This is needed (in case that the input top volume != master volume)
539  // because ROOT always sets the coordinate system origin at the centre of
540  // the specified top volume (whereas GENIE assumes that the global reference
541  // frame is that of the master volume)
542 
543  TGeoIterator next(fGeometry->GetMasterVolume());
544  TGeoNode *node;
545  TString nodeName, volNameStr;
546  const char* volName = fTopVolumeName.c_str();
547  while ((node = next())) {
548  nodeName = node->GetVolume()->GetName();
549  if (nodeName == volName) {
550  if (fMasterToTop) delete fMasterToTop;
551  fMasterToTop = new TGeoHMatrix(*next.GetCurrentMatrix());
552  fMasterToTopIsIdentity = fMasterToTop->IsIdentity();
553  break;
554  }
555  }
556 
557  // set volume name
558  fTopVolume = gvol;
559  fGeometry->SetTopVolume(fTopVolume);
560 }
561 
562 //===========================================================================
563 // Geometry/Unit transforms:
564 
565 //___________________________________________________________________________
567 {
568 /// convert path lengths from current geometry units to SI units
569 ///
570 
571  double scaling_factor = this->LengthUnits();
572  if (this->WeightWithDensity()) { scaling_factor *= this->DensityUnits(); }
573 
574 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
575  LOG("GROOTGeom", pDEBUG)
576  << "Scaling path-lengths from local units -> meters "
577  << ((this->WeightWithDensity()) ? "* kgr/m^3" : "")
578  << " - scale = " << scaling_factor;
579 #endif
580 
581  PathLengthList::iterator pliter;
582  for(pliter = pl.begin(); pliter != pl.end(); ++pliter)
583  {
584  int pdgc = pliter->first;
585  pl.ScalePathLength(pdgc, scaling_factor);
586  }
587 }
588 
589 //___________________________________________________________________________
590 void ROOTGeomAnalyzer::Local2SI(TVector3 & vec) const
591 {
592 /// convert position vector from current geometry units to SI units
593 ///
594 
595 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
596  LOG("GROOTGeom", pDEBUG)
597  << "Position (loc): " << utils::print::Vec3AsString(&vec);
598 #endif
599 
600  vec *= (this->LengthUnits());
601 
602 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
603  LOG("GROOTGeom", pDEBUG)
604  << "Position (SI): " << utils::print::Vec3AsString(&vec);
605 #endif
606 }
607 
608 //___________________________________________________________________________
609 void ROOTGeomAnalyzer::SI2Local(TVector3 & vec) const
610 {
611 /// convert position vector from SI units to current geometry units
612 ///
613 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
614  LOG("GROOTGeom", pDEBUG)
615  << "Position (SI): " << utils::print::Vec3AsString(&vec);
616 #endif
617 
618  vec *= (1./this->LengthUnits());
619 
620 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
621  LOG("GROOTGeom", pDEBUG)
622  << "Position (loc): " << utils::print::Vec3AsString(&vec);
623 #endif
624 }
625 
626 //___________________________________________________________________________
627 void ROOTGeomAnalyzer::Master2Top(TVector3 & vec) const
628 {
629 /// transform the input position vector from the master volume coordinate
630 /// system to the specified top volume coordinate system, but not
631 /// change the units.
632 
633 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
634  LOG("GROOTGeom", pDEBUG)
635  << "Position (coord:master): " << utils::print::Vec3AsString(&vec);
636 #endif
637 
638  Double_t mast[3], top[3];
639  vec.GetXYZ(mast);
640  fMasterToTop->MasterToLocal(mast, top);
641  vec.SetXYZ(top[0], top[1], top[2]);
642 
643 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
644  LOG("GROOTGeom", pDEBUG)
645  << "Position (coord:top): " << utils::print::Vec3AsString(&vec);
646 #endif
647 }
648 
649 //___________________________________________________________________________
650 void ROOTGeomAnalyzer::Master2TopDir(TVector3 & vec) const
651 {
652 /// transform the input direction vector from the master volume coordinate
653 /// system to the specified top volume coordinate system, but not
654 /// change the units.
655 
656 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
657  LOG("GROOTGeom", pDEBUG)
658  << "Direction (coord:master): " << utils::print::Vec3AsString(&vec);
659 #endif
660 
661  Double_t mast[3], top[3];
662  vec.GetXYZ(mast);
663  fMasterToTop->MasterToLocalVect(mast, top);
664  vec.SetXYZ(top[0], top[1], top[2]);
665 
666 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
667  LOG("GROOTGeom", pDEBUG)
668  << "Direction (coord:top): " << utils::print::Vec3AsString(&vec);
669 #endif
670 }
671 
672 //___________________________________________________________________________
673 void ROOTGeomAnalyzer::Top2Master(TVector3 & vec) const
674 {
675 /// transform the input position vector from the specified top volume
676 /// coordinate system to the master volume coordinate system, but not
677 /// change the units.
678 
679 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
680  LOG("GROOTGeom", pDEBUG)
681  << "Position (coord:top): " << utils::print::Vec3AsString(&vec);
682 #endif
683 
684  Double_t mast[3], top[3];
685  vec.GetXYZ(top);
686  fMasterToTop->LocalToMaster(top, mast);
687  vec.SetXYZ(mast[0], mast[1], mast[2]);
688 
689 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
690  LOG("GROOTGeom", pDEBUG)
691  << "Position (coord:master): " << utils::print::Vec3AsString(&vec);
692 #endif
693 }
694 
695 //___________________________________________________________________________
696 void ROOTGeomAnalyzer::Top2MasterDir(TVector3 & vec) const
697 {
698 /// transform the input direction vector from the specified top volume
699 /// coordinate system to the master volume coordinate system.
700 
701 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
702  LOG("GROOTGeom", pDEBUG)
703  << "Direction (coord:top): " << utils::print::Vec3AsString(&vec);
704 #endif
705 
706  Double_t mast[3], top[3];
707  vec.GetXYZ(top);
708  fMasterToTop->LocalToMasterVect(top, mast);
709  vec.SetXYZ(mast[0], mast[1], mast[2]);
710 
711 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
712  LOG("GROOTGeom", pDEBUG)
713  << "Direction (coord:master): " << utils::print::Vec3AsString(&vec);
714 #endif
715 }
716 
717 //===========================================================================
718 // Private methods:
719 
720 //___________________________________________________________________________
722 {
723  LOG("GROOTGeom", pNOTICE)
724  << "Initializing ROOT geometry driver & setting defaults";
725 
729  fGeomVolSelector = 0;
730  fCurrPDGCodeList = 0;
731  fTopVolume = 0;
732  fTopVolumeName = "";
733  fKeepSegPath = false;
734 
735  // some defaults:
736  this -> SetScannerNPoints (200);
737  this -> SetScannerNRays (200);
738  this -> SetScannerNParticles (10000);
739  this -> SetScannerFlux (0);
740  this -> SetMaxPlSafetyFactor (1.1);
743  this -> SetWeightWithDensity (true);
744  this -> SetMixtureWeightsSum (-1.);
745 
746  fMasterToTopIsIdentity = true;
747 
748  fmxddist = 0;
749  fmxdstep = 0;
750  fDebugFlags = 0;
751 }
752 
753 //___________________________________________________________________________
755 {
756  LOG("GROOTGeom", pNOTICE) << "Cleaning up...";
757 
761  if ( fCurrPDGCodeList ) delete fCurrPDGCodeList;
762  if ( fMasterToTop ) delete fMasterToTop;
763 }
764 
765 //___________________________________________________________________________
767 {
768 /// Load the detector geometry from the input ROOT file
769 ///
770  LOG("GROOTGeom", pNOTICE) << "Loading geometry from: " << filename;
771 
772  bool is_accessible = ! (gSystem->AccessPathName( filename.c_str() ));
773  if (!is_accessible) {
774  LOG("GROOTGeom", pFATAL)
775  << "The ROOT geometry doesn't exist! Initialization failed!";
776  exit(1);
777  }
778  TGeoManager * gm = TGeoManager::Import(filename.c_str());
779 
780  this->Load(gm);
781 }
782 
783 //___________________________________________________________________________
784 void ROOTGeomAnalyzer::Load(TGeoManager * gm)
785 {
786 /// Load the detector geometry from the input TGeoManager
787 
788  LOG("GROOTGeom", pNOTICE)
789  << "A TGeoManager is being loaded to the geometry driver";
790  fGeometry = gm;
791 
792  if (!fGeometry) {
793  LOG("GROOTGeom", pFATAL) << "Null TGeoManager! Aborting";
794  }
795  assert(fGeometry);
796 
797  this->BuildListOfTargetNuclei();
798 
799  const PDGCodeList & pdglist = this->ListOfTargetNuclei();
800 
801  fTopVolume = 0;
803  fCurrPathLengthList = new PathLengthList(pdglist);
804  fCurrMaxPathLengthList = new PathLengthList(pdglist);
805  fCurrVertex = new TVector3(0.,0.,0.);
806 
807  // ask geometry manager for its top volume
808  fTopVolume = fGeometry->GetTopVolume();
809  if (!fTopVolume) {
810  LOG("GROOTGeom", pFATAL) << "Could not get top volume!!!";
811  }
813 
814  // load matrix (identity) of top volume
815  fMasterToTop = new TGeoHMatrix(*fGeometry->GetCurrentMatrix());
816  fMasterToTopIsIdentity = true;
817 
818 //#define PRINT_MATERIALS
819 #ifdef PRINT_MATERIALS
820  fGeometry->GetListOfMaterials()->Print();
821  fGeometry->GetListOfMedia()->Print();
822 #endif
823 
824 }
825 
826 //___________________________________________________________________________
828 {
829 /// Determine possible target PDG codes.
830 /// Note: If one is using a top volume other than the master level
831 /// then the final list might include PDG codes that will never
832 /// be seen during swimming through the volumes if those code are only
833 /// found in materials outside the top volume.
834 
836 
837  if (!fGeometry) {
838  LOG("GROOTGeom", pFATAL) << "No ROOT geometry is loaded!!";
839  exit(1);
840  }
841 
842  TObjArray * volume_list = fGeometry->GetListOfVolumes();
843  if (!volume_list) {
844  LOG("GROOTGeom", pERROR)
845  << "Null list of geometry volumes. Can not find build target list!";
846  return;
847  }
848 
849  std::set<Int_t> seen_mat; // list of materials we've already handled
850  std::vector<TGeoVolume*> volvec; //RWH
851 
852  int numVol = volume_list->GetEntries();
853 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
854  LOG("GROOTGeom", pDEBUG) << "Number of volumes found: " << numVol;
855 #endif
856 
857  for (int ivol = 0; ivol < numVol; ivol++) {
858  TGeoVolume * volume = dynamic_cast <TGeoVolume *>(volume_list->At(ivol));
859  if (!volume) {
860  LOG("GROOTGeom", pWARN)
861  << "Got a null geometry volume!! Skiping current list element";
862  continue;
863  }
864  TGeoMaterial * mat = volume->GetMedium()->GetMaterial();
865 
866  // shortcut if we've already seen this material
867  Int_t mat_indx = mat->GetIndex();
868  if ( seen_mat.find(mat_indx) != seen_mat.end() ) continue;
869  seen_mat.insert(mat_indx);
870  volvec.push_back(volume); //RWH
871 
872  if (mat->IsMixture()) {
873  TGeoMixture * mixt = dynamic_cast <TGeoMixture*> (mat);
874  int Nelements = mixt->GetNelements();
875  for (int i=0; i<Nelements; i++) {
876  int ion_pdgc = this->GetTargetPdgCode(mixt,i);
877  fCurrPDGCodeList->push_back(ion_pdgc);
878  }
879  } else {
880  int ion_pdgc = this->GetTargetPdgCode(mat);
881  fCurrPDGCodeList->push_back(ion_pdgc);
882  }
883  }
884  // sort the list
885  // we don't calculate this list but once per geometry and a sorted
886  // list is easier to read so this doesn't cost much
887  std::sort(fCurrPDGCodeList->begin(),fCurrPDGCodeList->end());
888 
889 }
890 
891 //___________________________________________________________________________
892 int ROOTGeomAnalyzer::GetTargetPdgCode(const TGeoMaterial * const m) const
893 {
894  int A = TMath::Nint(m->GetA());
895  int Z = TMath::Nint(m->GetZ());
896 
897  int pdgc = pdg::IonPdgCode(A,Z);
898 
899  return pdgc;
900 }
901 
902 //___________________________________________________________________________
904  const TGeoMixture * const m, int ielement) const
905 {
906  int A = TMath::Nint(m->GetAmixt()[ielement]);
907  int Z = TMath::Nint(m->GetZmixt()[ielement]);
908 
909  int pdgc = pdg::IonPdgCode(A,Z);
910 
911  return pdgc;
912 }
913 
914 //___________________________________________________________________________
915 double ROOTGeomAnalyzer::GetWeight(const TGeoMaterial * mat, int pdgc)
916 {
917 /// Get the weight of the input material.
918 /// Return the weight only if the material's pdg code matches the input code.
919 /// If the material is found to be a mixture, call the corresponding method
920 /// for mixtures.
921 /// Weight is in the curr geom density units.
922 
923  if (!mat) {
924  LOG("GROOTGeom", pERROR) << "Null input material. Return weight = 0.";
925  return 0;
926  }
927 
928  bool exists = fCurrPDGCodeList->ExistsInPDGCodeList(pdgc);
929  if (!exists) {
930  LOG("GROOTGeom", pERROR) << "Target doesn't exist. Return weight = 0.";
931  return 0;
932  }
933 
934 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
935  LOG("GROOTGeom",pDEBUG)
936  << "Curr. material: A/Z = " << mat->GetA() << " / " << mat->GetZ();
937 #endif
938 
939  // if the input material is a mixture, get a the sum of weights for
940  // all matching elements
941  double weight = 0.;
942  if (mat->IsMixture()) {
943  const TGeoMixture * mixt = dynamic_cast <const TGeoMixture*> (mat);
944  if (!mixt) {
945  LOG("GROOTGeom", pERROR) << "Null input mixture. Return weight = 0.";
946  return 0;
947  }
948 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
949  LOG("GROOTGeom", pDEBUG)
950  << "Material : " << mat->GetName()
951  << " is a mixture with " << mixt->GetNelements() << " elements";
952 #endif
953  // loop over elements & sum weights of matching elements
954  weight = this->GetWeight(mixt,pdgc);
955  return weight;
956  } // is mixture?
957 
958  // pure material
959  int ion_pdgc = this->GetTargetPdgCode(mat);
960  if (ion_pdgc != pdgc) return 0.;
961 
962  if (this->WeightWithDensity())
963  weight = mat->GetDensity(); // material density (curr geom units)
964  else
965  weight = 1.0;
966 
967 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
968  LOG("GROOTGeom", pDEBUG)
969  << "Weight[mat:" << mat->GetName() << "] = " << weight;
970 #endif
971 
972  return weight;
973 }
974 
975 //___________________________________________________________________________
976 double ROOTGeomAnalyzer::GetWeight(const TGeoMixture * mixt, int pdgc)
977 {
978 /// Loop over the mixture elements, find the one matching the input pdgc
979 /// and return its weight.
980 /// Weight is in the curr geom density units.
981 
982  double weight = 0;
983 
984  int nm = 0;
985  for (int i = 0; i < mixt->GetNelements(); i++) {
986  double dw = (this->GetWeight(mixt,i,pdgc));
987  if (dw>0) nm++;
988  weight += dw;
989  }
990 
991  if (nm>1) {
992  for (int j = 0; j < mixt->GetNelements(); j++) {
993  LOG("GROOTGeom", pWARN)
994  << "[" << j << "] Z = " << mixt->GetZmixt()[j]
995  << ", A = " << mixt->GetAmixt()[j]
996  << " (pdgc = " << this->GetTargetPdgCode(mixt,j)
997  << "), w = " << mixt->GetWmixt()[j];
998  }
999  LOG("GROOTGeom", pERROR)
1000  << "Material pdgc = " << pdgc << " appears " << nm
1001  << " times (>1) in mixture = " << mixt->GetName();
1002  LOG("GROOTGeom", pFATAL)
1003  << "Your geometry must be incorrect - Aborting";
1004  exit(1);
1005  }
1006 
1007  // if we are not weighting with the density then the weight=1 if the pdg
1008  // code was matched for any element of this mixture
1009  if ( !this->WeightWithDensity() && weight>0. ) weight=1.0;
1010 
1011  return weight;
1012 }
1013 
1014 //___________________________________________________________________________
1015 double ROOTGeomAnalyzer::GetWeight(const TGeoMixture* mixt, int ielement, int pdgc)
1016 {
1017 /// Get the weight of the input ith element of the input material.
1018 /// Return the weight only if the element's pdg code matches the input code.
1019 /// Weight is in the curr geom density units.
1020 
1021 // int ion_pdgc = this->GetTargetPdgCode(mixt->GetElement(ielement));
1022  int ion_pdgc = this->GetTargetPdgCode(mixt, ielement);
1023  if (ion_pdgc != pdgc) return 0.;
1024 
1025  double d = mixt->GetDensity(); // mixture density (curr geom units)
1026  double w = mixt->GetWmixt()[ielement]; // relative proportion by mass
1027 
1028  double wtot = this->MixtureWeightsSum();
1029 
1030  // <0 forces explicit calculation of relative proportion normalization
1031  if (wtot < 0) {
1032  wtot = 0;
1033  for (int i = 0; i < mixt->GetNelements(); i++) {
1034  wtot += (mixt->GetWmixt()[i]);
1035  }
1036  }
1037  assert(wtot>0);
1038 
1039  w /= wtot;
1040  double weight = d*w;
1041 
1042  return weight;
1043 }
1044 
1045 //___________________________________________________________________________
1047 {
1048 /// Use the input flux driver to generate "rays", and then follow them through
1049 /// the detector and figure out the maximum path length for each material
1050 
1051  LOG("GROOTGeom", pNOTICE)
1052  << "Computing the maximum path lengths using the FLUX method";
1053 
1054  int iparticle = 0;
1055  PathLengthList::const_iterator pl_iter;
1056 
1057  const int nparticles = abs(this->ScannerNParticles());
1058 
1059  // if # scanner particles is negative, this signals that the user
1060  // desires to force rays to have the maximum energy (useful if the
1061  // volume considered changes size with neutrino energy)
1062  bool rescale_e = (this->ScannerNParticles() < 0 );
1063  double emax = fFlux->MaxEnergy();
1064  if ( rescale_e ) {
1065  LOG("GROOTGeom", pNOTICE)
1066  << "max path lengths with FLUX method forcing Enu=" << emax;
1067  }
1068 
1069  while (iparticle < nparticles ) {
1070 
1071  bool ok = fFlux->GenerateNext();
1072  if (!ok) {
1073  LOG("GROOTGeom", pWARN) << "Couldn't generate a flux neutrino";
1074  continue;
1075  }
1076 
1077  TLorentzVector nup4 = fFlux->Momentum();
1078  if ( rescale_e ) {
1079  double ecurr = nup4.E();
1080  if ( ecurr > 0 ) nup4 *= (emax/ecurr);
1081  }
1082  const TLorentzVector & nux4 = fFlux->Position();
1083 
1084  //LOG("GMCJDriver", pNOTICE)
1085  // << "\n [-] Generated flux neutrino: "
1086  // << "\n |----o 4-momentum : " << utils::print::P4AsString(&nup4)
1087  // << "\n |----o 4-position : " << utils::print::X4AsString(&nux4);
1088 
1089  const PathLengthList & pl = this->ComputePathLengths(nux4, nup4);
1090 
1091  bool enters = false;
1092 
1093  for (pl_iter = pl.begin(); pl_iter != pl.end(); ++pl_iter) {
1094  int pdgc = pl_iter->first;
1095  double pathlength = pl_iter->second;
1096 
1097  if ( pathlength > 0 ) {
1098  pathlength *= (this->MaxPlSafetyFactor());
1099 
1100  pathlength = TMath::Max(pathlength, fCurrMaxPathLengthList->PathLength(pdgc));
1101  fCurrMaxPathLengthList->SetPathLength(pdgc,pathlength);
1102  enters = true;
1103  }
1104  }
1105  if (enters) iparticle++;
1106  }
1107 }
1108 
1109 //___________________________________________________________________________
1111 {
1112 /// Generate points in the geometry's bounding box and for each point
1113 /// generate random rays, follow them through the detector and figure out
1114 /// the maximum path length for each material
1115 
1116  LOG("GROOTGeom", pNOTICE)
1117  << "Computing the maximum path lengths using the BOX method";
1118 #ifdef RWH_COUNTVOLS
1119  accum_vol_stat = true;
1120 #endif
1121 
1122  int iparticle = 0;
1123  bool ok = true;
1124  TLorentzVector nux4;
1125  TLorentzVector nup4;
1126 
1127  PathLengthList::const_iterator pl_iter;
1128 
1129  while ( (ok = this->GenBoxRay(iparticle++,nux4,nup4)) ) {
1130 
1131  //LOG("GMCJDriver", pNOTICE)
1132  // << "\n [-] Generated flux neutrino: "
1133  // << "\n |----o 4-momentum : " << utils::print::P4AsString(&nup4)
1134  // << "\n |----o 4-position : " << utils::print::X4AsString(&nux4);
1135 
1136  const PathLengthList & pllst = this->ComputePathLengths(nux4, nup4);
1137 
1138  for (pl_iter = pllst.begin(); pl_iter != pllst.end(); ++pl_iter) {
1139  int pdgc = pl_iter->first;
1140  double pl = pl_iter->second;
1141 
1142  if (pl>0) {
1143  pl *= (this->MaxPlSafetyFactor());
1144 
1145  pl = TMath::Max(pl, fCurrMaxPathLengthList->PathLength(pdgc));
1147  }
1148  }
1149  }
1150 
1151  // print out the results
1152  LOG("GROOTGeom", pDEBUG)
1153  << "DensWeight \"" << (fDensWeight?"true":"false")
1154  << "\" MixtWghtSum " << fMixtWghtSum;
1155  LOG("GROOTGeom", pDEBUG) << "CurrMaxPathLengthList: "
1157 
1158 #ifdef RWH_COUNTVOLS
1159  // rwh
1160  // print statistics for average,rms of number of volumes seen for
1161  // various rays for each face
1162  for (int j = 0; j < 6; ++j ) {
1163  long int ns = nswims[j];
1164  double x = dnvols[j];
1165  double x2 = dnvols2[j];
1166  if ( ns == 0 ) ns = 1;
1167  double avg = x / (double)ns;
1168  double rms = TMath::Sqrt((x2/(double)ns) - avg*avg);
1169  LOG("GROOTGeom", pNOTICE)
1170  << "RWH: nswim after BOX face " << j << " is " << ns
1171  << " avg " << avg << " rms " << rms
1172  << " never " << nnever[j];
1173  }
1174  LOG("GROOTGeom", pNOTICE)
1175  << "RWH: Max PathSegmentList size " << mxsegments;
1176  accum_vol_stat = false;
1177 #endif
1178 
1179 }
1180 
1181 //___________________________________________________________________________
1182 bool ROOTGeomAnalyzer::GenBoxRay(int indx, TLorentzVector& x4, TLorentzVector& p4)
1183 {
1184 /// Generate points in the geometry's bounding box and for each point generate
1185 /// random rays -- a pseudo-flux -- in master coordinates and SI units
1186 
1187  firay++;
1188  fnewpnt = false;
1189 
1190  // first time through ... special case
1191  if ( indx == 0 ) { fiface = 0; fipoint = 0; firay = 0; fnewpnt = true; }
1192 
1193  if ( firay >= fNRays ) { fipoint++; firay = 0; fnewpnt = true; }
1194  if ( fipoint >= fNPoints ) { fiface++; fipoint = 0; firay = 0; fnewpnt = true; }
1195  if ( fiface >= 3 ) {
1196  LOG("GROOTGeom",pINFO) << "Box surface scanned: " << indx << " points/rays";
1197  return false; // signal end
1198  }
1199 
1200  if ( indx == 0 ) {
1201  // get geometry's bounding box
1202  //LOG("GROOTGeom", pNOTICE) << "Getting a TGeoBBox enclosing the detector";
1203  TGeoShape * TS = fTopVolume->GetShape();
1204  TGeoBBox * box = (TGeoBBox *)TS;
1205  //get box origin and dimensions (in the same units as the geometry)
1206  fdx = box->GetDX(); // half-length
1207  fdy = box->GetDY(); // half-length
1208  fdz = box->GetDZ(); // half-length
1209  fox = (box->GetOrigin())[0];
1210  foy = (box->GetOrigin())[1];
1211  foz = (box->GetOrigin())[2];
1212 
1213  LOG("GROOTGeom",pNOTICE)
1214  << "Box size (GU) :"
1215  << " x = " << 2*fdx << ", y = " << 2*fdy << ", z = " << 2*fdz;
1216  LOG("GROOTGeom",pNOTICE)
1217  << "Box origin (GU) :"
1218  << " x = " << fox << ", y = " << foy << ", z = " << foz;
1219  LOG("GROOTGeom",pNOTICE)
1220  << "Will generate [" << fNPoints << "] random points / box surface";
1221  LOG("GROOTGeom",pNOTICE)
1222  << "Will generate [" << fNRays << "] rays / point";
1223 
1224 #ifdef VALIDATE_CORNERS
1225  // RWH: test that we know the coordinate transforms for the box corners
1226  for (int sz = -1; sz <= +1; ++sz) {
1227  for (int sy = -1; sy <= +1; ++sy) {
1228  for (int sx = -1; sx <= +1; ++sx) {
1229  if (sx == 0 || sy == 0 || sz == 0 ) continue;
1230  TVector3& pos = fGenBoxRayPos;
1231  pos.SetXYZ(fox+(sx*fdx),foy+(sy*fdy),foz+(sz*fdz));
1232  TVector3 master(fGenBoxRayPos);
1233  this->Top2Master(master); // transform position (top -> master)
1234  this->Local2SI(master);
1235  TVector3 pos2(master);
1236  this->Master2Top(pos2);
1237  this->SI2Local(pos2);
1238  LOG("GROOTGeom", pNOTICE)
1239  << "TopVol corner "
1240  << " [" << pos[0] << "," << pos[1] << "," << pos[2] << "] "
1241  << "Master corner "
1242  << " [" << master[0] << "," << master[1] << "," << master[2] << "] "
1243  << " top again"
1244  << " [" << pos2[0] << "," << pos2[1] << "," << pos2[2] << "] ";
1245  }
1246  }
1247  }
1248 #endif
1249 
1250  }
1251 
1253 
1254  double eps = 0.0; //1.0e-12; // tweak to make sure we're inside the box
1255 
1256  switch ( fiface ) {
1257  case 0:
1258  {
1259  //top:
1260  if ( firay == fNRays-1 && fipoint == fNPoints-1 )
1261  LOG("GROOTGeom",pINFO) << "Box surface scanned: [TOP]";
1262  fGenBoxRayDir.SetXYZ(-0.5+rnd->RndGeom().Rndm(),
1263  -rnd->RndGeom().Rndm(),
1264  -0.5+rnd->RndGeom().Rndm());
1265  if ( fnewpnt )
1266  fGenBoxRayPos.SetXYZ(fox-fdx+eps+2*(fdx-eps)*rnd->RndGeom().Rndm(),
1267  foy+fdy-eps,
1268  foz-fdz+eps+2*(fdz-eps)*rnd->RndGeom().Rndm());
1269  break;
1270  }
1271  case 1:
1272  {
1273  //left:
1274  if ( firay == fNRays-1 && fipoint == fNPoints-1 )
1275  LOG("GROOTGeom",pINFO) << "Box surface scanned: [LEFT]";
1276  fGenBoxRayDir.SetXYZ(rnd->RndGeom().Rndm(),
1277  -0.5+rnd->RndGeom().Rndm(),
1278  -0.5+rnd->RndGeom().Rndm());
1279  if ( fnewpnt )
1280  fGenBoxRayPos.SetXYZ(fox-fdx+eps,
1281  foy-fdy+eps+2*(fdy-eps)*rnd->RndGeom().Rndm(),
1282  foz-fdz+eps+2*(fdz-eps)*rnd->RndGeom().Rndm());
1283  break;
1284  }
1285  case 2:
1286  {
1287  //front: (really what I, RWH, would call the back)
1288  if ( firay == fNRays-1 && fipoint == fNPoints-1 )
1289  LOG("GROOTGeom",pINFO) << "Box surface scanned: [FRONT]";
1290  fGenBoxRayDir.SetXYZ(-0.5+rnd->RndGeom().Rndm(),
1291  -0.5+rnd->RndGeom().Rndm(),
1292  -rnd->RndGeom().Rndm());
1293  if ( fnewpnt )
1294  fGenBoxRayPos.SetXYZ(fox-fdx+eps+2*(fdx-eps)*rnd->RndGeom().Rndm(),
1295  foy-fdy+eps+2*(fdy-eps)*rnd->RndGeom().Rndm(),
1296  foz+fdz+eps);
1297  break;
1298  }
1299  }
1300 /*
1301  //bottom:
1302  pos.SetXYZ(ox-dx+2*dx*rnd->RndGeom().Rndm(),
1303  oy-dy,
1304  oz-dz+2*dz*rnd->RndGeom().Rndm());
1305  dir.SetXYZ(-0.5+rnd->RndGeom().Rndm(),
1306  rnd->RndGeom().Rndm(),
1307  -0.5+rnd->RndGeom().Rndm());
1308  //right:
1309  pos.SetXYZ(ox+dx,
1310  oy-dy+2*dy*rnd->RndGeom().Rndm(),
1311  oz-dz+2*dz*rnd->RndGeom().Rndm());
1312  dir.SetXYZ(-rnd->RndGeom().Rndm(),
1313  -0.5+rnd->RndGeom().Rndm(),
1314  -0.5+rnd->RndGeom().Rndm());
1315  //back:
1316  pos.SetXYZ(ox-dx+2*dx*rnd->RndGeom().Rndm(),
1317  oy-dy+2*dy*rnd->RndGeom().Rndm(),
1318  oz-dz);
1319  dir.SetXYZ(-0.5+rnd->RndGeom().Rndm(),
1320  -0.5+rnd->RndGeom().Rndm(),
1321  rnd->RndGeom().Rndm());
1322 */
1323 
1324 #ifdef RWH_COUNTVOLS
1325  curface = fiface;
1326 #endif
1327 
1328 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
1329  if ( fnewpnt )
1330  LOG("GROOTGeom", pNOTICE)
1331  << "GenBoxRay(topvol) "
1332  << " iface " << fiface << " ipoint " << fipoint << " iray " << firay
1333  << " newpnt " << (fnewpnt?"true":"false")
1336 #endif
1337 
1338  if ( fnewpnt ) {
1339  if ( ! fMasterToTopIsIdentity) {
1340  this->Top2Master(fGenBoxRayPos); // transform position (top -> master)
1341  }
1342  this->Local2SI(fGenBoxRayPos);
1343  }
1344  this->Top2MasterDir(fGenBoxRayDir); // transform direction (top -> master)
1345 
1346  x4.SetVect(fGenBoxRayPos);
1347  p4.SetVect(fGenBoxRayDir.Unit());
1348 
1349 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
1350  LOG("GROOTGeom", pNOTICE)
1351  << "GenBoxRay(master) "
1352  << " iface " << fiface << " ipoint " << fipoint << " iray " << firay
1353  << " newpnt " << (fnewpnt?"true":"false")
1356 #endif
1357 
1358  return true;
1359 }
1360 
1361 //________________________________________________________________________
1363  const TVector3 & r0, const TVector3 & udir, int pdgc)
1364 {
1365 /// Compute the path length for the material with pdg-code = pdc, staring
1366 /// from the input position r (top vol coord & units) and moving along the
1367 /// direction of the unit vector udir (top vol coord).
1368 
1369  double pl = 0; // path-length (x density, if density-weighting is ON)
1370 
1371  this->SwimOnce(r0,udir);
1372 
1373  double step = 0;
1374  double weight = 0;
1375 
1376  // const TGeoVolume * vol = 0;
1377  // const TGeoMedium * med = 0;
1378  const TGeoMaterial * mat = 0;
1379 
1380  // loop over independent materials, which is shorter or equal to # of volumes
1385  for ( ; itr != itr_end; ++itr ) {
1386  mat = itr->first;
1387  if ( ! mat ) continue; // segment outside geometry has no material
1388  step = itr->second;
1389  weight = this->GetWeight(mat,pdgc);
1390  pl += (step*weight);
1391  }
1392 
1393 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
1394  LOG("GROOTGeom", pDEBUG)
1395  << "PathLength[" << pdgc << "] = " << pl << " in curr geom units";
1396 #endif
1397 
1398  return pl;
1399 }
1400 
1401 //________________________________________________________________________
1402 void ROOTGeomAnalyzer::SwimOnce(const TVector3 & r0, const TVector3 & udir)
1403 {
1404 /// Swim through the geometry from the from the input position
1405 /// r0 (top vol coord & units) and moving along the direction of the
1406 /// unit vector udir (topvol coord) to create a filled PathSegmentList
1407 
1408  int nvolswim = 0; //rwh
1409 
1411 
1412  // don't swim if the current PathSegmentList is up-to-date
1413  if ( fCurrPathSegmentList->IsSameStart(r0,udir) ) return;
1414 
1415  // start fresh
1417 
1418  // set start info so next time we don't swim for the same ray
1420 
1421  PathSegment ps_curr;
1422 
1423  bool found_vol (false);
1424  bool keep_on (true);
1425 
1426  double step = 0;
1427  double raydist = 0;
1428 
1429  const TGeoVolume * vol = 0;
1430  const TGeoMedium * med = 0;
1431  const TGeoMaterial * mat = 0;
1432 
1433  // determining the geometry path is expensive, do it only if necessary
1434  bool selneedspath = ( fGeomVolSelector && fGeomVolSelector->GetNeedPath() );
1435  const bool fill_path = fKeepSegPath || selneedspath;
1436 
1437 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
1438  LOG("GROOTGeom", pNOTICE)
1439  << "SwimOnce x [" << r0[0] << "," << r0[1] << "," << r0[2]
1440  << "] udir [" << udir[0] << "," << udir[1] << "," << udir[2];
1441 #endif
1442 
1443  fGeometry -> SetCurrentDirection (udir[0],udir[1],udir[2]);
1444  fGeometry -> SetCurrentPoint (r0[0], r0[1], r0[2] );
1445 
1446  while (!found_vol || keep_on) {
1447  keep_on = true;
1448 
1449  fGeometry->FindNode();
1450 
1451  ps_curr.SetEnter( fGeometry->GetCurrentPoint() , raydist );
1452  vol = fGeometry->GetCurrentVolume();
1453  med = vol->GetMedium();
1454  mat = med->GetMaterial();
1455  ps_curr.SetGeo(vol,med,mat);
1456 #ifdef PATHSEG_KEEP_PATH
1457  if (fill_path) ps_curr.SetPath(fGeometry->GetPath());
1458 #endif
1459 
1460 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
1461 #ifdef DUMP_SWIM
1462  LOG("GROOTGeom", pDEBUG) << "Current volume: " << vol->GetName()
1463  << " pos " << fGeometry->GetCurrentPoint()[0]
1464  << " " << fGeometry->GetCurrentPoint()[1]
1465  << " " << fGeometry->GetCurrentPoint()[2]
1466  << " dir " << fGeometry->GetCurrentDirection()[0]
1467  << " " << fGeometry->GetCurrentDirection()[1]
1468  << " " << fGeometry->GetCurrentDirection()[2]
1469  << "[path: " << fGeometry->GetPath() << "]";
1470 #endif
1471 #endif
1472 
1473  // find the start of top
1474  if (fGeometry->IsOutside() || !vol) {
1475  keep_on = false;
1476  if (found_vol) break;
1477  step = 0;
1478  this->StepToNextBoundary();
1479  //rwh//raydist += step; // STNB doesn't actually "step"
1480 
1481 #ifdef RWH_DEBUG
1482 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
1483  LOG("GROOTGeom", pDEBUG) << "Outside ToNextBoundary step: " << step
1484  << " raydist: " << raydist;
1485 #endif
1486 #endif
1487 
1488  while (!fGeometry->IsEntering()) {
1489  step = this->Step();
1490  raydist += step;
1491 #ifdef RWH_DEBUG
1492 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
1493  LOG("GROOTGeom", pDEBUG)
1494  << "Stepping... [step size = " << step << "]";
1495  LOG("GROOTGeom", pDEBUG) << "Outside step: " << step
1496  << " raydist: " << raydist;
1497 #endif
1498 #endif
1499  if (this->WillNeverEnter(step)) {
1500 #ifdef RWH_COUNTVOLS
1501  if ( accum_vol_stat ) {
1502  // this really shouldn't happen for the box exploration...
1503  // if coord transforms done right
1504  // could happen for neutrinos on a flux window
1505  nnever[curface]++; //rwh
1506  if ( nnever[curface]%21 == 0 )
1507  LOG("GROOTGeom", pNOTICE)
1508  << "curface " << curface << " " << nswims[curface]
1509  << " never " << nnever[curface]
1510  << " x [" << r0[0] << "," << r0[1] << "," << r0[2] << "] "
1511  << " p [" << udir[0] << "," << udir[1] << "," << udir[2] << "]";
1512  }
1513 #endif
1515  return;
1516  }
1517  } // finished while
1518 
1519  ps_curr.SetExit(fGeometry->GetCurrentPoint());
1520  ps_curr.SetStep(step);
1521  if ( ( fDebugFlags & 0x10 ) ) {
1522  // In general don't add the path segments from the start point to
1523  // the top volume (here for debug purposes)
1524  // Clear out the step range even if we keep it
1525  ps_curr.fStepRangeSet.clear();
1526  LOG("GROOTGeom", pNOTICE)
1527  << "debug: step towards top volume: " << ps_curr;
1528  fCurrPathSegmentList->AddSegment(ps_curr);
1529  }
1530 
1531  } // outside or !vol
1532 
1533  if (keep_on) {
1534  if (!found_vol) found_vol = true;
1535 
1536  step = this->StepUntilEntering();
1537  raydist += step;
1538 
1539  ps_curr.SetExit(fGeometry->GetCurrentPoint());
1540  ps_curr.SetStep(step);
1541  fCurrPathSegmentList->AddSegment(ps_curr);
1542 
1543  nvolswim++; //rwh
1544 
1545 #ifdef DUMP_SWIM
1546  LOG("GROOTGeom", pDEBUG) << "Current volume: " << vol->GetName()
1547  << " step " << step << " in " << mat->GetName();
1548 #endif
1549 
1550 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
1551  LOG("GROOTGeom", pDEBUG)
1552  << "Cur med.: " << med->GetName() << ", mat.: " << mat->GetName();
1553  LOG("GROOTGeom", pDEBUG)
1554  << "Step = " << step; // << ", weight = " << weight;
1555 #endif
1556  } //keep on
1557  }
1558 
1559 #ifdef RWH_COUNTVOLS
1560  if ( accum_vol_stat ) {
1561  nswims[curface]++; //rwh
1562  dnvols[curface] += (double)nvolswim;
1563  dnvols2[curface] += (double)nvolswim * (double)nvolswim;
1564  long int ns = fCurrPathSegmentList->size();
1565  if ( ns > mxsegments ) mxsegments = ns;
1566  }
1567 #endif
1568 
1569 //rwh:debug
1570 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
1571  LOG("GROOTGeom", pDEBUG)
1572  << "PathSegmentList size " << fCurrPathSegmentList->size();
1573 #endif
1574 
1575 #ifdef RWH_DEBUG_2
1576  if ( ( fDebugFlags & 0x20 ) ) {
1578  LOG("GROOTGeom", pNOTICE) << "Before trimming" << *fCurrPathSegmentList;
1579  double mxddist = 0, mxdstep = 0;
1580  fCurrPathSegmentList->CrossCheck(mxddist,mxdstep);
1581  fmxddist = TMath::Max(fmxddist,mxddist);
1582  fmxdstep = TMath::Max(fmxdstep,mxdstep);
1583  }
1584 #endif
1585 
1586  // PathSegmentList trimming occurs here!
1587  if ( fGeomVolSelector ) {
1588  PathSegmentList* altlist =
1591  delete altlist; // after swap delete original
1592  }
1593 
1595 
1596 #ifdef RWH_DEBUG_2
1597  if ( fGeomVolSelector) {
1598  // after FillMatStepSum() so one can see the summed mass
1599  if ( ( fDebugFlags & 0x40 ) ) {
1601  LOG("GROOTGeom", pNOTICE) << "After trimming" << *fCurrPathSegmentList;
1603  }
1604  }
1605 #endif
1606 
1607 
1608  return;
1609 }
1610 
1611 //___________________________________________________________________________
1613 {
1614  TGeoVolume * vol = fGeometry -> GetCurrentVolume();
1615  if(vol) {
1616  TGeoMaterial * mat = vol->GetMedium()->GetMaterial();
1617  if(mat->IsMixture()) {
1618  TGeoMixture * mixt = dynamic_cast <TGeoMixture*> (mat);
1619  for(int i = 0; i < mixt->GetNelements(); i++) {
1620  int pdg = this->GetTargetPdgCode(mixt, i);
1621  if(tgtpdg == pdg) return true;
1622  }
1623  } else {
1624  int pdg = this->GetTargetPdgCode(mat);
1625  if(tgtpdg == pdg) return true;
1626  }
1627  } else {
1628  LOG("GROOTGeom", pWARN) << "Current volume is null!";
1629  return false;
1630  }
1631  return false;
1632 }
1633 //___________________________________________________________________________
1635 {
1636  fGeometry->FindNextBoundary();
1637  double step=fGeometry->GetStep();
1638  return step;
1639 }
1640 //___________________________________________________________________________
1642 {
1643  fGeometry->Step();
1644  double step=fGeometry->GetStep();
1645  return step;
1646 }
1647 //___________________________________________________________________________
1649 {
1650  this->StepToNextBoundary(); // doesn't actually step, so don't include in sum
1651  double step = 0; //
1652 
1653  while(!fGeometry->IsEntering()) {
1654  step += this->Step();
1655  }
1656 
1657 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
1658 
1659  bool isen = fGeometry->IsEntering();
1660  bool isob = fGeometry->IsOnBoundary();
1661 
1662  LOG("GROOTGeom",pDEBUG)
1663  << "IsEntering = " << utils::print::BoolAsYNString(isen)
1664  << ", IsOnBoundary = " << utils::print::BoolAsYNString(isob)
1665  << ", Step = " << step;
1666 #endif
1667 
1668  return step;
1669 }
1670 //___________________________________________________________________________
1672 {
1673 // If the neutrino trajectory would never enter the detector, then the
1674 // TGeoManager::GetStep returns the maximum step (1E30).
1675 // Compare surrent step with max step and figure out whether the particle
1676 // would never enter the detector
1677 
1678  if (step > 9.99E29) {
1679 
1680 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
1681  LOG("GROOTGeom", pINFO) << "Wow! Current step is dr = " << step;
1682  LOG("GROOTGeom", pINFO) << "This trajectory isn't entering the detector";
1683 #endif
1684  return true;
1685 
1686  } else
1687  return false;
1688 }
1689 
1690 //___________________________________________________________________________
virtual double MaxEnergy(void)=0
declare the max flux neutrino energy that can be generated (for init. purposes)
void ScalePathLength(int pdgc, double scale)
virtual int ScannerNParticles(void) const
virtual void SetMaxPlSafetyFactor(double sf)
TGeoManager * fGeometry
input detector geometry
double PathLength(int pdgc) const
void SetEnter(const TVector3 &p3enter, double raydist)
point of entry to geometry element
const XML_Char * name
Definition: expat.h:151
virtual double GetWeight(const TGeoMaterial *mat, int pdgc)
void SetPrintVerbose(bool doit=true)
string P4AsShortString(const TLorentzVector *p)
Definition: PrintUtils.cxx:52
virtual bool WeightWithDensity(void) const
virtual const PathLengthList & ComputeMaxPathLengths(void)
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
#define pERROR
Definition: Messenger.h:60
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:79
void CrossCheck(double &mxddist, double &mxdstep) const
const Var weight
void SetPath(const char *path)
string BoolAsYNString(bool b)
Definition: PrintUtils.cxx:115
StepRangeSet fStepRangeSet
collection of {steplo,stephi} pairs
bool fDensWeight
if true pathlengths are weighted with density [def:true]
const char * p
Definition: xmltok.h:285
#define pFATAL
Definition: Messenger.h:57
virtual const TLorentzVector & Position(void)=0
returns the flux neutrino 4-position (note: expect SI rather than physical units) ...
void AddSegment(const PathSegment &ps)
GFluxI * fFlux
a flux objects that can be used to scan the max path lengths
static constexpr Double_t nm
Definition: Munits.h:133
string P3AsString(const TVector3 *vec)
Definition: PrintUtils.cxx:76
xmlNodePtr FindNode(xmlDocPtr xml_doc, string node_path)
virtual void Top2Master(TVector3 &v) const
GENIE geometry drivers.
double fLengthScale
conversion factor: input geometry length units -> meters
Object to be filled with the neutrino path-segments representing geometry volume steps (generally bou...
void abs(TH1 *hist)
void SetPathLength(int pdgc, double pl)
bool ExistsInPDGCodeList(int pdg_code) const
string filename
Definition: shutoffs.py:106
int fNRays
max path length scanner (box method): rays/point [def:200]
TVector3 * fCurrVertex
current generated vertex
bool IsSameStart(const TVector3 &pos, const TVector3 &dir) const
double fMaxPlSafetyFactor
factor that can multiply the computed max path lengths
virtual void SetScannerNParticles(int np)
string fTopVolumeName
input top vol [other than TGeoManager::GetTopVolume()]
void SetStep(Double_t step, bool setlimits=true)
step taken in the geometry element
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:30
virtual void SetDensityUnits(double du)
virtual const PDGCodeList & ListOfTargetNuclei(void)
implement the GeomAnalyzerI interface
virtual bool FindMaterialInCurrentVol(int pdgc)
A list of PDG codes.
Definition: PDGCodeList.h:33
virtual double MixtureWeightsSum(void) const
void swap(art::HLTGlobalStatus &lhs, art::HLTGlobalStatus &rhs)
Float_t Z
Definition: plot.C:38
virtual const PathLengthList & ComputePathLengths(const TLorentzVector &x, const TLorentzVector &p)
Object to be filled with the neutrino path-length, for all detector geometry materials, when starting from a position x and travelling along the direction of the neutrino 4-momentum.
const TGeoMaterial * fMaterial
ref only ptr to TGeoMaterial
virtual int GetTargetPdgCode(const TGeoMaterial *const m) const
double foz
top vol size/origin (top vol units)
double fDensityScale
conversion factor: input geometry density units -> kgr/meters^3
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
int fNPoints
max path length scanner (box method): points/surface [def:200]
PathSegmentV_t::const_iterator PathSegVCItr_t
const double emax
virtual bool GenBoxRay(int indx, TLorentzVector &x4, TLorentzVector &p4)
virtual double LengthUnits(void) const
virtual void SetScannerNRays(int nr)
const TGeoVolume * fVolume
ref only ptr to TGeoVolume
PDGCodeList * fCurrPDGCodeList
current list of target nuclei
virtual void SetWeightWithDensity(bool wt)
PathSegmentList * fCurrPathSegmentList
current list of path-segments
PathLengthList * fCurrMaxPathLengthList
current list of max path-lengths
Float_t d
Definition: plot.C:236
Definition: structs.h:1
static const double meter
Definition: Units.h:33
static const double kilogram
Definition: Units.h:34
void SetStartInfo(const TVector3 &pos=TVector3(0, 0, 1e37), const TVector3 &dir=TVector3(0, 0, 0))
const double j
Definition: BetheBloch.cxx:29
virtual void Load(string geometry_filename)
double frac(double x)
Fractional part.
#define pINFO
Definition: Messenger.h:63
Eigen::VectorXd vec
Double_t GetSummedStepRange() const
get the sum of all the step range (in case step has been trimmed or split)
TRandom3 & RndGeom(void) const
rnd number generator used by geometry drivers
Definition: RandomGen.h:69
virtual double ComputePathLengthPDG(const TVector3 &r, const TVector3 &udir, int pdgc)
virtual void Local2SI(PathLengthList &pl) const
access to geometry coordinate/unit transforms for validation/test purposes
GeomVolSelectorI * fGeomVolSelector
optional path seg trimmer (owned)
const PathSegmentV_t & GetPathSegmentV(void) const
static const double ns
Module that plots metrics from reconstructed cosmic ray data.
double fmxdstep
max errors in pathsegmentlist
const ana::Var wgt
Misc GENIE control constants.
#define pWARN
Definition: Messenger.h:61
virtual double DensityUnits(void) const
virtual void Master2TopDir(TVector3 &v) const
virtual void Clear(Option_t *opt)=0
reset state variables based on opt
virtual bool GenerateNext(void)=0
generate the next flux neutrino (return false in err)
virtual const TVector3 & GenerateVertex(const TLorentzVector &x, const TLorentzVector &p, int tgtpdg)
bool GetNeedPath() const
allow toggle on only
MaterialMap_t::const_iterator MaterialMapCItr_t
void SetCurrentRay(const TLorentzVector &x4, const TLorentzVector &p4)
configure for individual neutrino ray
static const double A
Definition: Units.h:82
PathLengthList * fCurrPathLengthList
current list of path-lengths
virtual PathSegmentList * GenerateTrimmedList(const PathSegmentList *untrimmed) const
virtual void SetTopVolName(string nm)
void SetDoCrossCheck(bool doit=true)
TGeoHMatrix * fMasterToTop
matrix connecting master coordinates to top volume coordinates
bool fMasterToTopIsIdentity
is fMasterToTop matrix the identity matrix?
void SetGeo(const TGeoVolume *gvol, const TGeoMedium *gmed, const TGeoMaterial *gmat)
info about the geometry element
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:326
std::map< const TGeoMaterial *, Double_t > MaterialMap_t
TVector3 GetPosition(Double_t frac) const
calculate position within allowed ranges passed on fraction of total
virtual void SetScannerNPoints(int np)
set geometry driver&#39;s configuration options
virtual void SetScannerFlux(GFluxI *f)
const MaterialMap_t & GetMatStepSumMap(void) const
virtual void SetMixtureWeightsSum(double sum)
static const double meter3
Definition: Units.h:59
int IonPdgCode(int A, int Z)
Definition: PDGUtils.cxx:69
exit(0)
virtual const TLorentzVector & Momentum(void)=0
returns the flux neutrino 4-momentum
void SetExit(const TVector3 &p3exit)
point of exit from geometry element
TGeoVolume * fTopVolume
top volume
assert(nhit_max >=nhit_nbins)
virtual double MaxPlSafetyFactor(void) const
virtual void SetLengthUnits(double lu)
void AddPathLength(int pdgc, double pl)
virtual void Top2MasterDir(TVector3 &v) const
double fMixtWghtSum
norm of relative weights (<0 if explicit summing required)
#define pNOTICE
Definition: Messenger.h:62
Defines the GENIE Geometry Analyzer Interface.
Definition: GeomAnalyzerI.h:30
string X4AsString(const TLorentzVector *x)
Definition: PrintUtils.cxx:64
Double_t sum
Definition: plot.C:31
virtual void SI2Local(TVector3 &v) const
bool fKeepSegPath
need to fill path segment "path"
Float_t w
Definition: plot.C:20
std::list< PathSegment > PathSegmentV_t
void next()
Definition: show_event.C:84
string Vec3AsString(const TVector3 *vec)
Definition: PrintUtils.cxx:87
TGeoManager * gm
Definition: make_fe_box.C:4
void SetSI2Local(double scale)
set scale factor for SI to "raydist" units of PathSegmentList
virtual void SwimOnce(const TVector3 &r, const TVector3 &udir)
TGeoVolume * top
Definition: make_fe_box.C:9
virtual void Master2Top(TVector3 &v) const
void push_back(int pdg_code)
Definition: PDGCodeList.cxx:67
Eigen::MatrixXd mat
#define pDEBUG
Definition: Messenger.h:64
virtual bool WillNeverEnter(double step)