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