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