G4Torus Class Reference

#include <G4Torus.hh>

Inheritance diagram for G4Torus:

Public Member Functions
	G4Torus (const G4String &pName, G4double pRmin, G4double pRmax, G4double pRtor, G4double pSPhi, G4double pDPhi)
	~G4Torus ()
G4double	GetRmin () const
G4double	GetRmax () const
G4double	GetRtor () const
G4double	GetSPhi () const
G4double	GetDPhi () const
G4double	GetCubicVolume ()
G4double	GetSurfaceArea ()
EInside	Inside (const G4ThreeVector &p) const
G4bool	CalculateExtent (const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pmin, G4double &pmax) const
void	ComputeDimensions (G4VPVParameterisation *p, const G4int n, const G4VPhysicalVolume *pRep)
G4ThreeVector	SurfaceNormal (const G4ThreeVector &p) const
G4double	DistanceToIn (const G4ThreeVector &p, const G4ThreeVector &v) const
G4double	DistanceToIn (const G4ThreeVector &p) const
G4double	DistanceToOut (const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=G4bool(false), G4bool *validNorm=0, G4ThreeVector *n=0) const
G4double	DistanceToOut (const G4ThreeVector &p) const
G4GeometryType	GetEntityType () const
G4ThreeVector	GetPointOnSurface () const
G4VSolid *	Clone () const
std::ostream &	StreamInfo (std::ostream &os) const
void	DescribeYourselfTo (G4VGraphicsScene &scene) const
G4Polyhedron *	CreatePolyhedron () const
G4NURBS *	CreateNURBS () const
void	SetAllParameters (G4double pRmin, G4double pRmax, G4double pRtor, G4double pSPhi, G4double pDPhi)
	G4Torus (__void__ &)
	G4Torus (const G4Torus &rhs)
G4Torus &	operator= (const G4Torus &rhs)

---

Detailed Description

Definition at line 102 of file G4Torus.hh.

---

Constructor & Destructor Documentation

G4Torus::G4Torus	(	const G4String &	pName,
		G4double	pRmin,
		G4double	pRmax,
		G4double	pRtor,
		G4double	pSPhi,
		G4double	pDPhi
	)

Definition at line 84 of file G4Torus.cc.

References SetAllParameters().

Referenced by Clone().

  : G4CSGSolid(pName)
{
  SetAllParameters(pRmin, pRmax, pRtor, pSPhi, pDPhi);
}

G4Torus::~G4Torus

(

)

Definition at line 193 of file G4Torus.cc.

{}

G4Torus::G4Torus

(

__void__ &

)

Definition at line 182 of file G4Torus.cc.

  : G4CSGSolid(a), fRmin(0.), fRmax(0.), fRtor(0.), fSPhi(0.),
    fDPhi(0.), fRminTolerance(0.), fRmaxTolerance(0. ),
    kRadTolerance(0.), kAngTolerance(0.)
{
}

G4Torus::G4Torus

(

const G4Torus &

rhs

)

Definition at line 200 of file G4Torus.cc.

  : G4CSGSolid(rhs), fRmin(rhs.fRmin),fRmax(rhs.fRmax),
    fRtor(rhs.fRtor),fSPhi(rhs.fSPhi),fDPhi(rhs.fDPhi),
    fRminTolerance(rhs.fRminTolerance), fRmaxTolerance(rhs.fRmaxTolerance),
    kRadTolerance(rhs.kRadTolerance), kAngTolerance(rhs.kAngTolerance)
{
}

---

Member Function Documentation

G4bool G4Torus::CalculateExtent	(	const EAxis	pAxis,
		const G4VoxelLimits &	pVoxelLimit,
		const G4AffineTransform &	pTransform,
		G4double &	pmin,
		G4double &	pmax
	)			const [virtual]

Implements G4VSolid.

Definition at line 412 of file G4Torus.cc.

References G4VSolid::ClipBetweenSections(), G4VSolid::ClipCrossSection(), G4VoxelLimits::GetMaxExtent(), G4VoxelLimits::GetMaxXExtent(), G4VoxelLimits::GetMaxYExtent(), G4VoxelLimits::GetMaxZExtent(), G4VoxelLimits::GetMinExtent(), G4VoxelLimits::GetMinXExtent(), G4VoxelLimits::GetMinYExtent(), G4VoxelLimits::GetMinZExtent(), Inside(), G4AffineTransform::Inverse(), G4AffineTransform::IsRotated(), G4VoxelLimits::IsXLimited(), G4VoxelLimits::IsYLimited(), G4VoxelLimits::IsZLimited(), G4VSolid::kCarTolerance, kOutside, kXAxis, kYAxis, kZAxis, G4AffineTransform::NetTranslation(), and G4AffineTransform::TransformPoint().

{
  if ((!pTransform.IsRotated()) && (fDPhi==twopi) && (fRmin==0))
  {
    // Special case handling for unrotated solid torus
    // Compute x/y/z mins and maxs for bounding box respecting limits,
    // with early returns if outside limits. Then switch() on pAxis,
    // and compute exact x and y limit for x/y case
      
    G4double xoffset,xMin,xMax;
    G4double yoffset,yMin,yMax;
    G4double zoffset,zMin,zMax;

    G4double RTorus,delta,diff1,diff2,maxDiff,newMin,newMax;
    G4double xoff1,xoff2,yoff1,yoff2;

    xoffset = pTransform.NetTranslation().x();
    xMin    = xoffset - fRmax - fRtor ;
    xMax    = xoffset + fRmax + fRtor ;

    if (pVoxelLimit.IsXLimited())
    {
      if ( (xMin > pVoxelLimit.GetMaxXExtent()+kCarTolerance)
        || (xMax < pVoxelLimit.GetMinXExtent()-kCarTolerance) )
        return false ;
      else
      {
        if (xMin < pVoxelLimit.GetMinXExtent())
        {
          xMin = pVoxelLimit.GetMinXExtent() ;
        }
        if (xMax > pVoxelLimit.GetMaxXExtent())
        {
          xMax = pVoxelLimit.GetMaxXExtent() ;
        }
      }
    }
    yoffset = pTransform.NetTranslation().y();
    yMin    = yoffset - fRmax - fRtor ;
    yMax    = yoffset + fRmax + fRtor ;

    if (pVoxelLimit.IsYLimited())
    {
      if ( (yMin > pVoxelLimit.GetMaxYExtent()+kCarTolerance)
        || (yMax < pVoxelLimit.GetMinYExtent()-kCarTolerance) )
      {
        return false ;
      }
      else
      {
        if (yMin < pVoxelLimit.GetMinYExtent() )
        { 
          yMin = pVoxelLimit.GetMinYExtent() ;
        }
        if (yMax > pVoxelLimit.GetMaxYExtent() )
        {
          yMax = pVoxelLimit.GetMaxYExtent() ;
        }
      }
    }
    zoffset = pTransform.NetTranslation().z() ;
    zMin    = zoffset - fRmax ;
    zMax    = zoffset + fRmax ;

    if (pVoxelLimit.IsZLimited())
    {
      if ( (zMin > pVoxelLimit.GetMaxZExtent()+kCarTolerance)
        || (zMax < pVoxelLimit.GetMinZExtent()-kCarTolerance) )
      {
        return false ;
      }
      else
      {
        if (zMin < pVoxelLimit.GetMinZExtent() )
        {
          zMin = pVoxelLimit.GetMinZExtent() ;
        }
        if (zMax > pVoxelLimit.GetMaxZExtent() )
        {
          zMax = pVoxelLimit.GetMaxZExtent() ;
        }
      }
    }

    // Known to cut cylinder
    
    switch (pAxis)
    {
      case kXAxis:
        yoff1=yoffset-yMin;
        yoff2=yMax-yoffset;
        if ( yoff1 >= 0 && yoff2 >= 0 )
        {
          // Y limits cross max/min x => no change
          //
          pMin = xMin ;
          pMax = xMax ;
        }
        else
        {
          // Y limits don't cross max/min x => compute max delta x,
          // hence new mins/maxs
          //
       
          RTorus=fRmax+fRtor;
          delta   = RTorus*RTorus - yoff1*yoff1;
          diff1   = (delta>0.) ? std::sqrt(delta) : 0.;
          delta   = RTorus*RTorus - yoff2*yoff2;
          diff2   = (delta>0.) ? std::sqrt(delta) : 0.;
          maxDiff = (diff1 > diff2) ? diff1:diff2 ;
          newMin  = xoffset - maxDiff ;
          newMax  = xoffset + maxDiff ;
          pMin    = (newMin < xMin) ? xMin : newMin ;
          pMax    = (newMax > xMax) ? xMax : newMax ;
        }
        break;

      case kYAxis:
        xoff1 = xoffset - xMin ;
        xoff2 = xMax - xoffset ;
        if (xoff1 >= 0 && xoff2 >= 0 )
        {
          // X limits cross max/min y => no change
          //
          pMin = yMin ;
          pMax = yMax ;
        } 
        else
        {
          // X limits don't cross max/min y => compute max delta y,
          // hence new mins/maxs
          //
          RTorus=fRmax+fRtor;
          delta   = RTorus*RTorus - xoff1*xoff1;
          diff1   = (delta>0.) ? std::sqrt(delta) : 0.;
          delta   = RTorus*RTorus - xoff2*xoff2;
          diff2   = (delta>0.) ? std::sqrt(delta) : 0.;
          maxDiff = (diff1 > diff2) ? diff1 : diff2 ;
          newMin  = yoffset - maxDiff ;
          newMax  = yoffset + maxDiff ;
          pMin    = (newMin < yMin) ? yMin : newMin ;
          pMax    = (newMax > yMax) ? yMax : newMax ;
        }
        break;

      case kZAxis:
        pMin=zMin;
        pMax=zMax;
        break;

      default:
        break;
    }
    pMin -= kCarTolerance ;
    pMax += kCarTolerance ;

    return true;
  }
  else
  {
    G4int i, noEntries, noBetweenSections4 ;
    G4bool existsAfterClip = false ;

    // Calculate rotated vertex coordinates

    G4ThreeVectorList *vertices ;
    G4int noPolygonVertices ;  // will be 4 
    vertices = CreateRotatedVertices(pTransform,noPolygonVertices) ;

    pMin = +kInfinity ;
    pMax = -kInfinity ;

    noEntries          = vertices->size() ;
    noBetweenSections4 = noEntries - noPolygonVertices ;
      
    for (i=0;i<noEntries;i+=noPolygonVertices)
    {
      ClipCrossSection(vertices,i,pVoxelLimit,pAxis,pMin,pMax);
    }    
    for (i=0;i<noBetweenSections4;i+=noPolygonVertices)
    {
      ClipBetweenSections(vertices,i,pVoxelLimit,pAxis,pMin,pMax);
    }
    if (pMin!=kInfinity||pMax!=-kInfinity)
    {
      existsAfterClip = true ; // Add 2*tolerance to avoid precision troubles
      pMin           -= kCarTolerance ;
      pMax           += kCarTolerance ;
    }
    else
    {
      // Check for case where completely enveloping clipping volume
      // If point inside then we are confident that the solid completely
      // envelopes the clipping volume. Hence set min/max extents according
      // to clipping volume extents along the specified axis.

      G4ThreeVector clipCentre(
          (pVoxelLimit.GetMinXExtent()+pVoxelLimit.GetMaxXExtent())*0.5,
          (pVoxelLimit.GetMinYExtent()+pVoxelLimit.GetMaxYExtent())*0.5,
          (pVoxelLimit.GetMinZExtent()+pVoxelLimit.GetMaxZExtent())*0.5  ) ;
        
      if (Inside(pTransform.Inverse().TransformPoint(clipCentre)) != kOutside )
      {
        existsAfterClip = true ;
        pMin            = pVoxelLimit.GetMinExtent(pAxis) ;
        pMax            = pVoxelLimit.GetMaxExtent(pAxis) ;
      }
    }
    delete vertices;
    return existsAfterClip;
  }
}

G4VSolid * G4Torus::Clone

(

)

const [virtual]

Reimplemented from G4VSolid.

Definition at line 1703 of file G4Torus.cc.

References G4Torus().

{
  return new G4Torus(*this);
}

void G4Torus::ComputeDimensions	(	G4VPVParameterisation *	p,
		const G4int	n,
		const G4VPhysicalVolume *	pRep
	)			[virtual]

Reimplemented from G4VSolid.

Definition at line 237 of file G4Torus.cc.

References G4VPVParameterisation::ComputeDimensions().

{
  p->ComputeDimensions(*this,n,pRep);
}

G4NURBS * G4Torus::CreateNURBS

(

)

const [virtual]

Reimplemented from G4VSolid.

Definition at line 1793 of file G4Torus.cc.

{
  G4NURBS* pNURBS;
  if (fRmin != 0) 
  {
    if (fDPhi >= twopi) 
    {
      pNURBS = new G4NURBStube(fRmin, fRmax, fRtor);
    }
    else 
    {
      pNURBS = new G4NURBStubesector(fRmin, fRmax, fRtor, fSPhi, fSPhi + fDPhi);
    }
  }
  else 
  {
    if (fDPhi >= twopi) 
    {
      pNURBS = new G4NURBScylinder (fRmax, fRtor);
    }
    else 
    {
      const G4double epsilon = 1.e-4; // Cylinder sector not yet available!
      pNURBS = new G4NURBStubesector (epsilon, fRmax, fRtor,
                                      fSPhi, fSPhi + fDPhi);
    }
  }
  return pNURBS;
}

G4Polyhedron * G4Torus::CreatePolyhedron

(

)

const [virtual]

Reimplemented from G4VSolid.

Definition at line 1788 of file G4Torus.cc.

{
  return new G4PolyhedronTorus (fRmin, fRmax, fRtor, fSPhi, fDPhi);
}

void G4Torus::DescribeYourselfTo

(

G4VGraphicsScene &

scene

)

const [virtual]

Implements G4VSolid.

Definition at line 1783 of file G4Torus.cc.

References G4VGraphicsScene::AddSolid().

{
  scene.AddSolid (*this);
}

G4double G4Torus::DistanceToIn

(

const G4ThreeVector &

p

)

const [virtual]

Implements G4VSolid.

Definition at line 1139 of file G4Torus.cc.

{
  G4double safe=0.0, safe1, safe2 ;
  G4double phiC, cosPhiC, sinPhiC, safePhi, ePhi, cosPsi ;
  G4double rho2, rho, pt2, pt ;
    
  rho2 = p.x()*p.x() + p.y()*p.y() ;
  rho  = std::sqrt(rho2) ;
  pt2  = std::fabs(rho2 + p.z()*p.z() + fRtor*fRtor - 2*fRtor*rho) ;
  pt   = std::sqrt(pt2) ;

  safe1 = fRmin - pt ;
  safe2 = pt - fRmax ;

  if (safe1 > safe2)  { safe = safe1; }
  else                { safe = safe2; }

  if ( fDPhi < twopi && rho )
  {
    phiC    = fSPhi + fDPhi*0.5 ;
    cosPhiC = std::cos(phiC) ;
    sinPhiC = std::sin(phiC) ;
    cosPsi  = (p.x()*cosPhiC + p.y()*sinPhiC)/rho ;

    if (cosPsi < std::cos(fDPhi*0.5) ) // Psi=angle from central phi to point
    {                                  // Point lies outside phi range
      if ((p.y()*cosPhiC - p.x()*sinPhiC) <= 0 )
      {
        safePhi = std::fabs(p.x()*std::sin(fSPhi) - p.y()*std::cos(fSPhi)) ;
      }
      else
      {
        ePhi    = fSPhi + fDPhi ;
        safePhi = std::fabs(p.x()*std::sin(ePhi) - p.y()*std::cos(ePhi)) ;
      }
      if (safePhi > safe)  { safe = safePhi ; }
    }
  }
  if (safe < 0 )  { safe = 0 ; }
  return safe;
}

G4double G4Torus::DistanceToIn	(	const G4ThreeVector &	p,
		const G4ThreeVector &	v
	)			const [virtual]

Implements G4VSolid.

Definition at line 987 of file G4Torus.cc.

References G4VSolid::kCarTolerance.

Referenced by SurfaceNormal().

{

  G4double snxt=kInfinity, sphi=kInfinity; // snxt = default return value

  G4double  sd[4] ;

  // Precalculated trig for phi intersections - used by r,z intersections to
  //                                            check validity

  G4bool seg;        // true if segmented
  G4double hDPhi;    // half dphi
  G4double cPhi,sinCPhi=0.,cosCPhi=0.;  // central phi

  G4double tolORMin2;  // `generous' radii squared
  G4double tolORMax2;

  static const G4double halfCarTolerance = 0.5*kCarTolerance;
 
  G4double Dist,xi,yi,zi,rhoi2,it2; // Intersection point variables

  G4double Comp;
  G4double cosSPhi,sinSPhi;       // Trig for phi start intersect
  G4double ePhi,cosEPhi,sinEPhi;  // for phi end intersect

  // Set phi divided flag and precalcs
  //
  if ( fDPhi < twopi )
  {
    seg        = true ;
    hDPhi      = 0.5*fDPhi ;    // half delta phi
    cPhi       = fSPhi + hDPhi ;
    sinCPhi    = std::sin(cPhi) ;
    cosCPhi    = std::cos(cPhi) ;
  }
  else
  {
    seg = false ;
  }

  if (fRmin > fRminTolerance) // Calculate tolerant rmin and rmax
  {
    tolORMin2 = (fRmin - fRminTolerance)*(fRmin - fRminTolerance) ;
  }
  else
  {
    tolORMin2 = 0 ;
  }
  tolORMax2 = (fRmax + fRmaxTolerance)*(fRmax + fRmaxTolerance) ;

  // Intersection with Rmax (possible return) and Rmin (must also check phi)

  G4double Rtor2 = fRtor*fRtor ;

  snxt = SolveNumericJT(p,v,fRmax,true);

  if (fRmin)  // Possible Rmin intersection
  {
    sd[0] = SolveNumericJT(p,v,fRmin,true);
    if ( sd[0] < snxt )  { snxt = sd[0] ; }
  }

  //
  // Phi segment intersection
  //
  // o Tolerant of points inside phi planes by up to kCarTolerance*0.5
  //
  // o NOTE: Large duplication of code between sphi & ephi checks
  //         -> only diffs: sphi -> ephi, Comp -> -Comp and half-plane
  //            intersection check <=0 -> >=0
  //         -> use some form of loop Construct ?

  if (seg)
  {
    sinSPhi = std::sin(fSPhi) ; // First phi surface ('S'tarting phi)
    cosSPhi = std::cos(fSPhi) ;
    Comp    = v.x()*sinSPhi - v.y()*cosSPhi ;  // Component in outwards
                                               // normal direction
    if (Comp < 0 )
    {
      Dist = (p.y()*cosSPhi - p.x()*sinSPhi) ;

      if (Dist < halfCarTolerance)
      {
        sphi = Dist/Comp ;
        if (sphi < snxt)
        {
          if ( sphi < 0 )  { sphi = 0 ; }

          xi    = p.x() + sphi*v.x() ;
          yi    = p.y() + sphi*v.y() ;
          zi    = p.z() + sphi*v.z() ;
          rhoi2 = xi*xi + yi*yi ;
          it2   = std::fabs(rhoi2 + zi*zi + Rtor2 - 2*fRtor*std::sqrt(rhoi2)) ;

          if ( it2 >= tolORMin2 && it2 <= tolORMax2 )
          {
            // r intersection is good - check intersecting
            // with correct half-plane
            //
            if ((yi*cosCPhi-xi*sinCPhi)<=0)  { snxt=sphi; }
          }
        }
      }
    }
    ePhi=fSPhi+fDPhi;    // Second phi surface ('E'nding phi)
    sinEPhi=std::sin(ePhi);
    cosEPhi=std::cos(ePhi);
    Comp=-(v.x()*sinEPhi-v.y()*cosEPhi);
        
    if ( Comp < 0 )   // Component in outwards normal dirn
    {
      Dist = -(p.y()*cosEPhi - p.x()*sinEPhi) ;

      if (Dist < halfCarTolerance )
      {
        sphi = Dist/Comp ;

        if (sphi < snxt )
        {
          if (sphi < 0 )  { sphi = 0 ; }
       
          xi    = p.x() + sphi*v.x() ;
          yi    = p.y() + sphi*v.y() ;
          zi    = p.z() + sphi*v.z() ;
          rhoi2 = xi*xi + yi*yi ;
          it2   = std::fabs(rhoi2 + zi*zi + Rtor2 - 2*fRtor*std::sqrt(rhoi2)) ;

          if (it2 >= tolORMin2 && it2 <= tolORMax2)
          {
            // z and r intersections good - check intersecting
            // with correct half-plane
            //
            if ((yi*cosCPhi-xi*sinCPhi)>=0)  { snxt=sphi; }
          }    
        }
      }
    }
  }
  if(snxt < halfCarTolerance)  { snxt = 0.0 ; }

  return snxt ;
}

G4double G4Torus::DistanceToOut

(

const G4ThreeVector &

p

)

const [virtual]

Implements G4VSolid.

Definition at line 1542 of file G4Torus.cc.

References G4VSolid::DumpInfo(), G4cout, G4Exception(), Inside(), JustWarning, and kOutside.

{
  G4double safe=0.0,safeR1,safeR2;
  G4double rho2,rho,pt2,pt ;
  G4double safePhi,phiC,cosPhiC,sinPhiC,ePhi;
  rho2 = p.x()*p.x() + p.y()*p.y() ;
  rho  = std::sqrt(rho2) ;
  pt2  = std::fabs(rho2 + p.z()*p.z() + fRtor*fRtor - 2*fRtor*rho) ;
  pt   = std::sqrt(pt2) ;

#ifdef G4CSGDEBUG
  if( Inside(p) == kOutside )
  {
     G4int oldprc = G4cout.precision(16) ;
     G4cout << G4endl ;
     DumpInfo();
     G4cout << "Position:"  << G4endl << G4endl ;
     G4cout << "p.x() = "   << p.x()/mm << " mm" << G4endl ;
     G4cout << "p.y() = "   << p.y()/mm << " mm" << G4endl ;
     G4cout << "p.z() = "   << p.z()/mm << " mm" << G4endl << G4endl ;
     G4cout.precision(oldprc);
     G4Exception("G4Torus::DistanceToOut(p)", "GeomSolids1002",
                 JustWarning, "Point p is outside !?" );
  }
#endif

  if (fRmin)
  {
    safeR1 = pt - fRmin ;
    safeR2 = fRmax - pt ;

    if (safeR1 < safeR2)  { safe = safeR1 ; }
    else                  { safe = safeR2 ; }
  }
  else
  {
    safe = fRmax - pt ;
  }  

  // Check if phi divided, Calc distances closest phi plane
  //
  if (fDPhi<twopi) // Above/below central phi of Torus?
  {
    phiC    = fSPhi + fDPhi*0.5 ;
    cosPhiC = std::cos(phiC) ;
    sinPhiC = std::sin(phiC) ;

    if ((p.y()*cosPhiC-p.x()*sinPhiC)<=0)
    {
      safePhi = -(p.x()*std::sin(fSPhi) - p.y()*std::cos(fSPhi)) ;
    }
    else
    {
      ePhi    = fSPhi + fDPhi ;
      safePhi = (p.x()*std::sin(ePhi) - p.y()*std::cos(ePhi)) ;
    }
    if (safePhi < safe)  { safe = safePhi ; }
  }
  if (safe < 0)  { safe = 0 ; }
  return safe ;  
}

G4double G4Torus::DistanceToOut	(	const G4ThreeVector &	p,
		const G4ThreeVector &	v,
		const G4bool	calcNorm = G4bool(false),
		G4bool *	validNorm = 0,
		G4ThreeVector *	n = 0
	)			const [virtual]

Implements G4VSolid.

Definition at line 1187 of file G4Torus.cc.

References G4VSolid::DumpInfo(), G4cout, G4Exception(), JustWarning, G4VSolid::kCarTolerance, and G4INCL::Math::pi.

Referenced by SurfaceNormal().

{
  ESide    side = kNull, sidephi = kNull ;
  G4double snxt = kInfinity, sphi, sd[4] ;

  static const G4double halfCarTolerance = 0.5*kCarTolerance;
  static const G4double halfAngTolerance = 0.5*kAngTolerance;

  // Vars for phi intersection
  //
  G4double sinSPhi, cosSPhi, ePhi, sinEPhi, cosEPhi;
  G4double cPhi, sinCPhi, cosCPhi ;
  G4double pDistS, compS, pDistE, compE, sphi2, xi, yi, zi, vphi ;

  // Radial Intersections Defenitions & General Precals


#if 1

  // This is the version with the calculation of CalcNorm = true 
  // To be done: Check the precision of this calculation.
  // If you want return always validNorm = false, then take the version below
  
  G4double rho2  = p.x()*p.x()+p.y()*p.y();
  G4double rho   = std::sqrt(rho2) ;

  G4double pt2   = rho2 + p.z()*p.z() + fRtor * (fRtor - 2.0*rho);
    // Regroup for slightly better FP accuracy

  if( pt2 < 0.0)
  {
     pt2= std::fabs( pt2 );
  }
     
  G4double pt    = std::sqrt(pt2) ;

  G4double pDotV = p.x()*v.x() + p.y()*v.y() + p.z()*v.z() ;

  G4double tolRMax = fRmax - fRmaxTolerance ;
   
  G4double vDotNmax   = pDotV - fRtor*(v.x()*p.x() + v.y()*p.y())/rho ;
  G4double pDotxyNmax = (1 - fRtor/rho) ;

  if( (pt2 > tolRMax*tolRMax) && (vDotNmax >= 0) )
  {
    // On tolerant boundary & heading outwards (or perpendicular to) outer
    // radial surface -> leaving immediately with *n for really convex part
    // only
      
    if ( calcNorm && (pDotxyNmax >= -2.*fRmaxTolerance) ) 
    {
      *n = G4ThreeVector( p.x()*(1 - fRtor/rho)/pt,
                          p.y()*(1 - fRtor/rho)/pt,
                          p.z()/pt                  ) ;
      *validNorm = true ;
    }
     
    return snxt = 0 ; // Leaving by Rmax immediately
  }
  
  snxt = SolveNumericJT(p,v,fRmax,false);  
  side = kRMax ;

  // rmin

  if ( fRmin )
  {
    G4double tolRMin = fRmin + fRminTolerance ;

    if ( (pt2 < tolRMin*tolRMin) && (vDotNmax < 0) )
    {
      if (calcNorm)  { *validNorm = false ; } // Concave surface of the torus
      return  snxt = 0 ;                      // Leaving by Rmin immediately
    }
    
    sd[0] = SolveNumericJT(p,v,fRmin,false);
    if ( sd[0] < snxt )
    {
      snxt = sd[0] ;
      side = kRMin ;
    }
  }

#else

  // this is the "conservative" version which return always validnorm = false
  // NOTE: using this version the unit test testG4Torus will break

  snxt = SolveNumericJT(p,v,fRmax,false);  
  side = kRMax ;

  if ( fRmin )
  {
    sd[0] = SolveNumericJT(p,v,fRmin,false);
    if ( sd[0] < snxt )
    {
      snxt = sd[0] ;
      side = kRMin ;
    }
  }

  if ( calcNorm && (snxt == 0.0) )
  {
    *validNorm = false ;    // Leaving solid, but possible re-intersection
    return snxt  ;
  }

#endif
  
  if (fDPhi < twopi)  // Phi Intersections
  {
    sinSPhi = std::sin(fSPhi) ;
    cosSPhi = std::cos(fSPhi) ;
    ePhi    = fSPhi + fDPhi ;
    sinEPhi = std::sin(ePhi) ;
    cosEPhi = std::cos(ePhi) ;
    cPhi    = fSPhi + fDPhi*0.5 ;
    sinCPhi = std::sin(cPhi) ;
    cosCPhi = std::cos(cPhi) ;
    
    // angle calculation with correction 
    // of difference in domain of atan2 and Sphi
    //
    vphi = std::atan2(v.y(),v.x()) ;
     
    if ( vphi < fSPhi - halfAngTolerance  )    { vphi += twopi; }
    else if ( vphi > ePhi + halfAngTolerance ) { vphi -= twopi; }

    if ( p.x() || p.y() ) // Check if on z axis (rho not needed later)
    {
      pDistS = p.x()*sinSPhi - p.y()*cosSPhi ; // pDist -ve when inside
      pDistE = -p.x()*sinEPhi + p.y()*cosEPhi ;

      // Comp -ve when in direction of outwards normal
      //
      compS   = -sinSPhi*v.x() + cosSPhi*v.y() ;
      compE   = sinEPhi*v.x() - cosEPhi*v.y() ;
      sidephi = kNull ;
     
      if( ( (fDPhi <= pi) && ( (pDistS <= halfCarTolerance)
                            && (pDistE <= halfCarTolerance) ) )
       || ( (fDPhi >  pi) && !((pDistS >  halfCarTolerance)
                            && (pDistE >  halfCarTolerance) ) )  )
      {
        // Inside both phi *full* planes

        if ( compS < 0 )
        {
          sphi = pDistS/compS ;
            
          if (sphi >= -halfCarTolerance)
          {
            xi = p.x() + sphi*v.x() ;
            yi = p.y() + sphi*v.y() ;
              
            // Check intersecting with correct half-plane
            // (if not -> no intersect)
            //
            if ( (std::fabs(xi)<=kCarTolerance)
              && (std::fabs(yi)<=kCarTolerance) )
            {
              sidephi = kSPhi;
              if ( ((fSPhi-halfAngTolerance)<=vphi)
                && ((ePhi+halfAngTolerance)>=vphi) )
              {
                sphi = kInfinity;
              }
            }
            else if ( yi*cosCPhi-xi*sinCPhi >=0 )
            {
              sphi = kInfinity ;
            }
            else
            {
              sidephi = kSPhi ;
            }       
          }
          else
          {
            sphi = kInfinity ;
          }
        }
        else
        {
          sphi = kInfinity ;
        }

        if ( compE < 0 )
        {
          sphi2 = pDistE/compE ;
            
          // Only check further if < starting phi intersection
          //
          if ( (sphi2 > -kCarTolerance) && (sphi2 < sphi) )
          {
            xi = p.x() + sphi2*v.x() ;
            yi = p.y() + sphi2*v.y() ;
              
            if ( (std::fabs(xi)<=kCarTolerance)
              && (std::fabs(yi)<=kCarTolerance) )
            {
              // Leaving via ending phi
              //
              if( !( (fSPhi-halfAngTolerance <= vphi)
                  && (ePhi+halfAngTolerance >= vphi) ) )
              {
                sidephi = kEPhi ;
                sphi = sphi2;
              }
            } 
            else    // Check intersecting with correct half-plane 
            {
              if ( (yi*cosCPhi-xi*sinCPhi) >= 0)
              {
                // Leaving via ending phi
                //
                sidephi = kEPhi ;
                sphi = sphi2;
               
              }
            }
          }
        }
      }
      else
      {
        sphi = kInfinity ;
      }
    } 
    else
    {
      // On z axis + travel not || to z axis -> if phi of vector direction
      // within phi of shape, Step limited by rmax, else Step =0

      vphi = std::atan2(v.y(),v.x());
 
      if ( ( fSPhi-halfAngTolerance <= vphi ) && 
           ( vphi <= ( ePhi+halfAngTolerance ) ) )
      {
        sphi = kInfinity;
      }
      else
      {
        sidephi = kSPhi ; // arbitrary 
        sphi=0;
      }
    }

    // Order intersections

    if (sphi<snxt)
    {
      snxt=sphi;
      side=sidephi;
    }     
  }

  G4double rhoi2,rhoi,it2,it,iDotxyNmax ;
  // Note: by numerical computation we know where the ray hits the torus
  // So I propose to return the side where the ray hits

  if (calcNorm)
  {
    switch(side)
    {
      case kRMax:                     // n is unit vector 
        xi    = p.x() + snxt*v.x() ;
        yi    =p.y() + snxt*v.y() ;
        zi    = p.z() + snxt*v.z() ;
        rhoi2 = xi*xi + yi*yi ;
        rhoi  = std::sqrt(rhoi2) ;
        it2   = std::fabs(rhoi2 + zi*zi + fRtor*fRtor - 2*fRtor*rhoi) ;
        it    = std::sqrt(it2) ;
        iDotxyNmax = (1-fRtor/rhoi) ;
        if(iDotxyNmax >= -2.*fRmaxTolerance) // really convex part of Rmax
        {                       
          *n = G4ThreeVector( xi*(1-fRtor/rhoi)/it,
                              yi*(1-fRtor/rhoi)/it,
                              zi/it                 ) ;
          *validNorm = true ;
        }
        else
        {
          *validNorm = false ; // concave-convex part of Rmax
        }
        break ;

      case kRMin:
        *validNorm = false ;  // Rmin is concave or concave-convex
        break;

      case kSPhi:
        if (fDPhi <= pi )
        {
          *n=G4ThreeVector(std::sin(fSPhi),-std::cos(fSPhi),0);
          *validNorm=true;
        }
        else
        {
          *validNorm = false ;
        }
        break ;

      case kEPhi:
        if (fDPhi <= pi)
        {
          *n=G4ThreeVector(-std::sin(fSPhi+fDPhi),std::cos(fSPhi+fDPhi),0);
          *validNorm=true;
        }
        else
        {
          *validNorm = false ;
        }
        break;

      default:

        // It seems we go here from time to time ...

        G4cout << G4endl;
        DumpInfo();
        std::ostringstream message;
        G4int oldprc = message.precision(16);
        message << "Undefined side for valid surface normal to solid."
                << G4endl
                << "Position:"  << G4endl << G4endl
                << "p.x() = "   << p.x()/mm << " mm" << G4endl
                << "p.y() = "   << p.y()/mm << " mm" << G4endl
                << "p.z() = "   << p.z()/mm << " mm" << G4endl << G4endl
                << "Direction:" << G4endl << G4endl
                << "v.x() = "   << v.x() << G4endl
                << "v.y() = "   << v.y() << G4endl
                << "v.z() = "   << v.z() << G4endl << G4endl
                << "Proposed distance :" << G4endl << G4endl
                << "snxt = " << snxt/mm << " mm" << G4endl;
        message.precision(oldprc);
        G4Exception("G4Torus::DistanceToOut(p,v,..)",
                    "GeomSolids1002",JustWarning, message);
        break;
    }
  }
  if ( snxt<halfCarTolerance )  { snxt=0 ; }

  return snxt;
}

G4double G4Torus::GetCubicVolume

(

)

[inline, virtual]

Reimplemented from G4VSolid.

Definition at line 68 of file G4Torus.icc.

References G4CSGSolid::fCubicVolume, and G4INCL::Math::pi.

{
  if(fCubicVolume != 0.) {;}
  else  { fCubicVolume = fDPhi*CLHEP::pi*fRtor*(fRmax*fRmax-fRmin*fRmin); }
  return fCubicVolume;
}

G4double G4Torus::GetDPhi

(

)

const [inline]

Definition at line 62 of file G4Torus.icc.

Referenced by G4tgbGeometryDumper::GetSolidParams(), G4GDMLWriteParamvol::Torus_dimensionsWrite(), and G4GDMLWriteSolids::TorusWrite().

{
  return fDPhi ;
}

G4GeometryType G4Torus::GetEntityType

(

)

const [virtual]

Implements G4VSolid.

Definition at line 1694 of file G4Torus.cc.

{
  return G4String("G4Torus");
}

G4ThreeVector G4Torus::GetPointOnSurface

(

)

const [virtual]

Reimplemented from G4VSolid.

Definition at line 1735 of file G4Torus.cc.

References G4CSGSolid::GetRadiusInRing(), and G4INCL::Math::pi.

{
  G4double cosu, sinu,cosv, sinv, aOut, aIn, aSide, chose, phi, theta, rRand;
   
  phi   = RandFlat::shoot(fSPhi,fSPhi+fDPhi);
  theta = RandFlat::shoot(0.,twopi);
  
  cosu   = std::cos(phi);    sinu = std::sin(phi);
  cosv   = std::cos(theta);  sinv = std::sin(theta); 

  // compute the areas

  aOut   = (fDPhi)*twopi*fRtor*fRmax;
  aIn    = (fDPhi)*twopi*fRtor*fRmin;
  aSide  = pi*(fRmax*fRmax-fRmin*fRmin);
  
  if ((fSPhi == 0) && (fDPhi == twopi)){ aSide = 0; }
  chose = RandFlat::shoot(0.,aOut + aIn + 2.*aSide);

  if(chose < aOut)
  {
    return G4ThreeVector ((fRtor+fRmax*cosv)*cosu,
                          (fRtor+fRmax*cosv)*sinu, fRmax*sinv);
  }
  else if( (chose >= aOut) && (chose < aOut + aIn) )
  {
    return G4ThreeVector ((fRtor+fRmin*cosv)*cosu,
                          (fRtor+fRmin*cosv)*sinu, fRmin*sinv);
  }
  else if( (chose >= aOut + aIn) && (chose < aOut + aIn + aSide) )
  {
    rRand = GetRadiusInRing(fRmin,fRmax);
    return G4ThreeVector ((fRtor+rRand*cosv)*std::cos(fSPhi),
                          (fRtor+rRand*cosv)*std::sin(fSPhi), rRand*sinv);
  }
  else
  {   
    rRand = GetRadiusInRing(fRmin,fRmax);
    return G4ThreeVector ((fRtor+rRand*cosv)*std::cos(fSPhi+fDPhi),
                          (fRtor+rRand*cosv)*std::sin(fSPhi+fDPhi), 
                          rRand*sinv);
   }
}

G4double G4Torus::GetRmax

(

)

const [inline]

Definition at line 44 of file G4Torus.icc.

Referenced by G4tgbGeometryDumper::GetSolidParams(), G4GDMLWriteParamvol::Torus_dimensionsWrite(), and G4GDMLWriteSolids::TorusWrite().

{
  return fRmax ;
}

G4double G4Torus::GetRmin

(

)

const [inline]

Definition at line 38 of file G4Torus.icc.

Referenced by G4tgbGeometryDumper::GetSolidParams(), G4GDMLWriteParamvol::Torus_dimensionsWrite(), and G4GDMLWriteSolids::TorusWrite().

{
  return fRmin ;
}

G4double G4Torus::GetRtor

(

)

const [inline]

Definition at line 50 of file G4Torus.icc.

Referenced by G4tgbGeometryDumper::GetSolidParams(), G4GDMLWriteParamvol::Torus_dimensionsWrite(), and G4GDMLWriteSolids::TorusWrite().

{
  return fRtor ;
}

G4double G4Torus::GetSPhi

(

)

const [inline]

Definition at line 56 of file G4Torus.icc.

Referenced by G4tgbGeometryDumper::GetSolidParams(), G4GDMLWriteParamvol::Torus_dimensionsWrite(), and G4GDMLWriteSolids::TorusWrite().

{
  return fSPhi ;
}

G4double G4Torus::GetSurfaceArea

(

)

[inline, virtual]

Reimplemented from G4VSolid.

Definition at line 76 of file G4Torus.icc.

References G4CSGSolid::fSurfaceArea.

{
  if(fSurfaceArea != 0.) {;}
  else
  { 
    fSurfaceArea = fDPhi*CLHEP::twopi*fRtor*(fRmax+fRmin);
    if(fDPhi < CLHEP::twopi)
    {
      fSurfaceArea = fSurfaceArea + CLHEP::twopi*(fRmax*fRmax-fRmin*fRmin);
    } 
  }
  return fSurfaceArea;
}

EInside G4Torus::Inside

(

const G4ThreeVector &

p

)

const [virtual]

Implements G4VSolid.

Definition at line 632 of file G4Torus.cc.

References kInside, kOutside, and kSurface.

Referenced by CalculateExtent(), DistanceToOut(), and SurfaceNormal().

{
  G4double r2, pt2, pPhi, tolRMin, tolRMax ;

  EInside in = kOutside ;
  static const G4double halfAngTolerance = 0.5*kAngTolerance;

                                               // General precals
  r2  = p.x()*p.x() + p.y()*p.y() ;
  pt2 = r2 + p.z()*p.z() + fRtor*fRtor - 2*fRtor*std::sqrt(r2) ;

  if (fRmin) tolRMin = fRmin + fRminTolerance ;
  else       tolRMin = 0 ;

  tolRMax = fRmax - fRmaxTolerance;
      
  if (pt2 >= tolRMin*tolRMin && pt2 <= tolRMax*tolRMax )
  {
    if ( fDPhi == twopi || pt2 == 0 )  // on torus swept axis
    {
      in = kInside ;
    }
    else
    {
      // Try inner tolerant phi boundaries (=>inside)
      // if not inside, try outer tolerant phi boundaries

      pPhi = std::atan2(p.y(),p.x()) ;

      if ( pPhi < -halfAngTolerance )  { pPhi += twopi ; }  // 0<=pPhi<2pi
      if ( fSPhi >= 0 )
      {
        if ( (std::fabs(pPhi) < halfAngTolerance)
            && (std::fabs(fSPhi + fDPhi - twopi) < halfAngTolerance) )
        { 
            pPhi += twopi ; // 0 <= pPhi < 2pi
        }
        if ( (pPhi >= fSPhi + halfAngTolerance)
            && (pPhi <= fSPhi + fDPhi - halfAngTolerance) )
        {
          in = kInside ;
        }
          else if ( (pPhi >= fSPhi - halfAngTolerance)
                 && (pPhi <= fSPhi + fDPhi + halfAngTolerance) )
        {
          in = kSurface ;
        }
      }
      else  // fSPhi < 0
      {
          if ( (pPhi <= fSPhi + twopi - halfAngTolerance)
            && (pPhi >= fSPhi + fDPhi  + halfAngTolerance) )  {;}
          else
          {
            in = kSurface ;
          }
      }
    }
  }
  else   // Try generous boundaries
  {
    tolRMin = fRmin - fRminTolerance ;
    tolRMax = fRmax + fRmaxTolerance ;

    if (tolRMin < 0 )  { tolRMin = 0 ; }

    if ( (pt2 >= tolRMin*tolRMin) && (pt2 <= tolRMax*tolRMax) )
    {
      if ( (fDPhi == twopi) || (pt2 == 0) ) // Continuous in phi or on z-axis
      {
        in = kSurface ;
      }
      else // Try outer tolerant phi boundaries only
      {
        pPhi = std::atan2(p.y(),p.x()) ;

        if ( pPhi < -halfAngTolerance )  { pPhi += twopi ; }  // 0<=pPhi<2pi
        if ( fSPhi >= 0 )
        {
          if ( (std::fabs(pPhi) < halfAngTolerance)
            && (std::fabs(fSPhi + fDPhi - twopi) < halfAngTolerance) )
          { 
            pPhi += twopi ; // 0 <= pPhi < 2pi
          }
          if ( (pPhi >= fSPhi - halfAngTolerance)
            && (pPhi <= fSPhi + fDPhi + halfAngTolerance) )
          {
            in = kSurface;
          }
        }
        else  // fSPhi < 0
        {
          if ( (pPhi <= fSPhi + twopi - halfAngTolerance)
            && (pPhi >= fSPhi + fDPhi  + halfAngTolerance) )  {;}
          else
          {
            in = kSurface ;
          }
        }
      }
    }
  }
  return in ;
}

G4Torus & G4Torus::operator=

(

const G4Torus &

rhs

)

Definition at line 212 of file G4Torus.cc.

References fDPhi, fRmax, fRmaxTolerance, fRmin, fRminTolerance, fRtor, fSPhi, kAngTolerance, kRadTolerance, and G4CSGSolid::operator=().

{
   // Check assignment to self
   //
   if (this == &rhs)  { return *this; }

   // Copy base class data
   //
   G4CSGSolid::operator=(rhs);

   // Copy data
   //
   fRmin = rhs.fRmin; fRmax = rhs.fRmax;
   fRtor = rhs.fRtor; fSPhi = rhs.fSPhi; fDPhi = rhs.fDPhi;
   fRminTolerance = rhs.fRminTolerance; fRmaxTolerance = rhs.fRmaxTolerance;
   kRadTolerance = rhs.kRadTolerance; kAngTolerance = rhs.kAngTolerance;

   return *this;
}

void G4Torus::SetAllParameters	(	G4double	pRmin,
		G4double	pRmax,
		G4double	pRtor,
		G4double	pSPhi,
		G4double	pDPhi
	)

Definition at line 100 of file G4Torus.cc.

References FatalException, G4CSGSolid::fCubicVolume, G4CSGSolid::fpPolyhedron, G4CSGSolid::fSurfaceArea, G4endl, G4Exception(), G4GeometryTolerance::GetAngularTolerance(), G4GeometryTolerance::GetInstance(), G4VSolid::GetName(), G4GeometryTolerance::GetRadialTolerance(), and G4VSolid::kCarTolerance.

Referenced by G4Torus().

{
  const G4double fEpsilon = 4.e-11;  // relative tolerance of radii

  fCubicVolume = 0.;
  fSurfaceArea = 0.;
  fpPolyhedron = 0;

  kRadTolerance = G4GeometryTolerance::GetInstance()->GetRadialTolerance();
  kAngTolerance = G4GeometryTolerance::GetInstance()->GetAngularTolerance();
  
  if ( pRtor >= pRmax+1.e3*kCarTolerance )  // Check swept radius, as in G4Cons
  {
    fRtor = pRtor ;
  }
  else
  {
    std::ostringstream message;
    message << "Invalid swept radius for Solid: " << GetName() << G4endl
            << "        pRtor = " << pRtor << ", pRmax = " << pRmax;
    G4Exception("G4Torus::SetAllParameters()",
                "GeomSolids0002", FatalException, message);
  }

  // Check radii, as in G4Cons
  //
  if ( pRmin < pRmax - 1.e2*kCarTolerance && pRmin >= 0 )
  {
    if (pRmin >= 1.e2*kCarTolerance) { fRmin = pRmin ; }
    else                             { fRmin = 0.0   ; }
    fRmax = pRmax ;
  }
  else
  {
    std::ostringstream message;
    message << "Invalid values of radii for Solid: " << GetName() << G4endl
            << "        pRmin = " << pRmin << ", pRmax = " << pRmax;
    G4Exception("G4Torus::SetAllParameters()",
                "GeomSolids0002", FatalException, message);
  }

  // Relative tolerances
  //
  fRminTolerance = (fRmin)
                 ? 0.5*std::max( kRadTolerance, fEpsilon*(fRtor-fRmin )) : 0;
  fRmaxTolerance = 0.5*std::max( kRadTolerance, fEpsilon*(fRtor+fRmax) );

  // Check angles
  //
  if ( pDPhi >= twopi )  { fDPhi = twopi ; }
  else
  {
    if (pDPhi > 0)       { fDPhi = pDPhi ; }
    else
    {
      std::ostringstream message;
      message << "Invalid Z delta-Phi for Solid: " << GetName() << G4endl
              << "        pDPhi = " << pDPhi;
      G4Exception("G4Torus::SetAllParameters()",
                  "GeomSolids0002", FatalException, message);
    }
  }
  
  // Ensure psphi in 0-2PI or -2PI-0 range if shape crosses 0
  //
  fSPhi = pSPhi;

  if (fSPhi < 0)  { fSPhi = twopi-std::fmod(std::fabs(fSPhi),twopi) ; }
  else            { fSPhi = std::fmod(fSPhi,twopi) ; }

  if (fSPhi+fDPhi > twopi)  { fSPhi-=twopi ; }
}

std::ostream & G4Torus::StreamInfo

(

std::ostream &

os

)

const [virtual]

Reimplemented from G4CSGSolid.

Definition at line 1712 of file G4Torus.cc.

References G4VSolid::GetName().

{
  G4int oldprc = os.precision(16);
  os << "-----------------------------------------------------------\n"
     << "    *** Dump for solid - " << GetName() << " ***\n"
     << "    ===================================================\n"
     << " Solid type: G4Torus\n"
     << " Parameters: \n"
     << "    inner radius: " << fRmin/mm << " mm \n"
     << "    outer radius: " << fRmax/mm << " mm \n"
     << "    swept radius: " << fRtor/mm << " mm \n"
     << "    starting phi: " << fSPhi/degree << " degrees \n"
     << "    delta phi   : " << fDPhi/degree << " degrees \n"
     << "-----------------------------------------------------------\n";
  os.precision(oldprc);

  return os;
}

G4ThreeVector G4Torus::SurfaceNormal

(

const G4ThreeVector &

p

)

const [virtual]

Implements G4VSolid.

Definition at line 743 of file G4Torus.cc.

References DistanceToIn(), DistanceToOut(), G4endl, G4Exception(), Inside(), JustWarning, G4VSolid::kCarTolerance, kInside, kOutside, and kSurface.

{
  G4int noSurfaces = 0;  
  G4double rho2, rho, pt2, pt, pPhi;
  G4double distRMin = kInfinity;
  G4double distSPhi = kInfinity, distEPhi = kInfinity;

  // To cope with precision loss
  //
  const G4double delta = std::max(10.0*kCarTolerance,
                                  1.0e-8*(fRtor+fRmax));
  const G4double dAngle = 10.0*kAngTolerance;

  G4ThreeVector nR, nPs, nPe;
  G4ThreeVector norm, sumnorm(0.,0.,0.);

  rho2 = p.x()*p.x() + p.y()*p.y();
  rho = std::sqrt(rho2);
  pt2 = rho2+p.z()*p.z() +fRtor * (fRtor-2*rho);
  pt2 = std::max(pt2, 0.0); // std::fabs(pt2);
  pt = std::sqrt(pt2) ;

  G4double  distRMax = std::fabs(pt - fRmax);
  if(fRmin) distRMin = std::fabs(pt - fRmin);

  if( rho > delta && pt != 0.0 )
  {
    G4double redFactor= (rho-fRtor)/rho;
    nR = G4ThreeVector( p.x()*redFactor,  // p.x()*(1.-fRtor/rho),
                        p.y()*redFactor,  // p.y()*(1.-fRtor/rho),
                        p.z()          );
    nR *= 1.0/pt;
  }

  if ( fDPhi < twopi ) // && rho ) // old limitation against (0,0,z)
  {
    if ( rho )
    {
      pPhi = std::atan2(p.y(),p.x());

      if(pPhi < fSPhi-delta)            { pPhi += twopi; }
      else if(pPhi > fSPhi+fDPhi+delta) { pPhi -= twopi; }

      distSPhi = std::fabs( pPhi - fSPhi );
      distEPhi = std::fabs(pPhi-fSPhi-fDPhi);
    }
    nPs = G4ThreeVector(std::sin(fSPhi),-std::cos(fSPhi),0);
    nPe = G4ThreeVector(-std::sin(fSPhi+fDPhi),std::cos(fSPhi+fDPhi),0);
  } 
  if( distRMax <= delta )
  {
    noSurfaces ++;
    sumnorm += nR;
  }
  else if( fRmin && (distRMin <= delta) ) // Must not be on both Outer and Inner
  {
    noSurfaces ++;
    sumnorm -= nR;
  }

  //  To be on one of the 'phi' surfaces,
  //  it must be within the 'tube' - with tolerance

  if( (fDPhi < twopi) && (fRmin-delta <= pt) && (pt <= (fRmax+delta)) )
  {
    if (distSPhi <= dAngle)
    {
      noSurfaces ++;
      sumnorm += nPs;
    }
    if (distEPhi <= dAngle) 
    {
      noSurfaces ++;
      sumnorm += nPe;
    }
  }
  if ( noSurfaces == 0 )
  {
#ifdef G4CSGDEBUG
     G4ExceptionDescription ed;
     ed.precision(16);

     EInside  inIt= Inside( p );
     
     if( inIt != kSurface )
     {
        ed << " ERROR>  Surface Normal was called for Torus,"
           << " with point not on surface." << G4endl;
     }
     else
     {
        ed << " ERROR>  Surface Normal has not found a surface, "
           << " despite the point being on the surface. " <<G4endl;
     }

     if( inIt != kInside)
     {
         ed << " Safety (Dist To In)  = " << DistanceToIn(p) << G4endl;
     }
     if( inIt != kOutside)
     {
         ed << " Safety (Dist to Out) = " << DistanceToOut(p) << G4endl;
     }
     ed << " Coordinates of point : " << p << G4endl;
     ed << " Parameters  of solid : " << G4endl << *this << G4endl;

     if( inIt == kSurface )
     {
        G4Exception("G4Torus::SurfaceNormal(p)", "GeomSolids1002",
                    JustWarning, ed,
                    "Failing to find normal, even though point is on surface!" );
     }
     else
     {
        static const char* NameInside[3]= { "Inside", "Surface", "Outside" };
        ed << "  The point is " << NameInside[inIt] << " the solid. "<< G4endl;
        G4Exception("G4Torus::SurfaceNormal(p)", "GeomSolids1002",
                    JustWarning, ed, "Point p is not on surface !?" );
     }
#endif
     norm = ApproxSurfaceNormal(p);
  }
  else if ( noSurfaces == 1 )  { norm = sumnorm; }
  else                         { norm = sumnorm.unit(); }

  // G4cout << "G4Torus::SurfaceNormal p= " << p << " returns norm= " << norm << G4endl;

  return norm ;
}

---

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

• G4Torus.hh
• G4Torus.icc
• G4Torus.cc

---