GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GeomVolSelectorRockBox.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  Robert Hatcher <rhatcher@fnal.gov>
7 */
8 //____________________________________________________________________________
9 
13 
14 using namespace genie;
15 using namespace genie::geometry;
16 
17 #include <TGeoVolume.h>
18 #include <TGeoMaterial.h>
19 #include <TGeoMedium.h>
20 
21 //____________________________________________________________________________
22 GeomVolSelectorRockBox::GeomVolSelectorRockBox()
23  : GeomVolSelectorFiducial(), fMinimumWall(0.), fDeDx(1.)
24  , fExpandInclusion(false)
25  , fRockBoxShape(0), fROOTGeom(0)
26 {
27  fName = "RockBox";
28  // base class' fiducial volume always treated as reverse (if even exists)
29  this->SetReverseFiducial(true);
30 
31  for (int i=0; i<3; ++i) {
32  fMinimalXYZMin[i] = 0;
33  fMinimalXYZMax[i] = 0;
34  fInclusionXYZMin[i] = 0;
35  fInclusionXYZMax[i] = 0;
36  }
37 }
38 
39 //___________________________________________________________________________
41 {
42  if ( fRockBoxShape ) delete fRockBoxShape;
43  fRockBoxShape = 0;
44  fROOTGeom = 0; // was reference only
45 }
46 
47 //___________________________________________________________________________
49 {
50  // First trim the segment based on the ray vs. cylinder or box
51  // Then trim futher according to the Basic parameters
52 
53  if ( ! fInterceptRock.fIsHit ) {
54  // want in rock box, ray misses => reject all segments
55  ps.fStepRangeSet.clear(); //
56  } else {
57  // ray hit rock box volume, some segments steps need rejection, some need splitting...
58  // check the steps in this segment
59  Double_t dist = ps.fRayDist;
60  StepRangeSet::iterator srs_itr = ps.fStepRangeSet.begin();
61  StepRangeSet::iterator srs_end = ps.fStepRangeSet.end();
62  StepRangeSet modifiedStepRangeSet;
63  Bool_t ismod = false;
64 
65  // loop over steps within this segement
66  for ( ; srs_itr != srs_end; ++srs_itr ) {
67  Double_t slo = srs_itr->first;
68  Double_t shi = srs_itr->second;
69  Bool_t split = false;
70  StepRange step1, step2;
71  // determine new trimmed or split steps
72  ismod |= NewStepPairs(false,dist,slo,shi,
73  fInterceptRock,split,step1,step2);
74  // build up new step list
75  bool nonzerostep = ( step1.first != step1.second );
76  if ( nonzerostep || ! fRemoveEntries ) {
77  modifiedStepRangeSet.push_back(step1);
78  if (split) {
79  modifiedStepRangeSet.push_back(step2);
80  }
81  }
82  } // loop over step range set elements
83  if ( ismod ) ps.fStepRangeSet = modifiedStepRangeSet;
84 
85  } // fIsHit
86 
88 }
89 
90 //___________________________________________________________________________
92 {
93  // A new neutrino ray has been set, calculate the entrance/exit distances.
94 
96 
97  fCurrPathSegmentList = untrimmed;
98 
99  MakeRockBox();
100 
101  if ( ! fRockBoxShape ) {
102  LOG("GeomVolSel", pFATAL) << "no shape defined";
104  } else {
108  }
109 
110  //cout << "BeginPSList: " << endl
111  // << " fid: " << fIntercept << endl
112  // << " rock: " << fInterceptRock << endl;
113 
114 }
115 
116 //___________________________________________________________________________
118 {
119  // Completed current path segment list processsing
120 }
121 
122 //___________________________________________________________________________
124 {
126  fROOTGeom = rgeom; // stash away a copy
127 }
128 //___________________________________________________________________________
130  Double_t* xyzmax)
131 {
132  // This sets parameters for a minimal box
133 
134  for ( int j = 0; j < 3; ++j ) {
135  fMinimalXYZMin[j] = TMath::Min(xyzmin[j],xyzmax[j]);
136  fMinimalXYZMax[j] = TMath::Max(xyzmax[j],xyzmax[j]);
137  }
138  // create the default inner (exclusion) box
140 }
141 //___________________________________________________________________________
143  Double_t* xyzmax)
144 {
145  // This sets parameters for the inclusion (outer) box
146 
147  for ( int j = 0; j < 3; ++j ) {
148  fInclusionXYZMin[j] = TMath::Min(xyzmin[j],xyzmax[j]);
149  fInclusionXYZMax[j] = TMath::Max(xyzmax[j],xyzmax[j]);
150  }
151 }
152 //___________________________________________________________________________
154 {
155  fMinimumWall = w;
156  for ( int j = 0; j < 3; ++j ) {
159  }
160 }
161 //___________________________________________________________________________
163 {
164  // This sets parameters for a box
165 
166  // expanded box
167  double energy = fP4.Energy();
168  double boxXYZMin[3], boxXYZMax[3];
169  for ( int j = 0; j < 3; ++j ) {
170  double dmin = 0, dmax = 0;
171  double dircos = fCurrPathSegmentList->GetDirection()[j];
172  if ( dircos > 0 ) dmin = dircos*energy/fDeDx; // pad upstream
173  else dmax = -dircos*energy/fDeDx;
174 
175 //#define RWH_DEBUG
176 #ifdef RWH_DEBUG
177  cout << "MakeRockBox [" << j << "] wall " << fMinimumWall
178  << " dm " << dmin << " " << dmax << " dircos " << dircos
179  << " e " << energy << " dedx " << fDeDx << endl;
180 #endif
181 
182  if ( fExpandInclusion ) {
183  boxXYZMin[j] = fInclusionXYZMin[j] - dmin;
184  boxXYZMax[j] = fInclusionXYZMax[j] + dmax;
185  } else {
186  boxXYZMin[j] = TMath::Min(fMinimalXYZMin[j]-dmin,fInclusionXYZMin[j]);
187  boxXYZMax[j] = TMath::Max(fMinimalXYZMax[j]+dmin,fInclusionXYZMax[j]);
188  }
189  }
190 
191  FidPolyhedron* poly = new FidPolyhedron();
192  // careful about sign of "d" vs. direction normal
193  PlaneParam pln0(-1,0,0, boxXYZMin[0]); poly->push_back(pln0);
194  PlaneParam pln1(0,-1,0, boxXYZMin[1]); poly->push_back(pln1);
195  PlaneParam pln2(0,0,-1, boxXYZMin[2]); poly->push_back(pln2);
196  PlaneParam pln3(+1,0,0,-boxXYZMax[0]); poly->push_back(pln3);
197  PlaneParam pln4(0,+1,0,-boxXYZMax[1]); poly->push_back(pln4);
198  PlaneParam pln5(0,0,+1,-boxXYZMax[2]); poly->push_back(pln5);
199 
200  if ( fRockBoxShape ) delete fRockBoxShape;
201  fRockBoxShape = poly;
202 
203 #ifdef RWH_DEBUG
204  static bool first = true;
205  if ( first ) {
206  cout << "MakeRockBox first Minimal min ["
207  << fMinimalXYZMin[0] << ","
208  << fMinimalXYZMin[1] << ","
209  << fMinimalXYZMin[2] << "]"
210  << " max ["
211  << fMinimalXYZMax[0] << ","
212  << fMinimalXYZMax[1] << ","
213  << fMinimalXYZMax[2] << "]" << endl;
214  cout << "MakeRockBox first Inclusion min ["
215  << fInclusionXYZMin[0] << ","
216  << fInclusionXYZMin[1] << ","
217  << fInclusionXYZMin[2] << "]"
218  << " max ["
219  << fInclusionXYZMax[0] << ","
220  << fInclusionXYZMax[1] << ","
221  << fInclusionXYZMax[2] << "]" << endl;
222  first = false;
223  }
224  cout << "MakeRockBox this ray using min ["
225  << boxXYZMin[0] << ","
226  << boxXYZMin[1] << ","
227  << boxXYZMin[2] << "]"
228  << " max ["
229  << boxXYZMax[0] << ","
230  << boxXYZMax[1] << ","
231  << boxXYZMax[2] << "]" << endl;
232  cout << "rock before:" << *fRockBoxShape << endl;
233 #endif
234 
236 
237 #ifdef RWH_DEBUG
238  cout << "rock after: " << *fRockBoxShape << endl;
239  cout << "fid after: " << *fShape << endl;
240 #endif
241 
242 }
243 
244 //___________________________________________________________________________
GENIE Interface for user-defined volume selector functors Trim path segments based on the intersectio...
void TrimSegment(PathSegment &segment) const
const TVector3 & GetStartPos() const
void TrimSegment(PathSegment &segment) const
TLorentzVector fP4
current neutrino ray&#39;s momentum (global)
std::vector< StepRange > StepRangeSet
void SetRockBoxInclusion(Double_t *xyzmin, Double_t *xyzmax)
Bool_t fIsHit
distance along ray to exit fid volume
Definition: FidShape.h:52
StepRangeSet fStepRangeSet
collection of {steplo,stephi} pairs
#define pFATAL
Definition: Messenger.h:56
virtual void ConvertShapeMaster2Top(const ROOTGeomAnalyzer *rgeom)
Object to be filled with the neutrino path-segments representing geometry volume steps (generally bou...
const TVector3 & GetDirection() const
Double_t fInclusionXYZMin[3]
minimum distance around (XYZmin,XYZmax)
Double_t fMinimumWall
interior box upper corner
std::string fName
volume selector name
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
bool fRemoveEntries
whether selector should remove entries or set hi=lo
FidShape * fRockBoxShape
expand from minimal or inclusion box?
Double_t fInclusionXYZMax[3]
box within which events are always
A ROOT/GEANT4 geometry driver.
void BeginPSList(const PathSegmentList *untrimmed) const
std::pair< Double_t, Double_t > StepRange
static constexpr double ps
Definition: Units.h:99
const ROOTGeomAnalyzer * fROOTGeom
shape changes for every nu ray
Double_t fRayDist
distance from start of ray
virtual RayIntercept Intercept(const TVector3 &start, const TVector3 &dir) const =0
virtual void ConvertShapeMaster2Top(const ROOTGeomAnalyzer *rgeom)
virtual void ConvertMaster2Top(const ROOTGeomAnalyzer *rgeom)=0
static Bool_t NewStepPairs(Bool_t selectReverse, Double_t raydist, Double_t slo, Double_t shi, const RayIntercept &intercept, Bool_t &split, StepRange &step1, StepRange &step2)
void BeginPSList(const PathSegmentList *untrimmed) const
Bool_t fExpandInclusion
how to scale from energy to distance
void SetRockBoxMinimal(Double_t *xyzmin, Double_t *xyzmax)
void MakeBox(Double_t *xyzmin, Double_t *xyzmax)
FidShape * fShape
select for &quot;outside&quot; fiducial?
void push_back(const PlaneParam &pln)
Definition: FidShape.h:140
Double_t fMinimalXYZMax[3]
interior box lower corner
const PathSegmentList * fCurrPathSegmentList
shape