GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
FidShape.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 
11 #include <iomanip>
12 using namespace std;
13 
16 
17 using namespace genie;
18 using namespace genie::geometry;
19 
20 //___________________________________________________________________________
21 std::ostream&
22 genie::geometry::operator<< (std::ostream& stream,
23  const genie::geometry::PlaneParam& pparam)
24 {
25  pparam.Print(stream);
26  return stream;
27 }
28 
29 std::ostream&
30 genie::geometry::operator<< (std::ostream& stream,
32 {
33  stream << "RayIntercept: dist in/out " << ri.fDistIn << "/" << ri.fDistOut
34  << " hit=" << ((ri.fIsHit)?"true":"false")
35  << " surf " << ri.fSurfIn << "/" << ri.fSurfOut;
36  return stream;
37 }
38 
39 std::ostream&
40 genie::geometry::operator<< (std::ostream& stream,
41  const genie::geometry::FidShape& shape)
42 {
43  shape.Print(stream);
44  return stream;
45 }
46 
47 //___________________________________________________________________________
48 void PlaneParam::ConvertMaster2Top(const ROOTGeomAnalyzer* rgeom)
49 {
50  // convert a plane equation from master coordinates to "top vol" coordinates
51  TVector3 dir(a,b,c);
52  rgeom->Master2TopDir(dir); // new a,b,c
53  TVector3 zero(0,0,0);
54  rgeom->Master2Top(zero);
55  a = dir.X();
56  b = dir.Y();
57  c = dir.Z();
58  d = d - ( a*zero.X() + b*zero.Y() + c*zero.Z() );
59 
60 }
61 
62 //___________________________________________________________________________
63 void PlaneParam::Print(std::ostream& stream) const
64 {
65  stream << "PlaneParam=[" << a << "," << b << "," << c << "," << d << "]";
66 }
67 
68 //___________________________________________________________________________
69 RayIntercept FidSphere::Intercept(const TVector3& start, const TVector3& dir) const
70 {
71  // A new neutrino ray has been set, calculate the entrance/exit distances.
72  // This sets fDistIn/fDistOut for an sphere
73  RayIntercept intercept;
74 
75  TVector3 oc = fCenter - start;
76  Double_t loc2 = oc.Mag2();
77  Double_t r2 = fSRadius*fSRadius;
78  //LOG("GeomVolSel", pNOTICE) << " loc2 = " << loc2 << " r2 " << r2;
79  // if ( loc2 > r2 ) ray originates outside the sphere
80  const TVector3& d = dir;
81  Double_t d2 = d.Mag2();
82  Double_t tca = oc.Dot(d)/d2;
83  //if ( tca < 0.0 ) sphere _center_ behind the ray orgin
84  Double_t lhc2 = ( r2 -loc2 )/d2 + tca*tca;
85  if ( lhc2 < 0.0 ) return intercept; // ray misses the sphere
86  intercept.fIsHit = true;
87  Double_t lhc = TMath::Sqrt(lhc2);
88 
89  intercept.fDistIn = tca - lhc;
90  intercept.fSurfIn = 1;
91  intercept.fDistOut = tca + lhc;
92  intercept.fSurfOut = 1;
93 
94  return intercept;
95 }
96 
97 //___________________________________________________________________________
98 void FidSphere::ConvertMaster2Top(const ROOTGeomAnalyzer* rgeom)
99 {
100  rgeom->Master2Top(fCenter);
101 }
102 
103 //___________________________________________________________________________
104 void FidSphere::Print(std::ostream& stream) const
105 {
106  stream << "FidSphere @ ["
107  << fCenter.X() << ","
108  << fCenter.Y() << ","
109  << fCenter.Z() << "] r = " << fSRadius;
110 }
111 
112 //___________________________________________________________________________
113 RayIntercept FidCylinder::InterceptUncapped(const TVector3& start, const TVector3& dir) const
114 {
115  // A new neutrino ray has been set, calculate the entrance/exit distances.
116  // This sets fDistIn/fDistOut for an infinite cylinder
117  // Take as "hit" if the ray is parallel to the axis but inside the radius
118  RayIntercept intercept;
119 
120  TVector3 rc = start - fCylBase;
121  TVector3 n = dir.Cross(fCylAxis);
122  Double_t len = n.Mag();
123  Double_t dist = 0;
124  if ( len == 0.0 ) {
125  // ray is parallel to axis
126  dist = rc.Dot(fCylAxis);
127  TVector3 d = rc - dist*fCylAxis;
128  dist = d.Mag();
129  //LOG("GeomVolSel", pNOTICE) << " len = " << len << " is parallel, dist " << dist << ", " << fCylRadius;
130  if ( dist <= fCylRadius ) {
131  intercept.fIsHit = true; // inside is considered a hit
132  intercept.fSurfIn = 0;
133  intercept.fSurfOut = 0;
134  }
135  return intercept;
136  }
137  // ray is not parallel
138  if ( len != 1.0 ) n.SetMag(1.); // normalize if it isn't already
139  dist = TMath::Abs(rc.Dot(n)); // closest approach distance
140  //LOG("GeomVolSel", pNOTICE) << " len = " << len << " not parallel, dist " << dist << ", " << fCylRadius;
141  if ( dist <= fCylRadius ) {
142  intercept.fIsHit = true; // yes, it hits
143  intercept.fSurfIn = 0;
144  intercept.fSurfOut = 0;
145  TVector3 o = rc.Cross(fCylAxis);
146  Double_t t = - o.Dot(n)/len;
147  o = n.Cross(fCylAxis);
148  o.SetMag(1.);
149  Double_t s = TMath::Abs( TMath::Sqrt(fCylRadius*fCylRadius-dist*dist) /
150  dir.Dot(o) );
151  intercept.fDistIn = t - s;
152  intercept.fDistOut = t + s;
153  //LOG("GeomVolSel", pNOTICE) << " hits, t = " << t << " s = " << s;
154  }
155  return intercept;
156 }
157 
158 //___________________________________________________________________________
159 RayIntercept FidCylinder::Intercept(const TVector3& start, const TVector3& dir) const
160 {
161  // A new neutrino ray has been set, calculate the entrance/exit distances.
162 
163  RayIntercept intercept = InterceptUncapped(start,dir); // find where ray hits the infinite cylinder
164  // trim this down by applying the end caps
165  if ( ! intercept.fIsHit ) return intercept;
166  for ( int icap=1; icap <= 2; ++icap ) {
167  const PlaneParam& cap = (icap==1) ? fCylCap1 : fCylCap2;
168  if ( ! cap.IsValid() ) continue;
169  Double_t vd = cap.Vd(dir);
170  Double_t vn = cap.Vn(start);
171  //std::cout << "FidCyl::Intercept cap " << icap
172  // << " vd " << vd << " vn " << vn;
173  if ( vd == 0.0 ) { // parallel to surface, is it on the right side?
174  //std::cout << " vd=0, vn " << ((vn>0)?"wrong":"right") << "side " << std::endl;
175  if ( vn > 0 ) { intercept.fIsHit = false; break; } // wrong side
176  } else {
177  Double_t t = -vn / vd;
178  //std::cout << " t " << t << " in/out "
179  // << intercept.fDistIn << "/" << intercept.fDistOut << std::endl;
180  if ( vd < 0.0 ) { // t is the entering point
181  if ( t > intercept.fDistIn )
182  { intercept.fDistIn = t; intercept.fSurfIn = 1; }
183  } else { // t is the exiting point
184  if ( t < intercept.fDistOut )
185  { intercept.fDistOut = t; intercept.fSurfOut = 1; }
186  }
187  }
188  }
189  if ( intercept.fDistIn > intercept.fDistOut ) intercept.fIsHit = false;
190  return intercept;
191 }
192 
193 //___________________________________________________________________________
194 void FidCylinder::ConvertMaster2Top(const ROOTGeomAnalyzer* rgeom)
195 {
196  rgeom->Master2Top(fCylBase);
197  rgeom->Master2TopDir(fCylAxis);
198  fCylCap1.ConvertMaster2Top(rgeom);
199  fCylCap2.ConvertMaster2Top(rgeom);
200 }
201 
202 //___________________________________________________________________________
203 void FidCylinder::Print(std::ostream& stream) const
204 {
205  stream << "FidCylinder @ ["
206  << fCylBase.X() << ","
207  << fCylBase.Y() << ","
208  << fCylBase.Z() << "] dir ["
209  << fCylAxis.X() << ","
210  << fCylAxis.Y() << ","
211  << fCylAxis.Z() << "] r = " << fCylRadius;
212  stream << " cap1=" << fCylCap1 << " cap2=" << fCylCap2;
213 }
214 
215 //___________________________________________________________________________
216 RayIntercept FidPolyhedron::Intercept(const TVector3& start, const TVector3& dir) const
217 {
218  // A new neutrino ray has been set, calculate the entrance/exit distances.
219  // Calculate the point of intersection of a ray (directed line) and a
220  // *convex* polyhedron constructed from the intersection of a list
221  // of planar equations (no check on convex condition).
222 
223  RayIntercept intercept;
224 
225  Double_t tnear = -DBL_MAX;
226  Double_t tfar = DBL_MAX;
227  Int_t surfNear = -1;
228  Int_t surfFar = -1;
229  Bool_t parallel = false;
230 
231  // test each plane in the polyhedron
232  for ( size_t iface=0; iface < fPolyFaces.size(); ++iface ) {
233  const PlaneParam& pln = fPolyFaces[iface];
234  if ( ! pln.IsValid() ) continue;
235 
236  // calculate numerator, denominator to "t" = distance along ray to intersection w/ pln
237  Double_t vd = pln.Vd(dir);
238  Double_t vn = pln.Vn(start);
239 
240  //LOG("GeomVolSel", pNOTICE)
241  // << " face " << iface << " [" << pln.a << "," << pln.b << "," << pln.c << "," << pln.d
242  // << "] vd=" << vd << " vn=" << vn;
243 
244  if ( vd == 0.0 ) {
245  // ray is parallel to plane - check if ray origin is inside plane's half-space
246  //LOG("GeomVolSel", pNOTICE)
247  // << " vd=0, " << " possibly parallel ";
248  if ( vn > 0.0 ) { parallel = true; break; } // wrong side ... complete miss
249  } else {
250  // ray is not parallel to plane -- get the distance to the plane
251  Double_t t = -vn / vd; // notice negative sign!
252  //LOG("GeomVolSel", pNOTICE) << " t=" << t << " tnear=" << tnear << " tfar=" << tfar;
253  if ( vd < 0.0 ) {
254  // front face: t is a near point
255  if ( t > tnear ) {
256  surfNear = iface;
257  tnear = t;
258  }
259  } else {
260  // back face: t is a far point
261  if ( t < tfar ) {
262  surfFar = iface;
263  tfar = t;
264  }
265  }
266  //LOG("GeomVolSel", pNOTICE) << " new surf " << surfNear << "," << surfFar
267  // << " tnear=" << tnear << " tfar=" << tfar;
268  }
269  }
270  if ( ! parallel ) {
271  // survived all the tests
272  if ( tnear > 0.0 ) {
273  if ( tnear < tfar ) {
274  //LOG("GeomVolSel", pNOTICE) << "is hit case1 ";
275  intercept.fIsHit = true;
276  intercept.fSurfIn = surfNear;
277  intercept.fSurfOut = surfFar;
278  }
279  } else {
280  if ( tfar > 0.0 ) {
281  //LOG("GeomVolSel", pNOTICE) << "is hit case2 ";
282  intercept.fIsHit = true;
283  intercept.fSurfIn = -1;
284  intercept.fSurfOut = surfFar;
285  }
286  }
287  }
288  intercept.fDistIn = tnear;
289  intercept.fDistOut = tfar;
290  //LOG("GeomVolSel", pNOTICE) << " hit? " << (fIsHit?"true":"false")
291  // << " dist in " << fDistIn << " out " << fDistOut;
292  return intercept;
293 }
294 
295 //___________________________________________________________________________
296 void FidPolyhedron::ConvertMaster2Top(const ROOTGeomAnalyzer* rgeom)
297 {
298  for (unsigned int i = 0; i < fPolyFaces.size(); ++i ) {
299  PlaneParam& aplane = fPolyFaces[i];
300  aplane.ConvertMaster2Top(rgeom);
301  }
302 }
303 void FidPolyhedron::Print(std::ostream& stream) const
304 {
305  size_t nfaces = fPolyFaces.size();
306  stream << "FidPolyhedron n=" << nfaces;
307  for ( size_t i=0; i<nfaces; ++i) {
308  const PlaneParam& aface = fPolyFaces[i];
309  stream << std::endl << "[" << setw(2) << i << "]" << aface;
310  }
311 }
312 
313 //___________________________________________________________________________
Some simple volumes that know how to calculate where a ray intercepts them.
Definition: FidShape.h:90
Int_t fSurfOut
what surface was hit on way in
Definition: FidShape.h:54
void ConvertMaster2Top(const ROOTGeomAnalyzer *rgeom)
Definition: FidShape.cxx:48
Double_t Vn(const TVector3 &raybase) const
Definition: FidShape.h:73
Bool_t fIsHit
distance along ray to exit fid volume
Definition: FidShape.h:52
static constexpr double s
Definition: Units.h:95
string dir
static constexpr double b
Definition: Units.h:78
const double a
Int_t fSurfIn
was the volume hit
Definition: FidShape.h:53
A ROOT/GEANT4 geometry driver.
virtual void Master2TopDir(TVector3 &v) const
void Print(std::ostream &stream) const
Definition: FidShape.cxx:63
Double_t fDistOut
distance along ray to enter fid volume
Definition: FidShape.h:51
Bool_t IsValid() const
Definition: FidShape.h:77
virtual void Print(std::ostream &stream) const =0
std::ostream & operator<<(std::ostream &stream, const genie::geometry::PlaneParam &pparam)
Definition: FidShape.cxx:22
virtual void Master2Top(TVector3 &v) const
Double_t Vd(const TVector3 &raycos) const
Definition: FidShape.h:75