#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

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 ()

virtual G4Polyhedron *	GetPolyhedron () const

	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)

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

G4double	EstimateCubicVolume (G4int nStat, G4double epsilon) const

G4double	EstimateSurfaceArea (G4int nStat, G4double ell) const

Protected Attributes inherited from G4CSGSolid
G4double	fCubicVolume

G4double	fSurfaceArea

G4Polyhedron *	fpPolyhedron

Protected Attributes inherited from G4VSolid
G4double	kCarTolerance

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 80 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.

194 {}

G4Torus::G4Torus ( __void__ & a )

Definition at line 181 of file G4Torus.cc.

   : G4CSGSolid(a), fRmin(0.), fRmax(0.), fRtor(0.), fSPhi(0.),
     fDPhi(0.), fRminTolerance(0.), fRmaxTolerance(0. ),
     kRadTolerance(0.), kAngTolerance(0.),
     halfCarTolerance(0.), halfAngTolerance(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),
     halfCarTolerance(rhs.halfCarTolerance),
     halfAngTolerance(rhs.halfAngTolerance)
 {
 }

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 414 of file G4Torus.cc.

 {
   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 1700 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 241 of file G4Torus.cc.

References G4VPVParameterisation::ComputeDimensions().

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

G4Polyhedron * G4Torus::CreatePolyhedron ( ) const

virtual

Reimplemented from G4VSolid.

Definition at line 1785 of file G4Torus.cc.

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

void G4Torus::DescribeYourselfTo ( G4VGraphicsScene & scene ) const

virtual

Implements G4VSolid.

Definition at line 1780 of file G4Torus.cc.

References G4VGraphicsScene::AddSolid().

 {
   scene.AddSolid (*this);
 }

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

virtual

Implements G4VSolid.

Definition at line 989 of file G4Torus.cc.

References python.hepunit::twopi, CLHEP::Hep3Vector::x(), CLHEP::Hep3Vector::y(), and CLHEP::Hep3Vector::z().

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;
 
   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::DistanceToIn ( const G4ThreeVector & p ) const

virtual

Implements G4VSolid.

Definition at line 1139 of file G4Torus.cc.

References python.hepunit::twopi, CLHEP::Hep3Vector::x(), CLHEP::Hep3Vector::y(), and CLHEP::Hep3Vector::z().

 {
   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::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, G4endl, G4Exception(), JustWarning, G4VSolid::kCarTolerance, python.hepunit::mm, python.hepunit::pi, python.hepunit::twopi, CLHEP::Hep3Vector::x(), CLHEP::Hep3Vector::y(), and CLHEP::Hep3Vector::z().

Referenced by SurfaceNormal().

 {
   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 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::DistanceToOut ( const G4ThreeVector & p ) const

virtual

Implements G4VSolid.

Definition at line 1539 of file G4Torus.cc.

References G4VSolid::DumpInfo(), G4cout, G4endl, G4Exception(), Inside(), JustWarning, kOutside, python.hepunit::mm, python.hepunit::twopi, CLHEP::Hep3Vector::x(), CLHEP::Hep3Vector::y(), and CLHEP::Hep3Vector::z().

 {
   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::GetCubicVolume ( )

inlinevirtual

Reimplemented from G4VSolid.

G4double G4Torus::GetDPhi ( ) const

inline

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

G4GeometryType G4Torus::GetEntityType ( ) const

virtual

Implements G4VSolid.

Definition at line 1691 of file G4Torus.cc.

 {
   return G4String("G4Torus");
 }

G4ThreeVector G4Torus::GetPointOnSurface ( ) const

virtual

Reimplemented from G4VSolid.

Definition at line 1732 of file G4Torus.cc.

References G4CSGSolid::GetRadiusInRing(), python.hepunit::pi, G4INCL::DeJongSpin::shoot(), and python.hepunit::twopi.

 {
   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

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

G4double G4Torus::GetRmin ( ) const

inline

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

G4double G4Torus::GetRtor ( ) const

inline

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

G4double G4Torus::GetSPhi ( ) const

inline

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

G4double G4Torus::GetSurfaceArea ( )

inlinevirtual

Reimplemented from G4VSolid.

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

virtual

Implements G4VSolid.

Definition at line 634 of file G4Torus.cc.

References kInside, kOutside, kSurface, python.hepunit::twopi, CLHEP::Hep3Vector::x(), CLHEP::Hep3Vector::y(), and CLHEP::Hep3Vector::z().

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

 {
   G4double r2, pt2, pPhi, tolRMin, tolRMax ;
 
   EInside in = kOutside ;
 
   // 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 214 of file G4Torus.cc.

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

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

Definition at line 96 of file G4Torus.cc.

References FatalException, G4CSGSolid::fCubicVolume, G4CSGSolid::fpPolyhedron, G4CSGSolid::fSurfaceArea, G4endl, G4Exception(), G4GeometryTolerance::GetAngularTolerance(), G4GeometryTolerance::GetInstance(), G4VSolid::GetName(), G4GeometryTolerance::GetRadialTolerance(), G4VSolid::kCarTolerance, G4INCL::Math::max(), and python.hepunit::twopi.

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();
 
   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.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 1709 of file G4Torus.cc.

References python.hepunit::degree, G4VSolid::GetName(), and python.hepunit::mm.

 {
   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 745 of file G4Torus.cc.

References DistanceToIn(), DistanceToOut(), G4endl, G4Exception(), Inside(), JustWarning, G4VSolid::kCarTolerance, kInside, kOutside, kSurface, G4INCL::Math::max(), python.hepunit::twopi, CLHEP::Hep3Vector::unit(), CLHEP::Hep3Vector::x(), CLHEP::Hep3Vector::y(), and CLHEP::Hep3Vector::z().

 {
   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:

Public Member Functions

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation