#include <UIntersectingCone.hh>

Public Member Functions
	UIntersectingCone (const double r[2], const double z[2])

virtual	~UIntersectingCone ()

int	LineHitsCone (const UVector3 &p, const UVector3 &v, double &s1, double &s2)

bool	HitOn (const double r, const double z)

double	RLo () const

double	RHi () const

double	ZLo () const

double	ZHi () const

Protected Member Functions
int	LineHitsCone1 (const UVector3 &p, const UVector3 &v, double &s1, double &s2)

int	LineHitsCone1Optimized (const UVector3 &p, const UVector3 &v, double &s1, double &s2)

int	LineHitsCone2 (const UVector3 &p, const UVector3 &v, double &s1, double &s2)

Protected Attributes
double	zLo

double	zHi

double	rLo

double	rHi

bool	type1

double	A

double	B

Static Protected Attributes
static const double	EpsilonQuad = 1.0 / 9.0E99

Detailed Description

Definition at line 27 of file UIntersectingCone.hh.

Constructor & Destructor Documentation

UIntersectingCone::UIntersectingCone	(	const double	r[2],
		const double	z[2]
	)

Definition at line 29 of file UIntersectingCone.cc.

References A, B, rHi, rLo, VUSolid::Tolerance(), type1, zHi, and zLo.

 {
   static const double halfCarTolerance
     = 0.5 * VUSolid::Tolerance(); // UGeometryTolerance::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;
   }
 }

UIntersectingCone::~UIntersectingCone ( )

virtual

Definition at line 92 of file UIntersectingCone.cc.

93 {

94 }

Member Function Documentation

bool UIntersectingCone::HitOn	(	const double	r,
		const double	z
	)

Definition at line 103 of file UIntersectingCone.cc.

References rHi, type1, and zHi.

Referenced by UPolyconeSide::PointOnCone().

 {
   //
   // 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;
 }

int UIntersectingCone::LineHitsCone	(	const UVector3 &	p,
		const UVector3 &	v,
		double &	s1,
		double &	s2
	)

Definition at line 129 of file UIntersectingCone.cc.

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

Referenced by UPolyhedraSide::Distance(), UPolyconeSide::Distance(), and UPolyhedraSide::LineHitsSegments().

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

int UIntersectingCone::LineHitsCone1	(	const UVector3 &	p,
		const UVector3 &	v,
		double &	s1,
		double &	s2
	)

protected

Definition at line 209 of file UIntersectingCone.cc.

References test::a, A, test::b, B, test::c, EpsilonQuad, UUtils::sqr(), UVector3::x, UVector3::y, UVector3::z, and G4InuclParticleNames::z0.

Referenced by LineHitsCone().

 {
   double x0 = p.x, y0 = p.y, z0 = p.z;
   double tx = v.x, ty = v.y, tz = v.z;
 
   double a = tx * tx + ty * ty - UUtils::sqr(B * tz);
   double b = 2 * (x0 * tx + y0 * ty - (A * B + B * B * z0) * tz);
   double c = x0 * x0 + y0 * y0 - UUtils::sqr(A + B * z0);
 
   double radical = b * b - 4 * a * c;
   double radicalSqrt;
 
   double minRadical = 1E-6 * std::fabs(b);
 
   if (radical < -minRadical)
   {
     return 0;  // No solution
   }
 
   if (radical < minRadical)
   {
     //
     // The radical is roughly zero: check for special, very rare, cases
     //
     if (std::fabs(a) > EpsilonQuad)
     {
       if (B == 0.)
       {
         return 0;
       }
       if (std::fabs(x0 * ty - y0 * tx) < std::fabs(1E-6 / B))
       {
         s1 = -0.5 * b / a;
         return 1;
       }
       return 0;
     }
     radicalSqrt = radical; //TODO: check this case
   }
   else
   {
     radicalSqrt = std::sqrt(radical);
   }
 
   if (a > EpsilonQuad)
   {
     double sa, sb, q = -0.5 * (b + (b < 0 ? -radicalSqrt : +radicalSqrt));
     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;
     }
     // if ((z0 + (s1)*tz - A)/B < 0)  { return 0; } // these lines are equivalent
     return 2;
   }
   else if (a < -EpsilonQuad)
   {
     double sa, sb, q = -0.5 * (b + (b < 0 ? -radicalSqrt : +radicalSqrt));
     sa = q / a;
     sb = c / q;
     s1 = (tz * B > 0) ^ (sa > sb) ? sb : sa;
     return 1;
   }
   else if (std::fabs(b) < EpsilonQuad)
   {
     return 0;
   }
   else
   {
     s1 = -c / b;
     if (A + B * (z0 + (s1)*tz) < 0)
     {
       return 0;
     }
     return 1;
   }
 }

int UIntersectingCone::LineHitsCone1Optimized	(	const UVector3 &	p,
		const UVector3 &	v,
		double &	s1,
		double &	s2
	)

protected

Definition at line 301 of file UIntersectingCone.cc.

References test::a, A, test::b, B, test::c, EpsilonQuad, UUtils::sqr(), UVector3::x, UVector3::y, UVector3::z, and G4InuclParticleNames::z0.

 {
   double x0 = p.x, y0 = p.y, z0 = p.z;
   double tx = v.x, ty = v.y, tz = v.z;
 
   double a = tx * tx + ty * ty - UUtils::sqr(B * tz);
   double b = 2 * (x0 * tx + y0 * ty - (A * B + B * B * z0) * tz);
   double c = x0 * x0 + y0 * y0 - UUtils::sqr(A + B * z0);
 
   double radical = b * b - 4 * a * c;
 
   double minRadical = 1E-6 * std::fabs(b);
 
   if (radical < -minRadical)
   {
     return 0;  // No solution
   }
 
   if (std::fabs(a) > EpsilonQuad)
   {
     if (radical < minRadical)
     {
       //
       // The radical is roughly zero: check for special, very rare, cases
       //
 
       if (B == 0.)
       {
         return 0;
       }
       if (std::fabs(x0 * ty - y0 * tx) < std::fabs(1E-6 / B))
       {
         s1 = -0.5 * b / a;
         return 1;
       }
       return 0;
     }
     else
     {
       double radicalSqrt = std::sqrt(radical);
 
       if (a > 0)
       {
         double sa, sb, q = -0.5 * (b + (b < 0 ? -radicalSqrt : +radicalSqrt));
         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;
         }
         // if ((z0 + (s1)*tz - A)/B < 0)  { return 0; } // these lines are equivalent
         return 2;
       }
       else
       {
         double sa, sb, q = -0.5 * (b + (b < 0 ? -radicalSqrt : +radicalSqrt));
         sa = q / a;
         sb = c / q;
         s1 = (tz * B > 0) ^ (sa > sb) ? sb : sa;
         return 1;
       }
 
     }
   }
 
   if (std::fabs(b) < EpsilonQuad)
   {
     return 0;
   }
   else
   {
     s1 = -c / b;
     if (A + B * (z0 + (s1)*tz) < 0)
     {
       return 0;
     }
     return 1;
   }
 }

int UIntersectingCone::LineHitsCone2	(	const UVector3 &	p,
		const UVector3 &	v,
		double &	s1,
		double &	s2
	)

protected

Definition at line 416 of file UIntersectingCone.cc.

References test::a, A, test::b, B, test::c, EpsilonQuad, UUtils::sqr(), UVector3::x, UVector3::y, UVector3::z, and G4InuclParticleNames::z0.

Referenced by LineHitsCone().

 {
   double x0 = p.x, y0 = p.y, z0 = p.z;
   double 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) < EpsilonQuad)
     {
       return 0;
     }
 
     s1 = (A - z0) / tz;
     return 1;
   }
 
   double B2 = B * B;
 
   double a = tz * tz - B2 * (tx * tx + ty * ty);
   double b = 2 * ((z0 - A) * tz - B2 * (x0 * tx + y0 * ty));
   double c = UUtils::sqr(z0 - A) - B2 * (x0 * x0 + y0 * y0);
 
   double radical = b * b - 4 * a * c;
 
   if (radical < -1E-6 * std::fabs(b))
   {
     return 0;  // No solution
   }
 
   if (radical < 1E-6 * std::fabs(b))
   {
     //
     // The radical is roughly zero: check for special, very rare, cases
     //
     if (std::fabs(a) > EpsilonQuad)
     {
       if (std::fabs(x0 * ty - y0 * tx) < std::fabs(1E-6 / B))
       {
         s1 = -0.5 * b / a;
         return 1;
       }
       return 0;
     }
   }
   else
   {
     radical = std::sqrt(radical);
   }
 
   if (a < -EpsilonQuad)
   {
     double 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 > EpsilonQuad)
   {
     double 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) < EpsilonQuad)
   {
     return 0;
   }
   else
   {
     s1 = -c / b;
     if ((z0 + (s1)*tz - A) / B < 0)
     {
       return 0;
     }
     return 1;
   }
 }