G4IntersectingCone Class Reference

#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 = false
G4double	A
G4double	B

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.

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

~G4IntersectingCone()

G4IntersectingCone::~G4IntersectingCone ( )

virtual

Definition at line 76 of file G4IntersectingCone.cc.

{
}

G4IntersectingCone() [2/2]

G4IntersectingCone::G4IntersectingCone

(

__void__ &

)

Definition at line 69 of file G4IntersectingCone.cc.

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

Member Function Documentation

HitOn()

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

Definition at line 85 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;
}

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.

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

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.

{
  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();
 
  // Value of radical can be inaccurate due to loss of precision
  // if to calculate the coefficiets a,b,c like the following:
  //     G4double a = tx*tx + ty*ty - sqr(B*tz);
  //     G4double b = 2*( x0*tx + y0*ty - B*(A + B*z0)*tz);
  //     G4double c = x0*x0 + y0*y0 - sqr(A + B*z0);
  //
  // For more accurate calculation of radical the coefficients
  // are splitted in two components, radial and along z-axis
  //
  G4double ar = tx*tx + ty*ty;
  G4double az = sqr(B*tz);
  G4double br = 2*(x0*tx + y0*ty);
  G4double bz = 2*B*(A + B*z0)*tz;
  G4double cr = x0*x0 + y0*y0;
  G4double cz = sqr(A + B*z0);
 
  // Instead radical = b*b - 4*a*c
  G4double arcz = 4*ar*cz;
  G4double azcr = 4*az*cr;
  G4double radical = (br*br - 4*ar*cr) + ((std::max(arcz,azcr) - 2*bz*br) + std::min(arcz,azcr));
 
  // Find the coefficients
  G4double a = ar - az;
  G4double b = br - bz;
  G4double c = cr - cz;
 
  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;
  }
}

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.

{
  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;
  }
 
  // Value of radical can be inaccurate due to loss of precision
  // if to calculate the coefficiets a,b,c like the following:
  //   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 );
  //
  // For more accurate calculation of radical the coefficients
  // are splitted in two components, radial and along z-axis
  //
  G4double B2 = B*B;
 
  G4double az = tz*tz;
  G4double ar = B2*(tx*tx + ty*ty);
  G4double bz = 2*(z0-A)*tz;
  G4double br = 2*B2*(x0*tx + y0*ty);
  G4double cz = sqr(z0-A);
  G4double cr = B2*(x0*x0 + y0*y0);
 
  // Instead radical = b*b - 4*a*c
  G4double arcz = 4*ar*cz;
  G4double azcr = 4*az*cr;
  G4double radical = (br*br - 4*ar*cr) + ((std::max(arcz,azcr) - 2*bz*br) + std::min(arcz,azcr));
 
  // Find the coefficients
  G4double a = az - ar;
  G4double b = bz - br;
  G4double c = cz - cr;
 
  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;
  }
}

Referenced by LineHitsCone().

RHi()

G4double G4IntersectingCone::RHi ( ) const

inline

Definition at line 55 of file G4IntersectingCone.hh.

{ return rHi; }

RLo()

G4double G4IntersectingCone::RLo ( ) const

inline

Definition at line 54 of file G4IntersectingCone.hh.

{ return rLo; }

ZHi()

G4double G4IntersectingCone::ZHi ( ) const

inline

Definition at line 57 of file G4IntersectingCone.hh.

{ return zHi; }

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

ZLo()

G4double G4IntersectingCone::ZLo ( ) const

inline

Definition at line 56 of file G4IntersectingCone.hh.

{ return zLo; }

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

Member Data Documentation

A

G4double G4IntersectingCone::A

protected

Definition at line 74 of file G4IntersectingCone.hh.

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

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:

• G4IntersectingCone.hh
• G4IntersectingCone.cc