Geant4-11
Public Member Functions | Protected Member Functions | Protected Attributes
G4IntersectingCone Class Reference

#include <G4IntersectingCone.hh>

Public Member Functions

 G4IntersectingCone (__void__ &)
 
 G4IntersectingCone (const G4double r[2], const G4double z[2])
 
G4bool HitOn (const G4double r, const G4double z)
 
G4int LineHitsCone (const G4ThreeVector &p, const G4ThreeVector &v, G4double *s1, G4double *s2)
 
G4double RHi () const
 
G4double RLo () const
 
G4double ZHi () const
 
G4double ZLo () const
 
virtual ~G4IntersectingCone ()
 

Protected Member Functions

G4int LineHitsCone1 (const G4ThreeVector &p, const G4ThreeVector &v, G4double *s1, G4double *s2)
 
G4int LineHitsCone2 (const G4ThreeVector &p, const G4ThreeVector &v, G4double *s1, G4double *s2)
 

Protected Attributes

G4double A
 
G4double B
 
G4double rHi
 
G4double rLo
 
G4bool type1 = false
 
G4double zHi
 
G4double zLo
 

Detailed Description

Definition at line 42 of file G4IntersectingCone.hh.

Constructor & Destructor Documentation

◆ G4IntersectingCone() [1/2]

G4IntersectingCone::G4IntersectingCone ( const G4double  r[2],
const G4double  z[2] 
)

Definition at line 37 of file G4IntersectingCone.cc.

39{
40 const G4double halfCarTolerance
42
43 // What type of cone are we?
44 //
45 type1 = (std::abs(z[1]-z[0]) > std::abs(r[1]-r[0]));
46
47 if (type1) // tube like
48 {
49 B = (r[1] - r[0]) / (z[1] - z[0]);
50 A = (r[0]*z[1] - r[1]*z[0]) / (z[1] -z[0]);
51 }
52 else // disk like
53 {
54 B = (z[1] - z[0]) / (r[1] - r[0]);
55 A = (z[0]*r[1] - z[1]*r[0]) / (r[1] - r[0]);
56 }
57
58 // Calculate extent
59 //
60 rLo = std::min(r[0], r[1]) - halfCarTolerance;
61 rHi = std::max(r[0], r[1]) + halfCarTolerance;
62 zLo = std::min(z[0], z[1]) - halfCarTolerance;
63 zHi = std::max(z[0], z[1]) + halfCarTolerance;
64}
double G4double
Definition: G4Types.hh:83
G4double GetSurfaceTolerance() const
static G4GeometryTolerance * GetInstance()
T max(const T t1, const T t2)
brief Return the largest of the two arguments
T min(const T t1, const T t2)
brief Return the smallest of the two arguments

References A, B, G4GeometryTolerance::GetInstance(), G4GeometryTolerance::GetSurfaceTolerance(), G4INCL::Math::max(), G4INCL::Math::min(), rHi, rLo, type1, zHi, and zLo.

◆ ~G4IntersectingCone()

G4IntersectingCone::~G4IntersectingCone ( )
virtual

Definition at line 76 of file G4IntersectingCone.cc.

77{
78}

◆ G4IntersectingCone() [2/2]

G4IntersectingCone::G4IntersectingCone ( __void__ &  )

Definition at line 69 of file G4IntersectingCone.cc.

70 : zLo(0.), zHi(0.), rLo(0.), rHi(0.), A(0.), B(0.)
71{
72}

Member Function Documentation

◆ HitOn()

G4bool G4IntersectingCone::HitOn ( const G4double  r,
const G4double  z 
)

Definition at line 85 of file G4IntersectingCone.cc.

87{
88 //
89 // Be careful! The inequalities cannot be "<=" and ">=" here without
90 // punching a tiny hole in our shape!
91 //
92 if (type1)
93 {
94 if (z < zLo || z > zHi) return false;
95 }
96 else
97 {
98 if (r < rLo || r > rHi) return false;
99 }
100
101 return true;
102}

References rHi, type1, and zHi.

Referenced by G4PolyconeSide::PointOnCone().

◆ LineHitsCone()

G4int G4IntersectingCone::LineHitsCone ( const G4ThreeVector p,
const G4ThreeVector v,
G4double s1,
G4double s2 
)

Definition at line 109 of file G4IntersectingCone.cc.

112{
113 if (type1)
114 {
115 return LineHitsCone1( p, v, s1, s2 );
116 }
117 else
118 {
119 return LineHitsCone2( p, v, s1, s2 );
120 }
121}
G4int LineHitsCone1(const G4ThreeVector &p, const G4ThreeVector &v, G4double *s1, G4double *s2)
G4int LineHitsCone2(const G4ThreeVector &p, const G4ThreeVector &v, G4double *s1, G4double *s2)

References LineHitsCone1(), LineHitsCone2(), and type1.

Referenced by G4PolyconeSide::Intersect(), and G4PolyhedraSide::LineHitsSegments().

◆ LineHitsCone1()

G4int G4IntersectingCone::LineHitsCone1 ( const G4ThreeVector p,
const G4ThreeVector v,
G4double s1,
G4double s2 
)
protected

Definition at line 181 of file G4IntersectingCone.cc.

184{
185 static const G4double EPS = DBL_EPSILON; // Precision constant,
186 // originally it was 1E-6
187 G4double x0 = p.x(), y0 = p.y(), z0 = p.z();
188 G4double tx = v.x(), ty = v.y(), tz = v.z();
189
190 // Value of radical can be inaccurate due to loss of precision
191 // if to calculate the coefficiets a,b,c like the following:
192 // G4double a = tx*tx + ty*ty - sqr(B*tz);
193 // G4double b = 2*( x0*tx + y0*ty - B*(A + B*z0)*tz);
194 // G4double c = x0*x0 + y0*y0 - sqr(A + B*z0);
195 //
196 // For more accurate calculation of radical the coefficients
197 // are splitted in two components, radial and along z-axis
198 //
199 G4double ar = tx*tx + ty*ty;
200 G4double az = sqr(B*tz);
201 G4double br = 2*(x0*tx + y0*ty);
202 G4double bz = 2*B*(A + B*z0)*tz;
203 G4double cr = x0*x0 + y0*y0;
204 G4double cz = sqr(A + B*z0);
205
206 // Instead radical = b*b - 4*a*c
207 G4double arcz = 4*ar*cz;
208 G4double azcr = 4*az*cr;
209 G4double radical = (br*br - 4*ar*cr) + ((std::max(arcz,azcr) - 2*bz*br) + std::min(arcz,azcr));
210
211 // Find the coefficients
212 G4double a = ar - az;
213 G4double b = br - bz;
214 G4double c = cr - cz;
215
216 if (radical < -EPS*std::fabs(b)) { return 0; } // No solution
217
218 if (radical < EPS*std::fabs(b))
219 {
220 //
221 // The radical is roughly zero: check for special, very rare, cases
222 //
223 if (std::fabs(a) > 1/kInfinity)
224 {
225 if(B==0.) { return 0; }
226 if ( std::fabs(x0*ty - y0*tx) < std::fabs(EPS/B) )
227 {
228 *s1 = -0.5*b/a;
229 return 1;
230 }
231 return 0;
232 }
233 }
234 else
235 {
236 radical = std::sqrt(radical);
237 }
238
239 if (a > 1/kInfinity)
240 {
241 G4double sa, sb, q = -0.5*( b + (b < 0 ? -radical : +radical) );
242 sa = q/a;
243 sb = c/q;
244 if (sa < sb) { *s1 = sa; *s2 = sb; } else { *s1 = sb; *s2 = sa; }
245 if (A + B*(z0+(*s1)*tz) < 0) { return 0; }
246 return 2;
247 }
248 else if (a < -1/kInfinity)
249 {
250 G4double sa, sb, q = -0.5*( b + (b < 0 ? -radical : +radical) );
251 sa = q/a;
252 sb = c/q;
253 *s1 = (B*tz > 0)^(sa > sb) ? sb : sa;
254 return 1;
255 }
256 else if (std::fabs(b) < 1/kInfinity)
257 {
258 return 0;
259 }
260 else
261 {
262 *s1 = -c/b;
263 if (A + B*(z0+(*s1)*tz) < 0) { return 0; }
264 return 1;
265 }
266}
double z() const
double x() const
double y() const
static const G4double kInfinity
Definition: geomdefs.hh:41
#define EPS
T sqr(const T &x)
Definition: templates.hh:128
#define DBL_EPSILON
Definition: templates.hh:66

References A, B, DBL_EPSILON, EPS, kInfinity, G4INCL::Math::max(), G4INCL::Math::min(), sqr(), CLHEP::Hep3Vector::x(), CLHEP::Hep3Vector::y(), CLHEP::Hep3Vector::z(), and G4InuclParticleNames::z0.

Referenced by LineHitsCone().

◆ LineHitsCone2()

G4int G4IntersectingCone::LineHitsCone2 ( const G4ThreeVector p,
const G4ThreeVector v,
G4double s1,
G4double s2 
)
protected

Definition at line 291 of file G4IntersectingCone.cc.

294{
295 static const G4double EPS = DBL_EPSILON; // Precision constant,
296 // originally it was 1E-6
297 G4double x0 = p.x(), y0 = p.y(), z0 = p.z();
298 G4double tx = v.x(), ty = v.y(), tz = v.z();
299
300 // Special case which might not be so rare: B = 0 (precisely)
301 //
302 if (B==0)
303 {
304 if (std::fabs(tz) < 1/kInfinity) { return 0; }
305
306 *s1 = (A-z0)/tz;
307 return 1;
308 }
309
310 // Value of radical can be inaccurate due to loss of precision
311 // if to calculate the coefficiets a,b,c like the following:
312 // G4double a = tz*tz - B2*(tx*tx + ty*ty);
313 // G4double b = 2*( (z0-A)*tz - B2*(x0*tx + y0*ty) );
314 // G4double c = sqr(z0-A) - B2*( x0*x0 + y0*y0 );
315 //
316 // For more accurate calculation of radical the coefficients
317 // are splitted in two components, radial and along z-axis
318 //
319 G4double B2 = B*B;
320
321 G4double az = tz*tz;
322 G4double ar = B2*(tx*tx + ty*ty);
323 G4double bz = 2*(z0-A)*tz;
324 G4double br = 2*B2*(x0*tx + y0*ty);
325 G4double cz = sqr(z0-A);
326 G4double cr = B2*(x0*x0 + y0*y0);
327
328 // Instead radical = b*b - 4*a*c
329 G4double arcz = 4*ar*cz;
330 G4double azcr = 4*az*cr;
331 G4double radical = (br*br - 4*ar*cr) + ((std::max(arcz,azcr) - 2*bz*br) + std::min(arcz,azcr));
332
333 // Find the coefficients
334 G4double a = az - ar;
335 G4double b = bz - br;
336 G4double c = cz - cr;
337
338 if (radical < -EPS*std::fabs(b)) { return 0; } // No solution
339
340 if (radical < EPS*std::fabs(b))
341 {
342 //
343 // The radical is roughly zero: check for special, very rare, cases
344 //
345 if (std::fabs(a) > 1/kInfinity)
346 {
347 if ( std::fabs(x0*ty - y0*tx) < std::fabs(EPS/B) )
348 {
349 *s1 = -0.5*b/a;
350 return 1;
351 }
352 return 0;
353 }
354 }
355 else
356 {
357 radical = std::sqrt(radical);
358 }
359
360 if (a < -1/kInfinity)
361 {
362 G4double sa, sb, q = -0.5*( b + (b < 0 ? -radical : +radical) );
363 sa = q/a;
364 sb = c/q;
365 if (sa < sb) { *s1 = sa; *s2 = sb; } else { *s1 = sb; *s2 = sa; }
366 if ((z0 + (*s1)*tz - A)/B < 0) { return 0; }
367 return 2;
368 }
369 else if (a > 1/kInfinity)
370 {
371 G4double sa, sb, q = -0.5*( b + (b < 0 ? -radical : +radical) );
372 sa = q/a;
373 sb = c/q;
374 *s1 = (tz*B > 0)^(sa > sb) ? sb : sa;
375 return 1;
376 }
377 else if (std::fabs(b) < 1/kInfinity)
378 {
379 return 0;
380 }
381 else
382 {
383 *s1 = -c/b;
384 if ((z0 + (*s1)*tz - A)/B < 0) { return 0; }
385 return 1;
386 }
387}

References A, B, DBL_EPSILON, EPS, kInfinity, G4INCL::Math::max(), G4INCL::Math::min(), sqr(), CLHEP::Hep3Vector::x(), CLHEP::Hep3Vector::y(), CLHEP::Hep3Vector::z(), and G4InuclParticleNames::z0.

Referenced by LineHitsCone().

◆ RHi()

G4double G4IntersectingCone::RHi ( ) const
inline

Definition at line 55 of file G4IntersectingCone.hh.

55{ return rHi; }

References rHi.

◆ RLo()

G4double G4IntersectingCone::RLo ( ) const
inline

Definition at line 54 of file G4IntersectingCone.hh.

54{ return rLo; }

References rLo.

◆ ZHi()

G4double G4IntersectingCone::ZHi ( ) const
inline

Definition at line 57 of file G4IntersectingCone.hh.

57{ return zHi; }

References zHi.

Referenced by G4PolyconeSide::Extent(), and G4PolyhedraSide::Extent().

◆ ZLo()

G4double G4IntersectingCone::ZLo ( ) const
inline

Definition at line 56 of file G4IntersectingCone.hh.

56{ return zLo; }

References zLo.

Referenced by G4PolyconeSide::Extent(), and G4PolyhedraSide::Extent().

Field Documentation

◆ A

G4double G4IntersectingCone::A
protected

◆ B

G4double G4IntersectingCone::B
protected

Definition at line 74 of file G4IntersectingCone.hh.

Referenced by G4IntersectingCone(), LineHitsCone1(), and LineHitsCone2().

◆ rHi

G4double G4IntersectingCone::rHi
protected

Definition at line 70 of file G4IntersectingCone.hh.

Referenced by G4IntersectingCone(), HitOn(), and RHi().

◆ rLo

G4double G4IntersectingCone::rLo
protected

Definition at line 70 of file G4IntersectingCone.hh.

Referenced by G4IntersectingCone(), and RLo().

◆ type1

G4bool G4IntersectingCone::type1 = false
protected

Definition at line 72 of file G4IntersectingCone.hh.

Referenced by G4IntersectingCone(), HitOn(), and LineHitsCone().

◆ zHi

G4double G4IntersectingCone::zHi
protected

Definition at line 69 of file G4IntersectingCone.hh.

Referenced by G4IntersectingCone(), HitOn(), and ZHi().

◆ zLo

G4double G4IntersectingCone::zLo
protected

Definition at line 69 of file G4IntersectingCone.hh.

Referenced by G4IntersectingCone(), and ZLo().


The documentation for this class was generated from the following files: