Geant4  10.00.p03
UIntersectingCone.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * This Software is part of the AIDA Unified Solids Library package *
4 // * See: https://aidasoft.web.cern.ch/USolids *
5 // ********************************************************************
6 //
7 // $Id:$
8 //
9 // --------------------------------------------------------------------
10 //
11 // UIntersectingCone
12 //
13 // 19.02.13 Marek Gayer
14 // Created from original implementation in Geant4
15 // --------------------------------------------------------------------
16 
17 #include "UUtils.hh"
18 #include <string>
19 #include <cmath>
20 #include <sstream>
21 #include "UIntersectingCone.hh"
22 #include "VUSolid.hh"
23 
24 const double UIntersectingCone::EpsilonQuad = 1.0 / 9.0E99;
25 
26 //
27 // Constructor
28 //
30  const double z[2])
31 {
32  static const double halfCarTolerance
33  = 0.5 * VUSolid::Tolerance(); // UGeometryTolerance::GetInstance()->GetSurfaceTolerance();
34 
35  //
36  // What type of cone are we?
37  //
38  type1 = (std::fabs(z[1] - z[0]) > std::fabs(r[1] - r[0]));
39 
40  if (type1)
41  {
42  B = (r[1] - r[0]) / (z[1] - z[0]); // tube like
43  A = 0.5 * (r[1] + r[0] - B * (z[1] + z[0]));
44  }
45  else
46  {
47  B = (z[1] - z[0]) / (r[1] - r[0]); // disk like
48  A = 0.5 * (z[1] + z[0] - B * (r[1] + r[0]));
49  }
50  //
51  // Calculate extent
52  //
53  if (r[0] < r[1])
54  {
55  rLo = r[0] - halfCarTolerance;
56  rHi = r[1] + halfCarTolerance;
57  }
58  else
59  {
60  rLo = r[1] - halfCarTolerance;
61  rHi = r[0] + halfCarTolerance;
62  }
63 
64  if (z[0] < z[1])
65  {
66  zLo = z[0] - halfCarTolerance;
67  zHi = z[1] + halfCarTolerance;
68  }
69  else
70  {
71  zLo = z[1] - halfCarTolerance;
72  zHi = z[0] + halfCarTolerance;
73  }
74 }
75 
76 
77 /*
78 //
79 // Fake default constructor - sets only member data and allocates memory
80 // for usage restricted to object persistency.
81 //
82 UIntersectingCone::UIntersectingCone( __void__& )
83 : zLo(0.), zHi(0.), rLo(0.), rHi(0.), type1(false), A(0.), B(0.)
84 {
85 }
86 */
87 
88 
89 //
90 // Destructor
91 //
93 {
94 }
95 
96 
97 //
98 // HitOn
99 //
100 // Check r or z extent, as appropriate, to see if the point is possibly
101 // on the cone.
102 //
103 bool UIntersectingCone::HitOn(const double r,
104  const double z)
105 {
106  //
107  // Be careful! The inequalities cannot be "<=" and ">=" here without
108  // punching a tiny hole in our shape!
109  //
110  if (type1)
111  {
112  if (z < zLo || z > zHi) return false;
113  }
114  else
115  {
116  if (r < rLo || r > rHi) return false;
117  }
118 
119  return true;
120 }
121 
122 
123 //
124 // LineHitsCone
125 //
126 // Calculate the intersection of a line with our conical surface, ignoring
127 // any phi division
128 //
129 int UIntersectingCone::LineHitsCone(const UVector3& p, const UVector3& v, double& s1, double& s2)
130 {
131  if (type1)
132  {
133  return LineHitsCone1(p, v, s1, s2);
134  }
135  else
136  {
137  return LineHitsCone2(p, v, s1, s2);
138  }
139 }
140 
141 
142 //
143 // LineHitsCone1
144 //
145 // Calculate the intersections of a line with a conical surface. Only
146 // suitable if zPlane[0] != zPlane[1].
147 //
148 // Equation of a line:
149 //
150 // x = x0 + s*tx y = y0 + s*ty z = z0 + s*tz
151 //
152 // Equation of a conical surface:
153 //
154 // x**2 + y**2 = (A + B*z)**2
155 //
156 // Solution is quadratic:
157 //
158 // a*s**2 + b*s + c = 0
159 //
160 // where:
161 //
162 // a = x0**2 + y0**2 - (A + B*z0)**2
163 //
164 // b = 2*( x0*tx + y0*ty - (A*B - B*B*z0)*tz)
165 //
166 // c = tx**2 + ty**2 - (B*tz)**2
167 //
168 // Notice, that if a < 0, this indicates that the two solutions (assuming
169 // they exist) are in opposite cones (that is, given z0 = -A/B, one z < z0
170 // and the other z > z0). For our shapes, the invalid solution is one
171 // which produces A + Bz < 0, or the one where Bz is smallest (most negative).
172 // Since Bz = B*s*tz, if B*tz > 0, we want the largest s, otherwise,
173 // the smaller.
174 //
175 // If there are two solutions on one side of the cone, we want to make
176 // sure that they are on the "correct" side, that is A + B*z0 + s*B*tz >= 0.
177 //
178 // If a = 0, we have a linear problem: s = c/b, which again gives one solution.
179 // This should be rare.
180 //
181 // For b*b - 4*a*c = 0, we also have one solution, which is almost always
182 // a line just grazing the surface of a the cone, which we want to ignore.
183 // However, there are two other, very rare, possibilities:
184 // a line intersecting the z axis and either:
185 // 1. At the same angle std::atan(B) to just miss one side of the cone, or
186 // 2. Intersecting the cone apex (0,0,-A/B)
187 // We *don't* want to miss these! How do we identify them? Well, since
188 // this case is rare, we can at least swallow a little more CPU than we would
189 // normally be comfortable with. Intersection with the z axis means
190 // x0*ty - y0*tx = 0. Case (1) means a==0, and we've already dealt with that
191 // above. Case (2) means a < 0.
192 //
193 // Now: x0*tx + y0*ty = 0 in terms of roundoff error. We can write:
194 // Delta = x0*tx + y0*ty
195 // b = 2*( Delta - (A*B + B*B*z0)*tz )
196 // For:
197 // b*b - 4*a*c = epsilon
198 // where epsilon is small, then:
199 // Delta = epsilon/2/B
200 //
201 
202 /*
203 int UIntersectingCone::Solution (const UVector3 &p, const UVector3 &v, double a, double b, double c, double &s1, double &s2)
204 {
205  return 0 || 1 || 2;
206 }
207 */
208 
209 int UIntersectingCone::LineHitsCone1(const UVector3& p, const UVector3& v, double& s1, double& s2)
210 {
211  double x0 = p.x, y0 = p.y, z0 = p.z;
212  double tx = v.x, ty = v.y, tz = v.z;
213 
214  double a = tx * tx + ty * ty - UUtils::sqr(B * tz);
215  double b = 2 * (x0 * tx + y0 * ty - (A * B + B * B * z0) * tz);
216  double c = x0 * x0 + y0 * y0 - UUtils::sqr(A + B * z0);
217 
218  double radical = b * b - 4 * a * c;
219  double radicalSqrt;
220 
221  double minRadical = 1E-6 * std::fabs(b);
222 
223  if (radical < -minRadical)
224  {
225  return 0; // No solution
226  }
227 
228  if (radical < minRadical)
229  {
230  //
231  // The radical is roughly zero: check for special, very rare, cases
232  //
233  if (std::fabs(a) > EpsilonQuad)
234  {
235  if (B == 0.)
236  {
237  return 0;
238  }
239  if (std::fabs(x0 * ty - y0 * tx) < std::fabs(1E-6 / B))
240  {
241  s1 = -0.5 * b / a;
242  return 1;
243  }
244  return 0;
245  }
246  radicalSqrt = radical; //TODO: check this case
247  }
248  else
249  {
250  radicalSqrt = std::sqrt(radical);
251  }
252 
253  if (a > EpsilonQuad)
254  {
255  double sa, sb, q = -0.5 * (b + (b < 0 ? -radicalSqrt : +radicalSqrt));
256  sa = q / a;
257  sb = c / q;
258  if (sa < sb)
259  {
260  s1 = sa;
261  s2 = sb;
262  }
263  else
264  {
265  s1 = sb;
266  s2 = sa;
267  }
268  if (A + B * (z0 + (s1)*tz) < 0)
269  {
270  return 0;
271  }
272  // if ((z0 + (s1)*tz - A)/B < 0) { return 0; } // these lines are equivalent
273  return 2;
274  }
275  else if (a < -EpsilonQuad)
276  {
277  double sa, sb, q = -0.5 * (b + (b < 0 ? -radicalSqrt : +radicalSqrt));
278  sa = q / a;
279  sb = c / q;
280  s1 = (tz * B > 0) ^ (sa > sb) ? sb : sa;
281  return 1;
282  }
283  else if (std::fabs(b) < EpsilonQuad)
284  {
285  return 0;
286  }
287  else
288  {
289  s1 = -c / b;
290  if (A + B * (z0 + (s1)*tz) < 0)
291  {
292  return 0;
293  }
294  return 1;
295  }
296 }
297 
298 
299 
300 
301 int UIntersectingCone::LineHitsCone1Optimized(const UVector3& p, const UVector3& v, double& s1, double& s2)
302 {
303  double x0 = p.x, y0 = p.y, z0 = p.z;
304  double tx = v.x, ty = v.y, tz = v.z;
305 
306  double a = tx * tx + ty * ty - UUtils::sqr(B * tz);
307  double b = 2 * (x0 * tx + y0 * ty - (A * B + B * B * z0) * tz);
308  double c = x0 * x0 + y0 * y0 - UUtils::sqr(A + B * z0);
309 
310  double radical = b * b - 4 * a * c;
311 
312  double minRadical = 1E-6 * std::fabs(b);
313 
314  if (radical < -minRadical)
315  {
316  return 0; // No solution
317  }
318 
319  if (std::fabs(a) > EpsilonQuad)
320  {
321  if (radical < minRadical)
322  {
323  //
324  // The radical is roughly zero: check for special, very rare, cases
325  //
326 
327  if (B == 0.)
328  {
329  return 0;
330  }
331  if (std::fabs(x0 * ty - y0 * tx) < std::fabs(1E-6 / B))
332  {
333  s1 = -0.5 * b / a;
334  return 1;
335  }
336  return 0;
337  }
338  else
339  {
340  double radicalSqrt = std::sqrt(radical);
341 
342  if (a > 0)
343  {
344  double sa, sb, q = -0.5 * (b + (b < 0 ? -radicalSqrt : +radicalSqrt));
345  sa = q / a;
346  sb = c / q;
347  if (sa < sb)
348  {
349  s1 = sa;
350  s2 = sb;
351  }
352  else
353  {
354  s1 = sb;
355  s2 = sa;
356  }
357  if (A + B * (z0 + (s1)*tz) < 0)
358  {
359  return 0;
360  }
361  // if ((z0 + (s1)*tz - A)/B < 0) { return 0; } // these lines are equivalent
362  return 2;
363  }
364  else
365  {
366  double sa, sb, q = -0.5 * (b + (b < 0 ? -radicalSqrt : +radicalSqrt));
367  sa = q / a;
368  sb = c / q;
369  s1 = (tz * B > 0) ^ (sa > sb) ? sb : sa;
370  return 1;
371  }
372 
373  }
374  }
375 
376  if (std::fabs(b) < EpsilonQuad)
377  {
378  return 0;
379  }
380  else
381  {
382  s1 = -c / b;
383  if (A + B * (z0 + (s1)*tz) < 0)
384  {
385  return 0;
386  }
387  return 1;
388  }
389 }
390 
391 
392 //
393 // LineHitsCone2
394 //
395 // See comments under LineHitsCone1. In this routine, case2, we have:
396 //
397 // Z = A + B*R
398 //
399 // The solution is still quadratic:
400 //
401 // a = tz**2 - B*B*(tx**2 + ty**2)
402 //
403 // b = 2*( (z0-A)*tz - B*B*(x0*tx+y0*ty) )
404 //
405 // c = ( (z0-A)**2 - B*B*(x0**2 + y0**2) )
406 //
407 // The rest is much the same, except some details.
408 //
409 // a > 0 now means we intersect only once in the correct hemisphere.
410 //
411 // a > 0 ? We only want solution which produces R > 0.
412 // since R = (z0+s*tz-A)/B, for tz/B > 0, this is the largest s
413 // for tz/B < 0, this is the smallest s
414 // thus, same as in case 1 ( since sign(tz/B) = sign(tz*B) )
415 //
417  const UVector3& v,
418  double& s1, double& s2)
419 {
420  double x0 = p.x, y0 = p.y, z0 = p.z;
421  double tx = v.x, ty = v.y, tz = v.z;
422 
423  // Special case which might not be so rare: B = 0 (precisely)
424  //
425  if (B == 0)
426  {
427  if (std::fabs(tz) < EpsilonQuad)
428  {
429  return 0;
430  }
431 
432  s1 = (A - z0) / tz;
433  return 1;
434  }
435 
436  double B2 = B * B;
437 
438  double a = tz * tz - B2 * (tx * tx + ty * ty);
439  double b = 2 * ((z0 - A) * tz - B2 * (x0 * tx + y0 * ty));
440  double c = UUtils::sqr(z0 - A) - B2 * (x0 * x0 + y0 * y0);
441 
442  double radical = b * b - 4 * a * c;
443 
444  if (radical < -1E-6 * std::fabs(b))
445  {
446  return 0; // No solution
447  }
448 
449  if (radical < 1E-6 * std::fabs(b))
450  {
451  //
452  // The radical is roughly zero: check for special, very rare, cases
453  //
454  if (std::fabs(a) > EpsilonQuad)
455  {
456  if (std::fabs(x0 * ty - y0 * tx) < std::fabs(1E-6 / B))
457  {
458  s1 = -0.5 * b / a;
459  return 1;
460  }
461  return 0;
462  }
463  }
464  else
465  {
466  radical = std::sqrt(radical);
467  }
468 
469  if (a < -EpsilonQuad)
470  {
471  double sa, sb, q = -0.5 * (b + (b < 0 ? -radical : +radical));
472  sa = q / a;
473  sb = c / q;
474  if (sa < sb)
475  {
476  s1 = sa;
477  s2 = sb;
478  }
479  else
480  {
481  s1 = sb;
482  s2 = sa;
483  }
484  if ((z0 + (s1)*tz - A) / B < 0)
485  {
486  return 0;
487  }
488  return 2;
489  }
490  else if (a > EpsilonQuad)
491  {
492  double sa, sb, q = -0.5 * (b + (b < 0 ? -radical : +radical));
493  sa = q / a;
494  sb = c / q;
495  s1 = (tz * B > 0) ^ (sa > sb) ? sb : sa;
496  return 1;
497  }
498  else if (std::fabs(b) < EpsilonQuad)
499  {
500  return 0;
501  }
502  else
503  {
504  s1 = -c / b;
505  if ((z0 + (s1)*tz - A) / B < 0)
506  {
507  return 0;
508  }
509  return 1;
510  }
511 }
512 
int LineHitsCone1(const UVector3 &p, const UVector3 &v, double &s1, double &s2)
int LineHitsCone1Optimized(const UVector3 &p, const UVector3 &v, double &s1, double &s2)
int LineHitsCone2(const UVector3 &p, const UVector3 &v, double &s1, double &s2)
G4double z
Definition: TRTMaterials.hh:39
virtual ~UIntersectingCone()
G4double a
Definition: TRTMaterials.hh:39
static double Tolerance()
Definition: VUSolid.hh:127
double x
Definition: UVector3.hh:136
bool HitOn(const double r, const double z)
static const double EpsilonQuad
T sqr(const T &x)
Definition: UUtils.hh:137
double z
Definition: UVector3.hh:138
UIntersectingCone(const double r[2], const double z[2])
double y
Definition: UVector3.hh:137
int LineHitsCone(const UVector3 &p, const UVector3 &v, double &s1, double &s2)