#include <G4IntersectingCone.hh>

Public Member Functions
	G4IntersectingCone (const G4double r[2], const G4double z[2])

virtual	~G4IntersectingCone ()

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

G4bool	HitOn (const G4double r, const G4double z)

G4double	RLo () const

G4double	RHi () const

G4double	ZLo () const

G4double	ZHi () const

	G4IntersectingCone (__void__ &)

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	zLo

G4double	zHi

G4double	rLo

G4double	rHi

G4bool	type1

G4double	A

G4double	B

Detailed Description

Definition at line 51 of file G4IntersectingCone.hh.

Constructor & Destructor Documentation

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

Definition at line 46 of file G4IntersectingCone.cc.

 { 
   const G4double halfCarTolerance
     = 0.5 * G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
 
   //
   // What type of cone are we?
   //
   type1 = (std::fabs(z[1]-z[0]) > std::fabs(r[1]-r[0]));
   
   if (type1)
   {
     B = (r[1]-r[0])/(z[1]-z[0]);      // tube like
     A = 0.5*( r[1]+r[0] - B*(z[1]+z[0]) );
   }
   else
   {
     B = (z[1]-z[0])/(r[1]-r[0]);      // disk like
     A = 0.5*( z[1]+z[0] - B*(r[1]+r[0]) );
   }
   //
   // Calculate extent
   //
   if (r[0] < r[1])
   {
     rLo = r[0]-halfCarTolerance; rHi = r[1]+halfCarTolerance;
   }
   else
   {
     rLo = r[1]-halfCarTolerance; rHi = r[0]+halfCarTolerance;
   }
   
   if (z[0] < z[1])
   {
     zLo = z[0]-halfCarTolerance; zHi = z[1]+halfCarTolerance;
   }
   else
   {
     zLo = z[1]-halfCarTolerance; zHi = z[0]+halfCarTolerance;
   }
 }

Here is the call graph for this function:

G4IntersectingCone::~G4IntersectingCone ( )

virtual

Definition at line 103 of file G4IntersectingCone.cc.

104 {

105 }

G4IntersectingCone::G4IntersectingCone ( __void__ & )

Definition at line 94 of file G4IntersectingCone.cc.

   : zLo(0.), zHi(0.), rLo(0.), rHi(0.), type1(false), A(0.), B(0.)
 {
 }

Member Function Documentation

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

Definition at line 114 of file G4IntersectingCone.cc.

 {
   //
   // Be careful! The inequalities cannot be "<=" and ">=" here without
   // punching a tiny hole in our shape!
   //
   if (type1)
   {
     if (z < zLo || z > zHi) return false;
   }
   else
   {
     if (r < rLo || r > rHi) return false;
   }
 
   return true;
 }

Here is the caller graph for this function:

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

Definition at line 140 of file G4IntersectingCone.cc.

 {
   if (type1)
   {
     return LineHitsCone1( p, v, s1, s2 );
   }
   else
   {
     return LineHitsCone2( p, v, s1, s2 );
   }
 }

Here is the call graph for this function:

Here is the caller graph for this function:

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

protected

Definition at line 214 of file G4IntersectingCone.cc.

 {
   static const G4double EPS = DBL_EPSILON; // Precision constant,
                                            // originally it was 1E-6
   G4double x0 = p.x(), y0 = p.y(), z0 = p.z();
   G4double tx = v.x(), ty = v.y(), tz = v.z();
 
   G4double a = tx*tx + ty*ty - sqr(B*tz);
   G4double b = 2*( x0*tx + y0*ty - (A*B + B*B*z0)*tz);
   G4double c = x0*x0 + y0*y0 - sqr(A + B*z0);
   
   G4double radical = b*b - 4*a*c;
  
   if (radical < -EPS*std::fabs(b))  { return 0; }    // No solution
   
   if (radical < EPS*std::fabs(b))
   {
     //
     // The radical is roughly zero: check for special, very rare, cases
     //
     if (std::fabs(a) > 1/kInfinity)
       {
       if(B==0.) { return 0; }
       if ( std::fabs(x0*ty - y0*tx) < std::fabs(EPS/B) )
       {
          *s1 = -0.5*b/a;
          return 1;
       }
       return 0;
     }
   }
   else
   {
     radical = std::sqrt(radical);
   }
   
   if (a > 1/kInfinity)
   {
     G4double sa, sb, q = -0.5*( b + (b < 0 ? -radical : +radical) );
     sa = q/a;
     sb = c/q;
     if (sa < sb) { *s1 = sa; *s2 = sb; } else { *s1 = sb; *s2 = sa; }
     if (A + B*(z0+(*s1)*tz) < 0)  { return 0; }
     return 2;
   }
   else if (a < -1/kInfinity)
   {
     G4double sa, sb, q = -0.5*( b + (b < 0 ? -radical : +radical) );
     sa = q/a;
     sb = c/q;
     *s1 = (B*tz > 0)^(sa > sb) ? sb : sa;
     return 1;
   }
   else if (std::fabs(b) < 1/kInfinity)
   {
     return 0;
   }
   else
   {
     *s1 = -c/b;
     if (A + B*(z0+(*s1)*tz) < 0)  { return 0; }
     return 1;
   }
 }

Here is the call graph for this function:

Here is the caller graph for this function:

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

protected

Definition at line 306 of file G4IntersectingCone.cc.

 {
   static const G4double EPS = DBL_EPSILON; // Precision constant,
                                            // originally it was 1E-6
   G4double x0 = p.x(), y0 = p.y(), z0 = p.z();
   G4double tx = v.x(), ty = v.y(), tz = v.z();
   
   // Special case which might not be so rare: B = 0 (precisely)
   //
   if (B==0)
   {
     if (std::fabs(tz) < 1/kInfinity)  { return 0; }
     
     *s1 = (A-z0)/tz;
     return 1;
   }
 
   G4double B2 = B*B;
 
   G4double a = tz*tz - B2*(tx*tx + ty*ty);
   G4double b = 2*( (z0-A)*tz - B2*(x0*tx + y0*ty) );
   G4double c = sqr(z0-A) - B2*( x0*x0 + y0*y0 );
   
   G4double radical = b*b - 4*a*c;
  
   if (radical < -EPS*std::fabs(b)) { return 0; }   // No solution
   
   if (radical < EPS*std::fabs(b))
   {
     //
     // The radical is roughly zero: check for special, very rare, cases
     //
     if (std::fabs(a) > 1/kInfinity)
     {
       if ( std::fabs(x0*ty - y0*tx) < std::fabs(EPS/B) )
       {
         *s1 = -0.5*b/a;
         return 1;
       }
       return 0;
     }
   }
   else
   {
     radical = std::sqrt(radical);
   }
   
   if (a < -1/kInfinity)
   {
     G4double sa, sb, q = -0.5*( b + (b < 0 ? -radical : +radical) );
     sa = q/a;
     sb = c/q;
     if (sa < sb) { *s1 = sa; *s2 = sb; } else { *s1 = sb; *s2 = sa; }
     if ((z0 + (*s1)*tz  - A)/B < 0)  { return 0; }
     return 2;
   }
   else if (a > 1/kInfinity)
   {
     G4double sa, sb, q = -0.5*( b + (b < 0 ? -radical : +radical) );
     sa = q/a;
     sb = c/q;
     *s1 = (tz*B > 0)^(sa > sb) ? sb : sa;
     return 1;
   }
   else if (std::fabs(b) < 1/kInfinity)
   {
     return 0;
   }
   else
   {
     *s1 = -c/b;
     if ((z0 + (*s1)*tz  - A)/B < 0)  { return 0; }
     return 1;
   }
 }