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