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)
Public Member Functions inherited from G4CSGSolid
	G4CSGSolid (const G4String &pName)
virtual	~G4CSGSolid ()
virtual std::ostream &	StreamInfo (std::ostream &os) const
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
virtual G4bool	CalculateExtent (const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pMin, G4double &pMax) const =0
virtual EInside	Inside (const G4ThreeVector &p) const =0
virtual G4ThreeVector	SurfaceNormal (const G4ThreeVector &p) const =0
virtual G4double	DistanceToIn (const G4ThreeVector &p, const G4ThreeVector &v) const =0
virtual G4double	DistanceToIn (const G4ThreeVector &p) const =0
virtual G4double	DistanceToOut (const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=false, G4bool *validNorm=0, G4ThreeVector *n=0) const =0
virtual G4double	DistanceToOut (const G4ThreeVector &p) const =0
virtual void	ComputeDimensions (G4VPVParameterisation *p, const G4int n, const G4VPhysicalVolume *pRep)
virtual G4double	GetCubicVolume ()
virtual G4double	GetSurfaceArea ()
virtual G4GeometryType	GetEntityType () const =0
virtual G4ThreeVector	GetPointOnSurface () const
virtual G4VSolid *	Clone () const
virtual std::ostream &	StreamInfo (std::ostream &os) const =0
void	DumpInfo () const
virtual void	DescribeYourselfTo (G4VGraphicsScene &scene) const =0
virtual G4VisExtent	GetExtent () const
virtual G4Polyhedron *	CreatePolyhedron () const
virtual G4NURBS *	CreateNURBS () const
virtual G4Polyhedron *	GetPolyhedron () 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() [1/3]

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

Definition at line 84 of file G4Torus.cc.

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

~G4Torus()

G4Torus::~G4Torus

(

)

Definition at line 193 of file G4Torus.cc.

{}

G4Torus() [2/3]

G4Torus::G4Torus

(

__void__ &

a

)

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() [3/3]

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

CalculateExtent()

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.

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

Clone()

G4VSolid * G4Torus::Clone ( ) const

virtual

Reimplemented from G4VSolid.

Definition at line 1703 of file G4Torus.cc.

{
  return new G4Torus(*this);
}

ComputeDimensions()

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

virtual

Reimplemented from G4VSolid.

Definition at line 237 of file G4Torus.cc.

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

CreateNURBS()

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

CreatePolyhedron()

G4Polyhedron * G4Torus::CreatePolyhedron ( ) const

virtual

Reimplemented from G4VSolid.

Definition at line 1788 of file G4Torus.cc.

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

DescribeYourselfTo()

void G4Torus::DescribeYourselfTo ( G4VGraphicsScene &  scene ) const

virtual

Implements G4VSolid.

Definition at line 1783 of file G4Torus.cc.

{
  scene.AddSolid (*this);
}

DistanceToIn() [1/2]

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

DistanceToIn() [2/2]

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

virtual

Implements G4VSolid.

Definition at line 987 of file G4Torus.cc.

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

Referenced by SurfaceNormal().

DistanceToOut() [1/2]

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

virtual

Implements G4VSolid.

Definition at line 1542 of file G4Torus.cc.

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

DistanceToOut() [2/2]

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.

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

Referenced by SurfaceNormal().

GetCubicVolume()

G4double G4Torus::GetCubicVolume ( )

inlinevirtual

Reimplemented from G4VSolid.

GetDPhi()

G4double G4Torus::GetDPhi ( ) const

inline

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

GetEntityType()

G4GeometryType G4Torus::GetEntityType ( ) const

virtual

Implements G4VSolid.

Definition at line 1694 of file G4Torus.cc.

{
  return G4String("G4Torus");
}

GetPointOnSurface()

G4ThreeVector G4Torus::GetPointOnSurface ( ) const

virtual

Reimplemented from G4VSolid.

Definition at line 1735 of file G4Torus.cc.

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

GetRmax()

G4double G4Torus::GetRmax ( ) const

inline

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

GetRmin()

G4double G4Torus::GetRmin ( ) const

inline

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

GetRtor()

G4double G4Torus::GetRtor ( ) const

inline

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

GetSPhi()

G4double G4Torus::GetSPhi ( ) const

inline

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

GetSurfaceArea()

G4double G4Torus::GetSurfaceArea ( )

inlinevirtual

Reimplemented from G4VSolid.

Inside()

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

virtual

Implements G4VSolid.

Definition at line 632 of file G4Torus.cc.

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

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

operator=()

G4Torus & G4Torus::operator=

(

const G4Torus &

rhs

)

Definition at line 212 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;
 
   return *this;
}

SetAllParameters()

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

Definition at line 100 of file G4Torus.cc.

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

Referenced by G4Torus().

StreamInfo()

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

virtual

Reimplemented from G4CSGSolid.

Definition at line 1712 of file G4Torus.cc.

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

SurfaceNormal()

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

virtual

Implements G4VSolid.

Definition at line 743 of file G4Torus.cc.

{
  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.cc