G4Torus Class Reference

#include <G4Torus.hh>

Inheritance diagram for G4Torus:

Collaboration 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)
G4double	EstimateCubicVolume (G4int nStat, G4double epsilon) const
G4double	EstimateSurfaceArea (G4int nStat, G4double ell) const

Private Types
enum	ESide {   kNull, kRMin, kRMax, kSPhi,   kEPhi }
enum	ENorm { kNRMin, kNRMax, kNSPhi, kNEPhi }

Private Member Functions
void	TorusRootsJT (const G4ThreeVector &p, const G4ThreeVector &v, G4double r, std::vector< G4double > &roots) const
G4double	SolveNumericJT (const G4ThreeVector &p, const G4ThreeVector &v, G4double r, G4bool IsDistanceToIn) const
G4ThreeVectorList *	CreateRotatedVertices (const G4AffineTransform &pTransform, G4int &noPolygonVertices) const
G4ThreeVector	ApproxSurfaceNormal (const G4ThreeVector &p) const

Private Attributes
G4double	fRmin
G4double	fRmax
G4double	fRtor
G4double	fSPhi
G4double	fDPhi
G4double	fRminTolerance
G4double	fRmaxTolerance
G4double	kRadTolerance
G4double	kAngTolerance
G4double	halfCarTolerance
G4double	halfAngTolerance

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
G4double	fSurfaceArea
G4bool	fRebuildPolyhedron
G4Polyhedron *	fpPolyhedron
Protected Attributes inherited from G4VSolid
G4double	kCarTolerance

Detailed Description

Definition at line 111 of file G4Torus.hh.

Member Enumeration Documentation

ENorm

enum G4Torus::ENorm

private

Enumerator
kNRMin
kNRMax
kNSPhi
kNEPhi

Definition at line 209 of file G4Torus.hh.

{kNRMin,kNRMax,kNSPhi,kNEPhi};

ESide

enum G4Torus::ESide

private

Enumerator
kNull
kRMin
kRMax
kSPhi
kEPhi

Definition at line 206 of file G4Torus.hh.

{kNull,kRMin,kRMax,kSPhi,kEPhi};

Constructor & Destructor Documentation

G4Torus() [1/3]

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

Definition at line 82 of file G4Torus.cc.

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

Here is the call graph for this function:

Here is the caller graph for this function:

~G4Torus()

G4Torus::~G4Torus

(

)

Definition at line 195 of file G4Torus.cc.

{}

G4Torus() [2/3]

G4Torus::G4Torus

(

__void__ &

a

)

Definition at line 183 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() [3/3]

G4Torus::G4Torus

(

const G4Torus &

rhs

)

Definition at line 202 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

ApproxSurfaceNormal()

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

private

Definition at line 882 of file G4Torus.cc.

{
  ENorm side ;
  G4ThreeVector norm;
  G4double rho2,rho,pt2,pt,phi;
  G4double distRMin,distRMax,distSPhi,distEPhi,distMin;

  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
  G4cout << " G4Torus::ApproximateSurfaceNormal called for point " << p
         << G4endl;
#endif
   
  distRMax = std::fabs(pt - fRmax) ;

  if(fRmin)  // First minimum radius
  {
    distRMin = std::fabs(pt - fRmin) ;

    if (distRMin < distRMax)
    {
      distMin = distRMin ;
      side    = kNRMin ;
    }
    else
    {
      distMin = distRMax ;
      side    = kNRMax ;
    }
  }
  else
  {
    distMin = distRMax ;
    side    = kNRMax ;
  }    
  if ( (fDPhi < twopi) && rho )
  {
    phi = std::atan2(p.y(),p.x()) ; // Protected against (0,0,z) (above rho!=0)

    if (phi < 0)  { phi += twopi ; }

    if (fSPhi < 0 )  { distSPhi = std::fabs(phi-(fSPhi+twopi))*rho ; }
    else             { distSPhi = std::fabs(phi-fSPhi)*rho ; }

    distEPhi = std::fabs(phi - fSPhi - fDPhi)*rho ;

    if (distSPhi < distEPhi) // Find new minimum
    {
      if (distSPhi<distMin) side = kNSPhi ;
    }
    else
    {
      if (distEPhi < distMin)  { side = kNEPhi ; }
    }
  }  
  switch (side)
  {
    case kNRMin:      // Inner radius
      norm = G4ThreeVector( -p.x()*(1-fRtor/rho)/pt,
                            -p.y()*(1-fRtor/rho)/pt,
                            -p.z()/pt                 ) ;
      break ;
    case kNRMax:      // Outer radius
      norm = G4ThreeVector( p.x()*(1-fRtor/rho)/pt,
                            p.y()*(1-fRtor/rho)/pt,
                            p.z()/pt                  ) ;
      break;
    case kNSPhi:
      norm = G4ThreeVector(std::sin(fSPhi),-std::cos(fSPhi),0) ;
      break;
    case kNEPhi:
      norm = G4ThreeVector(-std::sin(fSPhi+fDPhi),std::cos(fSPhi+fDPhi),0) ;
      break;
    default:          // Should never reach this case ...
      DumpInfo();
      G4Exception("G4Torus::ApproxSurfaceNormal()",
                  "GeomSolids1002", JustWarning,
                  "Undefined side for valid surface normal to solid.");
      break ;
  } 
  return norm ;
}

Here is the call graph for this function:

Here is the caller graph for this function:

CalculateExtent()

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

virtual

Implements G4VSolid.

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

Here is the call graph for this function:

Clone()

G4VSolid * G4Torus::Clone ( ) const

virtual

Reimplemented from G4VSolid.

Definition at line 1702 of file G4Torus.cc.

{
  return new G4Torus(*this);
}

Here is the call graph for this function:

ComputeDimensions()

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

virtual

Reimplemented from G4VSolid.

Definition at line 243 of file G4Torus.cc.

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

Here is the call graph for this function:

CreatePolyhedron()

G4Polyhedron * G4Torus::CreatePolyhedron ( ) const

virtual

Reimplemented from G4VSolid.

Definition at line 1787 of file G4Torus.cc.

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

CreateRotatedVertices()

G4ThreeVectorList * G4Torus::CreateRotatedVertices ( const G4AffineTransform &  pTransform, G4int &  noPolygonVertices  ) const

private

Definition at line 1615 of file G4Torus.cc.

{
  G4ThreeVectorList *vertices;
  G4ThreeVector vertex0,vertex1,vertex2,vertex3;
  G4double meshAngle,meshRMax,crossAngle,cosCrossAngle,sinCrossAngle,sAngle;
  G4double rMaxX,rMaxY,rMinX,rMinY;
  G4int crossSection,noCrossSections;

  // Compute no of cross-sections necessary to mesh tube
  //
  noCrossSections = G4int (fDPhi/kMeshAngleDefault) + 1 ;

  if (noCrossSections < kMinMeshSections)
  {
    noCrossSections = kMinMeshSections ;
  }
  else if (noCrossSections>kMaxMeshSections)
  {
    noCrossSections=kMaxMeshSections;
  }
  meshAngle = fDPhi/(noCrossSections - 1) ;
  meshRMax  = (fRtor + fRmax)/std::cos(meshAngle*0.5) ;

  // If complete in phi, set start angle such that mesh will be at fRmax
  // on the x axis. Will give better extent calculations when not rotated

  if ( (fDPhi == twopi) && (fSPhi == 0) )
  {
    sAngle = -meshAngle*0.5 ;
  }
  else
  {
    sAngle = fSPhi ;
  }
  vertices = new G4ThreeVectorList();
  
  if (vertices)
  {
    vertices->reserve(noCrossSections*4) ;
    for (crossSection=0;crossSection<noCrossSections;crossSection++)
    {
      // Compute coordinates of cross section at section crossSection

      crossAngle=sAngle+crossSection*meshAngle;
      cosCrossAngle=std::cos(crossAngle);
      sinCrossAngle=std::sin(crossAngle);

      rMaxX=meshRMax*cosCrossAngle;
      rMaxY=meshRMax*sinCrossAngle;
      rMinX=(fRtor-fRmax)*cosCrossAngle;
      rMinY=(fRtor-fRmax)*sinCrossAngle;
      vertex0=G4ThreeVector(rMinX,rMinY,-fRmax);
      vertex1=G4ThreeVector(rMaxX,rMaxY,-fRmax);
      vertex2=G4ThreeVector(rMaxX,rMaxY,+fRmax);
      vertex3=G4ThreeVector(rMinX,rMinY,+fRmax);

      vertices->push_back(pTransform.TransformPoint(vertex0));
      vertices->push_back(pTransform.TransformPoint(vertex1));
      vertices->push_back(pTransform.TransformPoint(vertex2));
      vertices->push_back(pTransform.TransformPoint(vertex3));
    }
    noPolygonVertices = 4 ;
  }
  else
  {
    DumpInfo();
    G4Exception("G4Torus::CreateRotatedVertices()",
                "GeomSolids0003", FatalException,
                "Error in allocation of vertices. Out of memory !");
  }
  return vertices;
}

Here is the call graph for this function:

Here is the caller graph for this function:

DescribeYourselfTo()

void G4Torus::DescribeYourselfTo ( G4VGraphicsScene &  scene ) const

virtual

Implements G4VSolid.

Definition at line 1782 of file G4Torus.cc.

{
  scene.AddSolid (*this);
}

Here is the call graph for this function:

DistanceToIn() [1/2]

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

virtual

Implements G4VSolid.

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

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

Here is the call graph for this function:

Here is the caller graph for this function:

DistanceToIn() [2/2]

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

virtual

Implements G4VSolid.

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

Here is the call graph for this function:

DistanceToOut() [1/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 1189 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


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

Here is the call graph for this function:

Here is the caller graph for this function:

DistanceToOut() [2/2]

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

virtual

Implements G4VSolid.

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

Here is the call graph for this function:

GetCubicVolume()

G4double G4Torus::GetCubicVolume ( )

inlinevirtual

Reimplemented from G4VSolid.

GetDPhi()

G4double G4Torus::GetDPhi ( ) const

inline

Here is the caller graph for this function:

GetEntityType()

G4GeometryType G4Torus::GetEntityType ( ) const

virtual

Implements G4VSolid.

Definition at line 1693 of file G4Torus.cc.

{
  return G4String("G4Torus");
}

GetPointOnSurface()

G4ThreeVector G4Torus::GetPointOnSurface ( ) const

virtual

Reimplemented from G4VSolid.

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

Here is the call graph for this function:

GetRmax()

G4double G4Torus::GetRmax ( ) const

inline

Here is the caller graph for this function:

GetRmin()

G4double G4Torus::GetRmin ( ) const

inline

Here is the caller graph for this function:

GetRtor()

G4double G4Torus::GetRtor ( ) const

inline

Here is the caller graph for this function:

GetSPhi()

G4double G4Torus::GetSPhi ( ) const

inline

Here is the caller graph for this function:

GetSurfaceArea()

G4double G4Torus::GetSurfaceArea ( )

inlinevirtual

Reimplemented from G4VSolid.

Inside()

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

virtual

Implements G4VSolid.

Definition at line 636 of file G4Torus.cc.

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

Here is the call graph for this function:

Here is the caller graph for this function:

operator=()

G4Torus & G4Torus::operator=

(

const G4Torus &

rhs

)

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

Here is the call graph for this function:

SetAllParameters()

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

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

Here is the call graph for this function:

Here is the caller graph for this function:

SolveNumericJT()

G4double G4Torus::SolveNumericJT ( const G4ThreeVector &  p, const G4ThreeVector &  v, G4double  r, G4bool  IsDistanceToIn  ) const

private

Definition at line 297 of file G4Torus.cc.

{
  G4double bigdist = 10*mm ;
  G4double tmin = kInfinity ;
  G4double t, scal ;

  // calculate the distances to the intersections with the Torus
  // from a given point p and direction v.
  //
  std::vector<G4double> roots ;
  std::vector<G4double> rootsrefined ;
  TorusRootsJT(p,v,r,roots) ;

  G4ThreeVector ptmp ;

  // determine the smallest non-negative solution
  //
  for ( size_t k = 0 ; k<roots.size() ; k++ )
  {
    t = roots[k] ;

    if ( t < -halfCarTolerance )  { continue ; }  // skip negative roots

    if ( t > bigdist && t<kInfinity )    // problem with big distances
    {
      ptmp = p + t*v ;
      TorusRootsJT(ptmp,v,r,rootsrefined) ;
      if ( rootsrefined.size()==roots.size() )
      {
        t = t + rootsrefined[k] ;
      }
    }

    ptmp = p + t*v ;   // calculate the position of the proposed intersection

    G4double theta = std::atan2(ptmp.y(),ptmp.x());
    
    if ( fSPhi >= 0 )
    {
      if ( theta < - halfAngTolerance )  { theta += twopi; }
      if ( (std::fabs(theta) < halfAngTolerance)
        && (std::fabs(fSPhi + fDPhi - twopi) < halfAngTolerance) )
      { 
        theta += twopi ; // 0 <= theta < 2pi
      }
    }
    if ((fSPhi <= -pi )&&(theta>halfAngTolerance)) { theta = theta-twopi; }
       
    // We have to verify if this root is inside the region between
    // fSPhi and fSPhi + fDPhi
    //
    if ( (theta - fSPhi >= - halfAngTolerance)
      && (theta - (fSPhi + fDPhi) <=  halfAngTolerance) )
    {
      // check if P is on the surface, and called from DistanceToIn
      // DistanceToIn has to return 0.0 if particle is going inside the solid

      if ( IsDistanceToIn == true )
      {
        if (std::fabs(t) < halfCarTolerance )
        {
          // compute scalar product at position p : v.n
          // ( n taken from SurfaceNormal, not normalized )

          scal = v* G4ThreeVector( p.x()*(1-fRtor/std::sqrt(p.x()*p.x()
                                          + p.y()*p.y())),
                                   p.y()*(1-fRtor/std::sqrt(p.x()*p.x()
                                          + p.y()*p.y())),
                                   p.z() );

          // change sign in case of inner radius
          //
          if ( r == GetRmin() )  { scal = -scal ; }
          if ( scal < 0 )  { return 0.0  ; }
        }
      }

      // check if P is on the surface, and called from DistanceToOut
      // DistanceToIn has to return 0.0 if particle is leaving the solid

      if ( IsDistanceToIn == false )
      {
        if (std::fabs(t) < halfCarTolerance )
        {
          // compute scalar product at position p : v.n   
          //
          scal = v* G4ThreeVector( p.x()*(1-fRtor/std::sqrt(p.x()*p.x()
                                          + p.y()*p.y())),
                                   p.y()*(1-fRtor/std::sqrt(p.x()*p.x()
                                          + p.y()*p.y())),
                                   p.z() );

          // change sign in case of inner radius
          //
          if ( r == GetRmin() )  { scal = -scal ; }
          if ( scal > 0 )  { return 0.0  ; }
        }
      }

      // check if distance is larger than 1/2 kCarTolerance
      //
      if(  t > halfCarTolerance  )
      {
        tmin = t  ;
        return tmin  ;
      }
    }
  }

  return tmin;
}

Here is the call graph for this function:

Here is the caller graph for this function:

StreamInfo()

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

virtual

Reimplemented from G4CSGSolid.

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

Here is the call graph for this function:

SurfaceNormal()

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

virtual

Implements G4VSolid.

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

Here is the call graph for this function:

TorusRootsJT()

void G4Torus::TorusRootsJT ( const G4ThreeVector &  p, const G4ThreeVector &  v, G4double  r, std::vector< G4double > &  roots  ) const

private

Definition at line 257 of file G4Torus.cc.

{

  G4int i, num ;
  G4double c[5], srd[4], si[4] ;

  G4double Rtor2 = fRtor*fRtor, r2 = r*r  ;

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

  c[0] = 1.0 ;
  c[1] = 4*pDotV ;
  c[2] = 2*(pRad2 + 2*pDotV*pDotV - Rtor2 - r2 + 2*Rtor2*v.z()*v.z()) ;
  c[3] = 4*(pDotV*(pRad2 - Rtor2 - r2) + 2*Rtor2*p.z()*v.z()) ;
  c[4] = pRad2*pRad2 - 2*pRad2*(Rtor2+r2) 
       + 4*Rtor2*p.z()*p.z() + (Rtor2-r2)*(Rtor2-r2) ;
  
  G4JTPolynomialSolver  torusEq;

  num = torusEq.FindRoots( c, 4, srd, si );
  
  for ( i = 0; i < num; i++ ) 
  {
    if( si[i] == 0. )  { roots.push_back(srd[i]) ; }  // store real roots
  }  

  std::sort(roots.begin() , roots.end() ) ;  // sorting  with <
}

Here is the call graph for this function:

Here is the caller graph for this function:

Member Data Documentation

fDPhi

G4double G4Torus::fDPhi

private

Definition at line 203 of file G4Torus.hh.

fRmax

G4double G4Torus::fRmax

private

Definition at line 203 of file G4Torus.hh.

fRmaxTolerance

G4double G4Torus::fRmaxTolerance

private

Definition at line 211 of file G4Torus.hh.

fRmin

G4double G4Torus::fRmin

private

Definition at line 203 of file G4Torus.hh.

fRminTolerance

G4double G4Torus::fRminTolerance

private

Definition at line 211 of file G4Torus.hh.

fRtor

G4double G4Torus::fRtor

private

Definition at line 203 of file G4Torus.hh.

fSPhi

G4double G4Torus::fSPhi

private

Definition at line 203 of file G4Torus.hh.

halfAngTolerance

G4double G4Torus::halfAngTolerance

private

Definition at line 214 of file G4Torus.hh.

halfCarTolerance

G4double G4Torus::halfCarTolerance

private

Definition at line 214 of file G4Torus.hh.

kAngTolerance

G4double G4Torus::kAngTolerance

private

Definition at line 211 of file G4Torus.hh.

kRadTolerance

G4double G4Torus::kRadTolerance

private

Definition at line 211 of file G4Torus.hh.

---

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

• G4Torus.hh
• G4Torus.cc