Geant4_10
G4IntersectingCone.cc
Go to the documentation of this file.
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 //
27 // $Id: G4IntersectingCone.cc 72937 2013-08-14 13:20:38Z gcosmo $
28 //
29 //
30 // --------------------------------------------------------------------
31 // GEANT 4 class source file
32 //
33 //
34 // G4IntersectingCone.cc
35 //
36 // Implementation of a utility class which calculates the intersection
37 // of an arbitrary line with a fixed cone
38 // --------------------------------------------------------------------
39 
40 #include "G4IntersectingCone.hh"
41 #include "G4GeometryTolerance.hh"
42 
43 //
44 // Constructor
45 //
47  const G4double z[2] )
48 {
49  const G4double halfCarTolerance
51 
52  //
53  // What type of cone are we?
54  //
55  type1 = (std::fabs(z[1]-z[0]) > std::fabs(r[1]-r[0]));
56 
57  if (type1)
58  {
59  B = (r[1]-r[0])/(z[1]-z[0]); // tube like
60  A = 0.5*( r[1]+r[0] - B*(z[1]+z[0]) );
61  }
62  else
63  {
64  B = (z[1]-z[0])/(r[1]-r[0]); // disk like
65  A = 0.5*( z[1]+z[0] - B*(r[1]+r[0]) );
66  }
67  //
68  // Calculate extent
69  //
70  if (r[0] < r[1])
71  {
72  rLo = r[0]-halfCarTolerance; rHi = r[1]+halfCarTolerance;
73  }
74  else
75  {
76  rLo = r[1]-halfCarTolerance; rHi = r[0]+halfCarTolerance;
77  }
78 
79  if (z[0] < z[1])
80  {
81  zLo = z[0]-halfCarTolerance; zHi = z[1]+halfCarTolerance;
82  }
83  else
84  {
85  zLo = z[1]-halfCarTolerance; zHi = z[0]+halfCarTolerance;
86  }
87 }
88 
89 
90 //
91 // Fake default constructor - sets only member data and allocates memory
92 // for usage restricted to object persistency.
93 //
95  : zLo(0.), zHi(0.), rLo(0.), rHi(0.), type1(false), A(0.), B(0.)
96 {
97 }
98 
99 
100 //
101 // Destructor
102 //
104 {
105 }
106 
107 
108 //
109 // HitOn
110 //
111 // Check r or z extent, as appropriate, to see if the point is possibly
112 // on the cone.
113 //
115  const G4double z )
116 {
117  //
118  // Be careful! The inequalities cannot be "<=" and ">=" here without
119  // punching a tiny hole in our shape!
120  //
121  if (type1)
122  {
123  if (z < zLo || z > zHi) return false;
124  }
125  else
126  {
127  if (r < rLo || r > rHi) return false;
128  }
129 
130  return true;
131 }
132 
133 
134 //
135 // LineHitsCone
136 //
137 // Calculate the intersection of a line with our conical surface, ignoring
138 // any phi division
139 //
141  const G4ThreeVector &v,
142  G4double *s1, G4double *s2 )
143 {
144  if (type1)
145  {
146  return LineHitsCone1( p, v, s1, s2 );
147  }
148  else
149  {
150  return LineHitsCone2( p, v, s1, s2 );
151  }
152 }
153 
154 
155 //
156 // LineHitsCone1
157 //
158 // Calculate the intersections of a line with a conical surface. Only
159 // suitable if zPlane[0] != zPlane[1].
160 //
161 // Equation of a line:
162 //
163 // x = x0 + s*tx y = y0 + s*ty z = z0 + s*tz
164 //
165 // Equation of a conical surface:
166 //
167 // x**2 + y**2 = (A + B*z)**2
168 //
169 // Solution is quadratic:
170 //
171 // a*s**2 + b*s + c = 0
172 //
173 // where:
174 //
175 // a = x0**2 + y0**2 - (A + B*z0)**2
176 //
177 // b = 2*( x0*tx + y0*ty - (A*B - B*B*z0)*tz)
178 //
179 // c = tx**2 + ty**2 - (B*tz)**2
180 //
181 // Notice, that if a < 0, this indicates that the two solutions (assuming
182 // they exist) are in opposite cones (that is, given z0 = -A/B, one z < z0
183 // and the other z > z0). For our shapes, the invalid solution is one
184 // which produces A + Bz < 0, or the one where Bz is smallest (most negative).
185 // Since Bz = B*s*tz, if B*tz > 0, we want the largest s, otherwise,
186 // the smaller.
187 //
188 // If there are two solutions on one side of the cone, we want to make
189 // sure that they are on the "correct" side, that is A + B*z0 + s*B*tz >= 0.
190 //
191 // If a = 0, we have a linear problem: s = c/b, which again gives one solution.
192 // This should be rare.
193 //
194 // For b*b - 4*a*c = 0, we also have one solution, which is almost always
195 // a line just grazing the surface of a the cone, which we want to ignore.
196 // However, there are two other, very rare, possibilities:
197 // a line intersecting the z axis and either:
198 // 1. At the same angle std::atan(B) to just miss one side of the cone, or
199 // 2. Intersecting the cone apex (0,0,-A/B)
200 // We *don't* want to miss these! How do we identify them? Well, since
201 // this case is rare, we can at least swallow a little more CPU than we would
202 // normally be comfortable with. Intersection with the z axis means
203 // x0*ty - y0*tx = 0. Case (1) means a==0, and we've already dealt with that
204 // above. Case (2) means a < 0.
205 //
206 // Now: x0*tx + y0*ty = 0 in terms of roundoff error. We can write:
207 // Delta = x0*tx + y0*ty
208 // b = 2*( Delta - (A*B + B*B*z0)*tz )
209 // For:
210 // b*b - 4*a*c = epsilon
211 // where epsilon is small, then:
212 // Delta = epsilon/2/B
213 //
215  const G4ThreeVector &v,
216  G4double *s1, G4double *s2 )
217 {
218  G4double x0 = p.x(), y0 = p.y(), z0 = p.z();
219  G4double tx = v.x(), ty = v.y(), tz = v.z();
220 
221  G4double a = tx*tx + ty*ty - sqr(B*tz);
222  G4double b = 2*( x0*tx + y0*ty - (A*B + B*B*z0)*tz);
223  G4double c = x0*x0 + y0*y0 - sqr(A + B*z0);
224 
225  G4double radical = b*b - 4*a*c;
226 
227  if (radical < -1E-6*std::fabs(b)) { return 0; } // No solution
228 
229  if (radical < 1E-6*std::fabs(b))
230  {
231  //
232  // The radical is roughly zero: check for special, very rare, cases
233  //
234  if (std::fabs(a) > 1/kInfinity)
235  {
236  if(B==0.) { return 0; }
237  if ( std::fabs(x0*ty - y0*tx) < std::fabs(1E-6/B) )
238  {
239  *s1 = -0.5*b/a;
240  return 1;
241  }
242  return 0;
243  }
244  }
245  else
246  {
247  radical = std::sqrt(radical);
248  }
249 
250  if (a > 1/kInfinity)
251  {
252  G4double sa, sb, q = -0.5*( b + (b < 0 ? -radical : +radical) );
253  sa = q/a;
254  sb = c/q;
255  if (sa < sb) { *s1 = sa; *s2 = sb; } else { *s1 = sb; *s2 = sa; }
256  if (A + B*(z0+(*s1)*tz) < 0) { return 0; }
257  return 2;
258  }
259  else if (a < -1/kInfinity)
260  {
261  G4double sa, sb, q = -0.5*( b + (b < 0 ? -radical : +radical) );
262  sa = q/a;
263  sb = c/q;
264  *s1 = (B*tz > 0)^(sa > sb) ? sb : sa;
265  return 1;
266  }
267  else if (std::fabs(b) < 1/kInfinity)
268  {
269  return 0;
270  }
271  else
272  {
273  *s1 = -c/b;
274  if (A + B*(z0+(*s1)*tz) < 0) { return 0; }
275  return 1;
276  }
277 }
278 
279 
280 //
281 // LineHitsCone2
282 //
283 // See comments under LineHitsCone1. In this routine, case2, we have:
284 //
285 // Z = A + B*R
286 //
287 // The solution is still quadratic:
288 //
289 // a = tz**2 - B*B*(tx**2 + ty**2)
290 //
291 // b = 2*( (z0-A)*tz - B*B*(x0*tx+y0*ty) )
292 //
293 // c = ( (z0-A)**2 - B*B*(x0**2 + y0**2) )
294 //
295 // The rest is much the same, except some details.
296 //
297 // a > 0 now means we intersect only once in the correct hemisphere.
298 //
299 // a > 0 ? We only want solution which produces R > 0.
300 // since R = (z0+s*tz-A)/B, for tz/B > 0, this is the largest s
301 // for tz/B < 0, this is the smallest s
302 // thus, same as in case 1 ( since sign(tz/B) = sign(tz*B) )
303 //
305  const G4ThreeVector &v,
306  G4double *s1, G4double *s2 )
307 {
308  G4double x0 = p.x(), y0 = p.y(), z0 = p.z();
309  G4double tx = v.x(), ty = v.y(), tz = v.z();
310 
311 
312  // Special case which might not be so rare: B = 0 (precisely)
313  //
314  if (B==0)
315  {
316  if (std::fabs(tz) < 1/kInfinity) { return 0; }
317 
318  *s1 = (A-z0)/tz;
319  return 1;
320  }
321 
322  G4double B2 = B*B;
323 
324  G4double a = tz*tz - B2*(tx*tx + ty*ty);
325  G4double b = 2*( (z0-A)*tz - B2*(x0*tx + y0*ty) );
326  G4double c = sqr(z0-A) - B2*( x0*x0 + y0*y0 );
327 
328  G4double radical = b*b - 4*a*c;
329 
330  if (radical < -1E-6*std::fabs(b)) { return 0; } // No solution
331 
332  if (radical < 1E-6*std::fabs(b))
333  {
334  //
335  // The radical is roughly zero: check for special, very rare, cases
336  //
337  if (std::fabs(a) > 1/kInfinity)
338  {
339  if ( std::fabs(x0*ty - y0*tx) < std::fabs(1E-6/B) )
340  {
341  *s1 = -0.5*b/a;
342  return 1;
343  }
344  return 0;
345  }
346  }
347  else
348  {
349  radical = std::sqrt(radical);
350  }
351 
352  if (a < -1/kInfinity)
353  {
354  G4double sa, sb, q = -0.5*( b + (b < 0 ? -radical : +radical) );
355  sa = q/a;
356  sb = c/q;
357  if (sa < sb) { *s1 = sa; *s2 = sb; } else { *s1 = sb; *s2 = sa; }
358  if ((z0 + (*s1)*tz - A)/B < 0) { return 0; }
359  return 2;
360  }
361  else if (a > 1/kInfinity)
362  {
363  G4double sa, sb, q = -0.5*( b + (b < 0 ? -radical : +radical) );
364  sa = q/a;
365  sb = c/q;
366  *s1 = (tz*B > 0)^(sa > sb) ? sb : sa;
367  return 1;
368  }
369  else if (std::fabs(b) < 1/kInfinity)
370  {
371  return 0;
372  }
373  else
374  {
375  *s1 = -c/b;
376  if ((z0 + (*s1)*tz - A)/B < 0) { return 0; }
377  return 1;
378  }
379 }
tuple a
Definition: test.py:11
double x() const
const char * p
Definition: xmltok.h:285
G4double GetSurfaceTolerance() const
G4int LineHitsCone1(const G4ThreeVector &p, const G4ThreeVector &v, G4double *s1, G4double *s2)
G4bool HitOn(const G4double r, const G4double z)
G4int LineHitsCone(const G4ThreeVector &p, const G4ThreeVector &v, G4double *s1, G4double *s2)
int G4int
Definition: G4Types.hh:78
double z() const
tuple b
Definition: test.py:12
bool G4bool
Definition: G4Types.hh:79
tuple v
Definition: test.py:18
G4int LineHitsCone2(const G4ThreeVector &p, const G4ThreeVector &v, G4double *s1, G4double *s2)
jump r
Definition: plot.C:36
double y() const
tuple z
Definition: test.py:28
G4IntersectingCone(const G4double r[2], const G4double z[2])
T sqr(const T &x)
Definition: templates.hh:145
double G4double
Definition: G4Types.hh:76
tuple c
Definition: test.py:13
static G4GeometryTolerance * GetInstance()