ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4IntersectingCone.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4IntersectingCone.cc
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // Implementation of G4IntersectingCone, a utility class which calculates
27 // the intersection of an arbitrary line with a fixed cone
28 //
29 // Author: David C. Williams (davidw@scipp.ucsc.edu)
30 // --------------------------------------------------------------------
31 
32 #include "G4IntersectingCone.hh"
33 #include "G4GeometryTolerance.hh"
34 
35 // Constructor
36 //
38  const G4double z[2] )
39 {
40  const G4double halfCarTolerance
42 
43  // What type of cone are we?
44  //
45  type1 = (std::abs(z[1]-z[0]) > std::abs(r[1]-r[0]));
46 
47  if (type1) // tube like
48  {
49  B = (r[1] - r[0]) / (z[1] - z[0]);
50  A = (r[0]*z[1] - r[1]*z[0]) / (z[1] -z[0]);
51  }
52  else // disk like
53  {
54  B = (z[1] - z[0]) / (r[1] - r[0]);
55  A = (z[0]*r[1] - z[1]*r[0]) / (r[1] - r[0]);
56  }
57 
58  // Calculate extent
59  //
60  rLo = std::min(r[0], r[1]) - halfCarTolerance;
61  rHi = std::max(r[0], r[1]) + halfCarTolerance;
62  zLo = std::min(z[0], z[1]) - halfCarTolerance;
63  zHi = std::max(z[0], z[1]) + halfCarTolerance;
64 }
65 
66 // Fake default constructor - sets only member data and allocates memory
67 // for usage restricted to object persistency.
68 //
70  : zLo(0.), zHi(0.), rLo(0.), rHi(0.), A(0.), B(0.)
71 {
72 }
73 
74 // Destructor
75 //
77 {
78 }
79 
80 // HitOn
81 //
82 // Check r or z extent, as appropriate, to see if the point is possibly
83 // on the cone.
84 //
86  const G4double z )
87 {
88  //
89  // Be careful! The inequalities cannot be "<=" and ">=" here without
90  // punching a tiny hole in our shape!
91  //
92  if (type1)
93  {
94  if (z < zLo || z > zHi) return false;
95  }
96  else
97  {
98  if (r < rLo || r > rHi) return false;
99  }
100 
101  return true;
102 }
103 
104 // LineHitsCone
105 //
106 // Calculate the intersection of a line with our conical surface, ignoring
107 // any phi division
108 //
110  const G4ThreeVector& v,
111  G4double* s1, G4double* s2 )
112 {
113  if (type1)
114  {
115  return LineHitsCone1( p, v, s1, s2 );
116  }
117  else
118  {
119  return LineHitsCone2( p, v, s1, s2 );
120  }
121 }
122 
123 // LineHitsCone1
124 //
125 // Calculate the intersections of a line with a conical surface. Only
126 // suitable if zPlane[0] != zPlane[1].
127 //
128 // Equation of a line:
129 //
130 // x = x0 + s*tx y = y0 + s*ty z = z0 + s*tz
131 //
132 // Equation of a conical surface:
133 //
134 // x**2 + y**2 = (A + B*z)**2
135 //
136 // Solution is quadratic:
137 //
138 // a*s**2 + b*s + c = 0
139 //
140 // where:
141 //
142 // a = tx**2 + ty**2 - (B*tz)**2
143 //
144 // b = 2*( px*vx + py*vy - B*(A + B*pz)*vz )
145 //
146 // c = x0**2 + y0**2 - (A + B*z0)**2
147 //
148 // Notice, that if a < 0, this indicates that the two solutions (assuming
149 // they exist) are in opposite cones (that is, given z0 = -A/B, one z < z0
150 // and the other z > z0). For our shapes, the invalid solution is one
151 // which produces A + Bz < 0, or the one where Bz is smallest (most negative).
152 // Since Bz = B*s*tz, if B*tz > 0, we want the largest s, otherwise,
153 // the smaller.
154 //
155 // If there are two solutions on one side of the cone, we want to make
156 // sure that they are on the "correct" side, that is A + B*z0 + s*B*tz >= 0.
157 //
158 // If a = 0, we have a linear problem: s = c/b, which again gives one solution.
159 // This should be rare.
160 //
161 // For b*b - 4*a*c = 0, we also have one solution, which is almost always
162 // a line just grazing the surface of a the cone, which we want to ignore.
163 // However, there are two other, very rare, possibilities:
164 // a line intersecting the z axis and either:
165 // 1. At the same angle std::atan(B) to just miss one side of the cone, or
166 // 2. Intersecting the cone apex (0,0,-A/B)
167 // We *don't* want to miss these! How do we identify them? Well, since
168 // this case is rare, we can at least swallow a little more CPU than we would
169 // normally be comfortable with. Intersection with the z axis means
170 // x0*ty - y0*tx = 0. Case (1) means a==0, and we've already dealt with that
171 // above. Case (2) means a < 0.
172 //
173 // Now: x0*tx + y0*ty = 0 in terms of roundoff error. We can write:
174 // Delta = x0*tx + y0*ty
175 // b = 2*( Delta - B*(A + B*z0)*tz )
176 // For:
177 // b*b - 4*a*c = epsilon
178 // where epsilon is small, then:
179 // Delta = epsilon/2/B
180 //
182  const G4ThreeVector& v,
183  G4double* s1, G4double* s2 )
184 {
185  static const G4double EPS = DBL_EPSILON; // Precision constant,
186  // originally it was 1E-6
187  G4double x0 = p.x(), y0 = p.y(), z0 = p.z();
188  G4double tx = v.x(), ty = v.y(), tz = v.z();
189 
190  // Value of radical can be inaccurate due to loss of precision
191  // if to calculate the coefficiets a,b,c like the following:
192  // G4double a = tx*tx + ty*ty - sqr(B*tz);
193  // G4double b = 2*( x0*tx + y0*ty - B*(A + B*z0)*tz);
194  // G4double c = x0*x0 + y0*y0 - sqr(A + B*z0);
195  //
196  // For more accurate calculation of radical the coefficients
197  // are splitted in two components, radial and along z-axis
198  //
199  G4double ar = tx*tx + ty*ty;
200  G4double az = sqr(B*tz);
201  G4double br = 2*(x0*tx + y0*ty);
202  G4double bz = 2*B*(A + B*z0)*tz;
203  G4double cr = x0*x0 + y0*y0;
204  G4double cz = sqr(A + B*z0);
205 
206  // Instead radical = b*b - 4*a*c
207  G4double arcz = 4*ar*cz;
208  G4double azcr = 4*az*cr;
209  G4double radical = (br*br - 4*ar*cr) + ((std::max(arcz,azcr) - 2*bz*br) + std::min(arcz,azcr));
210 
211  // Find the coefficients
212  G4double a = ar - az;
213  G4double b = br - bz;
214  G4double c = cr - cz;
215 
216  if (radical < -EPS*std::fabs(b)) { return 0; } // No solution
217 
218  if (radical < EPS*std::fabs(b))
219  {
220  //
221  // The radical is roughly zero: check for special, very rare, cases
222  //
223  if (std::fabs(a) > 1/kInfinity)
224  {
225  if(B==0.) { return 0; }
226  if ( std::fabs(x0*ty - y0*tx) < std::fabs(EPS/B) )
227  {
228  *s1 = -0.5*b/a;
229  return 1;
230  }
231  return 0;
232  }
233  }
234  else
235  {
236  radical = std::sqrt(radical);
237  }
238 
239  if (a > 1/kInfinity)
240  {
241  G4double sa, sb, q = -0.5*( b + (b < 0 ? -radical : +radical) );
242  sa = q/a;
243  sb = c/q;
244  if (sa < sb) { *s1 = sa; *s2 = sb; } else { *s1 = sb; *s2 = sa; }
245  if (A + B*(z0+(*s1)*tz) < 0) { return 0; }
246  return 2;
247  }
248  else if (a < -1/kInfinity)
249  {
250  G4double sa, sb, q = -0.5*( b + (b < 0 ? -radical : +radical) );
251  sa = q/a;
252  sb = c/q;
253  *s1 = (B*tz > 0)^(sa > sb) ? sb : sa;
254  return 1;
255  }
256  else if (std::fabs(b) < 1/kInfinity)
257  {
258  return 0;
259  }
260  else
261  {
262  *s1 = -c/b;
263  if (A + B*(z0+(*s1)*tz) < 0) { return 0; }
264  return 1;
265  }
266 }
267 
268 // LineHitsCone2
269 //
270 // See comments under LineHitsCone1. In this routine, case2, we have:
271 //
272 // Z = A + B*R
273 //
274 // The solution is still quadratic:
275 //
276 // a = tz**2 - B*B*(tx**2 + ty**2)
277 //
278 // b = 2*( (z0-A)*tz - B*B*(x0*tx+y0*ty) )
279 //
280 // c = ( (z0-A)**2 - B*B*(x0**2 + y0**2) )
281 //
282 // The rest is much the same, except some details.
283 //
284 // a > 0 now means we intersect only once in the correct hemisphere.
285 //
286 // a > 0 ? We only want solution which produces R > 0.
287 // since R = (z0+s*tz-A)/B, for tz/B > 0, this is the largest s
288 // for tz/B < 0, this is the smallest s
289 // thus, same as in case 1 ( since sign(tz/B) = sign(tz*B) )
290 //
292  const G4ThreeVector& v,
293  G4double* s1, G4double* s2 )
294 {
295  static const G4double EPS = DBL_EPSILON; // Precision constant,
296  // originally it was 1E-6
297  G4double x0 = p.x(), y0 = p.y(), z0 = p.z();
298  G4double tx = v.x(), ty = v.y(), tz = v.z();
299 
300  // Special case which might not be so rare: B = 0 (precisely)
301  //
302  if (B==0)
303  {
304  if (std::fabs(tz) < 1/kInfinity) { return 0; }
305 
306  *s1 = (A-z0)/tz;
307  return 1;
308  }
309 
310  // Value of radical can be inaccurate due to loss of precision
311  // if to calculate the coefficiets a,b,c like the following:
312  // G4double a = tz*tz - B2*(tx*tx + ty*ty);
313  // G4double b = 2*( (z0-A)*tz - B2*(x0*tx + y0*ty) );
314  // G4double c = sqr(z0-A) - B2*( x0*x0 + y0*y0 );
315  //
316  // For more accurate calculation of radical the coefficients
317  // are splitted in two components, radial and along z-axis
318  //
319  G4double B2 = B*B;
320 
321  G4double az = tz*tz;
322  G4double ar = B2*(tx*tx + ty*ty);
323  G4double bz = 2*(z0-A)*tz;
324  G4double br = 2*B2*(x0*tx + y0*ty);
325  G4double cz = sqr(z0-A);
326  G4double cr = B2*(x0*x0 + y0*y0);
327 
328  // Instead radical = b*b - 4*a*c
329  G4double arcz = 4*ar*cz;
330  G4double azcr = 4*az*cr;
331  G4double radical = (br*br - 4*ar*cr) + ((std::max(arcz,azcr) - 2*bz*br) + std::min(arcz,azcr));
332 
333  // Find the coefficients
334  G4double a = az - ar;
335  G4double b = bz - br;
336  G4double c = cz - cr;
337 
338  if (radical < -EPS*std::fabs(b)) { return 0; } // No solution
339 
340  if (radical < EPS*std::fabs(b))
341  {
342  //
343  // The radical is roughly zero: check for special, very rare, cases
344  //
345  if (std::fabs(a) > 1/kInfinity)
346  {
347  if ( std::fabs(x0*ty - y0*tx) < std::fabs(EPS/B) )
348  {
349  *s1 = -0.5*b/a;
350  return 1;
351  }
352  return 0;
353  }
354  }
355  else
356  {
357  radical = std::sqrt(radical);
358  }
359 
360  if (a < -1/kInfinity)
361  {
362  G4double sa, sb, q = -0.5*( b + (b < 0 ? -radical : +radical) );
363  sa = q/a;
364  sb = c/q;
365  if (sa < sb) { *s1 = sa; *s2 = sb; } else { *s1 = sb; *s2 = sa; }
366  if ((z0 + (*s1)*tz - A)/B < 0) { return 0; }
367  return 2;
368  }
369  else if (a > 1/kInfinity)
370  {
371  G4double sa, sb, q = -0.5*( b + (b < 0 ? -radical : +radical) );
372  sa = q/a;
373  sb = c/q;
374  *s1 = (tz*B > 0)^(sa > sb) ? sb : sa;
375  return 1;
376  }
377  else if (std::fabs(b) < 1/kInfinity)
378  {
379  return 0;
380  }
381  else
382  {
383  *s1 = -c/b;
384  if ((z0 + (*s1)*tz - A)/B < 0) { return 0; }
385  return 1;
386  }
387 }