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