GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GeomVolSelectorFiducial.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 GeomVolSelectorFiducial::GeomVolSelectorFiducial()
23  : GeomVolSelectorBasic(), fSelectReverse(false)
24  , fShape(0), fCurrPathSegmentList(0)
25 {
26  fName = "Fiducial";
27 }
28 
29 //___________________________________________________________________________
31 {
32  if ( fShape ) delete fShape;
33  fShape = 0;
34  fCurrPathSegmentList = 0; // was reference only
35 }
36 
37 //___________________________________________________________________________
39 {
40  // First trim the segment based on the ray vs. cylinder or box
41  // Then trim futher according to the Basic parameters
42 
43  if ( ! fIntercept.fIsHit ) {
44  // simple case when ray doesn't intersect the fiducial volume at all
45  if ( fSelectReverse ) {
46  // fiducial is reversed (ie. only want regions outside)
47  /// so a miss means blindly accept all segments
48  } else {
49  // want in fiducial, ray misses => reject all segments
50  ps.fStepRangeSet.clear(); //
51  }
52  } else {
53  // ray hit fiducial volume, some segments steps need rejection, some need splitting...
54  // check the steps in this segment
55  Double_t dist = ps.fRayDist;
56  StepRangeSet::iterator srs_itr = ps.fStepRangeSet.begin();
57  StepRangeSet::iterator srs_end = ps.fStepRangeSet.end();
58  StepRangeSet modifiedStepRangeSet;
59  Bool_t ismod = false;
60 
61  // loop over steps within this segement
62  for ( ; srs_itr != srs_end; ++srs_itr ) {
63  Double_t slo = srs_itr->first;
64  Double_t shi = srs_itr->second;
65  Bool_t split = false;
66  StepRange step1, step2;
67  // determine new trimmed or split steps
68  ismod |= NewStepPairs(fSelectReverse,dist,slo,shi,
69  fIntercept,split,step1,step2);
70  // build up new step list
71  bool nonzerostep = ( step1.first != step1.second );
72  if ( nonzerostep || ! fRemoveEntries ) {
73  modifiedStepRangeSet.push_back(step1);
74  if (split) {
75  modifiedStepRangeSet.push_back(step2);
76  }
77  }
78  } // loop over step range set elements
79  if ( ismod ) ps.fStepRangeSet = modifiedStepRangeSet;
80 
81  } // fIsHit
82 
84 }
85 
86 //___________________________________________________________________________
87 Bool_t GeomVolSelectorFiducial::NewStepPairs(Bool_t selectReverse, Double_t raydist,
88  Double_t slo, Double_t shi,
89  const RayIntercept& intercept,
90  Bool_t& split,
91  StepRange& step1, StepRange& step2)
92 {
93  // modifying a step based on a range
94  // there seem to be six possible cases:
95  // step: |===| {slo,shi}
96  //
97  // range: in out
98  // # # normal reverse
99  // 0 |===| # # {0,0} {slo,shi}
100  // 1 |==#==| # {in,shi} {slo,in}
101  // 2 |==#========#==| {in,out} {slo,in}+{out,shi}
102  // 3 # |==| # {slo,shi} {0,0}
103  // 4 # |==#==| {slo,out} {out,shi}
104  // 5 # # |====| {0,0} {slo,shi}
105  // # #
106  //
107  bool ismodified = true;
108  step1 = StepRange(slo,shi);
109  split = false;
110  // dist in/out are relative to the ray origin
111  Double_t sdistin = intercept.fDistIn - raydist;
112  Double_t sdistout = intercept.fDistOut - raydist;
113  // do remaining calculations relative to steps within this segment
114  if ( slo < sdistin ) {
115  if ( shi < sdistin ) {
116  // case 0
117  if ( selectReverse ) { ismodified = false; }
118  else { step1 = StepRange(0,0); }
119  } else if ( shi < sdistout ) {
120  // case 1
121  if ( selectReverse ) { step1.second = sdistin; }
122  else { step1.first = sdistin; }
123  } else {
124  // case 2
125  if ( selectReverse ) { split = true;
126  step1 = StepRange(slo,sdistin);
127  step2 = StepRange(sdistout,shi); }
128  else { step1 = StepRange(sdistin,sdistout); }
129  }
130  } else if ( slo < sdistout ) {
131  if ( shi < sdistout ) {
132  // case 3
133  if ( selectReverse ) { step1 = StepRange(0,0); }
134  else { ismodified = false; }
135  } else {
136  // case 4
137  if ( selectReverse ) { step1.first = sdistout; }
138  else { step1.second = sdistout; }
139  }
140  } else {
141  // case 5
142  if ( selectReverse ) { ismodified = false; }
143  else { step1 = StepRange(0,0); }
144  }
145 
146  return ismodified;
147 }
148 
149 //___________________________________________________________________________
151 {
152  // A new neutrino ray has been set, calculate the entrance/exit distances.
153 
154  GeomVolSelectorBasic::BeginPSList(untrimmed); // initialize base class in case it is needed
155 
156  fCurrPathSegmentList = untrimmed;
157 
158  if ( ! fShape ) {
159  LOG("GeomVolSel", pFATAL) << "no shape defined";
161  } else {
164  }
165 
166 }
167 
168 //___________________________________________________________________________
170 {
171  // Completed current path segment list processsing
174 }
175 
176 //___________________________________________________________________________
178 {
179  if ( fShape ) delete fShape;
180  fShape = shape;
181 }
182 //___________________________________________________________________________
184 {
185  if ( fShape ) fShape->ConvertMaster2Top(rgeom);
186 }
187 //___________________________________________________________________________
188 void GeomVolSelectorFiducial::MakeSphere(Double_t x0, Double_t y0, Double_t z0, Double_t radius)
189 {
190  AdoptFidShape(new FidSphere(TVector3(x0,y0,z0),radius));
191 }
192 
193 //___________________________________________________________________________
194 void GeomVolSelectorFiducial::MakeXCylinder(Double_t y0, Double_t z0, Double_t radius,
195  Double_t xmin, Double_t xmax)
196 {
197  // This sets parameters for a cylinder parallel to the x-axis
198  Double_t base[] = { 0, y0, z0 };
199  Double_t axis[] = { 1, 0, 0 };
200  Double_t cap1[] = { -1, 0, 0, xmin };
201  Double_t cap2[] = { +1, 0, 0, -xmax }; // note sign change
202  this->MakeCylinder(base,axis,radius,cap1,cap2);
203 }
204 
205 //___________________________________________________________________________
206 void GeomVolSelectorFiducial::MakeYCylinder(Double_t x0, Double_t z0, Double_t radius,
207  Double_t ymin, Double_t ymax)
208 {
209  // This sets parameters for a cylinder parallel to the y-axis
210  Double_t base[] = { x0, 0, z0 };
211  Double_t axis[] = { 0, 1, 0 };
212  Double_t cap1[] = { 0, -1, 0, ymin };
213  Double_t cap2[] = { 0, +1, 0, -ymax }; // note sign change
214  this->MakeCylinder(base,axis,radius,cap1,cap2);
215 }
216 
217 //___________________________________________________________________________
218 void GeomVolSelectorFiducial::MakeZCylinder(Double_t x0, Double_t y0, Double_t radius,
219  Double_t zmin, Double_t zmax)
220 {
221  // This sets parameters for a cylinder parallel to the z-axis
222  Double_t base[] = { x0, y0, 0 };
223  Double_t axis[] = { 0, 0, 1 };
224  Double_t cap1[] = { 0, 0, -1, zmin };
225  Double_t cap2[] = { 0, 0, +1, -zmax }; // note sign change
226  this->MakeCylinder(base,axis,radius,cap1,cap2);
227 }
228 
229 //___________________________________________________________________________
230 void GeomVolSelectorFiducial::MakeCylinder(Double_t* base, Double_t* axis, Double_t radius,
231  Double_t* cap1, Double_t* cap2)
232 {
233  // This sets parameters for a cylinder
234  AdoptFidShape(new FidCylinder(TVector3(base),TVector3(axis),radius,
235  PlaneParam(cap1),PlaneParam(cap2)));
236 }
237 
238 //___________________________________________________________________________
239 void GeomVolSelectorFiducial::MakeBox(Double_t* xyzmin, Double_t* xyzmax)
240 {
241  // This sets parameters for a box
242 
243  double boxXYZmin[3], boxXYZmax[3];
244  for ( int j = 0; j < 3; ++j ) {
245  boxXYZmin[j] = TMath::Min(xyzmin[j],xyzmax[j]);
246  boxXYZmax[j] = TMath::Max(xyzmin[j],xyzmax[j]);
247  }
248 
249  FidPolyhedron* poly = new FidPolyhedron();
250  // careful about sign of "d" vs. direction normal
251  PlaneParam pln0(-1,0,0, boxXYZmin[0]); poly->push_back(pln0);
252  PlaneParam pln1(0,-1,0, boxXYZmin[1]); poly->push_back(pln1);
253  PlaneParam pln2(0,0,-1, boxXYZmin[2]); poly->push_back(pln2);
254  PlaneParam pln3(+1,0,0,-boxXYZmax[0]); poly->push_back(pln3);
255  PlaneParam pln4(0,+1,0,-boxXYZmax[1]); poly->push_back(pln4);
256  PlaneParam pln5(0,0,+1,-boxXYZmax[2]); poly->push_back(pln5);
257 
258  AdoptFidShape(poly);
259 }
260 
261 //___________________________________________________________________________
262 void GeomVolSelectorFiducial::MakeZPolygon(Int_t n, Double_t x0, Double_t y0, Double_t inradius,
263  Double_t phi0deg, Double_t zmin, Double_t zmax)
264 {
265  // phi=0 will put flat face towards +x
266  // inscribed radius
267 
268  FidPolyhedron* poly = new FidPolyhedron();
269 
270  Double_t dphir = TMath::TwoPi()/n;
271  Double_t phi0r = phi0deg * TMath::DegToRad();
272  for ( int iface = 0; iface < n; ++iface ) {
273  Double_t dx = TMath::Cos(dphir*iface + phi0r);
274  Double_t dy = TMath::Sin(dphir*iface + phi0r);
275  PlaneParam face(dx,dy,0,-inradius-dx*x0-dy*y0); poly->push_back(face);
276  }
277  PlaneParam plnz1(0,0,-1, zmin); poly->push_back(plnz1);
278  PlaneParam plnz2(0,0,+1,-zmax); poly->push_back(plnz2);
279 
280  AdoptFidShape(poly);
281 }
282 
283 //___________________________________________________________________________
Some simple volumes that know how to calculate where a ray intercepts them.
Definition: FidShape.h:90
void MakeSphere(Double_t x0, Double_t y0, Double_t z0, Double_t radius)
const TVector3 & GetStartPos() const
void TrimSegment(PathSegment &segment) const
std::vector< StepRange > StepRangeSet
void MakeZCylinder(Double_t x0, Double_t y0, Double_t radius, Double_t zmin, Double_t zmax)
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
void MakeCylinder(Double_t *base, Double_t *axis, Double_t radius, Double_t *cap1, Double_t *cap2)
void MakeXCylinder(Double_t y0, Double_t z0, Double_t radius, Double_t xmin, Double_t xmax)
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
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
void MakeYCylinder(Double_t x0, Double_t z0, Double_t radius, Double_t ymin, Double_t ymax)
bool fRemoveEntries
whether selector should remove entries or set hi=lo
void MakeZPolygon(Int_t n, Double_t x0, Double_t y0, Double_t inradius, Double_t phi0deg, Double_t zmin, Double_t zmax)
A ROOT/GEANT4 geometry driver.
GENIE Interface for user-defined volume selector functors This basic version allows configurations th...
std::pair< Double_t, Double_t > StepRange
static constexpr double ps
Definition: Units.h:99
Double_t fRayDist
distance from start of ray
virtual RayIntercept Intercept(const TVector3 &start, const TVector3 &dir) const =0
virtual void ConvertMaster2Top(const ROOTGeomAnalyzer *rgeom)=0
Double_t fDistOut
distance along ray to enter fid volume
Definition: FidShape.h:51
void TrimSegment(PathSegment &segment) const
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
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
void BeginPSList(const PathSegmentList *untrimmed) const
const PathSegmentList * fCurrPathSegmentList
shape