Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
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 //
27 // $Id: G4IntersectingCone.cc 95997 2016-03-07 13:16:25Z 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  static const G4double EPS = DBL_EPSILON; // Precision constant,
219  // originally it was 1E-6
220  G4double x0 = p.x(), y0 = p.y(), z0 = p.z();
221  G4double tx = v.x(), ty = v.y(), tz = v.z();
222 
223  G4double a = tx*tx + ty*ty - sqr(B*tz);
224  G4double b = 2*( x0*tx + y0*ty - (A*B + B*B*z0)*tz);
225  G4double c = x0*x0 + y0*y0 - sqr(A + B*z0);
226 
227  G4double radical = b*b - 4*a*c;
228 
229  if (radical < -EPS*std::fabs(b)) { return 0; } // No solution
230 
231  if (radical < EPS*std::fabs(b))
232  {
233  //
234  // The radical is roughly zero: check for special, very rare, cases
235  //
236  if (std::fabs(a) > 1/kInfinity)
237  {
238  if(B==0.) { return 0; }
239  if ( std::fabs(x0*ty - y0*tx) < std::fabs(EPS/B) )
240  {
241  *s1 = -0.5*b/a;
242  return 1;
243  }
244  return 0;
245  }
246  }
247  else
248  {
249  radical = std::sqrt(radical);
250  }
251 
252  if (a > 1/kInfinity)
253  {
254  G4double sa, sb, q = -0.5*( b + (b < 0 ? -radical : +radical) );
255  sa = q/a;
256  sb = c/q;
257  if (sa < sb) { *s1 = sa; *s2 = sb; } else { *s1 = sb; *s2 = sa; }
258  if (A + B*(z0+(*s1)*tz) < 0) { return 0; }
259  return 2;
260  }
261  else if (a < -1/kInfinity)
262  {
263  G4double sa, sb, q = -0.5*( b + (b < 0 ? -radical : +radical) );
264  sa = q/a;
265  sb = c/q;
266  *s1 = (B*tz > 0)^(sa > sb) ? sb : sa;
267  return 1;
268  }
269  else if (std::fabs(b) < 1/kInfinity)
270  {
271  return 0;
272  }
273  else
274  {
275  *s1 = -c/b;
276  if (A + B*(z0+(*s1)*tz) < 0) { return 0; }
277  return 1;
278  }
279 }
280 
281 
282 //
283 // LineHitsCone2
284 //
285 // See comments under LineHitsCone1. In this routine, case2, we have:
286 //
287 // Z = A + B*R
288 //
289 // The solution is still quadratic:
290 //
291 // a = tz**2 - B*B*(tx**2 + ty**2)
292 //
293 // b = 2*( (z0-A)*tz - B*B*(x0*tx+y0*ty) )
294 //
295 // c = ( (z0-A)**2 - B*B*(x0**2 + y0**2) )
296 //
297 // The rest is much the same, except some details.
298 //
299 // a > 0 now means we intersect only once in the correct hemisphere.
300 //
301 // a > 0 ? We only want solution which produces R > 0.
302 // since R = (z0+s*tz-A)/B, for tz/B > 0, this is the largest s
303 // for tz/B < 0, this is the smallest s
304 // thus, same as in case 1 ( since sign(tz/B) = sign(tz*B) )
305 //
307  const G4ThreeVector &v,
308  G4double *s1, G4double *s2 )
309 {
310  static const G4double EPS = DBL_EPSILON; // Precision constant,
311  // originally it was 1E-6
312  G4double x0 = p.x(), y0 = p.y(), z0 = p.z();
313  G4double tx = v.x(), ty = v.y(), tz = v.z();
314 
315  // Special case which might not be so rare: B = 0 (precisely)
316  //
317  if (B==0)
318  {
319  if (std::fabs(tz) < 1/kInfinity) { return 0; }
320 
321  *s1 = (A-z0)/tz;
322  return 1;
323  }
324 
325  G4double B2 = B*B;
326 
327  G4double a = tz*tz - B2*(tx*tx + ty*ty);
328  G4double b = 2*( (z0-A)*tz - B2*(x0*tx + y0*ty) );
329  G4double c = sqr(z0-A) - B2*( x0*x0 + y0*y0 );
330 
331  G4double radical = b*b - 4*a*c;
332 
333  if (radical < -EPS*std::fabs(b)) { return 0; } // No solution
334 
335  if (radical < EPS*std::fabs(b))
336  {
337  //
338  // The radical is roughly zero: check for special, very rare, cases
339  //
340  if (std::fabs(a) > 1/kInfinity)
341  {
342  if ( std::fabs(x0*ty - y0*tx) < std::fabs(EPS/B) )
343  {
344  *s1 = -0.5*b/a;
345  return 1;
346  }
347  return 0;
348  }
349  }
350  else
351  {
352  radical = std::sqrt(radical);
353  }
354 
355  if (a < -1/kInfinity)
356  {
357  G4double sa, sb, q = -0.5*( b + (b < 0 ? -radical : +radical) );
358  sa = q/a;
359  sb = c/q;
360  if (sa < sb) { *s1 = sa; *s2 = sb; } else { *s1 = sb; *s2 = sa; }
361  if ((z0 + (*s1)*tz - A)/B < 0) { return 0; }
362  return 2;
363  }
364  else if (a > 1/kInfinity)
365  {
366  G4double sa, sb, q = -0.5*( b + (b < 0 ? -radical : +radical) );
367  sa = q/a;
368  sb = c/q;
369  *s1 = (tz*B > 0)^(sa > sb) ? sb : sa;
370  return 1;
371  }
372  else if (std::fabs(b) < 1/kInfinity)
373  {
374  return 0;
375  }
376  else
377  {
378  *s1 = -c/b;
379  if ((z0 + (*s1)*tz - A)/B < 0) { return 0; }
380  return 1;
381  }
382 }
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
static const G4double kInfinity
Definition: geomdefs.hh:42
const char * p
Definition: xmltok.h:285
Double_t z
double z() const
double G4double
Definition: G4Types.hh:76
bool G4bool
Definition: G4Types.hh:79
G4IntersectingCone(const G4double r[2], const G4double z[2])
G4int LineHitsCone(const G4ThreeVector &p, const G4ThreeVector &v, G4double *s1, G4double *s2)
double A(double temperature)
#define DBL_EPSILON
Definition: templates.hh:87
int G4int
Definition: G4Types.hh:78
#define EPS
double x() const
static G4GeometryTolerance * GetInstance()
T sqr(const T &x)
Definition: templates.hh:145
double y() const
G4int LineHitsCone1(const G4ThreeVector &p, const G4ThreeVector &v, G4double *s1, G4double *s2)
double B(double temperature)
G4int LineHitsCone2(const G4ThreeVector &p, const G4ThreeVector &v, G4double *s1, G4double *s2)
G4double GetSurfaceTolerance() const
G4bool HitOn(const G4double r, const G4double z)