00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040 #include "G4IntersectingCone.hh"
00041 #include "G4GeometryTolerance.hh"
00042
00043
00044
00045
00046 G4IntersectingCone::G4IntersectingCone( const G4double r[2],
00047 const G4double z[2] )
00048 {
00049 static const G4double halfCarTolerance
00050 = 0.5 * G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
00051
00052
00053
00054
00055 type1 = (std::fabs(z[1]-z[0]) > std::fabs(r[1]-r[0]));
00056
00057 if (type1)
00058 {
00059 B = (r[1]-r[0])/(z[1]-z[0]);
00060 A = 0.5*( r[1]+r[0] - B*(z[1]+z[0]) );
00061 }
00062 else
00063 {
00064 B = (z[1]-z[0])/(r[1]-r[0]);
00065 A = 0.5*( z[1]+z[0] - B*(r[1]+r[0]) );
00066 }
00067
00068
00069
00070 if (r[0] < r[1])
00071 {
00072 rLo = r[0]-halfCarTolerance; rHi = r[1]+halfCarTolerance;
00073 }
00074 else
00075 {
00076 rLo = r[1]-halfCarTolerance; rHi = r[0]+halfCarTolerance;
00077 }
00078
00079 if (z[0] < z[1])
00080 {
00081 zLo = z[0]-halfCarTolerance; zHi = z[1]+halfCarTolerance;
00082 }
00083 else
00084 {
00085 zLo = z[1]-halfCarTolerance; zHi = z[0]+halfCarTolerance;
00086 }
00087 }
00088
00089
00090
00091
00092
00093
00094 G4IntersectingCone::G4IntersectingCone( __void__& )
00095 : zLo(0.), zHi(0.), rLo(0.), rHi(0.), type1(false), A(0.), B(0.)
00096 {
00097 }
00098
00099
00100
00101
00102
00103 G4IntersectingCone::~G4IntersectingCone()
00104 {
00105 }
00106
00107
00108
00109
00110
00111
00112
00113
00114 G4bool G4IntersectingCone::HitOn( const G4double r,
00115 const G4double z )
00116 {
00117
00118
00119
00120
00121 if (type1)
00122 {
00123 if (z < zLo || z > zHi) return false;
00124 }
00125 else
00126 {
00127 if (r < rLo || r > rHi) return false;
00128 }
00129
00130 return true;
00131 }
00132
00133
00134
00135
00136
00137
00138
00139
00140 G4int G4IntersectingCone::LineHitsCone( const G4ThreeVector &p,
00141 const G4ThreeVector &v,
00142 G4double *s1, G4double *s2 )
00143 {
00144 if (type1)
00145 {
00146 return LineHitsCone1( p, v, s1, s2 );
00147 }
00148 else
00149 {
00150 return LineHitsCone2( p, v, s1, s2 );
00151 }
00152 }
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214 G4int G4IntersectingCone::LineHitsCone1( const G4ThreeVector &p,
00215 const G4ThreeVector &v,
00216 G4double *s1, G4double *s2 )
00217 {
00218 G4double x0 = p.x(), y0 = p.y(), z0 = p.z();
00219 G4double tx = v.x(), ty = v.y(), tz = v.z();
00220
00221 G4double a = tx*tx + ty*ty - sqr(B*tz);
00222 G4double b = 2*( x0*tx + y0*ty - (A*B + B*B*z0)*tz);
00223 G4double c = x0*x0 + y0*y0 - sqr(A + B*z0);
00224
00225 G4double radical = b*b - 4*a*c;
00226
00227 if (radical < -1E-6*std::fabs(b)) { return 0; }
00228
00229 if (radical < 1E-6*std::fabs(b))
00230 {
00231
00232
00233
00234 if (std::fabs(a) > 1/kInfinity)
00235 {
00236 if(B==0.) { return 0; }
00237 if ( std::fabs(x0*ty - y0*tx) < std::fabs(1E-6/B) )
00238 {
00239 *s1 = -0.5*b/a;
00240 return 1;
00241 }
00242 return 0;
00243 }
00244 }
00245 else
00246 {
00247 radical = std::sqrt(radical);
00248 }
00249
00250 if (a > 1/kInfinity)
00251 {
00252 G4double sa, sb, q = -0.5*( b + (b < 0 ? -radical : +radical) );
00253 sa = q/a;
00254 sb = c/q;
00255 if (sa < sb) { *s1 = sa; *s2 = sb; } else { *s1 = sb; *s2 = sa; }
00256 if (A + B*(z0+(*s1)*tz) < 0) { return 0; }
00257 return 2;
00258 }
00259 else if (a < -1/kInfinity)
00260 {
00261 G4double sa, sb, q = -0.5*( b + (b < 0 ? -radical : +radical) );
00262 sa = q/a;
00263 sb = c/q;
00264 *s1 = (B*tz > 0)^(sa > sb) ? sb : sa;
00265 return 1;
00266 }
00267 else if (std::fabs(b) < 1/kInfinity)
00268 {
00269 return 0;
00270 }
00271 else
00272 {
00273 *s1 = -c/b;
00274 if (A + B*(z0+(*s1)*tz) < 0) { return 0; }
00275 return 1;
00276 }
00277 }
00278
00279
00280
00281
00282
00283
00284
00285
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303
00304 G4int G4IntersectingCone::LineHitsCone2( const G4ThreeVector &p,
00305 const G4ThreeVector &v,
00306 G4double *s1, G4double *s2 )
00307 {
00308 G4double x0 = p.x(), y0 = p.y(), z0 = p.z();
00309 G4double tx = v.x(), ty = v.y(), tz = v.z();
00310
00311
00312
00313
00314 if (B==0)
00315 {
00316 if (std::fabs(tz) < 1/kInfinity) { return 0; }
00317
00318 *s1 = (A-z0)/tz;
00319 return 1;
00320 }
00321
00322 G4double B2 = B*B;
00323
00324 G4double a = tz*tz - B2*(tx*tx + ty*ty);
00325 G4double b = 2*( (z0-A)*tz - B2*(x0*tx + y0*ty) );
00326 G4double c = sqr(z0-A) - B2*( x0*x0 + y0*y0 );
00327
00328 G4double radical = b*b - 4*a*c;
00329
00330 if (radical < -1E-6*std::fabs(b)) { return 0; }
00331
00332 if (radical < 1E-6*std::fabs(b))
00333 {
00334
00335
00336
00337 if (std::fabs(a) > 1/kInfinity)
00338 {
00339 if ( std::fabs(x0*ty - y0*tx) < std::fabs(1E-6/B) )
00340 {
00341 *s1 = -0.5*b/a;
00342 return 1;
00343 }
00344 return 0;
00345 }
00346 }
00347 else
00348 {
00349 radical = std::sqrt(radical);
00350 }
00351
00352 if (a < -1/kInfinity)
00353 {
00354 G4double sa, sb, q = -0.5*( b + (b < 0 ? -radical : +radical) );
00355 sa = q/a;
00356 sb = c/q;
00357 if (sa < sb) { *s1 = sa; *s2 = sb; } else { *s1 = sb; *s2 = sa; }
00358 if ((z0 + (*s1)*tz - A)/B < 0) { return 0; }
00359 return 2;
00360 }
00361 else if (a > 1/kInfinity)
00362 {
00363 G4double sa, sb, q = -0.5*( b + (b < 0 ? -radical : +radical) );
00364 sa = q/a;
00365 sb = c/q;
00366 *s1 = (tz*B > 0)^(sa > sb) ? sb : sa;
00367 return 1;
00368 }
00369 else if (std::fabs(b) < 1/kInfinity)
00370 {
00371 return 0;
00372 }
00373 else
00374 {
00375 *s1 = -c/b;
00376 if ((z0 + (*s1)*tz - A)/B < 0) { return 0; }
00377 return 1;
00378 }
00379 }