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