G4Sphere Class Reference

#include <G4Sphere.hh>

Inheritance diagram for G4Sphere:

Collaboration diagram for G4Sphere:

Public Member Functions
	G4Sphere (const G4String &pName, G4double pRmin, G4double pRmax, G4double pSPhi, G4double pDPhi, G4double pSTheta, G4double pDTheta)
	~G4Sphere ()
G4double	GetInnerRadius () const
G4double	GetOuterRadius () const
G4double	GetStartPhiAngle () const
G4double	GetDeltaPhiAngle () const
G4double	GetStartThetaAngle () const
G4double	GetDeltaThetaAngle () const
void	SetInnerRadius (G4double newRMin)
void	SetOuterRadius (G4double newRmax)
void	SetStartPhiAngle (G4double newSphi, G4bool trig=true)
void	SetDeltaPhiAngle (G4double newDphi)
void	SetStartThetaAngle (G4double newSTheta)
void	SetDeltaThetaAngle (G4double newDTheta)
G4double	GetCubicVolume ()
G4double	GetSurfaceArea ()
void	ComputeDimensions (G4VPVParameterisation *p, const G4int n, const G4VPhysicalVolume *pRep)
G4bool	CalculateExtent (const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pmin, G4double &pmax) const
EInside	Inside (const G4ThreeVector &p) const
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
G4VisExtent	GetExtent () const
void	DescribeYourselfTo (G4VGraphicsScene &scene) const
G4Polyhedron *	CreatePolyhedron () const
	G4Sphere (__void__ &)
	G4Sphere (const G4Sphere &rhs)
G4Sphere &	operator= (const G4Sphere &rhs)
G4double	GetRmin () const
G4double	GetRmax () const
G4double	GetSPhi () const
G4double	GetDPhi () const
G4double	GetSTheta () const
G4double	GetDTheta () const
G4double	GetInsideRadius () const
void	SetInsideRadius (G4double newRmin)
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 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, kSTheta, kETheta }
enum	ENorm {   kNRMin, kNRMax, kNSPhi, kNEPhi,   kNSTheta, kNETheta }

Private Member Functions
G4ThreeVectorList *	CreateRotatedVertices (const G4AffineTransform &pTransform, G4int &noPolygonVertices) const
void	Initialize ()
void	CheckThetaAngles (G4double sTheta, G4double dTheta)
void	CheckSPhiAngle (G4double sPhi)
void	CheckDPhiAngle (G4double dPhi)
void	CheckPhiAngles (G4double sPhi, G4double dPhi)
void	InitializePhiTrigonometry ()
void	InitializeThetaTrigonometry ()
G4ThreeVector	ApproxSurfaceNormal (const G4ThreeVector &p) const

Private Attributes
G4double	fRminTolerance
G4double	fRmaxTolerance
G4double	kAngTolerance
G4double	kRadTolerance
G4double	fEpsilon
G4double	fRmin
G4double	fRmax
G4double	fSPhi
G4double	fDPhi
G4double	fSTheta
G4double	fDTheta
G4double	sinCPhi
G4double	cosCPhi
G4double	cosHDPhiOT
G4double	cosHDPhiIT
G4double	sinSPhi
G4double	cosSPhi
G4double	sinEPhi
G4double	cosEPhi
G4double	hDPhi
G4double	cPhi
G4double	ePhi
G4double	sinSTheta
G4double	cosSTheta
G4double	sinETheta
G4double	cosETheta
G4double	tanSTheta
G4double	tanSTheta2
G4double	tanETheta
G4double	tanETheta2
G4double	eTheta
G4bool	fFullPhiSphere
G4bool	fFullThetaSphere
G4bool	fFullSphere
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 87 of file G4Sphere.hh.

Member Enumeration Documentation

ENorm

enum G4Sphere::ENorm

private

Enumerator
kNRMin
kNRMax
kNSPhi
kNEPhi
kNSTheta
kNETheta

Definition at line 227 of file G4Sphere.hh.

{kNRMin,kNRMax,kNSPhi,kNEPhi,kNSTheta,kNETheta};

ESide

enum G4Sphere::ESide

private

Enumerator
kNull
kRMin
kRMax
kSPhi
kEPhi
kSTheta
kETheta

Definition at line 223 of file G4Sphere.hh.

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

Constructor & Destructor Documentation

G4Sphere() [1/3]

G4Sphere::G4Sphere	(	const G4String &	pName,
		G4double	pRmin,
		G4double	pRmax,
		G4double	pSPhi,
		G4double	pDPhi,
		G4double	pSTheta,
		G4double	pDTheta
	)

Definition at line 89 of file G4Sphere.cc.

  : G4CSGSolid(pName), fEpsilon(2.e-11), fSPhi(0.0),
    fFullPhiSphere(true), fFullThetaSphere(true)
{
  kAngTolerance = G4GeometryTolerance::GetInstance()->GetAngularTolerance();
  kRadTolerance = G4GeometryTolerance::GetInstance()->GetRadialTolerance();

  halfCarTolerance = 0.5*kCarTolerance;
  halfAngTolerance = 0.5*kAngTolerance;

  // Check radii and set radial tolerances

  if ( (pRmin >= pRmax) || (pRmax < 1.1*kRadTolerance) || (pRmin < 0) )
  {
    std::ostringstream message;
    message << "Invalid radii for Solid: " << GetName() << G4endl
            << "        pRmin = " << pRmin << ", pRmax = " << pRmax;
    G4Exception("G4Sphere::G4Sphere()", "GeomSolids0002",
                FatalException, message);
  }
  fRmin=pRmin; fRmax=pRmax;
  fRminTolerance = (fRmin) ? std::max( kRadTolerance, fEpsilon*fRmin ) : 0;
  fRmaxTolerance = std::max( kRadTolerance, fEpsilon*fRmax );

  // Check angles

  CheckPhiAngles(pSPhi, pDPhi);
  CheckThetaAngles(pSTheta, pDTheta);
}

Here is the call graph for this function:

Here is the caller graph for this function:

~G4Sphere()

G4Sphere::~G4Sphere

(

)

Definition at line 144 of file G4Sphere.cc.

{
}

G4Sphere() [2/3]

G4Sphere::G4Sphere

(

__void__ &

a

)

Definition at line 127 of file G4Sphere.cc.

  : G4CSGSolid(a), fRminTolerance(0.), fRmaxTolerance(0.),
    kAngTolerance(0.), kRadTolerance(0.), fEpsilon(0.),
    fRmin(0.), fRmax(0.), fSPhi(0.), fDPhi(0.), fSTheta(0.),
    fDTheta(0.), sinCPhi(0.), cosCPhi(0.), cosHDPhiOT(0.), cosHDPhiIT(0.),
    sinSPhi(0.), cosSPhi(0.), sinEPhi(0.), cosEPhi(0.), hDPhi(0.), cPhi(0.),
    ePhi(0.), sinSTheta(0.), cosSTheta(0.), sinETheta(0.), cosETheta(0.),
    tanSTheta(0.), tanSTheta2(0.), tanETheta(0.), tanETheta2(0.), eTheta(0.),
    fFullPhiSphere(false), fFullThetaSphere(false), fFullSphere(true),
    halfCarTolerance(0.), halfAngTolerance(0.)
{
}

G4Sphere() [3/3]

G4Sphere::G4Sphere

(

const G4Sphere &

rhs

)

Definition at line 152 of file G4Sphere.cc.

  : G4CSGSolid(rhs), fRminTolerance(rhs.fRminTolerance),
    fRmaxTolerance(rhs.fRmaxTolerance), kAngTolerance(rhs.kAngTolerance),
    kRadTolerance(rhs.kRadTolerance), fEpsilon(rhs.fEpsilon),
    fRmin(rhs.fRmin), fRmax(rhs.fRmax), fSPhi(rhs.fSPhi), fDPhi(rhs.fDPhi),
    fSTheta(rhs.fSTheta), fDTheta(rhs.fDTheta),
    sinCPhi(rhs.sinCPhi), cosCPhi(rhs.cosCPhi),
    cosHDPhiOT(rhs.cosHDPhiOT), cosHDPhiIT(rhs.cosHDPhiIT),
    sinSPhi(rhs.sinSPhi), cosSPhi(rhs.cosSPhi),
    sinEPhi(rhs.sinEPhi), cosEPhi(rhs.cosEPhi),
    hDPhi(rhs.hDPhi), cPhi(rhs.cPhi), ePhi(rhs.ePhi),
    sinSTheta(rhs.sinSTheta), cosSTheta(rhs.cosSTheta),
    sinETheta(rhs.sinETheta), cosETheta(rhs.cosETheta),
    tanSTheta(rhs.tanSTheta), tanSTheta2(rhs.tanSTheta2),
    tanETheta(rhs.tanETheta), tanETheta2(rhs.tanETheta2), eTheta(rhs.eTheta),
    fFullPhiSphere(rhs.fFullPhiSphere), fFullThetaSphere(rhs.fFullThetaSphere),
    fFullSphere(rhs.fFullSphere),
    halfCarTolerance(rhs.halfCarTolerance),
    halfAngTolerance(rhs.halfAngTolerance)
{
}

Member Function Documentation

ApproxSurfaceNormal()

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

private

Definition at line 715 of file G4Sphere.cc.

{
  ENorm side;
  G4ThreeVector norm;
  G4double rho,rho2,radius,pPhi,pTheta;
  G4double distRMin,distRMax,distSPhi,distEPhi,
           distSTheta,distETheta,distMin;

  rho2=p.x()*p.x()+p.y()*p.y();
  radius=std::sqrt(rho2+p.z()*p.z());
  rho=std::sqrt(rho2);

  //
  // Distance to r shells
  //

  distRMax=std::fabs(radius-fRmax);
  if (fRmin)
  {
    distRMin=std::fabs(radius-fRmin);
      
    if (distRMin<distRMax)
    {
      distMin=distRMin;
      side=kNRMin;
    }
    else
    {   
      distMin=distRMax;
      side=kNRMax;
    }
  }
  else
  {
    distMin=distRMax;
    side=kNRMax;
  }

  //
  // Distance to phi planes
  //
  // Protected against (0,0,z) 
    
  pPhi = std::atan2(p.y(),p.x());
  if (pPhi<0) { pPhi += twopi; }

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

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

    // Find new minimum
    //
    if (distSPhi<distEPhi)
    {
      if (distSPhi<distMin)
      {
        distMin=distSPhi;
        side=kNSPhi;
      }
    }
    else
    {
      if (distEPhi<distMin)
      {
        distMin=distEPhi;
        side=kNEPhi;
      }
    }
  }

  //
  // Distance to theta planes
  //

  if (!fFullThetaSphere && radius)
  {
    pTheta=std::atan2(rho,p.z());
    distSTheta=std::fabs(pTheta-fSTheta)*radius;
    distETheta=std::fabs(pTheta-fSTheta-fDTheta)*radius;

    // Find new minimum
    //
    if (distSTheta<distETheta)
    {
      if (distSTheta<distMin)
      {
        distMin = distSTheta ;
        side = kNSTheta ;
      }
    }
    else
    {
      if (distETheta<distMin)
      {
        distMin = distETheta ;
        side = kNETheta ;
      }
    }
  }

  switch (side)
  {
    case kNRMin:      // Inner radius
      norm=G4ThreeVector(-p.x()/radius,-p.y()/radius,-p.z()/radius);
      break;
    case kNRMax:      // Outer radius
      norm=G4ThreeVector(p.x()/radius,p.y()/radius,p.z()/radius);
      break;
    case kNSPhi:
      norm=G4ThreeVector(sinSPhi,-cosSPhi,0);
      break;
    case kNEPhi:
      norm=G4ThreeVector(-sinEPhi,cosEPhi,0);
      break;
    case kNSTheta:
      norm=G4ThreeVector(-cosSTheta*std::cos(pPhi),
                         -cosSTheta*std::sin(pPhi),
                          sinSTheta            );
      break;
    case kNETheta:
      norm=G4ThreeVector( cosETheta*std::cos(pPhi),
                          cosETheta*std::sin(pPhi),
                         -sinETheta              );
      break;
    default:          // Should never reach this case ...
      DumpInfo();
      G4Exception("G4Sphere::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 G4Sphere::CalculateExtent ( const EAxis  pAxis, const G4VoxelLimits &  pVoxelLimit, const G4AffineTransform &  pTransform, G4double &  pmin, G4double &  pmax  ) const

virtual

Implements G4VSolid.

Definition at line 227 of file G4Sphere.cc.

{
  if ( fFullSphere )
  {
    // Special case handling for solid spheres-shells
    // (rotation doesn't influence).
    // 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 diff1,diff2,delta,maxDiff,newMin,newMax;
    G4double xoff1,xoff2,yoff1,yoff2;

    xoffset=pTransform.NetTranslation().x();
    xMin=xoffset-fRmax;
    xMax=xoffset+fRmax;
    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;
    yMax=yoffset+fRmax;
    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 sphere

    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
          //
          delta=fRmax*fRmax-yoff1*yoff1;
          diff1=(delta>0.) ? std::sqrt(delta) : 0.;
          delta=fRmax*fRmax-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
          //
          delta=fRmax*fRmax-xoff1*xoff1;
          diff1=(delta>0.) ? std::sqrt(delta) : 0.;
          delta=fRmax*fRmax-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       // Transformed cutted sphere
  {
    G4int i,j,noEntries,noBetweenSections;
    G4bool existsAfterClip=false;

    // Calculate rotated vertex coordinates

    G4ThreeVectorList* vertices;
    G4int  noPolygonVertices ;
    vertices=CreateRotatedVertices(pTransform,noPolygonVertices);

    pMin=+kInfinity;
    pMax=-kInfinity;

    noEntries=vertices->size();  // noPolygonVertices*noPhiCrossSections
    noBetweenSections=noEntries-noPolygonVertices;

    G4ThreeVectorList ThetaPolygon ;
    for (i=0;i<noEntries;i+=noPolygonVertices)
    {
      for(j=0;j<(noPolygonVertices/2)-1;j++)
      {
        ThetaPolygon.push_back((*vertices)[i+j]) ;      
        ThetaPolygon.push_back((*vertices)[i+j+1]) ;      
        ThetaPolygon.push_back((*vertices)[i+noPolygonVertices-2-j]) ;      
        ThetaPolygon.push_back((*vertices)[i+noPolygonVertices-1-j]) ;      
        CalculateClippedPolygonExtent(ThetaPolygon,pVoxelLimit,pAxis,pMin,pMax);
        ThetaPolygon.clear() ;
      }
    }
    for (i=0;i<noBetweenSections;i+=noPolygonVertices)
    {
      for(j=0;j<noPolygonVertices-1;j++)
      {
        ThetaPolygon.push_back((*vertices)[i+j]) ;      
        ThetaPolygon.push_back((*vertices)[i+j+1]) ;      
        ThetaPolygon.push_back((*vertices)[i+noPolygonVertices+j+1]) ;      
        ThetaPolygon.push_back((*vertices)[i+noPolygonVertices+j]) ;      
        CalculateClippedPolygonExtent(ThetaPolygon,pVoxelLimit,pAxis,pMin,pMax);
        ThetaPolygon.clear() ;
      }
      ThetaPolygon.push_back((*vertices)[i+noPolygonVertices-1]) ;      
      ThetaPolygon.push_back((*vertices)[i]) ;  
      ThetaPolygon.push_back((*vertices)[i+noPolygonVertices]) ;      
      ThetaPolygon.push_back((*vertices)[i+2*noPolygonVertices-1]) ;      
      CalculateClippedPolygonExtent(ThetaPolygon,pVoxelLimit,pAxis,pMin,pMax);
      ThetaPolygon.clear() ;
    }
      
    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:

CheckDPhiAngle()

void G4Sphere::CheckDPhiAngle ( G4double  dPhi )

inlineprivate

CheckPhiAngles()

void G4Sphere::CheckPhiAngles ( G4double  sPhi, G4double  dPhi  )

inlineprivate

Here is the caller graph for this function:

CheckSPhiAngle()

void G4Sphere::CheckSPhiAngle ( G4double  sPhi )

inlineprivate

CheckThetaAngles()

void G4Sphere::CheckThetaAngles ( G4double  sTheta, G4double  dTheta  )

inlineprivate

Here is the caller graph for this function:

Clone()

G4VSolid * G4Sphere::Clone ( ) const

virtual

Reimplemented from G4VSolid.

Definition at line 3054 of file G4Sphere.cc.

{
  return new G4Sphere(*this);
}

Here is the call graph for this function:

ComputeDimensions()

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

virtual

Reimplemented from G4VSolid.

Definition at line 216 of file G4Sphere.cc.

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

Here is the call graph for this function:

CreatePolyhedron()

G4Polyhedron * G4Sphere::CreatePolyhedron ( ) const

virtual

Reimplemented from G4VSolid.

Definition at line 3227 of file G4Sphere.cc.

{
  return new G4PolyhedronSphere (fRmin, fRmax, fSPhi, fDPhi, fSTheta, fDTheta);
}

CreateRotatedVertices()

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

private

Definition at line 2918 of file G4Sphere.cc.

{
  G4ThreeVectorList *vertices;
  G4ThreeVector vertex;
  G4double meshAnglePhi,meshRMax,crossAnglePhi,
           coscrossAnglePhi,sincrossAnglePhi,sAnglePhi;
  G4double meshTheta,crossTheta,startTheta;
  G4double rMaxX,rMaxY,rMinX,rMinY,rMinZ,rMaxZ;
  G4int crossSectionPhi,noPhiCrossSections,crossSectionTheta,noThetaSections;

  // Phi cross sections
    
  noPhiCrossSections = G4int(fDPhi/kMeshAngleDefault)+1;
    
  if (noPhiCrossSections<kMinMeshSections)
  {
    noPhiCrossSections=kMinMeshSections;
  }
  else if (noPhiCrossSections>kMaxMeshSections)
  {
    noPhiCrossSections=kMaxMeshSections;
  }
  meshAnglePhi=fDPhi/(noPhiCrossSections-1);
    
  // 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 (fFullPhiSphere)
  {
    sAnglePhi = -meshAnglePhi*0.5;
  }
    else
  {
    sAnglePhi=fSPhi;
  }    

  // Theta cross sections
    
  noThetaSections = G4int(fDTheta/kMeshAngleDefault)+1;
    
  if (noThetaSections<kMinMeshSections)
  {
    noThetaSections=kMinMeshSections;
  }
  else if (noThetaSections>kMaxMeshSections)
  {
    noThetaSections=kMaxMeshSections;
  }
  meshTheta=fDTheta/(noThetaSections-1);
    
  // If complete in Theta, set start angle such that mesh will be at fRMax
  // on the z axis. Will give better extent calculations when not rotated.
    
  if (fFullThetaSphere)
  {
    startTheta = -meshTheta*0.5;
  }
  else
  {
    startTheta=fSTheta;
  }    

  meshRMax = (meshAnglePhi >= meshTheta) ?
             fRmax/std::cos(meshAnglePhi*0.5) : fRmax/std::cos(meshTheta*0.5);
  G4double* cosCrossTheta = new G4double[noThetaSections];
  G4double* sinCrossTheta = new G4double[noThetaSections];    
  vertices=new G4ThreeVectorList();
  if (vertices && cosCrossTheta && sinCrossTheta)
  {
    vertices->reserve(noPhiCrossSections*(noThetaSections*2));
    for (crossSectionPhi=0;
         crossSectionPhi<noPhiCrossSections; crossSectionPhi++)
    {
      crossAnglePhi=sAnglePhi+crossSectionPhi*meshAnglePhi;
      coscrossAnglePhi=std::cos(crossAnglePhi);
      sincrossAnglePhi=std::sin(crossAnglePhi);
      for (crossSectionTheta=0;
           crossSectionTheta<noThetaSections;crossSectionTheta++)
      {
        // Compute coordinates of cross section at section crossSectionPhi
        //
        crossTheta=startTheta+crossSectionTheta*meshTheta;
        cosCrossTheta[crossSectionTheta]=std::cos(crossTheta);
        sinCrossTheta[crossSectionTheta]=std::sin(crossTheta);

        rMinX=fRmin*sinCrossTheta[crossSectionTheta]*coscrossAnglePhi;
        rMinY=fRmin*sinCrossTheta[crossSectionTheta]*sincrossAnglePhi;
        rMinZ=fRmin*cosCrossTheta[crossSectionTheta];
        
        vertex=G4ThreeVector(rMinX,rMinY,rMinZ);
        vertices->push_back(pTransform.TransformPoint(vertex));
        
      }    // Theta forward 
    
      for (crossSectionTheta=noThetaSections-1;
           crossSectionTheta>=0; crossSectionTheta--)
      {
        rMaxX=meshRMax*sinCrossTheta[crossSectionTheta]*coscrossAnglePhi;
        rMaxY=meshRMax*sinCrossTheta[crossSectionTheta]*sincrossAnglePhi;
        rMaxZ=meshRMax*cosCrossTheta[crossSectionTheta];
        
        vertex=G4ThreeVector(rMaxX,rMaxY,rMaxZ);
        vertices->push_back(pTransform.TransformPoint(vertex));

      }   // Theta back 
    }   // Phi
    noPolygonVertices = noThetaSections*2 ;
  }
  else
  {
    DumpInfo();
    G4Exception("G4Sphere::CreateRotatedVertices()",
                "GeomSolids0003", FatalException,
                "Error in allocation of vertices. Out of memory !");
  }

  delete [] cosCrossTheta;
  delete [] sinCrossTheta;

  return vertices;
}

Here is the call graph for this function:

Here is the caller graph for this function:

DescribeYourselfTo()

void G4Sphere::DescribeYourselfTo ( G4VGraphicsScene &  scene ) const

virtual

Implements G4VSolid.

Definition at line 3222 of file G4Sphere.cc.

{
  scene.AddSolid (*this);
}

Here is the call graph for this function:

DistanceToIn() [1/2]

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

virtual

Implements G4VSolid.

Definition at line 888 of file G4Sphere.cc.

{
  G4double snxt = kInfinity ;      // snxt = default return value
  G4double rho2, rad2, pDotV2d, pDotV3d, pTheta ;
  G4double tolSTheta=0., tolETheta=0. ;
  const G4double dRmax = 100.*fRmax;

  const G4double halfRmaxTolerance = fRmaxTolerance*0.5;
  const G4double halfRminTolerance = fRminTolerance*0.5;
  const G4double tolORMin2 = (fRmin>halfRminTolerance)
               ? (fRmin-halfRminTolerance)*(fRmin-halfRminTolerance) : 0;
  const G4double tolIRMin2 =
               (fRmin+halfRminTolerance)*(fRmin+halfRminTolerance);
  const G4double tolORMax2 =
               (fRmax+halfRmaxTolerance)*(fRmax+halfRmaxTolerance);
  const G4double tolIRMax2 =
               (fRmax-halfRmaxTolerance)*(fRmax-halfRmaxTolerance);

  // Intersection point
  //
  G4double xi, yi, zi, rhoi, rhoi2, radi2, iTheta ;

  // Phi intersection
  //
  G4double Comp ; 

  // Phi precalcs
  //
  G4double Dist, cosPsi ;

  // Theta precalcs
  //
  G4double dist2STheta, dist2ETheta ;
  G4double t1, t2, b, c, d2, d, sd = kInfinity ;

  // General Precalcs
  //
  rho2 = p.x()*p.x() + p.y()*p.y() ;
  rad2 = rho2 + p.z()*p.z() ;
  pTheta = std::atan2(std::sqrt(rho2),p.z()) ;

  pDotV2d = p.x()*v.x() + p.y()*v.y() ;
  pDotV3d = pDotV2d + p.z()*v.z() ;

  // Theta precalcs
  //
  if (!fFullThetaSphere)
  {
    tolSTheta = fSTheta - halfAngTolerance ;
    tolETheta = eTheta + halfAngTolerance ;

    // Special case rad2 = 0 comparing with direction
    //
    if ((rad2!=0.0) || (fRmin!=0.0))
    {
      // Keep going for computation of distance...
    }
    else  // Positioned on the sphere's origin
    {
      G4double vTheta = std::atan2(std::sqrt(v.x()*v.x()+v.y()*v.y()),v.z()) ;
      if ( (vTheta < tolSTheta) || (vTheta > tolETheta) )
      {
        return snxt ; // kInfinity
      }
      return snxt = 0.0 ;
    }
  }

  // Outer spherical shell intersection
  // - Only if outside tolerant fRmax
  // - Check for if inside and outer G4Sphere heading through solid (-> 0)
  // - No intersect -> no intersection with G4Sphere
  //
  // Shell eqn: x^2+y^2+z^2=RSPH^2
  //
  // => (px+svx)^2+(py+svy)^2+(pz+svz)^2=R^2
  //
  // => (px^2+py^2+pz^2) +2sd(pxvx+pyvy+pzvz)+sd^2(vx^2+vy^2+vz^2)=R^2
  // =>      rad2        +2sd(pDotV3d)       +sd^2                =R^2
  //
  // => sd=-pDotV3d+-std::sqrt(pDotV3d^2-(rad2-R^2))

  c = rad2 - fRmax*fRmax ;

  if (c > fRmaxTolerance*fRmax)
  {
    // If outside tolerant boundary of outer G4Sphere
    // [should be std::sqrt(rad2)-fRmax > halfRmaxTolerance]

    d2 = pDotV3d*pDotV3d - c ;

    if ( d2 >= 0 )
    {
      sd = -pDotV3d - std::sqrt(d2) ;

      if (sd >= 0 )
      {
        if ( sd>dRmax ) // Avoid rounding errors due to precision issues seen on
        {               // 64 bits systems. Split long distances and recompute
          G4double fTerm = sd-std::fmod(sd,dRmax);
          sd = fTerm + DistanceToIn(p+fTerm*v,v);
        } 
        xi   = p.x() + sd*v.x() ;
        yi   = p.y() + sd*v.y() ;
        rhoi = std::sqrt(xi*xi + yi*yi) ;

        if (!fFullPhiSphere && rhoi)    // Check phi intersection
        {
          cosPsi = (xi*cosCPhi + yi*sinCPhi)/rhoi ;

          if (cosPsi >= cosHDPhiOT)
          {
            if (!fFullThetaSphere)   // Check theta intersection
            {
              zi = p.z() + sd*v.z() ;

              // rhoi & zi can never both be 0
              // (=>intersect at origin =>fRmax=0)
              //
              iTheta = std::atan2(rhoi,zi) ;
              if ( (iTheta >= tolSTheta) && (iTheta <= tolETheta) )
              {
                return snxt = sd ;
              }
            }
            else
            {
              return snxt=sd;
            }
          }
        }
        else
        {
          if (!fFullThetaSphere)    // Check theta intersection
          {
            zi = p.z() + sd*v.z() ;

            // rhoi & zi can never both be 0
            // (=>intersect at origin => fRmax=0 !)
            //
            iTheta = std::atan2(rhoi,zi) ;
            if ( (iTheta >= tolSTheta) && (iTheta <= tolETheta) )
            {
              return snxt=sd;
            }
          }
          else
          {
            return snxt = sd;
          }
        }          
      }
    }
    else    // No intersection with G4Sphere
    {
      return snxt=kInfinity;
    }
  }
  else
  {
    // Inside outer radius
    // check not inside, and heading through G4Sphere (-> 0 to in)

    d2 = pDotV3d*pDotV3d - c ;

    if ( (rad2 > tolIRMax2)
      && ( (d2 >= fRmaxTolerance*fRmax) && (pDotV3d < 0) ) )
    {
      if (!fFullPhiSphere)
      {
        // Use inner phi tolerant boundary -> if on tolerant
        // phi boundaries, phi intersect code handles leaving/entering checks

        cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/std::sqrt(rho2) ;

        if (cosPsi>=cosHDPhiIT)
        { 
          // inside radii, delta r -ve, inside phi

          if ( !fFullThetaSphere )
          {
            if ( (pTheta >= tolSTheta + kAngTolerance)
              && (pTheta <= tolETheta - kAngTolerance) )
            {
              return snxt=0;
            }
          }
          else    // strictly inside Theta in both cases
          {
            return snxt=0;
          }
        }
      }
      else
      {
        if ( !fFullThetaSphere )
        {
          if ( (pTheta >= tolSTheta + kAngTolerance)
            && (pTheta <= tolETheta - kAngTolerance) )
          {
            return snxt=0;
          }
        }
        else   // strictly inside Theta in both cases
        {
          return snxt=0;
        }
      }
    }
  }

  // Inner spherical shell intersection
  // - Always farthest root, because would have passed through outer
  //   surface first.
  // - Tolerant check if travelling through solid

  if (fRmin)
  {
    c  = rad2 - fRmin*fRmin ;
    d2 = pDotV3d*pDotV3d - c ;

    // Within tolerance inner radius of inner G4Sphere
    // Check for immediate entry/already inside and travelling outwards

    if ( (c > -halfRminTolerance) && (rad2 < tolIRMin2)
      && ( (d2 < fRmin*kCarTolerance) || (pDotV3d >= 0) ) )
    {
      if ( !fFullPhiSphere )
      {
        // Use inner phi tolerant boundary -> if on tolerant
        // phi boundaries, phi intersect code handles leaving/entering checks

        cosPsi = (p.x()*cosCPhi+p.y()*sinCPhi)/std::sqrt(rho2) ;
        if (cosPsi >= cosHDPhiIT)
        { 
          // inside radii, delta r -ve, inside phi
          //
          if ( !fFullThetaSphere )
          {
            if ( (pTheta >= tolSTheta + kAngTolerance)
              && (pTheta <= tolETheta - kAngTolerance) )
            {
              return snxt=0;
            }
          }
          else
          {
            return snxt = 0 ;
          }
        }
      }
      else
      {
        if ( !fFullThetaSphere )
        {
          if ( (pTheta >= tolSTheta + kAngTolerance)
            && (pTheta <= tolETheta - kAngTolerance) )
          {
            return snxt = 0 ;
          }
        }
        else
        {
          return snxt=0;
        }
      }
    }
    else   // Not special tolerant case
    {
      if (d2 >= 0)
      {
        sd = -pDotV3d + std::sqrt(d2) ;
        if ( sd >= halfRminTolerance )  // It was >= 0 ??
        {
          xi   = p.x() + sd*v.x() ;
          yi   = p.y() + sd*v.y() ;
          rhoi = std::sqrt(xi*xi+yi*yi) ;

          if ( !fFullPhiSphere && rhoi )   // Check phi intersection
          {
            cosPsi = (xi*cosCPhi + yi*sinCPhi)/rhoi ;

            if (cosPsi >= cosHDPhiOT)
            {
              if ( !fFullThetaSphere )  // Check theta intersection
              {
                zi = p.z() + sd*v.z() ;

                // rhoi & zi can never both be 0
                // (=>intersect at origin =>fRmax=0)
                //
                iTheta = std::atan2(rhoi,zi) ;
                if ( (iTheta >= tolSTheta) && (iTheta<=tolETheta) )
                {
                  snxt = sd;
                }
              }
              else
              {
                snxt=sd;
              }
            }
          }
          else
          {
            if ( !fFullThetaSphere )   // Check theta intersection
            {
              zi = p.z() + sd*v.z() ;

              // rhoi & zi can never both be 0
              // (=>intersect at origin => fRmax=0 !)
              //
              iTheta = std::atan2(rhoi,zi) ;
              if ( (iTheta >= tolSTheta) && (iTheta <= tolETheta) )
              {
                snxt = sd;
              }
            }
            else
            {
              snxt = sd;
            }
          }
        }
      }
    }
  }

  // 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
  //         -> Should use some form of loop Construct
  //
  if ( !fFullPhiSphere )
  {
    // First phi surface ('S'tarting phi)
    // Comp = Component in outwards normal dirn
    //
    Comp = v.x()*sinSPhi - v.y()*cosSPhi ;
                    
    if ( Comp < 0 )
    {
      Dist = p.y()*cosSPhi - p.x()*sinSPhi ;

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

        if (sd < snxt)
        {
          if ( sd > 0 )
          {
            xi    = p.x() + sd*v.x() ;
            yi    = p.y() + sd*v.y() ;
            zi    = p.z() + sd*v.z() ;
            rhoi2 = xi*xi + yi*yi   ;
            radi2 = rhoi2 + zi*zi   ;
          }
          else
          {
            sd    = 0     ;
            xi    = p.x() ;
            yi    = p.y() ;
            zi    = p.z() ;
            rhoi2 = rho2  ;
            radi2 = rad2  ;
          }
          if ( (radi2 <= tolORMax2)
            && (radi2 >= tolORMin2)
            && ((yi*cosCPhi-xi*sinCPhi) <= 0) )
          {
            // Check theta intersection
            // rhoi & zi can never both be 0
            // (=>intersect at origin =>fRmax=0)
            //
            if ( !fFullThetaSphere )
            {
              iTheta = std::atan2(std::sqrt(rhoi2),zi) ;
              if ( (iTheta >= tolSTheta) && (iTheta <= tolETheta) )
              {
                // r and theta intersections good
                // - check intersecting with correct half-plane

                if ((yi*cosCPhi-xi*sinCPhi) <= 0)
                {
                  snxt = sd;
                }
              }
            }
            else
            {
              snxt = sd;
            }
          }
        }
      }
    }

    // Second phi surface ('E'nding phi)
    // Component in outwards normal dirn

    Comp = -( v.x()*sinEPhi-v.y()*cosEPhi ) ;
        
    if (Comp < 0)
    {
      Dist = -(p.y()*cosEPhi-p.x()*sinEPhi) ;
      if ( Dist < halfCarTolerance )
      {
        sd = Dist/Comp ;

        if ( sd < snxt )
        {
          if (sd > 0)
          {
            xi    = p.x() + sd*v.x() ;
            yi    = p.y() + sd*v.y() ;
            zi    = p.z() + sd*v.z() ;
            rhoi2 = xi*xi + yi*yi   ;
            radi2 = rhoi2 + zi*zi   ;
          }
          else
          {
            sd    = 0     ;
            xi    = p.x() ;
            yi    = p.y() ;
            zi    = p.z() ;
            rhoi2 = rho2  ;
            radi2 = rad2  ;
          }
          if ( (radi2 <= tolORMax2)
            && (radi2 >= tolORMin2)
            && ((yi*cosCPhi-xi*sinCPhi) >= 0) )
          {
            // Check theta intersection
            // rhoi & zi can never both be 0
            // (=>intersect at origin =>fRmax=0)
            //
            if ( !fFullThetaSphere )
            {
              iTheta = std::atan2(std::sqrt(rhoi2),zi) ;
              if ( (iTheta >= tolSTheta) && (iTheta <= tolETheta) )
              {
                // r and theta intersections good
                // - check intersecting with correct half-plane

                if ((yi*cosCPhi-xi*sinCPhi) >= 0)
                {
                  snxt = sd;
                }
              }
            }
            else
            {
              snxt = sd;
            }
          }
        }
      }
    }
  }

  // Theta segment intersection

  if ( !fFullThetaSphere )
  {

    // Intersection with theta surfaces
    // Known failure cases:
    // o  Inside tolerance of stheta surface, skim
    //    ~parallel to cone and Hit & enter etheta surface [& visa versa]
    //
    //    To solve: Check 2nd root of etheta surface in addition to stheta
    //
    // o  start/end theta is exactly pi/2 
    // Intersections with cones
    //
    // Cone equation: x^2+y^2=z^2tan^2(t)
    //
    // => (px+svx)^2+(py+svy)^2=(pz+svz)^2tan^2(t)
    //
    // => (px^2+py^2-pz^2tan^2(t))+2sd(pxvx+pyvy-pzvztan^2(t))
    //       + sd^2(vx^2+vy^2-vz^2tan^2(t)) = 0
    //
    // => sd^2(1-vz^2(1+tan^2(t))+2sd(pdotv2d-pzvztan^2(t))+(rho2-pz^2tan^2(t))=0

    if (fSTheta)
    {
      dist2STheta = rho2 - p.z()*p.z()*tanSTheta2 ;
    }
    else
    {
      dist2STheta = kInfinity ;
    }
    if ( eTheta < pi )
    {
      dist2ETheta=rho2-p.z()*p.z()*tanETheta2;
    }
    else
    {
      dist2ETheta=kInfinity;
    }      
    if ( pTheta < tolSTheta )
    {
      // Inside (theta<stheta-tol) stheta cone
      // First root of stheta cone, second if first root -ve

      t1 = 1 - v.z()*v.z()*(1 + tanSTheta2) ;
      t2 = pDotV2d - p.z()*v.z()*tanSTheta2 ;
      if (t1)
      {   
        b  = t2/t1 ;
        c  = dist2STheta/t1 ;
        d2 = b*b - c ;

        if ( d2 >= 0 )
        {
          d  = std::sqrt(d2) ;
          sd = -b - d ;    // First root
          zi = p.z() + sd*v.z();

          if ( (sd < 0) || (zi*(fSTheta - halfpi) > 0) )
          {
            sd = -b+d;    // Second root
          }
          if ((sd >= 0) && (sd < snxt))
          {
            xi    = p.x() + sd*v.x();
            yi    = p.y() + sd*v.y();
            zi    = p.z() + sd*v.z();
            rhoi2 = xi*xi + yi*yi;
            radi2 = rhoi2 + zi*zi;
            if ( (radi2 <= tolORMax2)
              && (radi2 >= tolORMin2)
              && (zi*(fSTheta - halfpi) <= 0) )
            {
              if ( !fFullPhiSphere && rhoi2 )  // Check phi intersection
              {
                cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
                if (cosPsi >= cosHDPhiOT)
                {
                  snxt = sd;
                }
              }
              else
              {
                snxt = sd;
              }
            }
          }
        }
      }

      // Possible intersection with ETheta cone. 
      // Second >= 0 root should be considered
        
      if ( eTheta < pi )
      {
        t1 = 1 - v.z()*v.z()*(1 + tanETheta2) ;
        t2 = pDotV2d - p.z()*v.z()*tanETheta2 ;
        if (t1)
        { 
          b  = t2/t1 ;
          c  = dist2ETheta/t1 ;
          d2 = b*b - c ;

          if (d2 >= 0)
          {
            d  = std::sqrt(d2) ;
            sd = -b + d ;    // Second root

            if ( (sd >= 0) && (sd < snxt) )
            {
              xi    = p.x() + sd*v.x() ;
              yi    = p.y() + sd*v.y() ;
              zi    = p.z() + sd*v.z() ;
              rhoi2 = xi*xi + yi*yi   ;
              radi2 = rhoi2 + zi*zi   ;

              if ( (radi2 <= tolORMax2)
                && (radi2 >= tolORMin2)
                && (zi*(eTheta - halfpi) <= 0) )
              {
                if (!fFullPhiSphere && rhoi2)   // Check phi intersection
                {
                  cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
                  if (cosPsi >= cosHDPhiOT)
                  {
                    snxt = sd;
                  }
                }
                else
                {
                  snxt = sd;
                }
              }
            }
          }
        }
      }
    }  
    else if ( pTheta > tolETheta ) 
    { 
      // dist2ETheta<-kRadTolerance*0.5 && dist2STheta>0)
      // Inside (theta > etheta+tol) e-theta cone
      // First root of etheta cone, second if first root 'imaginary'

      t1 = 1 - v.z()*v.z()*(1 + tanETheta2) ;
      t2 = pDotV2d - p.z()*v.z()*tanETheta2 ;
      if (t1)
      {  
        b  = t2/t1 ;
        c  = dist2ETheta/t1 ;
        d2 = b*b - c ;

        if (d2 >= 0)
        {
          d  = std::sqrt(d2) ;
          sd = -b - d ;    // First root
          zi = p.z() + sd*v.z();

          if ( (sd < 0) || (zi*(eTheta - halfpi) > 0) )
          {
            sd = -b + d ;           // second root
          }
          if ( (sd >= 0) && (sd < snxt) )
          {
            xi    = p.x() + sd*v.x() ;
            yi    = p.y() + sd*v.y() ;
            zi    = p.z() + sd*v.z() ;
            rhoi2 = xi*xi + yi*yi   ;
            radi2 = rhoi2 + zi*zi   ;

            if ( (radi2 <= tolORMax2)
              && (radi2 >= tolORMin2) 
              && (zi*(eTheta - halfpi) <= 0) )
            {
              if (!fFullPhiSphere && rhoi2)  // Check phi intersection
              {
                cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
                if (cosPsi >= cosHDPhiOT)
                {
                  snxt = sd;
                }
              }
              else
              {
                snxt = sd;
              }
            }
          }
        }
      }

      // Possible intersection with STheta cone. 
      // Second >= 0 root should be considered
        
      if ( fSTheta )
      {
        t1 = 1 - v.z()*v.z()*(1 + tanSTheta2) ;
        t2 = pDotV2d - p.z()*v.z()*tanSTheta2 ;
        if (t1)
        { 
          b  = t2/t1 ;
          c  = dist2STheta/t1 ;
          d2 = b*b - c ;

          if (d2 >= 0)
          {
            d  = std::sqrt(d2) ;
            sd = -b + d ;    // Second root

            if ( (sd >= 0) && (sd < snxt) )
            {
              xi    = p.x() + sd*v.x() ;
              yi    = p.y() + sd*v.y() ;
              zi    = p.z() + sd*v.z() ;
              rhoi2 = xi*xi + yi*yi   ;
              radi2 = rhoi2 + zi*zi   ;

              if ( (radi2 <= tolORMax2)
                && (radi2 >= tolORMin2)
                && (zi*(fSTheta - halfpi) <= 0) )
              {
                if (!fFullPhiSphere && rhoi2)   // Check phi intersection
                {
                  cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
                  if (cosPsi >= cosHDPhiOT)
                  {
                    snxt = sd;
                  }
                }
                else
                {
                  snxt = sd;
                }
              }
            }
          }
        }
      }  
    }     
    else if ( (pTheta < tolSTheta + kAngTolerance)
           && (fSTheta > halfAngTolerance) )
    {
      // In tolerance of stheta
      // If entering through solid [r,phi] => 0 to in
      // else try 2nd root

      t2 = pDotV2d - p.z()*v.z()*tanSTheta2 ;
      if ( (t2>=0 && tolIRMin2<rad2 && rad2<tolIRMax2 && fSTheta<halfpi)
        || (t2<0  && tolIRMin2<rad2 && rad2<tolIRMax2 && fSTheta>halfpi)
        || (v.z()<0 && tolIRMin2<rad2 && rad2<tolIRMax2 && fSTheta==halfpi) )
      {
        if (!fFullPhiSphere && rho2)  // Check phi intersection
        {
          cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/std::sqrt(rho2) ;
          if (cosPsi >= cosHDPhiIT)
          {
            return 0 ;
          }
        }
        else
        {
          return 0 ;
        }
      }

      // Not entering immediately/travelling through

      t1 = 1 - v.z()*v.z()*(1 + tanSTheta2) ;
      if (t1)
      { 
        b  = t2/t1 ;
        c  = dist2STheta/t1 ;
        d2 = b*b - c ;

        if (d2 >= 0)
        {
          d  = std::sqrt(d2) ;
          sd = -b + d ;
          if ( (sd >= halfCarTolerance) && (sd < snxt) && (fSTheta < halfpi) )
          {   // ^^^^^^^^^^^^^^^^^^^^^  shouldn't it be >=0 instead ?
            xi    = p.x() + sd*v.x() ;
            yi    = p.y() + sd*v.y() ;
            zi    = p.z() + sd*v.z() ;
            rhoi2 = xi*xi + yi*yi   ;
            radi2 = rhoi2 + zi*zi   ;

            if ( (radi2 <= tolORMax2)
              && (radi2 >= tolORMin2)
              && (zi*(fSTheta - halfpi) <= 0) )
            {
              if ( !fFullPhiSphere && rhoi2 )    // Check phi intersection
              {
                cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
                if ( cosPsi >= cosHDPhiOT )
                {
                  snxt = sd;
                }
              }
              else
              {
                snxt = sd;
              }
            }
          }
        }
      }
    }   
    else if ((pTheta > tolETheta-kAngTolerance) && (eTheta < pi-kAngTolerance))
    {

      // In tolerance of etheta
      // If entering through solid [r,phi] => 0 to in
      // else try 2nd root

      t2 = pDotV2d - p.z()*v.z()*tanETheta2 ;

      if (   ((t2<0) && (eTheta < halfpi)
          && (tolIRMin2 < rad2) && (rad2 < tolIRMax2))
        ||   ((t2>=0) && (eTheta > halfpi)
          && (tolIRMin2 < rad2) && (rad2 < tolIRMax2))
        ||   ((v.z()>0) && (eTheta == halfpi)
          && (tolIRMin2 < rad2) && (rad2 < tolIRMax2))  )
      {
        if (!fFullPhiSphere && rho2)   // Check phi intersection
        {
          cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/std::sqrt(rho2) ;
          if (cosPsi >= cosHDPhiIT)
          {
            return 0 ;
          }
        }
        else
        {
          return 0 ;
        }
      }

      // Not entering immediately/travelling through

      t1 = 1 - v.z()*v.z()*(1 + tanETheta2) ;
      if (t1)
      { 
        b  = t2/t1 ;
        c  = dist2ETheta/t1 ;
        d2 = b*b - c ;

        if (d2 >= 0)
        {
          d  = std::sqrt(d2) ;
          sd = -b + d ;
        
          if ( (sd >= halfCarTolerance)
            && (sd < snxt) && (eTheta > halfpi) )
          {
            xi    = p.x() + sd*v.x() ;
            yi    = p.y() + sd*v.y() ;
            zi    = p.z() + sd*v.z() ;
            rhoi2 = xi*xi + yi*yi   ;
            radi2 = rhoi2 + zi*zi   ;

            if ( (radi2 <= tolORMax2)
              && (radi2 >= tolORMin2)
              && (zi*(eTheta - halfpi) <= 0) )
            {
              if (!fFullPhiSphere && rhoi2)   // Check phi intersection
              {
                cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
                if (cosPsi >= cosHDPhiOT)
                {
                  snxt = sd;
                }
              }
              else
              {
                snxt = sd;
              }
            }
          }
        } 
      }   
    }  
    else
    {
      // stheta+tol<theta<etheta-tol
      // For BOTH stheta & etheta check 2nd root for validity [r,phi]

      t1 = 1 - v.z()*v.z()*(1 + tanSTheta2) ;
      t2 = pDotV2d - p.z()*v.z()*tanSTheta2 ;
      if (t1)
      { 
        b  = t2/t1;
        c  = dist2STheta/t1 ;
        d2 = b*b - c ;

        if (d2 >= 0)
        {
          d  = std::sqrt(d2) ;
          sd = -b + d ;    // second root

          if ((sd >= 0) && (sd < snxt))
          {
            xi    = p.x() + sd*v.x() ;
            yi    = p.y() + sd*v.y() ;
            zi    = p.z() + sd*v.z() ;
            rhoi2 = xi*xi + yi*yi   ;
            radi2 = rhoi2 + zi*zi   ;

            if ( (radi2 <= tolORMax2)
              && (radi2 >= tolORMin2)
              && (zi*(fSTheta - halfpi) <= 0) )
            {
              if (!fFullPhiSphere && rhoi2)   // Check phi intersection
              {
                cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
                if (cosPsi >= cosHDPhiOT)
                {
                  snxt = sd;
                }
              }
              else
              {
                snxt = sd;
              }
            }
          }
        }
      }        
      t1 = 1 - v.z()*v.z()*(1 + tanETheta2) ;
      t2 = pDotV2d - p.z()*v.z()*tanETheta2 ;
      if (t1)
      {   
        b  = t2/t1 ;
        c  = dist2ETheta/t1 ;
        d2 = b*b - c ;

        if (d2 >= 0)
        {
          d  = std::sqrt(d2) ;
          sd = -b + d;    // second root

          if ((sd >= 0) && (sd < snxt))
          {
            xi    = p.x() + sd*v.x() ;
            yi    = p.y() + sd*v.y() ;
            zi    = p.z() + sd*v.z() ;
            rhoi2 = xi*xi + yi*yi   ;
            radi2 = rhoi2 + zi*zi   ;

            if ( (radi2 <= tolORMax2)
              && (radi2 >= tolORMin2)
              && (zi*(eTheta - halfpi) <= 0) )
            {
              if (!fFullPhiSphere && rhoi2)   // Check phi intersection
              {
                cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
                if ( cosPsi >= cosHDPhiOT )
                {
                  snxt = sd;
                }
              }
              else
              {
                snxt = sd;
              }
            }
          }
        }
      }
    }  
  }
  return snxt;
}

Here is the call graph for this function:

DistanceToIn() [2/2]

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

virtual

Implements G4VSolid.

Definition at line 1836 of file G4Sphere.cc.

{
  G4double safe=0.0,safeRMin,safeRMax,safePhi,safeTheta;
  G4double rho2,rds,rho;
  G4double cosPsi;
  G4double pTheta,dTheta1,dTheta2;
  rho2=p.x()*p.x()+p.y()*p.y();
  rds=std::sqrt(rho2+p.z()*p.z());
  rho=std::sqrt(rho2);

  //
  // Distance to r shells
  //    
  if (fRmin)
  {
    safeRMin=fRmin-rds;
    safeRMax=rds-fRmax;
    if (safeRMin>safeRMax)
    {
      safe=safeRMin;
    }
    else
    {
      safe=safeRMax;
    }
  }
  else
  {
    safe=rds-fRmax;
  }

  //
  // Distance to phi extent
  //
  if (!fFullPhiSphere && rho)
  {
    // Psi=angle from central phi to point
    //
    cosPsi=(p.x()*cosCPhi+p.y()*sinCPhi)/rho;
    if (cosPsi<std::cos(hDPhi))
    {
      // Point lies outside phi range
      //
      if ((p.y()*cosCPhi-p.x()*sinCPhi)<=0)
      {
        safePhi=std::fabs(p.x()*sinSPhi-p.y()*cosSPhi);
      }
      else
      {
        safePhi=std::fabs(p.x()*sinEPhi-p.y()*cosEPhi);
      }
      if (safePhi>safe)  { safe=safePhi; }
    }
  }
  //
  // Distance to Theta extent
  //    
  if ((rds!=0.0) && (!fFullThetaSphere))
  {
    pTheta=std::acos(p.z()/rds);
    if (pTheta<0)  { pTheta+=pi; }
    dTheta1=fSTheta-pTheta;
    dTheta2=pTheta-eTheta;
    if (dTheta1>dTheta2)
    {
      if (dTheta1>=0)             // WHY ???????????
      {
        safeTheta=rds*std::sin(dTheta1);
        if (safe<=safeTheta)
        {
          safe=safeTheta;
        }
      }
    }
    else
    {
      if (dTheta2>=0)
      {
        safeTheta=rds*std::sin(dTheta2);
        if (safe<=safeTheta)
        {
          safe=safeTheta;
        }
      }
    }
  }

  if (safe<0)  { safe=0; }
  return safe;
}

Here is the call graph for this function:

DistanceToOut() [1/2]

G4double G4Sphere::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 1932 of file G4Sphere.cc.

{
  G4double snxt = kInfinity;     // snxt is default return value
  G4double sphi= kInfinity,stheta= kInfinity;
  ESide side=kNull,sidephi=kNull,sidetheta=kNull;  

  const G4double halfRmaxTolerance = fRmaxTolerance*0.5;
  const G4double halfRminTolerance = fRminTolerance*0.5;
  const G4double Rmax_plus  = fRmax + halfRmaxTolerance;
  const G4double Rmin_minus = (fRmin) ? fRmin-halfRminTolerance : 0;
  G4double t1,t2;
  G4double b,c,d;

  // Variables for phi intersection:

  G4double pDistS,compS,pDistE,compE,sphi2,vphi;
    
  G4double rho2,rad2,pDotV2d,pDotV3d;

  G4double xi,yi,zi;      // Intersection point

  // Theta precals
  //
  G4double rhoSecTheta;
  G4double dist2STheta, dist2ETheta, distTheta;
  G4double d2,sd;

  // General Precalcs
  //
  rho2 = p.x()*p.x()+p.y()*p.y();
  rad2 = rho2+p.z()*p.z();

  pDotV2d = p.x()*v.x()+p.y()*v.y();
  pDotV3d = pDotV2d+p.z()*v.z();

  // Radial Intersections from G4Sphere::DistanceToIn
  //
  // Outer spherical shell intersection
  // - Only if outside tolerant fRmax
  // - Check for if inside and outer G4Sphere heading through solid (-> 0)
  // - No intersect -> no intersection with G4Sphere
  //
  // Shell eqn: x^2+y^2+z^2=RSPH^2
  //
  // => (px+svx)^2+(py+svy)^2+(pz+svz)^2=R^2
  //
  // => (px^2+py^2+pz^2) +2sd(pxvx+pyvy+pzvz)+sd^2(vx^2+vy^2+vz^2)=R^2
  // =>      rad2        +2sd(pDotV3d)       +sd^2                =R^2
  //
  // => sd=-pDotV3d+-std::sqrt(pDotV3d^2-(rad2-R^2))

  if( (rad2 <= Rmax_plus*Rmax_plus) && (rad2 >= Rmin_minus*Rmin_minus) )
  {
    c = rad2 - fRmax*fRmax;

    if (c < fRmaxTolerance*fRmax) 
    {
      // Within tolerant Outer radius 
      // 
      // The test is
      //     rad  - fRmax < 0.5*kRadTolerance
      // =>  rad  < fRmax + 0.5*kRadTol
      // =>  rad2 < (fRmax + 0.5*kRadTol)^2
      // =>  rad2 < fRmax^2 + 2.*0.5*fRmax*kRadTol + 0.25*kRadTol*kRadTol
      // =>  rad2 - fRmax^2    <~    fRmax*kRadTol 

      d2 = pDotV3d*pDotV3d - c;

      if( (c >- fRmaxTolerance*fRmax)       // on tolerant surface
       && ((pDotV3d >=0) || (d2 < 0)) )     // leaving outside from Rmax 
                                            // not re-entering
      {
        if(calcNorm)
        {
          *validNorm = true ;
          *n         = G4ThreeVector(p.x()/fRmax,p.y()/fRmax,p.z()/fRmax) ;
        }
        return snxt = 0;
      }
      else 
      {
        snxt = -pDotV3d+std::sqrt(d2);    // second root since inside Rmax
        side =  kRMax ; 
      }
    }

    // Inner spherical shell intersection:
    // Always first >=0 root, because would have passed
    // from outside of Rmin surface .

    if (fRmin)
    {
      c  = rad2 - fRmin*fRmin;
      d2 = pDotV3d*pDotV3d - c;

      if (c >- fRminTolerance*fRmin) // 2.0 * (0.5*kRadTolerance) * fRmin
      {
        if ( (c < fRminTolerance*fRmin)              // leaving from Rmin
          && (d2 >= fRminTolerance*fRmin) && (pDotV3d < 0) )
        {
          if(calcNorm)  { *validNorm = false; }  // Rmin surface is concave         
          return snxt = 0 ;
        }
        else
        {  
          if ( d2 >= 0. )
          {
            sd = -pDotV3d-std::sqrt(d2);

            if ( sd >= 0. )     // Always intersect Rmin first
            {
              snxt = sd ;
              side = kRMin ;
            }
          }
        }
      }
    }
  }

  // Theta segment intersection

  if ( !fFullThetaSphere )
  {
    // Intersection with theta surfaces
    //
    // Known failure cases:
    // o  Inside tolerance of stheta surface, skim
    //    ~parallel to cone and Hit & enter etheta surface [& visa versa]
    //
    //    To solve: Check 2nd root of etheta surface in addition to stheta
    //
    // o  start/end theta is exactly pi/2 
    //
    // Intersections with cones
    //
    // Cone equation: x^2+y^2=z^2tan^2(t)
    //
    // => (px+svx)^2+(py+svy)^2=(pz+svz)^2tan^2(t)
    //
    // => (px^2+py^2-pz^2tan^2(t))+2sd(pxvx+pyvy-pzvztan^2(t))
    //       + sd^2(vx^2+vy^2-vz^2tan^2(t)) = 0
    //
    // => sd^2(1-vz^2(1+tan^2(t))+2sd(pdotv2d-pzvztan^2(t))+(rho2-pz^2tan^2(t))=0
    //
  
    if(fSTheta) // intersection with first cons
    {
      if( std::fabs(tanSTheta) > 5./kAngTolerance ) // kons is plane z=0
      {
        if( v.z() > 0. ) 
        {
          if ( std::fabs( p.z() ) <= halfRmaxTolerance )
          {
            if(calcNorm)
            {
              *validNorm = true;
              *n = G4ThreeVector(0.,0.,1.);
            }
            return snxt = 0 ;
          }  
          stheta    = -p.z()/v.z();
          sidetheta = kSTheta;
        }
      }
      else // kons is not plane 
      {
        t1          = 1-v.z()*v.z()*(1+tanSTheta2);
        t2          = pDotV2d-p.z()*v.z()*tanSTheta2;  // ~vDotN if p on cons
        dist2STheta = rho2-p.z()*p.z()*tanSTheta2;     // t3

        distTheta = std::sqrt(rho2)-p.z()*tanSTheta;

        if( std::fabs(t1) < halfAngTolerance ) // 1st order equation,
        {                                      // v parallel to kons
          if( v.z() > 0. )
          {
            if(std::fabs(distTheta) < halfRmaxTolerance) // p on surface
            {
              if( (fSTheta < halfpi) && (p.z() > 0.) )
              {
                if( calcNorm )  { *validNorm = false; }
                return snxt = 0.;
              }
              else if( (fSTheta > halfpi) && (p.z() <= 0) )
              {
                if( calcNorm ) 
                {
                  *validNorm = true;
                  if (rho2)
                  {
                    rhoSecTheta = std::sqrt(rho2*(1+tanSTheta2));
                   
                    *n = G4ThreeVector( p.x()/rhoSecTheta,   
                                        p.y()/rhoSecTheta,
                                        std::sin(fSTheta)  );
                  }
                  else *n = G4ThreeVector(0.,0.,1.);
                }
                return snxt = 0.;               
              }
            }
            stheta    = -0.5*dist2STheta/t2;
            sidetheta = kSTheta;
          }  
        }      // 2nd order equation, 1st root of fSTheta cone,
        else   // 2nd if 1st root -ve
        {
          if( std::fabs(distTheta) < halfRmaxTolerance )
          {
            if( (fSTheta > halfpi) && (t2 >= 0.) ) // leave
            {
              if( calcNorm ) 
              {
                *validNorm = true;
                if (rho2)
                {
                  rhoSecTheta = std::sqrt(rho2*(1+tanSTheta2));
                   
                  *n = G4ThreeVector( p.x()/rhoSecTheta,   
                                      p.y()/rhoSecTheta,
                                      std::sin(fSTheta)  );
                }
                else  { *n = G4ThreeVector(0.,0.,1.); }
              }
              return snxt = 0.;
            }
            else if( (fSTheta < halfpi) && (t2 < 0.) && (p.z() >=0.) ) // leave
            {
              if( calcNorm )  { *validNorm = false; }
              return snxt = 0.;
            }                               
          }
          b  = t2/t1;
          c  = dist2STheta/t1;
          d2 = b*b - c ;

          if ( d2 >= 0. )
          {
            d = std::sqrt(d2);

            if( fSTheta > halfpi )
            {
              sd = -b - d;         // First root

              if ( ((std::fabs(s) < halfRmaxTolerance) && (t2 < 0.))
               ||  (sd < 0.)  || ( (sd > 0.) && (p.z() + sd*v.z() > 0.) )     ) 
              {
                sd = -b + d ; // 2nd root
              }
              if( (sd > halfRmaxTolerance) && (p.z() + sd*v.z() <= 0.) )  
              {
                stheta    = sd;
                sidetheta = kSTheta;
              }
            }
            else // sTheta < pi/2, concave surface, no normal
            {
              sd = -b - d;         // First root

              if ( ( (std::fabs(sd) < halfRmaxTolerance) && (t2 >= 0.) )
                || (sd < 0.) || ( (sd > 0.) && (p.z() + sd*v.z() < 0.) )   )
              {
                sd = -b + d ; // 2nd root
              }
              if( (sd > halfRmaxTolerance) && (p.z() + sd*v.z() >= 0.) )  
              {
                stheta    = sd;
                sidetheta = kSTheta;
              }            
            }
          }
        }
      }
    }
    if (eTheta < pi) // intersection with second cons
    {
      if( std::fabs(tanETheta) > 5./kAngTolerance ) // kons is plane z=0
      {
        if( v.z() < 0. ) 
        {
          if ( std::fabs( p.z() ) <= halfRmaxTolerance )
          {
            if(calcNorm)
            {
              *validNorm = true;
              *n = G4ThreeVector(0.,0.,-1.);
            }
            return snxt = 0 ;
          }  
          sd = -p.z()/v.z();

          if( sd < stheta )
          {
            stheta    = sd;
            sidetheta = kETheta;
          }
        }
      }
      else // kons is not plane 
      {
        t1          = 1-v.z()*v.z()*(1+tanETheta2);
        t2          = pDotV2d-p.z()*v.z()*tanETheta2;  // ~vDotN if p on cons
        dist2ETheta = rho2-p.z()*p.z()*tanETheta2;     // t3

        distTheta = std::sqrt(rho2)-p.z()*tanETheta;

        if( std::fabs(t1) < halfAngTolerance ) // 1st order equation,
        {                                      // v parallel to kons
          if( v.z() < 0. )
          {
            if(std::fabs(distTheta) < halfRmaxTolerance) // p on surface
            {
              if( (eTheta > halfpi) && (p.z() < 0.) )
              {
                if( calcNorm )  { *validNorm = false; }
                return snxt = 0.;
              }
              else if ( (eTheta < halfpi) && (p.z() >= 0) )
              {
                if( calcNorm ) 
                {
                  *validNorm = true;
                  if (rho2)
                  {
                    rhoSecTheta = std::sqrt(rho2*(1+tanETheta2));
                    *n = G4ThreeVector( p.x()/rhoSecTheta,   
                                        p.y()/rhoSecTheta,
                                        -sinETheta  );
                  }
                  else  { *n = G4ThreeVector(0.,0.,-1.); }
                }
                return snxt = 0.;               
              }
            }
            sd = -0.5*dist2ETheta/t2;

            if( sd < stheta )
            {
              stheta    = sd;
              sidetheta = kETheta;
            }
          }  
        }      // 2nd order equation, 1st root of fSTheta cone
        else   // 2nd if 1st root -ve
        {
          if ( std::fabs(distTheta) < halfRmaxTolerance )
          {
            if( (eTheta < halfpi) && (t2 >= 0.) ) // leave
            {
              if( calcNorm ) 
              {
                *validNorm = true;
                if (rho2)
                {
                    rhoSecTheta = std::sqrt(rho2*(1+tanETheta2));
                    *n = G4ThreeVector( p.x()/rhoSecTheta,   
                                        p.y()/rhoSecTheta,
                                        -sinETheta  );
                }
                else *n = G4ThreeVector(0.,0.,-1.);
              }                           
              return snxt = 0.;
            }
            else if ( (eTheta > halfpi)
                   && (t2 < 0.) && (p.z() <=0.) ) // leave
            {
              if( calcNorm )  { *validNorm = false; }
              return snxt = 0.;
            }                               
          }
          b  = t2/t1;
          c  = dist2ETheta/t1;
          d2 = b*b - c ;
          if ( (d2 <halfRmaxTolerance) && (d2 > -halfRmaxTolerance) )
          {
            d2 = 0.;
          }
          if ( d2 >= 0. )
          {
            d = std::sqrt(d2);

            if( eTheta < halfpi )
            {
              sd = -b - d;         // First root

              if( ((std::fabs(sd) < halfRmaxTolerance) && (t2 < 0.))
               || (sd < 0.) ) 
              {
                sd = -b + d ; // 2nd root
              }
              if( sd > halfRmaxTolerance )  
              {
                if( sd < stheta )
                {
                  stheta    = sd;
                  sidetheta = kETheta;
                }
              }
            }
            else // sTheta+fDTheta > pi/2, concave surface, no normal
            {
              sd = -b - d;         // First root

              if ( ((std::fabs(sd) < halfRmaxTolerance) && (t2 >= 0.))
                || (sd < 0.) || ( (sd > 0.) && (p.z() + sd*v.z() > halfRmaxTolerance) ) )
              {
                sd = -b + d ; // 2nd root
              }
              if( (sd > halfRmaxTolerance) && (p.z() + sd*v.z() <= halfRmaxTolerance) )  
              {
                if( sd < stheta )
                {
                  stheta    = sd;
                  sidetheta = kETheta;
                }
              }            
            }
          }
        }
      }
    }

  } // end theta intersections

  // Phi Intersection
    
  if ( !fFullPhiSphere )
  {
    if ( p.x() || p.y() ) // Check if on z axis (rho not needed later)
    {
      // pDist -ve when inside

      pDistS=p.x()*sinSPhi-p.y()*cosSPhi;
      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 ( (pDistS <= 0) && (pDistE <= 0) )
      {
        // Inside both phi *full* planes

        if ( compS < 0 )
        {
          sphi = pDistS/compS ;
          xi   = p.x()+sphi*v.x() ;
          yi   = p.y()+sphi*v.y() ;

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

        if ( compE < 0 )
        {
          sphi2=pDistE/compE ;
          if (sphi2 < sphi) // Only check further if < starting phi intersection
          {
            xi = p.x()+sphi2*v.x() ;
            yi = p.y()+sphi2*v.y() ;

            // Check intersection with correct half-plane
            //
            if ((std::fabs(xi)<=kCarTolerance) && (std::fabs(yi)<=kCarTolerance))
            {
              // Leaving via ending phi
              //
              vphi = std::atan2(v.y(),v.x()) ;
               
              if( !((fSPhi-halfAngTolerance <= vphi)
                  &&(fSPhi+fDPhi+halfAngTolerance >= vphi)) )
              { 
                sidephi = kEPhi;
                if ( pDistE <= -halfCarTolerance )  { sphi = sphi2; }
                else                                { sphi = 0.0;   }
              }
            }
            else if ((yi*cosCPhi-xi*sinCPhi)>=0) // Leaving via ending phi
            {
              sidephi = kEPhi ;
              if ( pDistE <= -halfCarTolerance )
              {
                sphi=sphi2;
              }
              else 
              {
                sphi = 0 ;
              }
            }
          }
        }        
      }
      else if ((pDistS >= 0) && (pDistE >= 0)) // Outside both *full* phi planes
      {
        if ( pDistS <= pDistE )
        {
          sidephi = kSPhi ;
        }
        else
        {
          sidephi = kEPhi ;
        }
        if ( fDPhi > pi )
        {
          if ( (compS < 0) && (compE < 0) )  { sphi = 0; }
          else                               { sphi = kInfinity; }
        }
        else
        {
          // if towards both >=0 then once inside (after error)
          // will remain inside

          if ( (compS >= 0) && (compE >= 0) ) { sphi = kInfinity; }
          else                                { sphi = 0; }
        }    
      }
      else if ( (pDistS > 0) && (pDistE < 0) )
      {
        // Outside full starting plane, inside full ending plane

        if ( fDPhi > pi )
        {
          if ( compE < 0 )
          {
            sphi = pDistE/compE ;
            xi   = p.x() + sphi*v.x() ;
            yi   = p.y() + sphi*v.y() ;

            // Check intersection in correct half-plane
            // (if not -> not leaving phi extent)
            //
            if( (std::fabs(xi)<=kCarTolerance)&&(std::fabs(yi)<=kCarTolerance) )
            {
              vphi = std::atan2(v.y(),v.x());
              sidephi = kSPhi;
              if ( ( (fSPhi-halfAngTolerance) <= vphi)
                && ( (ePhi+halfAngTolerance)  >= vphi) )
              {
                sphi = kInfinity;
              }
            }
            else if ( ( yi*cosCPhi - xi*sinCPhi ) <= 0 )
            {
              sphi = kInfinity ;
            }
            else // Leaving via Ending phi
            {
              sidephi = kEPhi ;
              if ( pDistE > -halfCarTolerance )  { sphi = 0.; }
            }
          }
          else
          {
            sphi = kInfinity ;
          }
        }
        else
        {
          if ( compS >= 0 )
          {
            if ( compE < 0 )
            {            
              sphi = pDistE/compE ;
              xi   = p.x() + sphi*v.x() ;
              yi   = p.y() + sphi*v.y() ;

              // Check intersection in correct half-plane
              // (if not -> remain in extent)
              //
              if( (std::fabs(xi)<=kCarTolerance)&&(std::fabs(yi)<=kCarTolerance) )
              {
                vphi = std::atan2(v.y(),v.x());
                sidephi = kSPhi;
                if ( ( (fSPhi-halfAngTolerance) <= vphi)
                  && ( (ePhi+halfAngTolerance)  >= vphi) )
                {
                  sphi = kInfinity;
                }
              }
              else if ( ( yi*cosCPhi - xi*sinCPhi) <= 0 )
              {
                sphi=kInfinity;
              }
              else // otherwise leaving via Ending phi
              {
                sidephi = kEPhi ;
              }
            }
            else sphi=kInfinity;
          }
          else // leaving immediately by starting phi
          {
            sidephi = kSPhi ;
            sphi    = 0 ;
          }
        }
      }
      else
      {
        // Must be pDistS < 0 && pDistE > 0
        // Inside full starting plane, outside full ending plane

        if ( fDPhi > pi )
        {
          if ( compS < 0 )
          {
            sphi=pDistS/compS;
            xi=p.x()+sphi*v.x();
            yi=p.y()+sphi*v.y();
            
            // Check intersection in correct half-plane
            // (if not -> not leaving phi extent)
            //
            if( (std::fabs(xi)<=kCarTolerance)&&(std::fabs(yi)<=kCarTolerance) )
            {
              vphi = std::atan2(v.y(),v.x()) ;
              sidephi = kSPhi;
              if ( ( (fSPhi-halfAngTolerance) <= vphi)
                && ( (ePhi+halfAngTolerance)  >= vphi) )
              {
              sphi = kInfinity;
              }
            }
            else if ( ( yi*cosCPhi - xi*sinCPhi ) >= 0 )
            {
              sphi = kInfinity ;
            }
            else  // Leaving via Starting phi
            {
              sidephi = kSPhi ;
              if ( pDistS > -halfCarTolerance )  { sphi = 0; }
            }
          }
          else
          {
            sphi = kInfinity ;
          }
        }
        else
        {
          if ( compE >= 0 )
          {
            if ( compS < 0 )
            {
              sphi = pDistS/compS ;
              xi   = p.x()+sphi*v.x() ;
              yi   = p.y()+sphi*v.y() ;
              
              // Check intersection in correct half-plane
              // (if not -> remain in extent)
              //
              if((std::fabs(xi)<=kCarTolerance) && (std::fabs(yi)<=kCarTolerance))
              {
                vphi = std::atan2(v.y(),v.x()) ;
                sidephi = kSPhi;
                if ( ( (fSPhi-halfAngTolerance) <= vphi)
                  && ( (ePhi+halfAngTolerance)  >= vphi) )
                {
                  sphi = kInfinity;
                }
              }
              else if ( ( yi*cosCPhi - xi*sinCPhi ) >= 0 )
              {
                sphi = kInfinity ;
              }
              else // otherwise leaving via Starting phi
              {
                sidephi = kSPhi ;
              }
            }
            else
            {
              sphi = kInfinity ;
            }
          }
          else // leaving immediately by ending
          {
            sidephi = kEPhi ;
            sphi    = 0     ;
          }
        }
      }      
    }
    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

      if ( v.x() || v.y() )
      {
        vphi = std::atan2(v.y(),v.x()) ;
        if ((fSPhi-halfAngTolerance < vphi) && (vphi < ePhi+halfAngTolerance))
        {
          sphi = kInfinity;
        }
        else
        {
          sidephi = kSPhi ; // arbitrary 
          sphi    = 0     ;
        }
      }
      else  // travel along z - no phi intersection
      {
        sphi = kInfinity ;
      }
    }
    if ( sphi < snxt )  // Order intersecttions
    {
      snxt = sphi ;
      side = sidephi ;
    }
  }
  if (stheta < snxt ) // Order intersections
  {
    snxt = stheta ;
    side = sidetheta ;
  }

  if (calcNorm)    // Output switch operator
  {
    switch( side )
    {
      case kRMax:
        xi=p.x()+snxt*v.x();
        yi=p.y()+snxt*v.y();
        zi=p.z()+snxt*v.z();
        *n=G4ThreeVector(xi/fRmax,yi/fRmax,zi/fRmax);
        *validNorm=true;
        break;

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

      case kSPhi:
        if ( fDPhi <= pi )     // Normal to Phi-
        {
          *n=G4ThreeVector(sinSPhi,-cosSPhi,0);
          *validNorm=true;
        }
        else  { *validNorm=false; }
        break ;

      case kEPhi:
        if ( fDPhi <= pi )      // Normal to Phi+
        {
          *n=G4ThreeVector(-sinEPhi,cosEPhi,0);
          *validNorm=true;
        }
        else  { *validNorm=false; }
        break;

      case kSTheta:
        if( fSTheta == halfpi )
        {
          *n=G4ThreeVector(0.,0.,1.);
          *validNorm=true;
        }
        else if ( fSTheta > halfpi )
        {
          xi = p.x() + snxt*v.x();
          yi = p.y() + snxt*v.y();
          rho2=xi*xi+yi*yi;
          if (rho2)
          { 
            rhoSecTheta = std::sqrt(rho2*(1+tanSTheta2));
            *n = G4ThreeVector( xi/rhoSecTheta, yi/rhoSecTheta,
                               -tanSTheta/std::sqrt(1+tanSTheta2));
          }
          else
          {
            *n = G4ThreeVector(0.,0.,1.);
          }
          *validNorm=true;
        }
        else  { *validNorm=false; }  // Concave STheta cone
        break;

      case kETheta:
        if( eTheta == halfpi )
        {
          *n         = G4ThreeVector(0.,0.,-1.);
          *validNorm = true;
        }
        else if ( eTheta < halfpi )
        {
          xi=p.x()+snxt*v.x();
          yi=p.y()+snxt*v.y();
          rho2=xi*xi+yi*yi;
          if (rho2)
          { 
            rhoSecTheta = std::sqrt(rho2*(1+tanETheta2));
            *n = G4ThreeVector( xi/rhoSecTheta, yi/rhoSecTheta,
                               -tanETheta/std::sqrt(1+tanETheta2) );
          }
          else
          {
            *n = G4ThreeVector(0.,0.,-1.);
          }
          *validNorm=true;
        }
        else  { *validNorm=false; }   // Concave ETheta cone
        break;

      default:
        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("G4Sphere::DistanceToOut(p,v,..)",
                    "GeomSolids1002", JustWarning, message);
        break;
    }
  }
  if (snxt == kInfinity)
  {
    G4cout << G4endl;
    DumpInfo();
    std::ostringstream message;
    G4int oldprc = message.precision(16);
    message << "Logic error: snxt = kInfinity  ???" << 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
            << "Rp = "<< std::sqrt( p.x()*p.x()+p.y()*p.y()+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("G4Sphere::DistanceToOut(p,v,..)",
                "GeomSolids1002", JustWarning, message);
  }

  return snxt;
}

Here is the call graph for this function:

DistanceToOut() [2/2]

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

virtual

Implements G4VSolid.

Definition at line 2817 of file G4Sphere.cc.

{
  G4double safe=0.0,safeRMin,safeRMax,safePhi,safeTheta;
  G4double rho2,rds,rho;
  G4double pTheta,dTheta1 = kInfinity,dTheta2 = kInfinity;
  rho2=p.x()*p.x()+p.y()*p.y();
  rds=std::sqrt(rho2+p.z()*p.z());
  rho=std::sqrt(rho2);

#ifdef G4CSGDEBUG
  if( Inside(p) == kOutside )
  {
     G4int old_prc = 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(old_prc) ;
     G4Exception("G4Sphere::DistanceToOut(p)",
                 "GeomSolids1002", JustWarning, "Point p is outside !?" );
  }
#endif

  // Distance to r shells
  //
  safeRMax = fRmax-rds;
  safe = safeRMax;  
  if (fRmin)
  {
     safeRMin = rds-fRmin;
     safe = std::min( safeRMin, safeRMax ); 
  }

  // Distance to phi extent
  //
  if ( !fFullPhiSphere )
  {
     if (rho>0.0)
     {
        if ((p.y()*cosCPhi-p.x()*sinCPhi)<=0)
        {
           safePhi=-(p.x()*sinSPhi-p.y()*cosSPhi);
        }
        else
        {
           safePhi=(p.x()*sinEPhi-p.y()*cosEPhi);
        }
     }
     else
     {
        safePhi = 0.0;  // Distance to both Phi surfaces (extended)
     }
     // Both cases above can be improved - in case fRMin > 0.0
     //  although it may be costlier (good for precise, not fast version)
     
     safe= std::min(safe, safePhi);
  }

  // Distance to Theta extent
  //
  if ( !fFullThetaSphere )
  {
    if( rds > 0.0 )
    {
       pTheta=std::acos(p.z()/rds);
       if (pTheta<0) { pTheta+=pi; }
       if(fSTheta>0.)
       { dTheta1=pTheta-fSTheta;}
       if(eTheta<pi)
       { dTheta2=eTheta-pTheta;}
      
       safeTheta=rds*std::sin(std::min(dTheta1, dTheta2) );
    }
    else
    {
       safeTheta= 0.0;
         // An improvement will be to return negative answer if outside (TODO)
    }
    safe = std::min( safe, safeTheta );
  }

  if (safe<0.0) { safe=0; }
    // An improvement to return negative answer if outside (TODO)
  
  return safe;
}

Here is the call graph for this function:

GetCubicVolume()

G4double G4Sphere::GetCubicVolume ( )

inlinevirtual

Reimplemented from G4VSolid.

GetDeltaPhiAngle()

G4double G4Sphere::GetDeltaPhiAngle ( ) const

inline

Here is the caller graph for this function:

GetDeltaThetaAngle()

G4double G4Sphere::GetDeltaThetaAngle ( ) const

inline

Here is the caller graph for this function:

GetDPhi()

G4double G4Sphere::GetDPhi ( ) const

inline

GetDTheta()

G4double G4Sphere::GetDTheta ( ) const

inline

GetEntityType()

G4GeometryType G4Sphere::GetEntityType ( ) const

virtual

Implements G4VSolid.

Definition at line 3045 of file G4Sphere.cc.

{
  return G4String("G4Sphere");
}

GetExtent()

G4VisExtent G4Sphere::GetExtent ( ) const

virtual

Reimplemented from G4VSolid.

Definition at line 3216 of file G4Sphere.cc.

{
  return G4VisExtent(-fRmax, fRmax,-fRmax, fRmax,-fRmax, fRmax );
}

GetInnerRadius()

G4double G4Sphere::GetInnerRadius ( ) const

inline

Here is the caller graph for this function:

GetInsideRadius()

G4double G4Sphere::GetInsideRadius ( ) const

inline

Here is the caller graph for this function:

GetOuterRadius()

G4double G4Sphere::GetOuterRadius ( ) const

inline

Here is the caller graph for this function:

GetPointOnSurface()

G4ThreeVector G4Sphere::GetPointOnSurface ( ) const

virtual

Reimplemented from G4VSolid.

Definition at line 3087 of file G4Sphere.cc.

{
  G4double zRand, aOne, aTwo, aThr, aFou, aFiv, chose, phi, sinphi, cosphi;
  G4double height1, height2, slant1, slant2, costheta, sintheta, rRand;

  height1 = (fRmax-fRmin)*cosSTheta;
  height2 = (fRmax-fRmin)*cosETheta;
  slant1  = std::sqrt(sqr((fRmax - fRmin)*sinSTheta) + height1*height1);
  slant2  = std::sqrt(sqr((fRmax - fRmin)*sinETheta) + height2*height2);
  rRand   = GetRadiusInRing(fRmin,fRmax);
  
  aOne = fRmax*fRmax*fDPhi*(cosSTheta-cosETheta);
  aTwo = fRmin*fRmin*fDPhi*(cosSTheta-cosETheta);
  aThr = fDPhi*((fRmax + fRmin)*sinSTheta)*slant1;
  aFou = fDPhi*((fRmax + fRmin)*sinETheta)*slant2;
  aFiv = 0.5*fDTheta*(fRmax*fRmax-fRmin*fRmin);
  
  phi = G4RandFlat::shoot(fSPhi, ePhi); 
  cosphi = std::cos(phi); 
  sinphi = std::sin(phi);
  costheta = G4RandFlat::shoot(cosETheta,cosSTheta);
  sintheta = std::sqrt(1.-sqr(costheta));

  if(fFullPhiSphere) { aFiv = 0; }
  if(fSTheta == 0)   { aThr=0; }
  if(eTheta == pi) { aFou = 0; }
  if(fSTheta == halfpi) { aThr = pi*(fRmax*fRmax-fRmin*fRmin); }
  if(eTheta == halfpi)  { aFou = pi*(fRmax*fRmax-fRmin*fRmin); }

  chose = G4RandFlat::shoot(0.,aOne+aTwo+aThr+aFou+2.*aFiv);
  if( (chose>=0.) && (chose<aOne) )
  {
    return G4ThreeVector(fRmax*sintheta*cosphi,
                         fRmax*sintheta*sinphi, fRmax*costheta);
  }
  else if( (chose>=aOne) && (chose<aOne+aTwo) )
  {
    return G4ThreeVector(fRmin*sintheta*cosphi,
                         fRmin*sintheta*sinphi, fRmin*costheta);
  }
  else if( (chose>=aOne+aTwo) && (chose<aOne+aTwo+aThr) )
  {
    if (fSTheta != halfpi)
    {
      zRand = G4RandFlat::shoot(fRmin*cosSTheta,fRmax*cosSTheta);
      return G4ThreeVector(tanSTheta*zRand*cosphi,
                           tanSTheta*zRand*sinphi,zRand);
    }
    else
    {
      return G4ThreeVector(rRand*cosphi, rRand*sinphi, 0.);
    }    
  }
  else if( (chose>=aOne+aTwo+aThr) && (chose<aOne+aTwo+aThr+aFou) )
  {
    if(eTheta != halfpi)
    {
      zRand = G4RandFlat::shoot(fRmin*cosETheta, fRmax*cosETheta);
      return G4ThreeVector  (tanETheta*zRand*cosphi,
                             tanETheta*zRand*sinphi,zRand);
    }
    else
    {
      return G4ThreeVector(rRand*cosphi, rRand*sinphi, 0.);
    }
  }
  else if( (chose>=aOne+aTwo+aThr+aFou) && (chose<aOne+aTwo+aThr+aFou+aFiv) )
  {
    return G4ThreeVector(rRand*sintheta*cosSPhi,
                         rRand*sintheta*sinSPhi,rRand*costheta);
  }
  else
  {
    return G4ThreeVector(rRand*sintheta*cosEPhi,
                         rRand*sintheta*sinEPhi,rRand*costheta);
  }
}

Here is the call graph for this function:

GetRmax()

G4double G4Sphere::GetRmax ( ) const

inline

GetRmin()

G4double G4Sphere::GetRmin ( ) const

inline

GetSPhi()

G4double G4Sphere::GetSPhi ( ) const

inline

GetStartPhiAngle()

G4double G4Sphere::GetStartPhiAngle ( ) const

inline

Here is the caller graph for this function:

GetStartThetaAngle()

G4double G4Sphere::GetStartThetaAngle ( ) const

inline

Here is the caller graph for this function:

GetSTheta()

G4double G4Sphere::GetSTheta ( ) const

inline

GetSurfaceArea()

G4double G4Sphere::GetSurfaceArea ( )

virtual

Reimplemented from G4VSolid.

Definition at line 3169 of file G4Sphere.cc.

{
  if(fSurfaceArea != 0.) {;}
  else
  {   
    G4double Rsq=fRmax*fRmax;
    G4double rsq=fRmin*fRmin;
         
    fSurfaceArea = fDPhi*(rsq+Rsq)*(cosSTheta - cosETheta);
    if(!fFullPhiSphere)
    {
      fSurfaceArea = fSurfaceArea + fDTheta*(Rsq-rsq);
    }
    if(fSTheta >0)
    {
      G4double acos1=std::acos( std::pow(sinSTheta,2) * std::cos(fDPhi)
                              + std::pow(cosSTheta,2));
      if(fDPhi>pi)
      { 
        fSurfaceArea = fSurfaceArea + 0.5*(Rsq-rsq)*(twopi-acos1);
      }
      else
      {
        fSurfaceArea = fSurfaceArea + 0.5*(Rsq-rsq)*acos1;
      }
    }
    if(eTheta < pi)
    {
      G4double acos2=std::acos( std::pow(sinETheta,2) * std::cos(fDPhi)
                              + std::pow(cosETheta,2));
      if(fDPhi>pi)
      { 
        fSurfaceArea = fSurfaceArea + 0.5*(Rsq-rsq)*(twopi-acos2);
      }
      else
      {
        fSurfaceArea = fSurfaceArea + 0.5*(Rsq-rsq)*acos2;
      }
    }
  }
  return fSurfaceArea;
}

Initialize()

void G4Sphere::Initialize ( )

inlineprivate

InitializePhiTrigonometry()

void G4Sphere::InitializePhiTrigonometry ( )

inlineprivate

InitializeThetaTrigonometry()

void G4Sphere::InitializeThetaTrigonometry ( )

inlineprivate

Inside()

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

virtual

Implements G4VSolid.

Definition at line 472 of file G4Sphere.cc.

{
  G4double rho,rho2,rad2,tolRMin,tolRMax;
  G4double pPhi,pTheta;
  EInside in = kOutside;

  const G4double halfRmaxTolerance = fRmaxTolerance*0.5;
  const G4double halfRminTolerance = fRminTolerance*0.5;
  const G4double Rmax_minus = fRmax - halfRmaxTolerance;
  const G4double Rmin_plus  = (fRmin > 0) ? fRmin+halfRminTolerance : 0;

  rho2 = p.x()*p.x() + p.y()*p.y() ;
  rad2 = rho2 + p.z()*p.z() ;

  // Check radial surfaces. Sets 'in'

  tolRMin = Rmin_plus;
  tolRMax = Rmax_minus;

  if(rad2 == 0.0)
  { 
    if (fRmin > 0.0)
    {
      return in = kOutside;
    }
    if ( (!fFullPhiSphere) || (!fFullThetaSphere) )
    {
      return in = kSurface;
    }
    else
    {
      return in = kInside; 
    }
  }

  if ( (rad2 <= Rmax_minus*Rmax_minus) && (rad2 >= Rmin_plus*Rmin_plus) )
  {
    in = kInside;
  }
  else
  {
    tolRMax = fRmax + halfRmaxTolerance;                  // outside case
    tolRMin = std::max(fRmin-halfRminTolerance, 0.);      // outside case
    if ( (rad2 <= tolRMax*tolRMax) && (rad2 >= tolRMin*tolRMin) )
    {
      in = kSurface;
    }
    else
    {
      return in = kOutside;
    }
  }

  // Phi boundaries   : Do not check if it has no phi boundary!

  if ( !fFullPhiSphere && rho2 )  // [fDPhi < twopi] and [p.x or p.y]
  {
    pPhi = std::atan2(p.y(),p.x()) ;

    if      ( pPhi < fSPhi - halfAngTolerance  ) { pPhi += twopi; }
    else if ( pPhi > ePhi + halfAngTolerance )   { pPhi -= twopi; }
    
    if ( (pPhi < fSPhi - halfAngTolerance)
      || (pPhi > ePhi + halfAngTolerance) )      { return in = kOutside; }
    
    else if (in == kInside)  // else it's kSurface anyway already
    {
      if ( (pPhi < fSPhi + halfAngTolerance)
        || (pPhi > ePhi - halfAngTolerance) )    { in = kSurface; }     
    }
  }

  // Theta bondaries
  
  if ( (rho2 || p.z()) && (!fFullThetaSphere) )
  {
    rho    = std::sqrt(rho2);
    pTheta = std::atan2(rho,p.z());

    if ( in == kInside )
    {
      if ( ((fSTheta > 0.0) && (pTheta < fSTheta + halfAngTolerance))
    || ((eTheta < pi) && (pTheta > eTheta - halfAngTolerance)) )
      {
    if ( (( (fSTheta>0.0)&&(pTheta>=fSTheta-halfAngTolerance) )
             || (fSTheta == 0.0) )
          && ((eTheta==pi)||(pTheta <= eTheta + halfAngTolerance) ) )
        {
          in = kSurface;
        }
        else
        {
          in = kOutside;
        }
      }
    }
    else
    {
        if ( ((fSTheta > 0.0)&&(pTheta < fSTheta - halfAngTolerance))
       ||((eTheta < pi  )&&(pTheta > eTheta + halfAngTolerance)) )
      {
        in = kOutside;
      }
    }
  }
  return in;
}

Here is the call graph for this function:

Here is the caller graph for this function:

operator=()

G4Sphere & G4Sphere::operator=

(

const G4Sphere &

rhs

)

Definition at line 178 of file G4Sphere.cc.

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

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

   // Copy data
   //
   fRminTolerance = rhs.fRminTolerance; fRmaxTolerance = rhs.fRmaxTolerance;
   kAngTolerance = rhs.kAngTolerance; kRadTolerance = rhs.kRadTolerance;
   fEpsilon = rhs.fEpsilon; fRmin = rhs.fRmin; fRmax = rhs.fRmax;
   fSPhi = rhs.fSPhi; fDPhi = rhs.fDPhi; fSTheta = rhs.fSTheta;
   fDTheta = rhs.fDTheta; sinCPhi = rhs.sinCPhi; cosCPhi = rhs.cosCPhi;
   cosHDPhiOT = rhs.cosHDPhiOT; cosHDPhiIT = rhs.cosHDPhiIT;
   sinSPhi = rhs.sinSPhi; cosSPhi = rhs.cosSPhi;
   sinEPhi = rhs.sinEPhi; cosEPhi = rhs.cosEPhi;
   hDPhi = rhs.hDPhi; cPhi = rhs.cPhi; ePhi = rhs.ePhi;
   sinSTheta = rhs.sinSTheta; cosSTheta = rhs.cosSTheta;
   sinETheta = rhs.sinETheta; cosETheta = rhs.cosETheta;
   tanSTheta = rhs.tanSTheta; tanSTheta2 = rhs.tanSTheta2;
   tanETheta = rhs.tanETheta; tanETheta2 = rhs.tanETheta2;
   eTheta = rhs.eTheta; fFullPhiSphere = rhs.fFullPhiSphere;
   fFullThetaSphere = rhs.fFullThetaSphere; fFullSphere = rhs.fFullSphere;
   halfCarTolerance = rhs.halfCarTolerance;
   halfAngTolerance = rhs.halfAngTolerance;

   return *this;
}

Here is the call graph for this function:

SetDeltaPhiAngle()

void G4Sphere::SetDeltaPhiAngle ( G4double  newDphi )

inline

Here is the caller graph for this function:

SetDeltaThetaAngle()

void G4Sphere::SetDeltaThetaAngle ( G4double  newDTheta )

inline

Here is the caller graph for this function:

SetInnerRadius()

void G4Sphere::SetInnerRadius ( G4double  newRMin )

inline

Here is the caller graph for this function:

SetInsideRadius()

void G4Sphere::SetInsideRadius ( G4double  newRmin )

inline

Here is the caller graph for this function:

SetOuterRadius()

void G4Sphere::SetOuterRadius ( G4double  newRmax )

inline

Here is the caller graph for this function:

SetStartPhiAngle()

void G4Sphere::SetStartPhiAngle ( G4double  newSphi, G4bool  trig = true  )

inline

Here is the caller graph for this function:

SetStartThetaAngle()

void G4Sphere::SetStartThetaAngle ( G4double  newSTheta )

inline

Here is the caller graph for this function:

StreamInfo()

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

virtual

Reimplemented from G4CSGSolid.

Definition at line 3063 of file G4Sphere.cc.

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

  return os;
}

Here is the call graph for this function:

SurfaceNormal()

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

virtual

Implements G4VSolid.

Definition at line 586 of file G4Sphere.cc.

{
  G4int noSurfaces = 0;  
  G4double rho, rho2, radius, pTheta, pPhi=0.;
  G4double distRMin = kInfinity;
  G4double distSPhi = kInfinity, distEPhi = kInfinity;
  G4double distSTheta = kInfinity, distETheta = kInfinity;
  G4ThreeVector nR, nPs, nPe, nTs, nTe, nZ(0.,0.,1.);
  G4ThreeVector norm, sumnorm(0.,0.,0.);

  rho2 = p.x()*p.x()+p.y()*p.y();
  radius = std::sqrt(rho2+p.z()*p.z());
  rho  = std::sqrt(rho2);

  G4double    distRMax = std::fabs(radius-fRmax);
  if (fRmin)  distRMin = std::fabs(radius-fRmin);
    
  if ( rho && !fFullSphere )
  {
    pPhi = std::atan2(p.y(),p.x());

    if (pPhi < fSPhi-halfAngTolerance)     { pPhi += twopi; }
    else if (pPhi > ePhi+halfAngTolerance) { pPhi -= twopi; }
  }
  if ( !fFullPhiSphere )
  {
    if ( rho )
    {
      distSPhi = std::fabs( pPhi-fSPhi ); 
      distEPhi = std::fabs( pPhi-ePhi ); 
    }
    else if( !fRmin )
    {
      distSPhi = 0.; 
      distEPhi = 0.; 
    }
    nPs = G4ThreeVector(sinSPhi,-cosSPhi,0);
    nPe = G4ThreeVector(-sinEPhi,cosEPhi,0);
  }        
  if ( !fFullThetaSphere )
  {
    if ( rho )
    {
      pTheta     = std::atan2(rho,p.z());
      distSTheta = std::fabs(pTheta-fSTheta); 
      distETheta = std::fabs(pTheta-eTheta);
 
      nTs = G4ThreeVector(-cosSTheta*p.x()/rho,
                          -cosSTheta*p.y()/rho,
                           sinSTheta          );

      nTe = G4ThreeVector( cosETheta*p.x()/rho,
                           cosETheta*p.y()/rho,
                          -sinETheta          );    
    }
    else if( !fRmin )
    {
      if ( fSTheta )  
      {              
        distSTheta = 0.;
        nTs = G4ThreeVector(0.,0.,-1.);
      }
      if ( eTheta < pi )
      {              
        distETheta = 0.;
        nTe = G4ThreeVector(0.,0.,1.);
      }
    }    
  }
  if( radius )  { nR = G4ThreeVector(p.x()/radius,p.y()/radius,p.z()/radius); }

  if( distRMax <= halfCarTolerance )
  {
    noSurfaces ++;
    sumnorm += nR;
  }
  if( fRmin && (distRMin <= halfCarTolerance) )
  {
    noSurfaces ++;
    sumnorm -= nR;
  }
  if( !fFullPhiSphere )   
  {
    if (distSPhi <= halfAngTolerance)
    {
      noSurfaces ++;
      sumnorm += nPs;
    }
    if (distEPhi <= halfAngTolerance) 
    {
      noSurfaces ++;
      sumnorm += nPe;
    }
  }
  if ( !fFullThetaSphere )
  {
    if ((distSTheta <= halfAngTolerance) && (fSTheta > 0.))
    {
      noSurfaces ++;
      if ((radius <= halfCarTolerance) && fFullPhiSphere)  { sumnorm += nZ;  }
      else                                                 { sumnorm += nTs; }
    }
    if ((distETheta <= halfAngTolerance) && (eTheta < pi)) 
    {
      noSurfaces ++;
      if ((radius <= halfCarTolerance) && fFullPhiSphere)  { sumnorm -= nZ;  }
      else                                                 { sumnorm += nTe; }
      if(sumnorm.z() == 0.)  { sumnorm += nZ; }
    }
  }
  if ( noSurfaces == 0 )
  {
#ifdef G4CSGDEBUG
    G4Exception("G4Sphere::SurfaceNormal(p)", "GeomSolids1002",
                JustWarning, "Point p is not on surface !?" ); 
#endif
     norm = ApproxSurfaceNormal(p);
  }
  else if ( noSurfaces == 1 )  { norm = sumnorm; }
  else                         { norm = sumnorm.unit(); }
  return norm;
}

Here is the call graph for this function:

Member Data Documentation

cosCPhi

G4double G4Sphere::cosCPhi

private

Definition at line 238 of file G4Sphere.hh.

cosEPhi

G4double G4Sphere::cosEPhi

private

Definition at line 238 of file G4Sphere.hh.

cosETheta

G4double G4Sphere::cosETheta

private

Definition at line 243 of file G4Sphere.hh.

cosHDPhiIT

G4double G4Sphere::cosHDPhiIT

private

Definition at line 238 of file G4Sphere.hh.

cosHDPhiOT

G4double G4Sphere::cosHDPhiOT

private

Definition at line 238 of file G4Sphere.hh.

cosSPhi

G4double G4Sphere::cosSPhi

private

Definition at line 238 of file G4Sphere.hh.

cosSTheta

G4double G4Sphere::cosSTheta

private

Definition at line 243 of file G4Sphere.hh.

cPhi

G4double G4Sphere::cPhi

private

Definition at line 238 of file G4Sphere.hh.

ePhi

G4double G4Sphere::ePhi

private

Definition at line 238 of file G4Sphere.hh.

eTheta

G4double G4Sphere::eTheta

private

Definition at line 243 of file G4Sphere.hh.

fDPhi

G4double G4Sphere::fDPhi

private

Definition at line 234 of file G4Sphere.hh.

fDTheta

G4double G4Sphere::fDTheta

private

Definition at line 234 of file G4Sphere.hh.

fEpsilon

G4double G4Sphere::fEpsilon

private

Definition at line 229 of file G4Sphere.hh.

fFullPhiSphere

G4bool G4Sphere::fFullPhiSphere

private

Definition at line 248 of file G4Sphere.hh.

fFullSphere

G4bool G4Sphere::fFullSphere

private

Definition at line 248 of file G4Sphere.hh.

fFullThetaSphere

G4bool G4Sphere::fFullThetaSphere

private

Definition at line 248 of file G4Sphere.hh.

fRmax

G4double G4Sphere::fRmax

private

Definition at line 234 of file G4Sphere.hh.

fRmaxTolerance

G4double G4Sphere::fRmaxTolerance

private

Definition at line 229 of file G4Sphere.hh.

fRmin

G4double G4Sphere::fRmin

private

Definition at line 234 of file G4Sphere.hh.

fRminTolerance

G4double G4Sphere::fRminTolerance

private

Definition at line 229 of file G4Sphere.hh.

fSPhi

G4double G4Sphere::fSPhi

private

Definition at line 234 of file G4Sphere.hh.

fSTheta

G4double G4Sphere::fSTheta

private

Definition at line 234 of file G4Sphere.hh.

halfAngTolerance

G4double G4Sphere::halfAngTolerance

private

Definition at line 252 of file G4Sphere.hh.

halfCarTolerance

G4double G4Sphere::halfCarTolerance

private

Definition at line 252 of file G4Sphere.hh.

hDPhi

G4double G4Sphere::hDPhi

private

Definition at line 238 of file G4Sphere.hh.

kAngTolerance

G4double G4Sphere::kAngTolerance

private

Definition at line 229 of file G4Sphere.hh.

kRadTolerance

G4double G4Sphere::kRadTolerance

private

Definition at line 229 of file G4Sphere.hh.

sinCPhi

G4double G4Sphere::sinCPhi

private

Definition at line 238 of file G4Sphere.hh.

sinEPhi

G4double G4Sphere::sinEPhi

private

Definition at line 238 of file G4Sphere.hh.

sinETheta

G4double G4Sphere::sinETheta

private

Definition at line 243 of file G4Sphere.hh.

sinSPhi

G4double G4Sphere::sinSPhi

private

Definition at line 238 of file G4Sphere.hh.

sinSTheta

G4double G4Sphere::sinSTheta

private

Definition at line 243 of file G4Sphere.hh.

tanETheta

G4double G4Sphere::tanETheta

private

Definition at line 243 of file G4Sphere.hh.

tanETheta2

G4double G4Sphere::tanETheta2

private

Definition at line 243 of file G4Sphere.hh.

tanSTheta

G4double G4Sphere::tanSTheta

private

Definition at line 243 of file G4Sphere.hh.

tanSTheta2

G4double G4Sphere::tanSTheta2

private

Definition at line 243 of file G4Sphere.hh.

---

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

• G4Sphere.hh
• G4Sphere.cc