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 () override
G4double	GetRmin () const
G4double	GetRmax () const
G4double	GetRtor () const
G4double	GetSPhi () const
G4double	GetDPhi () const
G4double	GetSinStartPhi () const
G4double	GetCosStartPhi () const
G4double	GetSinEndPhi () const
G4double	GetCosEndPhi () const
G4double	GetCubicVolume () override
G4double	GetSurfaceArea () override
EInside	Inside (const G4ThreeVector &p) const override
void	BoundingLimits (G4ThreeVector &pMin, G4ThreeVector &pMax) const override
G4bool	CalculateExtent (const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pmin, G4double &pmax) const override
void	ComputeDimensions (G4VPVParameterisation *p, const G4int n, const G4VPhysicalVolume *pRep) override
G4ThreeVector	SurfaceNormal (const G4ThreeVector &p) const override
G4double	DistanceToIn (const G4ThreeVector &p, const G4ThreeVector &v) const override
G4double	DistanceToIn (const G4ThreeVector &p) const override
G4double	DistanceToOut (const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=false, G4bool *validNorm=nullptr, G4ThreeVector *n=nullptr) const override
G4double	DistanceToOut (const G4ThreeVector &p) const override
G4GeometryType	GetEntityType () const override
G4ThreeVector	GetPointOnSurface () const override
G4VSolid *	Clone () const override
std::ostream &	StreamInfo (std::ostream &os) const override
void	DescribeYourselfTo (G4VGraphicsScene &scene) const override
G4Polyhedron *	CreatePolyhedron () const override
void	SetAllParameters (G4double pRmin, G4double pRmax, G4double pRtor, G4double pSPhi, G4double pDPhi)
	G4Torus (__void__ &)
	G4Torus (const G4Torus &rhs)
G4Torus &	operator= (const G4Torus &rhs)
Public Member Functions inherited from G4CSGSolid
	G4CSGSolid (const G4String &pName)
virtual	~G4CSGSolid ()
G4Polyhedron *	GetPolyhedron () const override
	G4CSGSolid (__void__ &)
	G4CSGSolid (const G4CSGSolid &rhs)
G4CSGSolid &	operator= (const G4CSGSolid &rhs)
Public Member Functions inherited from G4VSolid
	G4VSolid (const G4String &name)
virtual	~G4VSolid ()
G4bool	operator== (const G4VSolid &s) const
G4String	GetName () const
void	SetName (const G4String &name)
G4double	GetTolerance () const
void	DumpInfo () const
virtual G4VisExtent	GetExtent () const
virtual const G4VSolid *	GetConstituentSolid (G4int no) const
virtual G4VSolid *	GetConstituentSolid (G4int no)
virtual const G4DisplacedSolid *	GetDisplacedSolidPtr () const
virtual G4DisplacedSolid *	GetDisplacedSolidPtr ()
	G4VSolid (__void__ &)
	G4VSolid (const G4VSolid &rhs)
G4VSolid &	operator= (const G4VSolid &rhs)
G4double	EstimateCubicVolume (G4int nStat, G4double epsilon) const
G4double	EstimateSurfaceArea (G4int nStat, G4double ell) const

Additional Inherited Members
Protected Member Functions inherited from G4CSGSolid
G4double	GetRadiusInRing (G4double rmin, G4double rmax) const
Protected Member Functions inherited from G4VSolid
void	CalculateClippedPolygonExtent (G4ThreeVectorList &pPolygon, const G4VoxelLimits &pVoxelLimit, const EAxis pAxis, G4double &pMin, G4double &pMax) const
void	ClipCrossSection (G4ThreeVectorList *pVertices, const G4int pSectionIndex, const G4VoxelLimits &pVoxelLimit, const EAxis pAxis, G4double &pMin, G4double &pMax) const
void	ClipBetweenSections (G4ThreeVectorList *pVertices, const G4int pSectionIndex, const G4VoxelLimits &pVoxelLimit, const EAxis pAxis, G4double &pMin, G4double &pMax) const
void	ClipPolygon (G4ThreeVectorList &pPolygon, const G4VoxelLimits &pVoxelLimit, const EAxis pAxis) const
Protected Attributes inherited from G4CSGSolid
G4double	fCubicVolume = 0.0
G4double	fSurfaceArea = 0.0
G4bool	fRebuildPolyhedron = false
G4Polyhedron *	fpPolyhedron = nullptr
Protected Attributes inherited from G4VSolid
G4double	kCarTolerance

Detailed Description

Definition at line 91 of file G4Torus.hh.

Constructor & Destructor Documentation

G4Torus() [1/3]

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

Definition at line 65 of file G4Torus.cc.

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

Referenced by Clone().

~G4Torus()

G4Torus::~G4Torus ( )

overridedefault

G4Torus() [2/3]

G4Torus::G4Torus

(

__void__ &

a

)

Definition at line 166 of file G4Torus.cc.

  : G4CSGSolid(a)
{
}

G4Torus() [3/3]

G4Torus::G4Torus ( const G4Torus & rhs )

default

Member Function Documentation

BoundingLimits()

void G4Torus::BoundingLimits ( G4ThreeVector & pMin, G4ThreeVector & pMax ) const

overridevirtual

Reimplemented from G4VSolid.

Definition at line 383 of file G4Torus.cc.

{
  G4double rmax = GetRmax();
  G4double rtor = GetRtor();
  G4double rint = rtor - rmax;
  G4double rext = rtor + rmax;
  G4double dz   = rmax;
 
  // Find bounding box
  //
  if (GetDPhi() >= twopi)
  {
    pMin.set(-rext,-rext,-dz);
    pMax.set( rext, rext, dz);
  }
  else
  {
    G4TwoVector vmin,vmax;
    G4GeomTools::DiskExtent(rint,rext,
                            GetSinStartPhi(),GetCosStartPhi(),
                            GetSinEndPhi(),GetCosEndPhi(),
                            vmin,vmax);
    pMin.set(vmin.x(),vmin.y(),-dz);
    pMax.set(vmax.x(),vmax.y(), dz);
  }
 
  // Check correctness of the bounding box
  //
  if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
  {
    std::ostringstream message;
    message << "Bad bounding box (min >= max) for solid: "
            << GetName() << " !"
            << "\npMin = " << pMin
            << "\npMax = " << pMax;
    G4Exception("G4Torus::BoundingLimits()", "GeomMgt0001",
                JustWarning, message);
    DumpInfo();
  }
}

Referenced by CalculateExtent().

CalculateExtent()

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

overridevirtual

Implements G4VSolid.

Definition at line 428 of file G4Torus.cc.

{
  G4ThreeVector bmin, bmax;
  G4bool exist;
 
  // Get bounding box
  BoundingLimits(bmin,bmax);
 
  // Check bounding box
  G4BoundingEnvelope bbox(bmin,bmax);
#ifdef G4BBOX_EXTENT
  return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
#endif
  if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVoxelLimit,pTransform,pMin,pMax))
  {
    return exist = pMin < pMax;
  }
 
  // Get parameters of the solid
  G4double rmin = GetRmin();
  G4double rmax = GetRmax();
  G4double rtor = GetRtor();
  G4double dphi = GetDPhi();
  G4double sinStart = GetSinStartPhi();
  G4double cosStart = GetCosStartPhi();
  G4double sinEnd   = GetSinEndPhi();
  G4double cosEnd   = GetCosEndPhi();
  G4double rint = rtor - rmax;
  G4double rext = rtor + rmax;
 
  // Find bounding envelope and calculate extent
  //
  static const G4int NPHI  = 24; // number of steps for whole torus
  static const G4int NDISK = 16; // number of steps for disk
  static const G4double sinHalfDisk = std::sin(pi/NDISK);
  static const G4double cosHalfDisk = std::cos(pi/NDISK);
  static const G4double sinStepDisk = 2.*sinHalfDisk*cosHalfDisk;
  static const G4double cosStepDisk = 1. - 2.*sinHalfDisk*sinHalfDisk;
 
  G4double astep = (360/NPHI)*deg; // max angle for one slice in phi
  G4int    kphi  = (dphi <= astep) ? 1 : (G4int)((dphi-deg)/astep) + 1;
  G4double ang   = dphi/kphi;
 
  G4double sinHalf = std::sin(0.5*ang);
  G4double cosHalf = std::cos(0.5*ang);
  G4double sinStep = 2.*sinHalf*cosHalf;
  G4double cosStep = 1. - 2.*sinHalf*sinHalf;
 
  // define vectors for bounding envelope
  G4ThreeVectorList pols[NDISK+1];
  for (auto & pol : pols) pol.resize(4);
 
  std::vector<const G4ThreeVectorList *> polygons;
  polygons.resize(NDISK+1);
  for (G4int k=0; k<NDISK+1; ++k) polygons[k] = &pols[k];
 
  // set internal and external reference circles
  G4TwoVector rzmin[NDISK];
  G4TwoVector rzmax[NDISK];
 
  if ((rtor-rmin*sinHalfDisk)/cosHalf > (rtor+rmin*sinHalfDisk)) rmin = 0;
  rmax /= cosHalfDisk;
  G4double sinCurDisk = sinHalfDisk;
  G4double cosCurDisk = cosHalfDisk;
  for (G4int k=0; k<NDISK; ++k)
  {
    G4double rmincur = rtor + rmin*cosCurDisk;
    if (cosCurDisk < 0 && rmin > 0) rmincur /= cosHalf;
    rzmin[k].set(rmincur,rmin*sinCurDisk);
 
    G4double rmaxcur = rtor + rmax*cosCurDisk;
    if (cosCurDisk > 0) rmaxcur /= cosHalf;
    rzmax[k].set(rmaxcur,rmax*sinCurDisk);
 
    G4double sinTmpDisk = sinCurDisk;
    sinCurDisk = sinCurDisk*cosStepDisk + cosCurDisk*sinStepDisk;
    cosCurDisk = cosCurDisk*cosStepDisk - sinTmpDisk*sinStepDisk;
  }
 
  // Loop along slices in Phi. The extent is calculated as cumulative
  // extent of the slices
  pMin =  kInfinity;
  pMax = -kInfinity;
  G4double eminlim = pVoxelLimit.GetMinExtent(pAxis);
  G4double emaxlim = pVoxelLimit.GetMaxExtent(pAxis);
  G4double sinCur1 = 0, cosCur1 = 0, sinCur2 = 0, cosCur2 = 0;
  for (G4int i=0; i<kphi+1; ++i)
  {
    if (i == 0)
    {
      sinCur1 = sinStart;
      cosCur1 = cosStart;
      sinCur2 = sinCur1*cosHalf + cosCur1*sinHalf;
      cosCur2 = cosCur1*cosHalf - sinCur1*sinHalf;
    }
    else
    {
      sinCur1 = sinCur2;
      cosCur1 = cosCur2;
      sinCur2 = (i == kphi) ? sinEnd : sinCur1*cosStep + cosCur1*sinStep;
      cosCur2 = (i == kphi) ? cosEnd : cosCur1*cosStep - sinCur1*sinStep;
    }
    for (G4int k=0; k<NDISK; ++k)
    {
      G4double r1 = rzmin[k].x(), r2 = rzmax[k].x();
      G4double z1 = rzmin[k].y(), z2 = rzmax[k].y();
      pols[k][0].set(r1*cosCur1,r1*sinCur1,z1);
      pols[k][1].set(r2*cosCur1,r2*sinCur1,z2);
      pols[k][2].set(r2*cosCur2,r2*sinCur2,z2);
      pols[k][3].set(r1*cosCur2,r1*sinCur2,z1);
    }
    pols[NDISK] = pols[0];
 
    // get bounding box of current slice
    G4TwoVector vmin,vmax;
    G4GeomTools::
      DiskExtent(rint,rext,sinCur1,cosCur1,sinCur2,cosCur2,vmin,vmax);
    bmin.setX(vmin.x()); bmin.setY(vmin.y());
    bmax.setX(vmax.x()); bmax.setY(vmax.y());
 
    // set bounding envelope for current slice and adjust extent
    G4double emin,emax;
    G4BoundingEnvelope benv(bmin,bmax,polygons);
    if (!benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,emin,emax)) continue;
    if (emin < pMin) pMin = emin;
    if (emax > pMax) pMax = emax;
    if (eminlim > pMin && emaxlim < pMax) break; // max possible extent
  }
  return (pMin < pMax);
}

Clone()

G4VSolid * G4Torus::Clone ( ) const

overridevirtual

Reimplemented from G4VSolid.

Definition at line 1553 of file G4Torus.cc.

{
  return new G4Torus(*this);
}

ComputeDimensions()

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

overridevirtual

Reimplemented from G4VSolid.

Definition at line 214 of file G4Torus.cc.

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

CreatePolyhedron()

G4Polyhedron * G4Torus::CreatePolyhedron ( ) const

overridevirtual

Reimplemented from G4VSolid.

Definition at line 1635 of file G4Torus.cc.

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

DescribeYourselfTo()

void G4Torus::DescribeYourselfTo ( G4VGraphicsScene & scene ) const

overridevirtual

Implements G4VSolid.

Definition at line 1630 of file G4Torus.cc.

{
  scene.AddSolid (*this);
}

DistanceToIn() [1/2]

G4double G4Torus::DistanceToIn ( const G4ThreeVector & p ) const

overridevirtual

Implements G4VSolid.

Definition at line 1092 of file G4Torus.cc.

{
  G4double safe=0.0, safe1, safe2 ;
  G4double phiC, cosPhiC, sinPhiC, safePhi, ePhi, cosPsi ;
  G4double rho, pt ;
  
  rho = std::hypot(p.x(),p.y());
  pt  = std::hypot(p.z(),rho-fRtor);
  safe1 = fRmin - pt ;
  safe2 = pt - fRmax ;
 
  if (safe1 > safe2)  { safe = safe1; }
  else                { safe = safe2; }
 
  if ( fDPhi < twopi && (rho != 0.0) )
  {
    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;
}

DistanceToIn() [2/2]

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

overridevirtual

Implements G4VSolid.

Definition at line 914 of file G4Torus.cc.

{
  // Get bounding box of full torus
  //
  G4double boxDx  = fRtor + fRmax;
  G4double boxDy  = boxDx;
  G4double boxDz  = fRmax;
  G4double boxMax = boxDx;
  G4double boxMin = boxDz;
 
  // Check if point is traveling away
  //
  G4double distX = std::abs(p.x()) - boxDx;
  G4double distY = std::abs(p.y()) - boxDy;
  G4double distZ = std::abs(p.z()) - boxDz;
  if (distX >= -halfCarTolerance && p.x()*v.x() >= 0) return kInfinity;
  if (distY >= -halfCarTolerance && p.y()*v.y() >= 0) return kInfinity;
  if (distZ >= -halfCarTolerance && p.z()*v.z() >= 0) return kInfinity;
 
  // Calculate safety distance to bounding box
  // If point is too far, move it closer and calculate distance
  //
  G4double Dmax = 32*boxMax; 
  G4double safe = std::max(std::max(distX,distY),distZ);
  if (safe > Dmax)
  {
    G4double dist = safe - 1.e-8*safe - boxMin; // stay outside after the move
    dist += DistanceToIn(p + dist*v, v);
    return (dist >= kInfinity) ? kInfinity : dist;
  }
 
  // Find intersection with torus
  //
  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;
 
  G4double Dist,xi,yi,zi,rhoi,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)
 
  snxt = SolveNumericJT(p,v,fRmax,true);
 
  if (fRmin != 0.0)  // 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() ;
          rhoi = std::hypot(xi,yi);
          it2 = zi*zi + (rhoi-fRtor)*(rhoi-fRtor);
 
          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() ;
          rhoi = std::hypot(xi,yi);
          it2 = zi*zi + (rhoi-fRtor)*(rhoi-fRtor);
 
          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 ;
}

Referenced by DistanceToIn(), and SurfaceNormal().

DistanceToOut() [1/2]

G4double G4Torus::DistanceToOut ( const G4ThreeVector & p ) const

overridevirtual

Implements G4VSolid.

Definition at line 1479 of file G4Torus.cc.

{
  G4double safe=0.0,safeR1,safeR2;
  G4double rho,pt ;
  G4double safePhi,phiC,cosPhiC,sinPhiC,ePhi;
  
  rho = std::hypot(p.x(),p.y());
  pt  = std::hypot(p.z(),rho-fRtor);
  
#ifdef G4CSGDEBUG
  if( Inside(p) == kOutside )
  {
     G4long 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 != 0.0)
  {
    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 ;  
}

DistanceToOut() [2/2]

G4double G4Torus::DistanceToOut ( const G4ThreeVector & p, const G4ThreeVector & v, const G4bool calcNorm = false, G4bool * validNorm = nullptr, G4ThreeVector * n = nullptr ) const

overridevirtual

Implements G4VSolid.

Definition at line 1137 of file G4Torus.cc.

{
  ESide    side = kNull, sidephi = kNull ;
  G4double snxt = kInfinity, sphi, sd[4] ;
 
  // 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
 
  //////////////////////// new calculation //////////////////////
 
#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 rho = std::hypot(p.x(),p.y());
  G4double pt = hypot(p.z(),rho-fRtor);
 
  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( (pt*pt > 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 != 0.0 )
  {
    G4double tolRMin = fRmin + fRminTolerance ;
 
    if ( (pt*pt < 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() != 0.0) || (p.y() != 0.0) ) // 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 rhoi,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() ;
        rhoi = std::hypot(xi,yi);
        it = hypot(zi,rhoi-fRtor);
 
        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;
        G4long 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;
}

Referenced by SurfaceNormal().

GetCosEndPhi()

G4double G4Torus::GetCosEndPhi ( ) const

inline

Referenced by BoundingLimits(), and CalculateExtent().

GetCosStartPhi()

G4double G4Torus::GetCosStartPhi ( ) const

inline

Referenced by BoundingLimits(), and CalculateExtent().

GetCubicVolume()

G4double G4Torus::GetCubicVolume ( )

inlineoverridevirtual

Reimplemented from G4VSolid.

GetDPhi()

G4double G4Torus::GetDPhi ( ) const

inline

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

GetEntityType()

G4GeometryType G4Torus::GetEntityType ( ) const

overridevirtual

Implements G4VSolid.

Definition at line 1544 of file G4Torus.cc.

{
  return {"G4Torus"};
}

GetPointOnSurface()

G4ThreeVector G4Torus::GetPointOnSurface ( ) const

overridevirtual

Reimplemented from G4VSolid.

Definition at line 1585 of file G4Torus.cc.

{
  G4double cosu, sinu,cosv, sinv, aOut, aIn, aSide, chose, phi, theta, rRand;
   
  phi   = G4RandFlat::shoot(fSPhi,fSPhi+fDPhi);
  theta = G4RandFlat::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 = G4RandFlat::shoot(0.,aOut + aIn + 2.*aSide);
 
  if(chose < aOut)
  {
    return { (fRtor+fRmax*cosv)*cosu, (fRtor+fRmax*cosv)*sinu, fRmax*sinv };
  }
  else if( (chose >= aOut) && (chose < aOut + aIn) )
  {
    return { (fRtor+fRmin*cosv)*cosu, (fRtor+fRmin*cosv)*sinu, fRmin*sinv };
  }
  else if( (chose >= aOut + aIn) && (chose < aOut + aIn + aSide) )
  {
    rRand = GetRadiusInRing(fRmin,fRmax);
    return { (fRtor+rRand*cosv)*std::cos(fSPhi),
             (fRtor+rRand*cosv)*std::sin(fSPhi), rRand*sinv };
  }
  else
  {   
    rRand = GetRadiusInRing(fRmin,fRmax);
    return { (fRtor+rRand*cosv)*std::cos(fSPhi+fDPhi),
             (fRtor+rRand*cosv)*std::sin(fSPhi+fDPhi), rRand*sinv };
   }
}

GetRmax()

G4double G4Torus::GetRmax ( ) const

inline

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

GetRmin()

G4double G4Torus::GetRmin ( ) const

inline

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

GetRtor()

G4double G4Torus::GetRtor ( ) const

inline

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

GetSinEndPhi()

G4double G4Torus::GetSinEndPhi ( ) const

inline

Referenced by BoundingLimits(), and CalculateExtent().

GetSinStartPhi()

G4double G4Torus::GetSinStartPhi ( ) const

inline

Referenced by BoundingLimits(), and CalculateExtent().

GetSPhi()

G4double G4Torus::GetSPhi ( ) const

inline

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

GetSurfaceArea()

G4double G4Torus::GetSurfaceArea ( )

inlineoverridevirtual

Reimplemented from G4VSolid.

Inside()

EInside G4Torus::Inside ( const G4ThreeVector & p ) const

overridevirtual

Implements G4VSolid.

Definition at line 566 of file G4Torus.cc.

{
  G4double r, pt2, pPhi, tolRMin, tolRMax ;
 
  EInside in = kOutside ;
 
  // General precals
  //
  r   = std::hypot(p.x(),p.y());
  pt2 = p.z()*p.z() + (r-fRtor)*(r-fRtor);
 
  if (fRmin != 0.0) 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 ;
}

Referenced by DistanceToOut(), and SurfaceNormal().

operator=()

G4Torus & G4Torus::operator=

(

const G4Torus &

rhs

)

Definition at line 187 of file G4Torus.cc.

{
   // 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;
   halfCarTolerance = rhs.halfCarTolerance;
   halfAngTolerance = rhs.halfAngTolerance;
 
   return *this;
}

SetAllParameters()

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

Definition at line 81 of file G4Torus.cc.

{
  const G4double fEpsilon = 4.e-11;  // relative tolerance of radii
 
  fCubicVolume = 0.;
  fSurfaceArea = 0.;
  fRebuildPolyhedron = true;
 
  kRadTolerance = G4GeometryTolerance::GetInstance()->GetRadialTolerance();
  kAngTolerance = G4GeometryTolerance::GetInstance()->GetAngularTolerance();
 
  halfCarTolerance = 0.5*kCarTolerance;
  halfAngTolerance = 0.5*kAngTolerance;
 
  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.0
                 ? 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 ; }
}

Referenced by G4Torus().

StreamInfo()

std::ostream & G4Torus::StreamInfo ( std::ostream & os ) const

overridevirtual

Reimplemented from G4CSGSolid.

Definition at line 1562 of file G4Torus.cc.

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

SurfaceNormal()

G4ThreeVector G4Torus::SurfaceNormal ( const G4ThreeVector & p ) const

overridevirtual

Implements G4VSolid.

Definition at line 677 of file G4Torus.cc.

{
  G4int noSurfaces = 0;  
  G4double rho, 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.);
 
  rho = std::hypot(p.x(),p.y());
  pt  = std::hypot(p.z(),rho-fRtor);
 
  G4double  distRMax = std::fabs(pt - fRmax);
  if(fRmin != 0.0) 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 != 0.0 )
    {
      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 != 0.0) && (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(); }
 
  return norm ;
}

---

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

• G4Torus.hh
• G4Torus.cc