Geant4 11.4.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4Sphere Class Reference

G4Sphere is, in the general case, a section of a spherical shell, between specified phi and theta angles. The phi and theta segments are described by a starting angle and the +ve delta angle for the shape. If the delta angle is >=2*pi, or >=pi the shape is treated as continuous in phi or theta respectively. Theta must lie between [0..pi]. More...

#include <G4Sphere.hh>

Inheritance diagram for G4Sphere:

Public Member Functions

 G4Sphere (const G4String &pName, G4double pRmin, G4double pRmax, G4double pSPhi, G4double pDPhi, G4double pSTheta, G4double pDTheta)
 ~G4Sphere () override=default
G4double GetInnerRadius () const
G4double GetOuterRadius () const
G4double GetStartPhiAngle () const
G4double GetDeltaPhiAngle () const
G4double GetStartThetaAngle () const
G4double GetDeltaThetaAngle () const
G4double GetSinStartPhi () const
G4double GetCosStartPhi () const
G4double GetSinEndPhi () const
G4double GetCosEndPhi () const
G4double GetSinStartTheta () const
G4double GetCosStartTheta () const
G4double GetSinEndTheta () const
G4double GetCosEndTheta () 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 () override
G4double GetSurfaceArea () override
void ComputeDimensions (G4VPVParameterisation *p, const G4int n, const G4VPhysicalVolume *pRep) override
void BoundingLimits (G4ThreeVector &pMin, G4ThreeVector &pMax) const override
G4bool CalculateExtent (const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pmin, G4double &pmax) const override
EInside Inside (const G4ThreeVector &p) const override
G4ThreeVector SurfaceNormal (const G4ThreeVector &p) const override
G4double DistanceToIn (const G4ThreeVector &p, const G4ThreeVector &v) const override
G4double DistanceToIn (const G4ThreeVector &p) const override
G4double DistanceToOut (const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=false, G4bool *validNorm=nullptr, G4ThreeVector *n=nullptr) const override
G4double DistanceToOut (const G4ThreeVector &p) const override
G4GeometryType GetEntityType () const override
G4ThreeVector GetPointOnSurface () const override
G4VSolid * Clone () const override
std::ostream & StreamInfo (std::ostream &os) const override
G4VisExtent GetExtent () const override
void DescribeYourselfTo (G4VGraphicsScene &scene) const override
G4Polyhedron * CreatePolyhedron () const override
 G4Sphere (__void__ &)
 G4Sphere (const G4Sphere &rhs)=default
G4Sphere & operator= (const G4Sphere &rhs)
Public Member Functions inherited from G4CSGSolid
 G4CSGSolid (const G4String &pName)
 ~G4CSGSolid () override
G4Polyhedron * GetPolyhedron () const override
 G4CSGSolid (__void__ &)
 G4CSGSolid (const G4CSGSolid &rhs)
G4CSGSolid & operator= (const G4CSGSolid &rhs)
Public Member Functions inherited from G4VSolid
 G4VSolid (const G4String &name)
virtual ~G4VSolid ()
 G4VSolid (const G4VSolid &rhs)
G4VSolid & operator= (const G4VSolid &rhs)
G4bool operator== (const G4VSolid &s) const
G4String GetName () const
void SetName (const G4String &name)
G4double GetTolerance () const
virtual G4int GetNumOfConstituents () const
virtual G4bool IsFaceted () 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__ &)
G4double EstimateCubicVolume (G4int nStat, G4double epsilon) const
G4double EstimateSurfaceArea (G4int nStat, G4double epsilon) const

Additional Inherited Members

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

Detailed Description

G4Sphere is, in the general case, a section of a spherical shell, between specified phi and theta angles. The phi and theta segments are described by a starting angle and the +ve delta angle for the shape. If the delta angle is >=2*pi, or >=pi the shape is treated as continuous in phi or theta respectively. Theta must lie between [0..pi].

Definition at line 88 of file G4Sphere.hh.

Constructor & Destructor Documentation

◆ G4Sphere() [1/3]

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

Constructs a sphere or sphere shell section with the given name and dimensions.

Parameters
[in]pNameThe name of the solid.
[in]pRminInner radius.
[in]pRmaxOuter radius.
[in]pSPhiStarting Phi angle of the segment in radians.
[in]pDPhiDelta Phi angle of the segment in radians.
[in]pSThetaStarting Theta angle of the segment in radians.
[in]pDThetaDelta Theta angle of the segment in radians.

Definition at line 81 of file G4Sphere.cc.

85 : G4CSGSolid(pName), fSPhi(0.0), fFullPhiSphere(true), fFullThetaSphere(true)
86{
87 kAngTolerance = G4GeometryTolerance::GetInstance()->GetAngularTolerance();
88 kRadTolerance = G4GeometryTolerance::GetInstance()->GetRadialTolerance();
89
90 halfCarTolerance = 0.5*kCarTolerance;
91 halfAngTolerance = 0.5*kAngTolerance;
92
93 // Check radii and set radial tolerances
94
95 if ( (pRmin >= pRmax) || (pRmax < 1.1*kRadTolerance) || (pRmin < 0) )
96 {
97 std::ostringstream message;
98 message << "Invalid radii for Solid: " << GetName() << G4endl
99 << " pRmin = " << pRmin << ", pRmax = " << pRmax;
100 G4Exception("G4Sphere::G4Sphere()", "GeomSolids0002",
101 FatalException, message);
102 }
103 fRmin=pRmin; fRmax=pRmax;
104 fRminTolerance = (fRmin) != 0.0 ? std::max( kRadTolerance, fEpsilon*fRmin ) : 0;
105 fRmaxTolerance = std::max( kRadTolerance, fEpsilon*fRmax );
106
107 // Check angles
108
109 CheckPhiAngles(pSPhi, pDPhi);
110 CheckThetaAngles(pSTheta, pDTheta);
111}
FatalException
@ FatalException
Definition G4ExceptionSeverity.hh:67
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition G4Exception.cc:59
G4endl
#define G4endl
Definition G4ios.hh:67
G4CSGSolid::G4CSGSolid
G4CSGSolid(const G4String &pName)
Definition G4CSGSolid.cc:49
G4GeometryTolerance::GetRadialTolerance
G4double GetRadialTolerance() const
Definition G4GeometryTolerance.cc:82
G4GeometryTolerance::GetInstance
static G4GeometryTolerance * GetInstance()
Definition G4GeometryTolerance.cc:58
G4GeometryTolerance::GetAngularTolerance
G4double GetAngularTolerance() const
Definition G4GeometryTolerance.cc:77
G4VSolid::GetName
G4String GetName() const
G4VSolid::kCarTolerance
G4double kCarTolerance
Definition G4VSolid.hh:418

Referenced by Clone(), G4Sphere(), and operator=().

◆ ~G4Sphere()

G4Sphere::~G4Sphere ( )
overridedefault

Default destructor.

◆ G4Sphere() [2/3]

G4Sphere::G4Sphere ( __void__ & a)

Fake default constructor for usage restricted to direct object persistency for clients requiring preallocation of memory for persistifiable objects.

Definition at line 118 of file G4Sphere.cc.

119 : G4CSGSolid(a)
120{
121}

◆ G4Sphere() [3/3]

G4Sphere::G4Sphere ( const G4Sphere & rhs)
default

Copy constructor and assignment operator.

Member Function Documentation

◆ BoundingLimits()

void G4Sphere::BoundingLimits ( G4ThreeVector & pMin,
G4ThreeVector & pMax ) const
overridevirtual

Computes the bounding limits of the solid.

Parameters
[out]pMinThe minimum bounding limit point.
[out]pMaxThe maximum bounding limit point.

Reimplemented from G4VSolid.

Definition at line 177 of file G4Sphere.cc.

178{
179 G4double rmin = GetInnerRadius();
180 G4double rmax = GetOuterRadius();
181
182 // Find bounding box
183 //
184 if (GetDeltaThetaAngle() >= pi && GetDeltaPhiAngle() >= twopi)
185 {
186 pMin.set(-rmax,-rmax,-rmax);
187 pMax.set( rmax, rmax, rmax);
188 }
189 else
190 {
191 G4double sinStart = GetSinStartTheta();
192 G4double cosStart = GetCosStartTheta();
193 G4double sinEnd = GetSinEndTheta();
194 G4double cosEnd = GetCosEndTheta();
195
196 G4double stheta = GetStartThetaAngle();
197 G4double etheta = stheta + GetDeltaThetaAngle();
198 G4double rhomin = rmin*std::min(sinStart,sinEnd);
199 G4double rhomax = rmax;
200 if (stheta > halfpi) { rhomax = rmax*sinStart; }
201 if (etheta < halfpi) { rhomax = rmax*sinEnd; }
202
203 G4TwoVector xymin,xymax;
204 G4GeomTools::DiskExtent(rhomin,rhomax,
205 GetSinStartPhi(),GetCosStartPhi(),
206 GetSinEndPhi(),GetCosEndPhi(),
207 xymin,xymax);
208
209 G4double zmin = std::min(rmin*cosEnd,rmax*cosEnd);
210 G4double zmax = std::max(rmin*cosStart,rmax*cosStart);
211 pMin.set(xymin.x(),xymin.y(),zmin);
212 pMax.set(xymax.x(),xymax.y(),zmax);
213 }
214
215 // Check correctness of the bounding box
216 //
217 if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
218 {
219 std::ostringstream message;
220 message << "Bad bounding box (min >= max) for solid: "
221 << GetName() << " !"
222 << "\npMin = " << pMin
223 << "\npMax = " << pMax;
224 G4Exception("G4Sphere::BoundingLimits()", "GeomMgt0001",
225 JustWarning, message);
226 DumpInfo();
227 }
228}
JustWarning
@ JustWarning
Definition G4ExceptionSeverity.hh:71
G4TwoVector
CLHEP::Hep2Vector G4TwoVector
Definition G4TwoVector.hh:36
G4double
double G4double
Definition G4Types.hh:83
CLHEP::Hep2Vector::x
double x() const
CLHEP::Hep2Vector::y
double y() const
CLHEP::Hep3Vector::z
double z() const
CLHEP::Hep3Vector::x
double x() const
CLHEP::Hep3Vector::y
double y() const
CLHEP::Hep3Vector::set
void set(double x, double y, double z)
G4GeomTools::DiskExtent
static G4bool DiskExtent(G4double rmin, G4double rmax, G4double startPhi, G4double delPhi, G4TwoVector &pmin, G4TwoVector &pmax)
Definition G4GeomTools.cc:399
G4Sphere::GetSinStartTheta
G4double GetSinStartTheta() const
G4Sphere::GetCosStartPhi
G4double GetCosStartPhi() const
G4Sphere::GetDeltaPhiAngle
G4double GetDeltaPhiAngle() const
G4Sphere::GetCosEndTheta
G4double GetCosEndTheta() const
G4Sphere::GetInnerRadius
G4double GetInnerRadius() const
G4Sphere::GetOuterRadius
G4double GetOuterRadius() const
G4Sphere::GetCosEndPhi
G4double GetCosEndPhi() const
G4Sphere::GetSinEndTheta
G4double GetSinEndTheta() const
G4Sphere::GetDeltaThetaAngle
G4double GetDeltaThetaAngle() const
G4Sphere::GetSinEndPhi
G4double GetSinEndPhi() const
G4Sphere::GetSinStartPhi
G4double GetSinStartPhi() const
G4Sphere::GetStartThetaAngle
G4double GetStartThetaAngle() const
G4Sphere::GetCosStartTheta
G4double GetCosStartTheta() const
G4VSolid::DumpInfo
void DumpInfo() const

Referenced by CalculateExtent().

◆ CalculateExtent()

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

Calculates the minimum and maximum extent of the solid, when under the specified transform, and within the specified limits.

Parameters
[in]pAxisThe axis along which compute the extent.
[in]pVoxelLimitThe limiting space dictated by voxels.
[in]pTransformThe internal transformation applied to the solid.
[out]pMinThe minimum extent value.
[out]pMaxThe maximum extent value.
Returns
True if the solid is intersected by the extent region.

Implements G4VSolid.

Definition at line 234 of file G4Sphere.cc.

238{
239 G4ThreeVector bmin, bmax;
240
241 // Get bounding box
242 BoundingLimits(bmin,bmax);
243
244 // Find extent
245 G4BoundingEnvelope bbox(bmin,bmax);
246 return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
247}
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition G4ThreeVector.hh:36
G4Sphere::BoundingLimits
void BoundingLimits(G4ThreeVector &pMin, G4ThreeVector &pMax) const override
Definition G4Sphere.cc:177

◆ Clone()

G4VSolid * G4Sphere::Clone ( ) const
overridevirtual

Makes a clone of the object for use in multi-treading.

Returns
A pointer to the new cloned allocated solid.

Reimplemented from G4VSolid.

Definition at line 2711 of file G4Sphere.cc.

2712{
2713 return new G4Sphere(*this);
2714}
G4Sphere::G4Sphere
G4Sphere(const G4String &pName, G4double pRmin, G4double pRmax, G4double pSPhi, G4double pDPhi, G4double pSTheta, G4double pDTheta)
Definition G4Sphere.cc:81

◆ ComputeDimensions()

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

Dispatch method for parameterisation replication mechanism and dimension computation.

Reimplemented from G4VSolid.

Definition at line 166 of file G4Sphere.cc.

169{
170 p->ComputeDimensions(*this,n,pRep);
171}
G4VPVParameterisation::ComputeDimensions
virtual void ComputeDimensions(G4Box &, const G4int, const G4VPhysicalVolume *) const
Definition G4VPVParameterisation.hh:123

◆ CreatePolyhedron()

G4Polyhedron * G4Sphere::CreatePolyhedron ( ) const
overridevirtual

Creates a Polyhedron used for Visualisation. It is the caller's responsibility to delete it. A null pointer means "not created".

Reimplemented from G4VSolid.

Definition at line 2843 of file G4Sphere.cc.

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

◆ DescribeYourselfTo()

void G4Sphere::DescribeYourselfTo ( G4VGraphicsScene & scene) const
overridevirtual

A "double dispatch" function which identifies the solid to the graphics scene for visualization.

Implements G4VSolid.

Definition at line 2838 of file G4Sphere.cc.

2839{
2840 scene.AddSolid (*this);
2841}
G4VGraphicsScene::AddSolid
virtual void AddSolid(const G4Box &)=0

◆ DistanceToIn() [1/2]

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

Calculates the distance to the nearest surface of a shape from an outside point. The distance can be an underestimate.

Parameters
[in]pThe point at offset p.
Returns
The safety distance to enter the shape.

Implements G4VSolid.

Definition at line 1618 of file G4Sphere.cc.

1619{
1620 G4double safe=0.0,safeRMin,safeRMax,safePhi,safeTheta;
1621 G4double rho2,rds,rho;
1622 G4double cosPsi;
1623 G4double pTheta,dTheta1,dTheta2;
1624 rho2=p.x()*p.x()+p.y()*p.y();
1625 rds=std::sqrt(rho2+p.z()*p.z());
1626 rho=std::sqrt(rho2);
1627
1628 //
1629 // Distance to r shells
1630 //
1631 if (fRmin != 0.0)
1632 {
1633 safeRMin=fRmin-rds;
1634 safeRMax=rds-fRmax;
1635 if (safeRMin>safeRMax)
1636 {
1637 safe=safeRMin;
1638 }
1639 else
1640 {
1641 safe=safeRMax;
1642 }
1643 }
1644 else
1645 {
1646 safe=rds-fRmax;
1647 }
1648
1649 //
1650 // Distance to phi extent
1651 //
1652 if (!fFullPhiSphere && (rho != 0.0))
1653 {
1654 // Psi=angle from central phi to point
1655 //
1656 cosPsi=(p.x()*cosCPhi+p.y()*sinCPhi)/rho;
1657 if (cosPsi<cosHDPhi)
1658 {
1659 // Point lies outside phi range
1660 //
1661 if ((p.y()*cosCPhi-p.x()*sinCPhi)<=0)
1662 {
1663 safePhi=std::fabs(p.x()*sinSPhi-p.y()*cosSPhi);
1664 }
1665 else
1666 {
1667 safePhi=std::fabs(p.x()*sinEPhi-p.y()*cosEPhi);
1668 }
1669 if (safePhi>safe) { safe=safePhi; }
1670 }
1671 }
1672 //
1673 // Distance to Theta extent
1674 //
1675 if ((rds!=0.0) && (!fFullThetaSphere))
1676 {
1677 pTheta=std::acos(p.z()/rds);
1678 if (pTheta<0) { pTheta+=pi; }
1679 dTheta1=fSTheta-pTheta;
1680 dTheta2=pTheta-eTheta;
1681 if (dTheta1>dTheta2)
1682 {
1683 if (dTheta1>=0) // WHY ???????????
1684 {
1685 safeTheta=rds*std::sin(dTheta1);
1686 if (safe<=safeTheta)
1687 {
1688 safe=safeTheta;
1689 }
1690 }
1691 }
1692 else
1693 {
1694 if (dTheta2>=0)
1695 {
1696 safeTheta=rds*std::sin(dTheta2);
1697 if (safe<=safeTheta)
1698 {
1699 safe=safeTheta;
1700 }
1701 }
1702 }
1703 }
1704
1705 if (safe<0) { safe=0; }
1706 return safe;
1707}
G4INCL::Math::pi
const G4double pi
Definition G4INCLGlobals.hh:69

◆ DistanceToIn() [2/2]

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

Returns the distance along the normalised vector 'v' to the shape, from the point at offset 'p'. If there is no intersection, returns kInfinity. The first intersection resulting from 'leaving' a surface/volume is discarded. Hence, it is tolerant of points on the surface of the shape.

Parameters
[in]pThe point at offset p.
[in]vThe normalised direction vector.
Returns
The distance to enter the shape.

Implements G4VSolid.

Definition at line 669 of file G4Sphere.cc.

671{
672 G4double snxt = kInfinity ; // snxt = default return value
673 G4double rho2, rad2, pDotV2d, pDotV3d, pTheta ;
674 G4double tolSTheta=0., tolETheta=0. ;
675 const G4double dRmax = 100.*fRmax;
676
677 const G4double halfRmaxTolerance = fRmaxTolerance*0.5;
678 const G4double halfRminTolerance = fRminTolerance*0.5;
679 const G4double tolORMin2 = (fRmin>halfRminTolerance)
680 ? (fRmin-halfRminTolerance)*(fRmin-halfRminTolerance) : 0;
681 const G4double tolIRMin2 =
682 (fRmin+halfRminTolerance)*(fRmin+halfRminTolerance);
683 const G4double tolORMax2 =
684 (fRmax+halfRmaxTolerance)*(fRmax+halfRmaxTolerance);
685 const G4double tolIRMax2 =
686 (fRmax-halfRmaxTolerance)*(fRmax-halfRmaxTolerance);
687
688 // Intersection point
689 //
690 G4double xi, yi, zi, rhoi, rhoi2, radi2, iTheta ;
691
692 // Phi intersection
693 //
694 G4double Comp ;
695
696 // Phi precalcs
697 //
698 G4double Dist, cosPsi ;
699
700 // Theta precalcs
701 //
702 G4double dist2STheta, dist2ETheta ;
703 G4double t1, t2, b, c, d2, d, sd = kInfinity ;
704
705 // General Precalcs
706 //
707 rho2 = p.x()*p.x() + p.y()*p.y() ;
708 rad2 = rho2 + p.z()*p.z() ;
709 pTheta = std::atan2(std::sqrt(rho2),p.z()) ;
710
711 pDotV2d = p.x()*v.x() + p.y()*v.y() ;
712 pDotV3d = pDotV2d + p.z()*v.z() ;
713
714 // Theta precalcs
715 //
716 if (!fFullThetaSphere)
717 {
718 tolSTheta = fSTheta - halfAngTolerance ;
719 tolETheta = eTheta + halfAngTolerance ;
720
721 // Special case rad2 = 0 comparing with direction
722 //
723 if ((rad2!=0.0) || (fRmin!=0.0))
724 {
725 // Keep going for computation of distance...
726 }
727 else // Positioned on the sphere's origin
728 {
729 G4double vTheta = std::atan2(std::sqrt(v.x()*v.x()+v.y()*v.y()),v.z()) ;
730 if ( (vTheta < tolSTheta) || (vTheta > tolETheta) )
731 {
732 return snxt ; // kInfinity
733 }
734 return snxt = 0.0 ;
735 }
736 }
737
738 // Outer spherical shell intersection
739 // - Only if outside tolerant fRmax
740 // - Check for if inside and outer G4Sphere heading through solid (-> 0)
741 // - No intersect -> no intersection with G4Sphere
742 //
743 // Shell eqn: x^2+y^2+z^2=RSPH^2
744 //
745 // => (px+svx)^2+(py+svy)^2+(pz+svz)^2=R^2
746 //
747 // => (px^2+py^2+pz^2) +2sd(pxvx+pyvy+pzvz)+sd^2(vx^2+vy^2+vz^2)=R^2
748 // => rad2 +2sd(pDotV3d) +sd^2 =R^2
749 //
750 // => sd=-pDotV3d+-std::sqrt(pDotV3d^2-(rad2-R^2))
751
752 c = rad2 - fRmax*fRmax ;
753
754 if (c > fRmaxTolerance*fRmax)
755 {
756 // If outside tolerant boundary of outer G4Sphere
757 // [should be std::sqrt(rad2)-fRmax > halfRmaxTolerance]
758
759 d2 = pDotV3d*pDotV3d - c ;
760
761 if ( d2 >= 0 )
762 {
763 sd = -pDotV3d - std::sqrt(d2) ;
764
765 if (sd >= 0 )
766 {
767 if ( sd>dRmax ) // Avoid rounding errors due to precision issues seen on
768 { // 64 bits systems. Split long distances and recompute
769 G4double fTerm = sd-std::fmod(sd,dRmax);
770 sd = fTerm + DistanceToIn(p+fTerm*v,v);
771 }
772 xi = p.x() + sd*v.x() ;
773 yi = p.y() + sd*v.y() ;
774 rhoi = std::sqrt(xi*xi + yi*yi) ;
775
776 if (!fFullPhiSphere && (rhoi != 0.0)) // Check phi intersection
777 {
778 cosPsi = (xi*cosCPhi + yi*sinCPhi)/rhoi ;
779
780 if (cosPsi >= cosHDPhiOT)
781 {
782 if (!fFullThetaSphere) // Check theta intersection
783 {
784 zi = p.z() + sd*v.z() ;
785
786 // rhoi & zi can never both be 0
787 // (=>intersect at origin =>fRmax=0)
788 //
789 iTheta = std::atan2(rhoi,zi) ;
790 if ( (iTheta >= tolSTheta) && (iTheta <= tolETheta) )
791 {
792 return snxt = sd ;
793 }
794 }
795 else
796 {
797 return snxt=sd;
798 }
799 }
800 }
801 else
802 {
803 if (!fFullThetaSphere) // Check theta intersection
804 {
805 zi = p.z() + sd*v.z() ;
806
807 // rhoi & zi can never both be 0
808 // (=>intersect at origin => fRmax=0 !)
809 //
810 iTheta = std::atan2(rhoi,zi) ;
811 if ( (iTheta >= tolSTheta) && (iTheta <= tolETheta) )
812 {
813 return snxt=sd;
814 }
815 }
816 else
817 {
818 return snxt = sd;
819 }
820 }
821 }
822 }
823 else // No intersection with G4Sphere
824 {
825 return snxt=kInfinity;
826 }
827 }
828 else
829 {
830 // Inside outer radius
831 // check not inside, and heading through G4Sphere (-> 0 to in)
832
833 d2 = pDotV3d*pDotV3d - c ;
834
835 if ( (rad2 > tolIRMax2)
836 && ( (d2 >= fRmaxTolerance*fRmax) && (pDotV3d < 0) ) )
837 {
838 if (!fFullPhiSphere)
839 {
840 // Use inner phi tolerant boundary -> if on tolerant
841 // phi boundaries, phi intersect code handles leaving/entering checks
842
843 cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/std::sqrt(rho2) ;
844
845 if (cosPsi>=cosHDPhiIT)
846 {
847 // inside radii, delta r -ve, inside phi
848
849 if ( !fFullThetaSphere )
850 {
851 if ( (pTheta >= tolSTheta + kAngTolerance)
852 && (pTheta <= tolETheta - kAngTolerance) )
853 {
854 return snxt=0;
855 }
856 }
857 else // strictly inside Theta in both cases
858 {
859 return snxt=0;
860 }
861 }
862 }
863 else
864 {
865 if ( !fFullThetaSphere )
866 {
867 if ( (pTheta >= tolSTheta + kAngTolerance)
868 && (pTheta <= tolETheta - kAngTolerance) )
869 {
870 return snxt=0;
871 }
872 }
873 else // strictly inside Theta in both cases
874 {
875 return snxt=0;
876 }
877 }
878 }
879 }
880
881 // Inner spherical shell intersection
882 // - Always farthest root, because would have passed through outer
883 // surface first.
884 // - Tolerant check if travelling through solid
885
886 if (fRmin != 0.0)
887 {
888 c = rad2 - fRmin*fRmin ;
889 d2 = pDotV3d*pDotV3d - c ;
890
891 // Within tolerance inner radius of inner G4Sphere
892 // Check for immediate entry/already inside and travelling outwards
893
894 if ( (c > -halfRminTolerance) && (rad2 < tolIRMin2)
895 && ( (d2 < fRmin*kCarTolerance) || (pDotV3d >= 0) ) )
896 {
897 if ( !fFullPhiSphere )
898 {
899 // Use inner phi tolerant boundary -> if on tolerant
900 // phi boundaries, phi intersect code handles leaving/entering checks
901
902 cosPsi = (p.x()*cosCPhi+p.y()*sinCPhi)/std::sqrt(rho2) ;
903 if (cosPsi >= cosHDPhiIT)
904 {
905 // inside radii, delta r -ve, inside phi
906 //
907 if ( !fFullThetaSphere )
908 {
909 if ( (pTheta >= tolSTheta + kAngTolerance)
910 && (pTheta <= tolETheta - kAngTolerance) )
911 {
912 return snxt=0;
913 }
914 }
915 else
916 {
917 return snxt = 0 ;
918 }
919 }
920 }
921 else
922 {
923 if ( !fFullThetaSphere )
924 {
925 if ( (pTheta >= tolSTheta + kAngTolerance)
926 && (pTheta <= tolETheta - kAngTolerance) )
927 {
928 return snxt = 0 ;
929 }
930 }
931 else
932 {
933 return snxt=0;
934 }
935 }
936 }
937 else // Not special tolerant case
938 {
939 if (d2 >= 0)
940 {
941 sd = -pDotV3d + std::sqrt(d2) ;
942 if ( sd >= halfRminTolerance ) // It was >= 0 ??
943 {
944 xi = p.x() + sd*v.x() ;
945 yi = p.y() + sd*v.y() ;
946 rhoi = std::sqrt(xi*xi+yi*yi) ;
947
948 if ( !fFullPhiSphere && (rhoi != 0.0) ) // Check phi intersection
949 {
950 cosPsi = (xi*cosCPhi + yi*sinCPhi)/rhoi ;
951
952 if (cosPsi >= cosHDPhiOT)
953 {
954 if ( !fFullThetaSphere ) // Check theta intersection
955 {
956 zi = p.z() + sd*v.z() ;
957
958 // rhoi & zi can never both be 0
959 // (=>intersect at origin =>fRmax=0)
960 //
961 iTheta = std::atan2(rhoi,zi) ;
962 if ( (iTheta >= tolSTheta) && (iTheta<=tolETheta) )
963 {
964 snxt = sd;
965 }
966 }
967 else
968 {
969 snxt=sd;
970 }
971 }
972 }
973 else
974 {
975 if ( !fFullThetaSphere ) // Check theta intersection
976 {
977 zi = p.z() + sd*v.z() ;
978
979 // rhoi & zi can never both be 0
980 // (=>intersect at origin => fRmax=0 !)
981 //
982 iTheta = std::atan2(rhoi,zi) ;
983 if ( (iTheta >= tolSTheta) && (iTheta <= tolETheta) )
984 {
985 snxt = sd;
986 }
987 }
988 else
989 {
990 snxt = sd;
991 }
992 }
993 }
994 }
995 }
996 }
997
998 // Phi segment intersection
999 //
1000 // o Tolerant of points inside phi planes by up to kCarTolerance*0.5
1001 //
1002 // o NOTE: Large duplication of code between sphi & ephi checks
1003 // -> only diffs: sphi -> ephi, Comp -> -Comp and half-plane
1004 // intersection check <=0 -> >=0
1005 // -> Should use some form of loop Construct
1006 //
1007 if ( !fFullPhiSphere )
1008 {
1009 // First phi surface ('S'tarting phi)
1010 // Comp = Component in outwards normal dirn
1011 //
1012 Comp = v.x()*sinSPhi - v.y()*cosSPhi ;
1013
1014 if ( Comp < 0 )
1015 {
1016 Dist = p.y()*cosSPhi - p.x()*sinSPhi ;
1017
1018 if (Dist < halfCarTolerance)
1019 {
1020 sd = Dist/Comp ;
1021
1022 if (sd < snxt)
1023 {
1024 if ( sd > 0 )
1025 {
1026 xi = p.x() + sd*v.x() ;
1027 yi = p.y() + sd*v.y() ;
1028 zi = p.z() + sd*v.z() ;
1029 rhoi2 = xi*xi + yi*yi ;
1030 radi2 = rhoi2 + zi*zi ;
1031 }
1032 else
1033 {
1034 sd = 0 ;
1035 xi = p.x() ;
1036 yi = p.y() ;
1037 zi = p.z() ;
1038 rhoi2 = rho2 ;
1039 radi2 = rad2 ;
1040 }
1041 if ( (radi2 <= tolORMax2)
1042 && (radi2 >= tolORMin2)
1043 && ((yi*cosCPhi-xi*sinCPhi) <= 0) )
1044 {
1045 // Check theta intersection
1046 // rhoi & zi can never both be 0
1047 // (=>intersect at origin =>fRmax=0)
1048 //
1049 if ( !fFullThetaSphere )
1050 {
1051 iTheta = std::atan2(std::sqrt(rhoi2),zi) ;
1052 if ( (iTheta >= tolSTheta) && (iTheta <= tolETheta) )
1053 {
1054 // r and theta intersections good
1055 // - check intersecting with correct half-plane
1056
1057 if ((yi*cosCPhi-xi*sinCPhi) <= 0)
1058 {
1059 snxt = sd;
1060 }
1061 }
1062 }
1063 else
1064 {
1065 snxt = sd;
1066 }
1067 }
1068 }
1069 }
1070 }
1071
1072 // Second phi surface ('E'nding phi)
1073 // Component in outwards normal dirn
1074
1075 Comp = -( v.x()*sinEPhi-v.y()*cosEPhi ) ;
1076
1077 if (Comp < 0)
1078 {
1079 Dist = -(p.y()*cosEPhi-p.x()*sinEPhi) ;
1080 if ( Dist < halfCarTolerance )
1081 {
1082 sd = Dist/Comp ;
1083
1084 if ( sd < snxt )
1085 {
1086 if (sd > 0)
1087 {
1088 xi = p.x() + sd*v.x() ;
1089 yi = p.y() + sd*v.y() ;
1090 zi = p.z() + sd*v.z() ;
1091 rhoi2 = xi*xi + yi*yi ;
1092 radi2 = rhoi2 + zi*zi ;
1093 }
1094 else
1095 {
1096 sd = 0 ;
1097 xi = p.x() ;
1098 yi = p.y() ;
1099 zi = p.z() ;
1100 rhoi2 = rho2 ;
1101 radi2 = rad2 ;
1102 }
1103 if ( (radi2 <= tolORMax2)
1104 && (radi2 >= tolORMin2)
1105 && ((yi*cosCPhi-xi*sinCPhi) >= 0) )
1106 {
1107 // Check theta intersection
1108 // rhoi & zi can never both be 0
1109 // (=>intersect at origin =>fRmax=0)
1110 //
1111 if ( !fFullThetaSphere )
1112 {
1113 iTheta = std::atan2(std::sqrt(rhoi2),zi) ;
1114 if ( (iTheta >= tolSTheta) && (iTheta <= tolETheta) )
1115 {
1116 // r and theta intersections good
1117 // - check intersecting with correct half-plane
1118
1119 if ((yi*cosCPhi-xi*sinCPhi) >= 0)
1120 {
1121 snxt = sd;
1122 }
1123 }
1124 }
1125 else
1126 {
1127 snxt = sd;
1128 }
1129 }
1130 }
1131 }
1132 }
1133 }
1134
1135 // Theta segment intersection
1136
1137 if ( !fFullThetaSphere )
1138 {
1139
1140 // Intersection with theta surfaces
1141 // Known failure cases:
1142 // o Inside tolerance of stheta surface, skim
1143 // ~parallel to cone and Hit & enter etheta surface [& visa versa]
1144 //
1145 // To solve: Check 2nd root of etheta surface in addition to stheta
1146 //
1147 // o start/end theta is exactly pi/2
1148 // Intersections with cones
1149 //
1150 // Cone equation: x^2+y^2=z^2tan^2(t)
1151 //
1152 // => (px+svx)^2+(py+svy)^2=(pz+svz)^2tan^2(t)
1153 //
1154 // => (px^2+py^2-pz^2tan^2(t))+2sd(pxvx+pyvy-pzvztan^2(t))
1155 // + sd^2(vx^2+vy^2-vz^2tan^2(t)) = 0
1156 //
1157 // => sd^2(1-vz^2(1+tan^2(t))+2sd(pdotv2d-pzvztan^2(t))
1158 // + (rho2-pz^2tan^2(t)) = 0
1159
1160 if (fSTheta != 0.0)
1161 {
1162 dist2STheta = rho2 - p.z()*p.z()*tanSTheta2 ;
1163 }
1164 else
1165 {
1166 dist2STheta = kInfinity ;
1167 }
1168 if ( eTheta < pi )
1169 {
1170 dist2ETheta=rho2-p.z()*p.z()*tanETheta2;
1171 }
1172 else
1173 {
1174 dist2ETheta=kInfinity;
1175 }
1176 if ( pTheta < tolSTheta )
1177 {
1178 // Inside (theta<stheta-tol) stheta cone
1179 // First root of stheta cone, second if first root -ve
1180
1181 t1 = 1 - v.z()*v.z()*(1 + tanSTheta2) ;
1182 t2 = pDotV2d - p.z()*v.z()*tanSTheta2 ;
1183 if (t1 != 0.0)
1184 {
1185 b = t2/t1 ;
1186 c = dist2STheta/t1 ;
1187 d2 = b*b - c ;
1188
1189 if ( d2 >= 0 )
1190 {
1191 d = std::sqrt(d2) ;
1192 sd = -b - d ; // First root
1193 zi = p.z() + sd*v.z();
1194
1195 if ( (sd < 0) || (zi*(fSTheta - halfpi) > 0) )
1196 {
1197 sd = -b+d; // Second root
1198 }
1199 if ((sd >= 0) && (sd < snxt))
1200 {
1201 xi = p.x() + sd*v.x();
1202 yi = p.y() + sd*v.y();
1203 zi = p.z() + sd*v.z();
1204 rhoi2 = xi*xi + yi*yi;
1205 radi2 = rhoi2 + zi*zi;
1206 if ( (radi2 <= tolORMax2)
1207 && (radi2 >= tolORMin2)
1208 && (zi*(fSTheta - halfpi) <= 0) )
1209 {
1210 if ( !fFullPhiSphere && (rhoi2 != 0.0) ) // Check phi intersection
1211 {
1212 cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
1213 if (cosPsi >= cosHDPhiOT)
1214 {
1215 snxt = sd;
1216 }
1217 }
1218 else
1219 {
1220 snxt = sd;
1221 }
1222 }
1223 }
1224 }
1225 }
1226
1227 // Possible intersection with ETheta cone.
1228 // Second >= 0 root should be considered
1229
1230 if ( eTheta < pi )
1231 {
1232 t1 = 1 - v.z()*v.z()*(1 + tanETheta2) ;
1233 t2 = pDotV2d - p.z()*v.z()*tanETheta2 ;
1234 if (t1 != 0.0)
1235 {
1236 b = t2/t1 ;
1237 c = dist2ETheta/t1 ;
1238 d2 = b*b - c ;
1239
1240 if (d2 >= 0)
1241 {
1242 d = std::sqrt(d2) ;
1243 sd = -b + d ; // Second root
1244
1245 if ( (sd >= 0) && (sd < snxt) )
1246 {
1247 xi = p.x() + sd*v.x() ;
1248 yi = p.y() + sd*v.y() ;
1249 zi = p.z() + sd*v.z() ;
1250 rhoi2 = xi*xi + yi*yi ;
1251 radi2 = rhoi2 + zi*zi ;
1252
1253 if ( (radi2 <= tolORMax2)
1254 && (radi2 >= tolORMin2)
1255 && (zi*(eTheta - halfpi) <= 0) )
1256 {
1257 if (!fFullPhiSphere && (rhoi2 != 0.0)) // Check phi intersection
1258 {
1259 cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
1260 if (cosPsi >= cosHDPhiOT)
1261 {
1262 snxt = sd;
1263 }
1264 }
1265 else
1266 {
1267 snxt = sd;
1268 }
1269 }
1270 }
1271 }
1272 }
1273 }
1274 }
1275 else if ( pTheta > tolETheta )
1276 {
1277 // dist2ETheta<-kRadTolerance*0.5 && dist2STheta>0)
1278 // Inside (theta > etheta+tol) e-theta cone
1279 // First root of etheta cone, second if first root 'imaginary'
1280
1281 t1 = 1 - v.z()*v.z()*(1 + tanETheta2) ;
1282 t2 = pDotV2d - p.z()*v.z()*tanETheta2 ;
1283 if (t1 != 0.0)
1284 {
1285 b = t2/t1 ;
1286 c = dist2ETheta/t1 ;
1287 d2 = b*b - c ;
1288
1289 if (d2 >= 0)
1290 {
1291 d = std::sqrt(d2) ;
1292 sd = -b - d ; // First root
1293 zi = p.z() + sd*v.z();
1294
1295 if ( (sd < 0) || (zi*(eTheta - halfpi) > 0) )
1296 {
1297 sd = -b + d ; // second root
1298 }
1299 if ( (sd >= 0) && (sd < snxt) )
1300 {
1301 xi = p.x() + sd*v.x() ;
1302 yi = p.y() + sd*v.y() ;
1303 zi = p.z() + sd*v.z() ;
1304 rhoi2 = xi*xi + yi*yi ;
1305 radi2 = rhoi2 + zi*zi ;
1306
1307 if ( (radi2 <= tolORMax2)
1308 && (radi2 >= tolORMin2)
1309 && (zi*(eTheta - halfpi) <= 0) )
1310 {
1311 if (!fFullPhiSphere && (rhoi2 != 0.0)) // Check phi intersection
1312 {
1313 cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
1314 if (cosPsi >= cosHDPhiOT)
1315 {
1316 snxt = sd;
1317 }
1318 }
1319 else
1320 {
1321 snxt = sd;
1322 }
1323 }
1324 }
1325 }
1326 }
1327
1328 // Possible intersection with STheta cone.
1329 // Second >= 0 root should be considered
1330
1331 if ( fSTheta != 0.0 )
1332 {
1333 t1 = 1 - v.z()*v.z()*(1 + tanSTheta2) ;
1334 t2 = pDotV2d - p.z()*v.z()*tanSTheta2 ;
1335 if (t1 != 0.0)
1336 {
1337 b = t2/t1 ;
1338 c = dist2STheta/t1 ;
1339 d2 = b*b - c ;
1340
1341 if (d2 >= 0)
1342 {
1343 d = std::sqrt(d2) ;
1344 sd = -b + d ; // Second root
1345
1346 if ( (sd >= 0) && (sd < snxt) )
1347 {
1348 xi = p.x() + sd*v.x() ;
1349 yi = p.y() + sd*v.y() ;
1350 zi = p.z() + sd*v.z() ;
1351 rhoi2 = xi*xi + yi*yi ;
1352 radi2 = rhoi2 + zi*zi ;
1353
1354 if ( (radi2 <= tolORMax2)
1355 && (radi2 >= tolORMin2)
1356 && (zi*(fSTheta - halfpi) <= 0) )
1357 {
1358 if (!fFullPhiSphere && (rhoi2 != 0.0)) // Check phi intersection
1359 {
1360 cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
1361 if (cosPsi >= cosHDPhiOT)
1362 {
1363 snxt = sd;
1364 }
1365 }
1366 else
1367 {
1368 snxt = sd;
1369 }
1370 }
1371 }
1372 }
1373 }
1374 }
1375 }
1376 else if ( (pTheta < tolSTheta + kAngTolerance)
1377 && (fSTheta > halfAngTolerance) )
1378 {
1379 // In tolerance of stheta
1380 // If entering through solid [r,phi] => 0 to in
1381 // else try 2nd root
1382
1383 t2 = pDotV2d - p.z()*v.z()*tanSTheta2 ;
1384 if ( (t2>=0 && tolIRMin2<rad2 && rad2<tolIRMax2 && fSTheta<halfpi)
1385 || (t2<0 && tolIRMin2<rad2 && rad2<tolIRMax2 && fSTheta>halfpi)
1386 || (v.z()<0 && tolIRMin2<rad2 && rad2<tolIRMax2 && fSTheta==halfpi) )
1387 {
1388 if (!fFullPhiSphere && (rho2 != 0.0)) // Check phi intersection
1389 {
1390 cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/std::sqrt(rho2) ;
1391 if (cosPsi >= cosHDPhiIT)
1392 {
1393 return 0 ;
1394 }
1395 }
1396 else
1397 {
1398 return 0 ;
1399 }
1400 }
1401
1402 // Not entering immediately/travelling through
1403
1404 t1 = 1 - v.z()*v.z()*(1 + tanSTheta2) ;
1405 if (t1 != 0.0)
1406 {
1407 b = t2/t1 ;
1408 c = dist2STheta/t1 ;
1409 d2 = b*b - c ;
1410
1411 if (d2 >= 0)
1412 {
1413 d = std::sqrt(d2) ;
1414 sd = -b + d ;
1415 if ( (sd >= halfCarTolerance) && (sd < snxt) && (fSTheta < halfpi) )
1416 { // ^^^^^^^^^^^^^^^^^^^^^ shouldn't it be >=0 instead ?
1417 xi = p.x() + sd*v.x() ;
1418 yi = p.y() + sd*v.y() ;
1419 zi = p.z() + sd*v.z() ;
1420 rhoi2 = xi*xi + yi*yi ;
1421 radi2 = rhoi2 + zi*zi ;
1422
1423 if ( (radi2 <= tolORMax2)
1424 && (radi2 >= tolORMin2)
1425 && (zi*(fSTheta - halfpi) <= 0) )
1426 {
1427 if ( !fFullPhiSphere && (rhoi2 != 0.0) ) // Check phi intersection
1428 {
1429 cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
1430 if ( cosPsi >= cosHDPhiOT )
1431 {
1432 snxt = sd;
1433 }
1434 }
1435 else
1436 {
1437 snxt = sd;
1438 }
1439 }
1440 }
1441 }
1442 }
1443 }
1444 else if ((pTheta > tolETheta-kAngTolerance) && (eTheta < pi-kAngTolerance))
1445 {
1446
1447 // In tolerance of etheta
1448 // If entering through solid [r,phi] => 0 to in
1449 // else try 2nd root
1450
1451 t2 = pDotV2d - p.z()*v.z()*tanETheta2 ;
1452
1453 if ( ((t2<0) && (eTheta < halfpi)
1454 && (tolIRMin2 < rad2) && (rad2 < tolIRMax2))
1455 || ((t2>=0) && (eTheta > halfpi)
1456 && (tolIRMin2 < rad2) && (rad2 < tolIRMax2))
1457 || ((v.z()>0) && (eTheta == halfpi)
1458 && (tolIRMin2 < rad2) && (rad2 < tolIRMax2)) )
1459 {
1460 if (!fFullPhiSphere && (rho2 != 0.0)) // Check phi intersection
1461 {
1462 cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/std::sqrt(rho2) ;
1463 if (cosPsi >= cosHDPhiIT)
1464 {
1465 return 0 ;
1466 }
1467 }
1468 else
1469 {
1470 return 0 ;
1471 }
1472 }
1473
1474 // Not entering immediately/travelling through
1475
1476 t1 = 1 - v.z()*v.z()*(1 + tanETheta2) ;
1477 if (t1 != 0.0)
1478 {
1479 b = t2/t1 ;
1480 c = dist2ETheta/t1 ;
1481 d2 = b*b - c ;
1482
1483 if (d2 >= 0)
1484 {
1485 d = std::sqrt(d2) ;
1486 sd = -b + d ;
1487
1488 if ( (sd >= halfCarTolerance)
1489 && (sd < snxt) && (eTheta > halfpi) )
1490 {
1491 xi = p.x() + sd*v.x() ;
1492 yi = p.y() + sd*v.y() ;
1493 zi = p.z() + sd*v.z() ;
1494 rhoi2 = xi*xi + yi*yi ;
1495 radi2 = rhoi2 + zi*zi ;
1496
1497 if ( (radi2 <= tolORMax2)
1498 && (radi2 >= tolORMin2)
1499 && (zi*(eTheta - halfpi) <= 0) )
1500 {
1501 if (!fFullPhiSphere && (rhoi2 != 0.0)) // Check phi intersection
1502 {
1503 cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
1504 if (cosPsi >= cosHDPhiOT)
1505 {
1506 snxt = sd;
1507 }
1508 }
1509 else
1510 {
1511 snxt = sd;
1512 }
1513 }
1514 }
1515 }
1516 }
1517 }
1518 else
1519 {
1520 // stheta+tol<theta<etheta-tol
1521 // For BOTH stheta & etheta check 2nd root for validity [r,phi]
1522
1523 t1 = 1 - v.z()*v.z()*(1 + tanSTheta2) ;
1524 t2 = pDotV2d - p.z()*v.z()*tanSTheta2 ;
1525 if (t1 != 0.0)
1526 {
1527 b = t2/t1;
1528 c = dist2STheta/t1 ;
1529 d2 = b*b - c ;
1530
1531 if (d2 >= 0)
1532 {
1533 d = std::sqrt(d2) ;
1534 sd = -b + d ; // second root
1535
1536 if ((sd >= 0) && (sd < snxt))
1537 {
1538 xi = p.x() + sd*v.x() ;
1539 yi = p.y() + sd*v.y() ;
1540 zi = p.z() + sd*v.z() ;
1541 rhoi2 = xi*xi + yi*yi ;
1542 radi2 = rhoi2 + zi*zi ;
1543
1544 if ( (radi2 <= tolORMax2)
1545 && (radi2 >= tolORMin2)
1546 && (zi*(fSTheta - halfpi) <= 0) )
1547 {
1548 if (!fFullPhiSphere && (rhoi2 != 0.0)) // Check phi intersection
1549 {
1550 cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
1551 if (cosPsi >= cosHDPhiOT)
1552 {
1553 snxt = sd;
1554 }
1555 }
1556 else
1557 {
1558 snxt = sd;
1559 }
1560 }
1561 }
1562 }
1563 }
1564 t1 = 1 - v.z()*v.z()*(1 + tanETheta2) ;
1565 t2 = pDotV2d - p.z()*v.z()*tanETheta2 ;
1566 if (t1 != 0.0)
1567 {
1568 b = t2/t1 ;
1569 c = dist2ETheta/t1 ;
1570 d2 = b*b - c ;
1571
1572 if (d2 >= 0)
1573 {
1574 d = std::sqrt(d2) ;
1575 sd = -b + d; // second root
1576
1577 if ((sd >= 0) && (sd < snxt))
1578 {
1579 xi = p.x() + sd*v.x() ;
1580 yi = p.y() + sd*v.y() ;
1581 zi = p.z() + sd*v.z() ;
1582 rhoi2 = xi*xi + yi*yi ;
1583 radi2 = rhoi2 + zi*zi ;
1584
1585 if ( (radi2 <= tolORMax2)
1586 && (radi2 >= tolORMin2)
1587 && (zi*(eTheta - halfpi) <= 0) )
1588 {
1589 if (!fFullPhiSphere && (rhoi2 != 0.0)) // Check phi intersection
1590 {
1591 cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
1592 if ( cosPsi >= cosHDPhiOT )
1593 {
1594 snxt = sd;
1595 }
1596 }
1597 else
1598 {
1599 snxt = sd;
1600 }
1601 }
1602 }
1603 }
1604 }
1605 }
1606 }
1607 return snxt;
1608}
G4Sphere::DistanceToIn
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const override
Definition G4Sphere.cc:669

Referenced by DistanceToIn().

◆ DistanceToOut() [1/2]

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

Calculates the distance to the nearest surface of a shape from an inside point 'p'. The distance can be an underestimate.

Parameters
[in]pThe point at offset p.
Returns
The safety distance to exit the shape.

Implements G4VSolid.

Definition at line 2609 of file G4Sphere.cc.

2610{
2611 G4double safe=0.0,safeRMin,safeRMax,safePhi,safeTheta;
2612 G4double rho2,rds,rho;
2613 G4double pTheta,dTheta1 = kInfinity,dTheta2 = kInfinity;
2614 rho2=p.x()*p.x()+p.y()*p.y();
2615 rds=std::sqrt(rho2+p.z()*p.z());
2616 rho=std::sqrt(rho2);
2617
2618#ifdef G4CSGDEBUG
2619 if( Inside(p) == kOutside )
2620 {
2621 G4long old_prc = G4cout.precision(16);
2622 G4cout << G4endl;
2623 DumpInfo();
2624 G4cout << "Position:" << G4endl << G4endl ;
2625 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl ;
2626 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ;
2627 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ;
2628 G4cout.precision(old_prc) ;
2629 G4Exception("G4Sphere::DistanceToOut(p)",
2630 "GeomSolids1002", JustWarning, "Point p is outside !?" );
2631 }
2632#endif
2633
2634 // Distance to r shells
2635 //
2636 safeRMax = fRmax-rds;
2637 safe = safeRMax;
2638 if (fRmin != 0.0)
2639 {
2640 safeRMin = rds-fRmin;
2641 safe = std::min( safeRMin, safeRMax );
2642 }
2643
2644 // Distance to phi extent
2645 //
2646 if ( !fFullPhiSphere )
2647 {
2648 if (rho>0.0)
2649 {
2650 if ((p.y()*cosCPhi-p.x()*sinCPhi)<=0)
2651 {
2652 safePhi=-(p.x()*sinSPhi-p.y()*cosSPhi);
2653 }
2654 else
2655 {
2656 safePhi=(p.x()*sinEPhi-p.y()*cosEPhi);
2657 }
2658 }
2659 else
2660 {
2661 safePhi = 0.0; // Distance to both Phi surfaces (extended)
2662 }
2663 // Both cases above can be improved - in case fRMin > 0.0
2664 // although it may be costlier (good for precise, not fast version)
2665
2666 safe= std::min(safe, safePhi);
2667 }
2668
2669 // Distance to Theta extent
2670 //
2671 if ( !fFullThetaSphere )
2672 {
2673 if( rds > 0.0 )
2674 {
2675 pTheta=std::acos(p.z()/rds);
2676 if (pTheta<0) { pTheta+=pi; }
2677 if(fSTheta>0.)
2678 { dTheta1=pTheta-fSTheta;}
2679 if(eTheta<pi)
2680 { dTheta2=eTheta-pTheta;}
2681
2682 safeTheta=rds*std::sin(std::min(dTheta1, dTheta2) );
2683 }
2684 else
2685 {
2686 safeTheta= 0.0;
2687 // An improvement will be to return negative answer if outside (TODO)
2688 }
2689 safe = std::min( safe, safeTheta );
2690 }
2691
2692 if (safe<0.0) { safe=0; }
2693 // An improvement to return negative answer if outside (TODO)
2694
2695 return safe;
2696}
G4long
long G4long
Definition G4Types.hh:87
G4cout
G4GLOB_DLL std::ostream G4cout
G4Sphere::Inside
EInside Inside(const G4ThreeVector &p) const override
Definition G4Sphere.cc:255
kOutside
@ kOutside
Definition geomdefs.hh:68

◆ DistanceToOut() [2/2]

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

Returns the distance along the normalised vector 'v' to the shape, from a point at an offset 'p' inside or on the surface of the shape. Intersections with surfaces, when the point is less than Tolerance/2 from a surface must be ignored.

Parameters
[in]pThe point at offset p.
[in]vThe normalised direction vector.
[in]calcNormFlag to indicate if to calculate the normal or not.
[out]validNormFlag set to true if the solid lies entirely behind or on the exiting surface. It is set false if the solid does not lie entirely behind or on the exiting surface. 'calcNorm' must be true, otherwise it is unused.
[out]nThe exiting outwards normal vector (undefined Magnitude). 'calcNorm' must be true, otherwise it is unused.
Returns
The distance to exit the shape.

Implements G4VSolid.

Definition at line 1714 of file G4Sphere.cc.

1719{
1720 G4double snxt = kInfinity; // snxt is default return value
1721 G4double sphi= kInfinity,stheta= kInfinity;
1722 ESide side=kNull,sidephi=kNull,sidetheta=kNull;
1723
1724 const G4double halfRmaxTolerance = fRmaxTolerance*0.5;
1725 const G4double halfRminTolerance = fRminTolerance*0.5;
1726 const G4double Rmax_plus = fRmax + halfRmaxTolerance;
1727 const G4double Rmin_minus = (fRmin) != 0.0 ? fRmin-halfRminTolerance : 0;
1728 G4double t1,t2;
1729 G4double b,c,d;
1730
1731 // Variables for phi intersection:
1732
1733 G4double pDistS,compS,pDistE,compE,sphi2,vphi;
1734
1735 G4double rho2,rad2,pDotV2d,pDotV3d;
1736
1737 G4double xi,yi,zi; // Intersection point
1738
1739 // Theta precals
1740 //
1741 G4double rhoSecTheta;
1742 G4double dist2STheta, dist2ETheta, distTheta;
1743 G4double d2,sd;
1744
1745 // General Precalcs
1746 //
1747 rho2 = p.x()*p.x()+p.y()*p.y();
1748 rad2 = rho2+p.z()*p.z();
1749
1750 pDotV2d = p.x()*v.x()+p.y()*v.y();
1751 pDotV3d = pDotV2d+p.z()*v.z();
1752
1753 // Radial Intersections from G4Sphere::DistanceToIn
1754 //
1755 // Outer spherical shell intersection
1756 // - Only if outside tolerant fRmax
1757 // - Check for if inside and outer G4Sphere heading through solid (-> 0)
1758 // - No intersect -> no intersection with G4Sphere
1759 //
1760 // Shell eqn: x^2+y^2+z^2=RSPH^2
1761 //
1762 // => (px+svx)^2+(py+svy)^2+(pz+svz)^2=R^2
1763 //
1764 // => (px^2+py^2+pz^2) +2sd(pxvx+pyvy+pzvz)+sd^2(vx^2+vy^2+vz^2)=R^2
1765 // => rad2 +2sd(pDotV3d) +sd^2 =R^2
1766 //
1767 // => sd=-pDotV3d+-std::sqrt(pDotV3d^2-(rad2-R^2))
1768
1769 if( (rad2 <= Rmax_plus*Rmax_plus) && (rad2 >= Rmin_minus*Rmin_minus) )
1770 {
1771 c = rad2 - fRmax*fRmax;
1772
1773 if (c < fRmaxTolerance*fRmax)
1774 {
1775 // Within tolerant Outer radius
1776 //
1777 // The test is
1778 // rad - fRmax < 0.5*kRadTolerance
1779 // => rad < fRmax + 0.5*kRadTol
1780 // => rad2 < (fRmax + 0.5*kRadTol)^2
1781 // => rad2 < fRmax^2 + 2.*0.5*fRmax*kRadTol + 0.25*kRadTol*kRadTol
1782 // => rad2 - fRmax^2 <~ fRmax*kRadTol
1783
1784 d2 = pDotV3d*pDotV3d - c;
1785
1786 if( (c >- fRmaxTolerance*fRmax) // on tolerant surface
1787 && ((pDotV3d >=0) || (d2 < 0)) ) // leaving outside from Rmax
1788 // not re-entering
1789 {
1790 if(calcNorm)
1791 {
1792 *validNorm = true ;
1793 *n = G4ThreeVector(p.x()/fRmax,p.y()/fRmax,p.z()/fRmax) ;
1794 }
1795 return snxt = 0;
1796 }
1797
1798 snxt = -pDotV3d+std::sqrt(d2); // second root since inside Rmax
1799 side = kRMax ;
1800 }
1801
1802 // Inner spherical shell intersection:
1803 // Always first >=0 root, because would have passed
1804 // from outside of Rmin surface .
1805
1806 if (fRmin != 0.0)
1807 {
1808 c = rad2 - fRmin*fRmin;
1809 d2 = pDotV3d*pDotV3d - c;
1810
1811 if (c >- fRminTolerance*fRmin) // 2.0 * (0.5*kRadTolerance) * fRmin
1812 {
1813 if ( (c < fRminTolerance*fRmin) // leaving from Rmin
1814 && (d2 >= fRminTolerance*fRmin) && (pDotV3d < 0) )
1815 {
1816 if(calcNorm) { *validNorm = false; } // Rmin surface is concave
1817 return snxt = 0 ;
1818 }
1819
1820 if ( d2 >= 0. )
1821 {
1822 sd = -pDotV3d-std::sqrt(d2);
1823
1824 if ( sd >= 0. ) // Always intersect Rmin first
1825 {
1826 snxt = sd ;
1827 side = kRMin ;
1828 }
1829 }
1830 }
1831 }
1832 }
1833
1834 // Theta segment intersection
1835
1836 if ( !fFullThetaSphere )
1837 {
1838 // Intersection with theta surfaces
1839 //
1840 // Known failure cases:
1841 // o Inside tolerance of stheta surface, skim
1842 // ~parallel to cone and Hit & enter etheta surface [& visa versa]
1843 //
1844 // To solve: Check 2nd root of etheta surface in addition to stheta
1845 //
1846 // o start/end theta is exactly pi/2
1847 //
1848 // Intersections with cones
1849 //
1850 // Cone equation: x^2+y^2=z^2tan^2(t)
1851 //
1852 // => (px+svx)^2+(py+svy)^2=(pz+svz)^2tan^2(t)
1853 //
1854 // => (px^2+py^2-pz^2tan^2(t))+2sd(pxvx+pyvy-pzvztan^2(t))
1855 // + sd^2(vx^2+vy^2-vz^2tan^2(t)) = 0
1856 //
1857 // => sd^2(1-vz^2(1+tan^2(t))+2sd(pdotv2d-pzvztan^2(t))
1858 // + (rho2-pz^2tan^2(t)) = 0
1859 //
1860
1861 if(fSTheta != 0.0) // intersection with first cons
1862 {
1863 if( std::fabs(tanSTheta) > 5./kAngTolerance ) // kons is plane z=0
1864 {
1865 if( v.z() > 0. )
1866 {
1867 if ( std::fabs( p.z() ) <= halfRmaxTolerance )
1868 {
1869 if(calcNorm)
1870 {
1871 *validNorm = true;
1872 *n = G4ThreeVector(0.,0.,1.);
1873 }
1874 return snxt = 0 ;
1875 }
1876 stheta = -p.z()/v.z();
1877 sidetheta = kSTheta;
1878 }
1879 }
1880 else // kons is not plane
1881 {
1882 t1 = 1-v.z()*v.z()*(1+tanSTheta2);
1883 t2 = pDotV2d-p.z()*v.z()*tanSTheta2; // ~vDotN if p on cons
1884 dist2STheta = rho2-p.z()*p.z()*tanSTheta2; // t3
1885
1886 distTheta = std::sqrt(rho2)-p.z()*tanSTheta;
1887
1888 if( std::fabs(t1) < halfAngTolerance ) // 1st order equation,
1889 { // v parallel to kons
1890 if( v.z() > 0. )
1891 {
1892 if(std::fabs(distTheta) < halfRmaxTolerance) // p on surface
1893 {
1894 if( (fSTheta < halfpi) && (p.z() > 0.) )
1895 {
1896 if( calcNorm ) { *validNorm = false; }
1897 return snxt = 0.;
1898 }
1899 if( (fSTheta > halfpi) && (p.z() <= 0) )
1900 {
1901 if( calcNorm )
1902 {
1903 *validNorm = true;
1904 if (rho2 != 0.0)
1905 {
1906 rhoSecTheta = std::sqrt(rho2*(1+tanSTheta2));
1907
1908 *n = G4ThreeVector( p.x()/rhoSecTheta,
1909 p.y()/rhoSecTheta,
1910 std::sin(fSTheta) );
1911 }
1912 else
1913 {
1914 *n = G4ThreeVector(0.,0.,1.);
1915 }
1916 }
1917 return snxt = 0.;
1918 }
1919 }
1920 stheta = -0.5*dist2STheta/t2;
1921 sidetheta = kSTheta;
1922 }
1923 } // 2nd order equation, 1st root of fSTheta cone,
1924 else // 2nd if 1st root -ve
1925 {
1926 if( std::fabs(distTheta) < halfRmaxTolerance )
1927 {
1928 if( (fSTheta > halfpi) && (t2 >= 0.) ) // leave
1929 {
1930 if( calcNorm )
1931 {
1932 *validNorm = true;
1933 if (rho2 != 0.0)
1934 {
1935 rhoSecTheta = std::sqrt(rho2*(1+tanSTheta2));
1936
1937 *n = G4ThreeVector( p.x()/rhoSecTheta,
1938 p.y()/rhoSecTheta,
1939 std::sin(fSTheta) );
1940 }
1941 else { *n = G4ThreeVector(0.,0.,1.); }
1942 }
1943 return snxt = 0.;
1944 }
1945 if( (fSTheta < halfpi) && (t2 < 0.) && (p.z() >=0.) ) // leave
1946 {
1947 if( calcNorm ) { *validNorm = false; }
1948 return snxt = 0.;
1949 }
1950 }
1951 b = t2/t1;
1952 c = dist2STheta/t1;
1953 d2 = b*b - c ;
1954
1955 if ( d2 >= 0. )
1956 {
1957 d = std::sqrt(d2);
1958
1959 if( fSTheta > halfpi )
1960 {
1961 sd = -b - d; // First root
1962
1963 if ( ((std::fabs(s) < halfRmaxTolerance) && (t2 < 0.))
1964 || (sd < 0.) || ( (sd > 0.) && (p.z() + sd*v.z() > 0.) ) )
1965 {
1966 sd = -b + d ; // 2nd root
1967 }
1968 if( (sd > halfRmaxTolerance) && (p.z() + sd*v.z() <= 0.) )
1969 {
1970 stheta = sd;
1971 sidetheta = kSTheta;
1972 }
1973 }
1974 else // sTheta < pi/2, concave surface, no normal
1975 {
1976 sd = -b - d; // First root
1977
1978 if ( ( (std::fabs(sd) < halfRmaxTolerance) && (t2 >= 0.) )
1979 || (sd < 0.) || ( (sd > 0.) && (p.z() + sd*v.z() < 0.) ) )
1980 {
1981 sd = -b + d ; // 2nd root
1982 }
1983 if( (sd > halfRmaxTolerance) && (p.z() + sd*v.z() >= 0.) )
1984 {
1985 stheta = sd;
1986 sidetheta = kSTheta;
1987 }
1988 }
1989 }
1990 }
1991 }
1992 }
1993 if (eTheta < pi) // intersection with second cons
1994 {
1995 if( std::fabs(tanETheta) > 5./kAngTolerance ) // kons is plane z=0
1996 {
1997 if( v.z() < 0. )
1998 {
1999 if ( std::fabs( p.z() ) <= halfRmaxTolerance )
2000 {
2001 if(calcNorm)
2002 {
2003 *validNorm = true;
2004 *n = G4ThreeVector(0.,0.,-1.);
2005 }
2006 return snxt = 0 ;
2007 }
2008 sd = -p.z()/v.z();
2009
2010 if( sd < stheta )
2011 {
2012 stheta = sd;
2013 sidetheta = kETheta;
2014 }
2015 }
2016 }
2017 else // kons is not plane
2018 {
2019 t1 = 1-v.z()*v.z()*(1+tanETheta2);
2020 t2 = pDotV2d-p.z()*v.z()*tanETheta2; // ~vDotN if p on cons
2021 dist2ETheta = rho2-p.z()*p.z()*tanETheta2; // t3
2022
2023 distTheta = std::sqrt(rho2)-p.z()*tanETheta;
2024
2025 if( std::fabs(t1) < halfAngTolerance ) // 1st order equation,
2026 { // v parallel to kons
2027 if( v.z() < 0. )
2028 {
2029 if(std::fabs(distTheta) < halfRmaxTolerance) // p on surface
2030 {
2031 if( (eTheta > halfpi) && (p.z() < 0.) )
2032 {
2033 if( calcNorm ) { *validNorm = false; }
2034 return snxt = 0.;
2035 }
2036 if ( (eTheta < halfpi) && (p.z() >= 0) )
2037 {
2038 if( calcNorm )
2039 {
2040 *validNorm = true;
2041 if (rho2 != 0.0)
2042 {
2043 rhoSecTheta = std::sqrt(rho2*(1+tanETheta2));
2044 *n = G4ThreeVector( p.x()/rhoSecTheta,
2045 p.y()/rhoSecTheta,
2046 -sinETheta );
2047 }
2048 else { *n = G4ThreeVector(0.,0.,-1.); }
2049 }
2050 return snxt = 0.;
2051 }
2052 }
2053 sd = -0.5*dist2ETheta/t2;
2054
2055 if( sd < stheta )
2056 {
2057 stheta = sd;
2058 sidetheta = kETheta;
2059 }
2060 }
2061 } // 2nd order equation, 1st root of fSTheta cone
2062 else // 2nd if 1st root -ve
2063 {
2064 if ( std::fabs(distTheta) < halfRmaxTolerance )
2065 {
2066 if( (eTheta < halfpi) && (t2 >= 0.) ) // leave
2067 {
2068 if( calcNorm )
2069 {
2070 *validNorm = true;
2071 if (rho2 != 0.0)
2072 {
2073 rhoSecTheta = std::sqrt(rho2*(1+tanETheta2));
2074 *n = G4ThreeVector( p.x()/rhoSecTheta,
2075 p.y()/rhoSecTheta,
2076 -sinETheta );
2077 }
2078 else
2079 {
2080 *n = G4ThreeVector(0.,0.,-1.);
2081 }
2082 }
2083 return snxt = 0.;
2084 }
2085 if ( (eTheta > halfpi) && (t2 < 0.) && (p.z() <=0.) ) // leave
2086 {
2087 if( calcNorm ) { *validNorm = false; }
2088 return snxt = 0.;
2089 }
2090 }
2091 b = t2/t1;
2092 c = dist2ETheta/t1;
2093 d2 = b*b - c ;
2094 if ( (d2 <halfRmaxTolerance) && (d2 > -halfRmaxTolerance) )
2095 {
2096 d2 = 0.;
2097 }
2098 if ( d2 >= 0. )
2099 {
2100 d = std::sqrt(d2);
2101
2102 if( eTheta < halfpi )
2103 {
2104 sd = -b - d; // First root
2105
2106 if( ((std::fabs(sd) < halfRmaxTolerance) && (t2 < 0.))
2107 || (sd < 0.) )
2108 {
2109 sd = -b + d ; // 2nd root
2110 }
2111 if( sd > halfRmaxTolerance )
2112 {
2113 if( sd < stheta )
2114 {
2115 stheta = sd;
2116 sidetheta = kETheta;
2117 }
2118 }
2119 }
2120 else // sTheta+fDTheta > pi/2, concave surface, no normal
2121 {
2122 sd = -b - d; // First root
2123
2124 if ( ((std::fabs(sd) < halfRmaxTolerance) && (t2 >= 0.))
2125 || (sd < 0.)
2126 || ( (sd > 0.) && (p.z() + sd*v.z() > halfRmaxTolerance) ) )
2127 {
2128 sd = -b + d ; // 2nd root
2129 }
2130 if ( ( sd>halfRmaxTolerance )
2131 && ( p.z()+sd*v.z() <= halfRmaxTolerance ) )
2132 {
2133 if( sd < stheta )
2134 {
2135 stheta = sd;
2136 sidetheta = kETheta;
2137 }
2138 }
2139 }
2140 }
2141 }
2142 }
2143 }
2144
2145 } // end theta intersections
2146
2147 // Phi Intersection
2148
2149 if ( !fFullPhiSphere )
2150 {
2151 if ( (p.x() != 0.0) || (p.y() != 0.0) ) // Check if on z axis (rho not needed later)
2152 {
2153 // pDist -ve when inside
2154
2155 pDistS=p.x()*sinSPhi-p.y()*cosSPhi;
2156 pDistE=-p.x()*sinEPhi+p.y()*cosEPhi;
2157
2158 // Comp -ve when in direction of outwards normal
2159
2160 compS = -sinSPhi*v.x()+cosSPhi*v.y() ;
2161 compE = sinEPhi*v.x()-cosEPhi*v.y() ;
2162 sidephi = kNull ;
2163
2164 if ( (pDistS <= 0) && (pDistE <= 0) )
2165 {
2166 // Inside both phi *full* planes
2167
2168 if ( compS < 0 )
2169 {
2170 sphi = pDistS/compS ;
2171 xi = p.x()+sphi*v.x() ;
2172 yi = p.y()+sphi*v.y() ;
2173
2174 // Check intersection with correct half-plane (if not -> no intersect)
2175 //
2176 if( (std::fabs(xi)<=kCarTolerance) && (std::fabs(yi)<=kCarTolerance) )
2177 {
2178 vphi = std::atan2(v.y(),v.x());
2179 sidephi = kSPhi;
2180 if ( ( (fSPhi-halfAngTolerance) <= vphi)
2181 && ( (ePhi+halfAngTolerance) >= vphi) )
2182 {
2183 sphi = kInfinity;
2184 }
2185 }
2186 else if ( ( yi*cosCPhi - xi*sinCPhi ) >= 0 )
2187 {
2188 sphi=kInfinity;
2189 }
2190 else
2191 {
2192 sidephi = kSPhi ;
2193 if ( pDistS > -halfCarTolerance) { sphi = 0; } // Leave by sphi
2194 }
2195 }
2196 else { sphi = kInfinity; }
2197
2198 if ( compE < 0 )
2199 {
2200 sphi2=pDistE/compE ;
2201 if (sphi2 < sphi) // Only check further if < starting phi intersection
2202 {
2203 xi = p.x()+sphi2*v.x() ;
2204 yi = p.y()+sphi2*v.y() ;
2205
2206 // Check intersection with correct half-plane
2207 //
2208 if ( (std::fabs(xi)<=kCarTolerance)
2209 && (std::fabs(yi)<=kCarTolerance))
2210 {
2211 // Leaving via ending phi
2212 //
2213 vphi = std::atan2(v.y(),v.x()) ;
2214
2215 if( (fSPhi-halfAngTolerance > vphi)
2216 ||(fSPhi+fDPhi+halfAngTolerance < vphi) )
2217 {
2218 sidephi = kEPhi;
2219 if ( pDistE <= -halfCarTolerance ) { sphi = sphi2; }
2220 else { sphi = 0.0; }
2221 }
2222 }
2223 else if ((yi*cosCPhi-xi*sinCPhi)>=0) // Leaving via ending phi
2224 {
2225 sidephi = kEPhi ;
2226 if ( pDistE <= -halfCarTolerance )
2227 {
2228 sphi=sphi2;
2229 }
2230 else
2231 {
2232 sphi = 0 ;
2233 }
2234 }
2235 }
2236 }
2237 }
2238 else if ((pDistS >= 0) && (pDistE >= 0)) // Outside both *full* phi planes
2239 {
2240 if ( pDistS <= pDistE )
2241 {
2242 sidephi = kSPhi ;
2243 }
2244 else
2245 {
2246 sidephi = kEPhi ;
2247 }
2248 if ( fDPhi > pi )
2249 {
2250 if ( (compS < 0) && (compE < 0) ) { sphi = 0; }
2251 else { sphi = kInfinity; }
2252 }
2253 else
2254 {
2255 // if towards both >=0 then once inside (after error)
2256 // will remain inside
2257
2258 if ( (compS >= 0) && (compE >= 0) ) { sphi = kInfinity; }
2259 else { sphi = 0; }
2260 }
2261 }
2262 else if ( (pDistS > 0) && (pDistE < 0) )
2263 {
2264 // Outside full starting plane, inside full ending plane
2265
2266 if ( fDPhi > pi )
2267 {
2268 if ( compE < 0 )
2269 {
2270 sphi = pDistE/compE ;
2271 xi = p.x() + sphi*v.x() ;
2272 yi = p.y() + sphi*v.y() ;
2273
2274 // Check intersection in correct half-plane
2275 // (if not -> not leaving phi extent)
2276 //
2277 if( (std::fabs(xi)<=kCarTolerance)&&(std::fabs(yi)<=kCarTolerance) )
2278 {
2279 vphi = std::atan2(v.y(),v.x());
2280 sidephi = kSPhi;
2281 if ( ( (fSPhi-halfAngTolerance) <= vphi)
2282 && ( (ePhi+halfAngTolerance) >= vphi) )
2283 {
2284 sphi = kInfinity;
2285 }
2286 }
2287 else if ( ( yi*cosCPhi - xi*sinCPhi ) <= 0 )
2288 {
2289 sphi = kInfinity ;
2290 }
2291 else // Leaving via Ending phi
2292 {
2293 sidephi = kEPhi ;
2294 if ( pDistE > -halfCarTolerance ) { sphi = 0.; }
2295 }
2296 }
2297 else
2298 {
2299 sphi = kInfinity ;
2300 }
2301 }
2302 else
2303 {
2304 if ( compS >= 0 )
2305 {
2306 if ( compE < 0 )
2307 {
2308 sphi = pDistE/compE ;
2309 xi = p.x() + sphi*v.x() ;
2310 yi = p.y() + sphi*v.y() ;
2311
2312 // Check intersection in correct half-plane
2313 // (if not -> remain in extent)
2314 //
2315 if( (std::fabs(xi)<=kCarTolerance)
2316 && (std::fabs(yi)<=kCarTolerance) )
2317 {
2318 vphi = std::atan2(v.y(),v.x());
2319 sidephi = kSPhi;
2320 if ( ( (fSPhi-halfAngTolerance) <= vphi)
2321 && ( (ePhi+halfAngTolerance) >= vphi) )
2322 {
2323 sphi = kInfinity;
2324 }
2325 }
2326 else if ( ( yi*cosCPhi - xi*sinCPhi) <= 0 )
2327 {
2328 sphi=kInfinity;
2329 }
2330 else // otherwise leaving via Ending phi
2331 {
2332 sidephi = kEPhi ;
2333 }
2334 }
2335 else
2336 {
2337 sphi=kInfinity;
2338 }
2339 }
2340 else // leaving immediately by starting phi
2341 {
2342 sidephi = kSPhi ;
2343 sphi = 0 ;
2344 }
2345 }
2346 }
2347 else
2348 {
2349 // Must be pDistS < 0 && pDistE > 0
2350 // Inside full starting plane, outside full ending plane
2351
2352 if ( fDPhi > pi )
2353 {
2354 if ( compS < 0 )
2355 {
2356 sphi=pDistS/compS;
2357 xi=p.x()+sphi*v.x();
2358 yi=p.y()+sphi*v.y();
2359
2360 // Check intersection in correct half-plane
2361 // (if not -> not leaving phi extent)
2362 //
2363 if( (std::fabs(xi)<=kCarTolerance)&&(std::fabs(yi)<=kCarTolerance) )
2364 {
2365 vphi = std::atan2(v.y(),v.x()) ;
2366 sidephi = kSPhi;
2367 if ( ( (fSPhi-halfAngTolerance) <= vphi)
2368 && ( (ePhi+halfAngTolerance) >= vphi) )
2369 {
2370 sphi = kInfinity;
2371 }
2372 }
2373 else if ( ( yi*cosCPhi - xi*sinCPhi ) >= 0 )
2374 {
2375 sphi = kInfinity ;
2376 }
2377 else // Leaving via Starting phi
2378 {
2379 sidephi = kSPhi ;
2380 if ( pDistS > -halfCarTolerance ) { sphi = 0; }
2381 }
2382 }
2383 else
2384 {
2385 sphi = kInfinity ;
2386 }
2387 }
2388 else
2389 {
2390 if ( compE >= 0 )
2391 {
2392 if ( compS < 0 )
2393 {
2394 sphi = pDistS/compS ;
2395 xi = p.x()+sphi*v.x() ;
2396 yi = p.y()+sphi*v.y() ;
2397
2398 // Check intersection in correct half-plane
2399 // (if not -> remain in extent)
2400 //
2401 if( (std::fabs(xi)<=kCarTolerance)
2402 && (std::fabs(yi)<=kCarTolerance))
2403 {
2404 vphi = std::atan2(v.y(),v.x()) ;
2405 sidephi = kSPhi;
2406 if ( ( (fSPhi-halfAngTolerance) <= vphi)
2407 && ( (ePhi+halfAngTolerance) >= vphi) )
2408 {
2409 sphi = kInfinity;
2410 }
2411 }
2412 else if ( ( yi*cosCPhi - xi*sinCPhi ) >= 0 )
2413 {
2414 sphi = kInfinity ;
2415 }
2416 else // otherwise leaving via Starting phi
2417 {
2418 sidephi = kSPhi ;
2419 }
2420 }
2421 else
2422 {
2423 sphi = kInfinity ;
2424 }
2425 }
2426 else // leaving immediately by ending
2427 {
2428 sidephi = kEPhi ;
2429 sphi = 0 ;
2430 }
2431 }
2432 }
2433 }
2434 else
2435 {
2436 // On z axis + travel not || to z axis -> if phi of vector direction
2437 // within phi of shape, Step limited by rmax, else Step =0
2438
2439 if ( (v.x() != 0.0) || (v.y() != 0.0) )
2440 {
2441 vphi = std::atan2(v.y(),v.x()) ;
2442 if ((fSPhi-halfAngTolerance < vphi) && (vphi < ePhi+halfAngTolerance))
2443 {
2444 sphi = kInfinity;
2445 }
2446 else
2447 {
2448 sidephi = kSPhi ; // arbitrary
2449 sphi = 0 ;
2450 }
2451 }
2452 else // travel along z - no phi intersection
2453 {
2454 sphi = kInfinity ;
2455 }
2456 }
2457 if ( sphi < snxt ) // Order intersecttions
2458 {
2459 snxt = sphi ;
2460 side = sidephi ;
2461 }
2462 }
2463 if (stheta < snxt ) // Order intersections
2464 {
2465 snxt = stheta ;
2466 side = sidetheta ;
2467 }
2468
2469 if (calcNorm) // Output switch operator
2470 {
2471 switch( side )
2472 {
2473 case kRMax:
2474 xi=p.x()+snxt*v.x();
2475 yi=p.y()+snxt*v.y();
2476 zi=p.z()+snxt*v.z();
2477 *n=G4ThreeVector(xi/fRmax,yi/fRmax,zi/fRmax);
2478 *validNorm=true;
2479 break;
2480
2481 case kRMin:
2482 *validNorm=false; // Rmin is concave
2483 break;
2484
2485 case kSPhi:
2486 if ( fDPhi <= pi ) // Normal to Phi-
2487 {
2488 *n=G4ThreeVector(sinSPhi,-cosSPhi,0);
2489 *validNorm=true;
2490 }
2491 else { *validNorm=false; }
2492 break ;
2493
2494 case kEPhi:
2495 if ( fDPhi <= pi ) // Normal to Phi+
2496 {
2497 *n=G4ThreeVector(-sinEPhi,cosEPhi,0);
2498 *validNorm=true;
2499 }
2500 else { *validNorm=false; }
2501 break;
2502
2503 case kSTheta:
2504 if( fSTheta == halfpi )
2505 {
2506 *n=G4ThreeVector(0.,0.,1.);
2507 *validNorm=true;
2508 }
2509 else if ( fSTheta > halfpi )
2510 {
2511 xi = p.x() + snxt*v.x();
2512 yi = p.y() + snxt*v.y();
2513 rho2=xi*xi+yi*yi;
2514 if (rho2 != 0.0)
2515 {
2516 rhoSecTheta = std::sqrt(rho2*(1+tanSTheta2));
2517 *n = G4ThreeVector( xi/rhoSecTheta, yi/rhoSecTheta,
2518 -tanSTheta/std::sqrt(1+tanSTheta2));
2519 }
2520 else
2521 {
2522 *n = G4ThreeVector(0.,0.,1.);
2523 }
2524 *validNorm=true;
2525 }
2526 else { *validNorm=false; } // Concave STheta cone
2527 break;
2528
2529 case kETheta:
2530 if( eTheta == halfpi )
2531 {
2532 *n = G4ThreeVector(0.,0.,-1.);
2533 *validNorm = true;
2534 }
2535 else if ( eTheta < halfpi )
2536 {
2537 xi=p.x()+snxt*v.x();
2538 yi=p.y()+snxt*v.y();
2539 rho2=xi*xi+yi*yi;
2540 if (rho2 != 0.0)
2541 {
2542 rhoSecTheta = std::sqrt(rho2*(1+tanETheta2));
2543 *n = G4ThreeVector( xi/rhoSecTheta, yi/rhoSecTheta,
2544 -tanETheta/std::sqrt(1+tanETheta2) );
2545 }
2546 else
2547 {
2548 *n = G4ThreeVector(0.,0.,-1.);
2549 }
2550 *validNorm=true;
2551 }
2552 else { *validNorm=false; } // Concave ETheta cone
2553 break;
2554
2555 default:
2556 G4cout << G4endl;
2557 DumpInfo();
2558 std::ostringstream message;
2559 G4long oldprc = message.precision(16);
2560 message << "Undefined side for valid surface normal to solid."
2561 << G4endl
2562 << "Position:" << G4endl << G4endl
2563 << "p.x() = " << p.x()/mm << " mm" << G4endl
2564 << "p.y() = " << p.y()/mm << " mm" << G4endl
2565 << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl
2566 << "Direction:" << G4endl << G4endl
2567 << "v.x() = " << v.x() << G4endl
2568 << "v.y() = " << v.y() << G4endl
2569 << "v.z() = " << v.z() << G4endl << G4endl
2570 << "Proposed distance :" << G4endl << G4endl
2571 << "snxt = " << snxt/mm << " mm" << G4endl;
2572 message.precision(oldprc);
2573 G4Exception("G4Sphere::DistanceToOut(p,v,..)",
2574 "GeomSolids1002", JustWarning, message);
2575 break;
2576 }
2577 }
2578 if (snxt == kInfinity)
2579 {
2580 G4cout << G4endl;
2581 DumpInfo();
2582 std::ostringstream message;
2583 G4long oldprc = message.precision(16);
2584 message << "Logic error: snxt = kInfinity ???" << G4endl
2585 << "Position:" << G4endl << G4endl
2586 << "p.x() = " << p.x()/mm << " mm" << G4endl
2587 << "p.y() = " << p.y()/mm << " mm" << G4endl
2588 << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl
2589 << "Rp = "<< std::sqrt( p.x()*p.x()+p.y()*p.y()+p.z()*p.z() )/mm
2590 << " mm" << G4endl << G4endl
2591 << "Direction:" << G4endl << G4endl
2592 << "v.x() = " << v.x() << G4endl
2593 << "v.y() = " << v.y() << G4endl
2594 << "v.z() = " << v.z() << G4endl << G4endl
2595 << "Proposed distance :" << G4endl << G4endl
2596 << "snxt = " << snxt/mm << " mm" << G4endl;
2597 message.precision(oldprc);
2598 G4Exception("G4Sphere::DistanceToOut(p,v,..)",
2599 "GeomSolids1002", JustWarning, message);
2600 }
2601
2602 return snxt;
2603}

◆ GetCosEndPhi()

G4double G4Sphere::GetCosEndPhi ( ) const
inline

Referenced by BoundingLimits().

◆ GetCosEndTheta()

G4double G4Sphere::GetCosEndTheta ( ) const
inline

Referenced by BoundingLimits().

◆ GetCosStartPhi()

G4double G4Sphere::GetCosStartPhi ( ) const
inline

Referenced by BoundingLimits().

◆ GetCosStartTheta()

G4double G4Sphere::GetCosStartTheta ( ) const
inline

Referenced by BoundingLimits().

◆ GetCubicVolume()

G4double G4Sphere::GetCubicVolume ( )
overridevirtual

Returning an estimation of the solid volume (capacity) and surface area, in internal units.

Reimplemented from G4VSolid.

Definition at line 2744 of file G4Sphere.cc.

2745{
2746 if (fCubicVolume == 0)
2747 {
2748 G4AutoLock l(&sphereMutex);
2749 G4double RRR = fRmax*fRmax*fRmax;
2750 G4double rrr = fRmin*fRmin*fRmin;
2751 fCubicVolume = fDPhi*(cosSTheta - cosETheta)*(RRR - rrr)/3.;
2752 l.unlock();
2753 }
2754 return fCubicVolume;
2755}
G4AutoLock
G4TemplateAutoLock< G4Mutex > G4AutoLock
Definition G4AutoLock.hh:582
G4CSGSolid::fCubicVolume
G4double fCubicVolume
Definition G4CSGSolid.hh:92

◆ GetDeltaPhiAngle()

G4double G4Sphere::GetDeltaPhiAngle ( ) const
inline

Referenced by BoundingLimits(), G4tgbGeometryDumper::GetSolidParams(), G4GDMLWriteParamvol::Sphere_dimensionsWrite(), and G4GDMLWriteSolids::SphereWrite().

◆ GetDeltaThetaAngle()

G4double G4Sphere::GetDeltaThetaAngle ( ) const
inline

Referenced by BoundingLimits(), G4tgbGeometryDumper::GetSolidParams(), G4GDMLWriteParamvol::Sphere_dimensionsWrite(), and G4GDMLWriteSolids::SphereWrite().

◆ GetEntityType()

G4GeometryType G4Sphere::GetEntityType ( ) const
overridevirtual

Returns the type ID, "G4Sphere" of the solid.

Implements G4VSolid.

Definition at line 2702 of file G4Sphere.cc.

2703{
2704 return {"G4Sphere"};
2705}

◆ GetExtent()

G4VisExtent G4Sphere::GetExtent ( ) const
overridevirtual

Methods for creating graphical representations (i.e. for visualisation).

Reimplemented from G4VSolid.

Definition at line 2832 of file G4Sphere.cc.

2833{
2834 return { -fRmax, fRmax,-fRmax, fRmax,-fRmax, fRmax };
2835}

◆ GetInnerRadius()

G4double G4Sphere::GetInnerRadius ( ) const
inline

Accessors.

Referenced by BoundingLimits(), G4tgbGeometryDumper::GetSolidParams(), G4PSSphereSurfaceCurrent::IsSelectedSurface(), G4PSSphereSurfaceFlux::IsSelectedSurface(), G4GDMLWriteParamvol::Sphere_dimensionsWrite(), and G4GDMLWriteSolids::SphereWrite().

◆ GetOuterRadius()

G4double G4Sphere::GetOuterRadius ( ) const
inline

Referenced by BoundingLimits(), G4tgbGeometryDumper::GetSolidParams(), G4GDMLWriteParamvol::Sphere_dimensionsWrite(), and G4GDMLWriteSolids::SphereWrite().

◆ GetPointOnSurface()

G4ThreeVector G4Sphere::GetPointOnSurface ( ) const
overridevirtual

Returns a random point located and uniformly distributed on the surface of the solid.

Reimplemented from G4VSolid.

Definition at line 2781 of file G4Sphere.cc.

2782{
2783 G4double RR = fRmax*fRmax;
2784 G4double rr = fRmin*fRmin;
2785
2786 // Find surface areas
2787 //
2788 G4double aInner = fDPhi*rr*(cosSTheta - cosETheta);
2789 G4double aOuter = fDPhi*RR*(cosSTheta - cosETheta);
2790 G4double aPhi = (!fFullPhiSphere) ? fDTheta*(RR - rr) : 0.;
2791 G4double aSTheta = (fSTheta > 0) ? 0.5*fDPhi*(RR - rr)*sinSTheta : 0.;
2792 G4double aETheta = (eTheta < pi) ? 0.5*fDPhi*(RR - rr)*sinETheta : 0.;
2793 G4double aTotal = aInner + aOuter + aPhi + aSTheta + aETheta;
2794
2795 // Select surface and generate a point
2796 //
2797 G4double select = aTotal*G4QuickRand();
2798 G4double u = G4QuickRand();
2799 G4double v = G4QuickRand();
2800 if (select < aInner + aOuter) // lateral surface
2801 {
2802 G4double r = (select < aInner) ? fRmin : fRmax;
2803 G4double z = cosSTheta + (cosETheta - cosSTheta)*u;
2804 G4double rho = std::sqrt(1. - z*z);
2805 G4double phi = fDPhi*v + fSPhi;
2806 return { r*rho*std::cos(phi), r*rho*std::sin(phi), r*z };
2807 }
2808 if (select < aInner + aOuter + aPhi) // cut in phi
2809 {
2810 G4double phi = (select < aInner + aOuter + 0.5*aPhi) ? fSPhi : fSPhi + fDPhi;
2811 G4double r = std::sqrt((RR - rr)*u + rr);
2812 G4double theta = fDTheta*v + fSTheta;
2813 G4double z = std::cos(theta);
2814 G4double rho = std::sin(theta);
2815 return { r*rho*std::cos(phi), r*rho*std::sin(phi), r*z };
2816 }
2817
2818 // cut in theta
2819
2820 G4double theta = (select < aTotal - aETheta) ? fSTheta : fSTheta + fDTheta;
2821 G4double r = std::sqrt((RR - rr)*u + rr);
2822 G4double phi = fDPhi*v + fSPhi;
2823 G4double z = std::cos(theta);
2824 G4double rho = std::sin(theta);
2825 return { r*rho*std::cos(phi), r*rho*std::sin(phi), r*z };
2826}
G4QuickRand
G4double G4QuickRand(uint32_t seed=0)
Definition G4QuickRand.hh:34

◆ GetSinEndPhi()

G4double G4Sphere::GetSinEndPhi ( ) const
inline

Referenced by BoundingLimits().

◆ GetSinEndTheta()

G4double G4Sphere::GetSinEndTheta ( ) const
inline

Referenced by BoundingLimits().

◆ GetSinStartPhi()

G4double G4Sphere::GetSinStartPhi ( ) const
inline

Referenced by BoundingLimits().

◆ GetSinStartTheta()

G4double G4Sphere::GetSinStartTheta ( ) const
inline

Referenced by BoundingLimits().

◆ GetStartPhiAngle()

G4double G4Sphere::GetStartPhiAngle ( ) const
inline

Referenced by G4tgbGeometryDumper::GetSolidParams(), G4GDMLWriteParamvol::Sphere_dimensionsWrite(), and G4GDMLWriteSolids::SphereWrite().

◆ GetStartThetaAngle()

G4double G4Sphere::GetStartThetaAngle ( ) const
inline

Referenced by BoundingLimits(), G4tgbGeometryDumper::GetSolidParams(), G4GDMLWriteParamvol::Sphere_dimensionsWrite(), and G4GDMLWriteSolids::SphereWrite().

◆ GetSurfaceArea()

G4double G4Sphere::GetSurfaceArea ( )
overridevirtual

Returns an estimation of the solid surface area in internal units. This method may be overloaded by derived classes to compute the exact geometrical quantity for solids where this is possible, or anyway to cache the computed value. Note: the computed value is NOT cached.

Reimplemented from G4VSolid.

Definition at line 2761 of file G4Sphere.cc.

2762{
2763 if (fSurfaceArea == 0)
2764 {
2765 G4AutoLock l(&sphereMutex);
2766 G4double RR = fRmax*fRmax;
2767 G4double rr = fRmin*fRmin;
2768 fSurfaceArea = fDPhi*(RR + rr)*(cosSTheta - cosETheta);
2769 if (!fFullPhiSphere) { fSurfaceArea += fDTheta*(RR - rr); }
2770 if (fSTheta > 0) { fSurfaceArea += 0.5*fDPhi*(RR - rr)*sinSTheta; }
2771 if (eTheta < CLHEP::pi) { fSurfaceArea += 0.5*fDPhi*(RR - rr)*sinETheta; }
2772 l.unlock();
2773 }
2774 return fSurfaceArea;
2775}
G4CSGSolid::fSurfaceArea
G4double fSurfaceArea
Definition G4CSGSolid.hh:93

◆ Inside()

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

Concrete implementations of the expected query interfaces for solids, as defined in the base class G4VSolid.

Implements G4VSolid.

Definition at line 255 of file G4Sphere.cc.

256{
257 G4double rho,rho2,rad2,tolRMin,tolRMax;
258 G4double pPhi,pTheta;
259 EInside in = kOutside;
260
261 const G4double halfRmaxTolerance = fRmaxTolerance*0.5;
262 const G4double halfRminTolerance = fRminTolerance*0.5;
263 const G4double Rmax_minus = fRmax - halfRmaxTolerance;
264 const G4double Rmin_plus = (fRmin > 0) ? fRmin+halfRminTolerance : 0;
265
266 rho2 = p.x()*p.x() + p.y()*p.y() ;
267 rad2 = rho2 + p.z()*p.z() ;
268
269 // Check radial surfaces. Sets 'in'
270
271 tolRMin = Rmin_plus;
272 tolRMax = Rmax_minus;
273
274 if(rad2 == 0.0)
275 {
276 if (fRmin > 0.0)
277 {
278 return in = kOutside;
279 }
280 if ( (!fFullPhiSphere) || (!fFullThetaSphere) )
281 {
282 return in = kSurface;
283 }
284
285 return in = kInside;
286 }
287
288 if ( (rad2 <= Rmax_minus*Rmax_minus) && (rad2 >= Rmin_plus*Rmin_plus) )
289 {
290 in = kInside;
291 }
292 else
293 {
294 tolRMax = fRmax + halfRmaxTolerance; // outside case
295 tolRMin = std::max(fRmin-halfRminTolerance, 0.); // outside case
296 if ( (rad2 <= tolRMax*tolRMax) && (rad2 >= tolRMin*tolRMin) )
297 {
298 in = kSurface;
299 }
300 else
301 {
302 return in = kOutside;
303 }
304 }
305
306 // Phi boundaries : Do not check if it has no phi boundary!
307
308 if ( !fFullPhiSphere && (rho2 != 0.0) ) // [fDPhi < twopi] and [p.x or p.y]
309 {
310 pPhi = std::atan2(p.y(),p.x()) ;
311
312 if ( pPhi < fSPhi - halfAngTolerance ) { pPhi += twopi; }
313 else if ( pPhi > ePhi + halfAngTolerance ) { pPhi -= twopi; }
314
315 if ( (pPhi < fSPhi - halfAngTolerance)
316 || (pPhi > ePhi + halfAngTolerance) ) { return in = kOutside; }
317
318 if (in == kInside) // else it's kSurface anyway already
319 {
320 if ( (pPhi < fSPhi + halfAngTolerance)
321 || (pPhi > ePhi - halfAngTolerance) ) { in = kSurface; }
322 }
323 }
324
325 // Theta bondaries
326
327 if ( ((rho2 != 0.0) || (p.z() != 0.0)) && (!fFullThetaSphere) )
328 {
329 rho = std::sqrt(rho2);
330 pTheta = std::atan2(rho,p.z());
331
332 if ( in == kInside )
333 {
334 if ( ((fSTheta > 0.0) && (pTheta < fSTheta + halfAngTolerance))
335 || ((eTheta < pi) && (pTheta > eTheta - halfAngTolerance)) )
336 {
337 if ( (( (fSTheta>0.0)&&(pTheta>=fSTheta-halfAngTolerance) )
338 || (fSTheta == 0.0) )
339 && ((eTheta==pi)||(pTheta <= eTheta + halfAngTolerance) ) )
340 {
341 in = kSurface;
342 }
343 else
344 {
345 in = kOutside;
346 }
347 }
348 }
349 else
350 {
351 if ( ((fSTheta > 0.0)&&(pTheta < fSTheta - halfAngTolerance))
352 ||((eTheta < pi )&&(pTheta > eTheta + halfAngTolerance)) )
353 {
354 in = kOutside;
355 }
356 }
357 }
358 return in;
359}
EInside
EInside
Definition geomdefs.hh:67
kInside
@ kInside
Definition geomdefs.hh:70
kSurface
@ kSurface
Definition geomdefs.hh:69

Referenced by DistanceToOut().

◆ operator=()

G4Sphere & G4Sphere::operator= ( const G4Sphere & rhs)

Definition at line 127 of file G4Sphere.cc.

128{
129 // Check assignment to self
130 //
131 if (this == &rhs) { return *this; }
132
133 // Copy base class data
134 //
135 G4CSGSolid::operator=(rhs);
136
137 // Copy data
138 //
139 fRminTolerance = rhs.fRminTolerance; fRmaxTolerance = rhs.fRmaxTolerance;
140 kAngTolerance = rhs.kAngTolerance; kRadTolerance = rhs.kRadTolerance;
141 fEpsilon = rhs.fEpsilon; fRmin = rhs.fRmin; fRmax = rhs.fRmax;
142 fSPhi = rhs.fSPhi; fDPhi = rhs.fDPhi; fSTheta = rhs.fSTheta;
143 fDTheta = rhs.fDTheta; sinCPhi = rhs.sinCPhi; cosCPhi = rhs.cosCPhi;
144 cosHDPhi = rhs.cosHDPhi;
145 cosHDPhiOT = rhs.cosHDPhiOT; cosHDPhiIT = rhs.cosHDPhiIT;
146 sinSPhi = rhs.sinSPhi; cosSPhi = rhs.cosSPhi;
147 sinEPhi = rhs.sinEPhi; cosEPhi = rhs.cosEPhi;
148 hDPhi = rhs.hDPhi; cPhi = rhs.cPhi; ePhi = rhs.ePhi;
149 sinSTheta = rhs.sinSTheta; cosSTheta = rhs.cosSTheta;
150 sinETheta = rhs.sinETheta; cosETheta = rhs.cosETheta;
151 tanSTheta = rhs.tanSTheta; tanSTheta2 = rhs.tanSTheta2;
152 tanETheta = rhs.tanETheta; tanETheta2 = rhs.tanETheta2;
153 eTheta = rhs.eTheta; fFullPhiSphere = rhs.fFullPhiSphere;
154 fFullThetaSphere = rhs.fFullThetaSphere; fFullSphere = rhs.fFullSphere;
155 halfCarTolerance = rhs.halfCarTolerance;
156 halfAngTolerance = rhs.halfAngTolerance;
157
158 return *this;
159}
G4CSGSolid::operator=
G4CSGSolid & operator=(const G4CSGSolid &rhs)
Definition G4CSGSolid.cc:89

◆ SetDeltaPhiAngle()

void G4Sphere::SetDeltaPhiAngle ( G4double newDphi)
inline

◆ SetDeltaThetaAngle()

void G4Sphere::SetDeltaThetaAngle ( G4double newDTheta)
inline

◆ SetInnerRadius()

void G4Sphere::SetInnerRadius ( G4double newRMin)
inline

Modifiers.

◆ SetOuterRadius()

void G4Sphere::SetOuterRadius ( G4double newRmax)
inline

◆ SetStartPhiAngle()

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

◆ SetStartThetaAngle()

void G4Sphere::SetStartThetaAngle ( G4double newSTheta)
inline

◆ StreamInfo()

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

Streams the object contents to an output stream.

Reimplemented from G4CSGSolid.

Definition at line 2720 of file G4Sphere.cc.

2721{
2722 G4long oldprc = os.precision(16);
2723 os << "-----------------------------------------------------------\n"
2724 << " *** Dump for solid - " << GetName() << " ***\n"
2725 << " ===================================================\n"
2726 << " Solid type: G4Sphere\n"
2727 << " Parameters: \n"
2728 << " inner radius: " << fRmin/mm << " mm \n"
2729 << " outer radius: " << fRmax/mm << " mm \n"
2730 << " starting phi of segment : " << fSPhi/degree << " degrees \n"
2731 << " delta phi of segment : " << fDPhi/degree << " degrees \n"
2732 << " starting theta of segment: " << fSTheta/degree << " degrees \n"
2733 << " delta theta of segment : " << fDTheta/degree << " degrees \n"
2734 << "-----------------------------------------------------------\n";
2735 os.precision(oldprc);
2736
2737 return os;
2738}

◆ SurfaceNormal()

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

Returns the outwards pointing unit normal of the shape for the surface closest to the point at offset 'p'.

Parameters
[in]pThe point at offset p.
Returns
The outwards pointing unit normal.

Implements G4VSolid.

Definition at line 367 of file G4Sphere.cc.

368{
369 G4int noSurfaces = 0;
370 G4double rho, rho2, radius, pTheta, pPhi=0.;
371 G4double distRMin = kInfinity;
372 G4double distSPhi = kInfinity, distEPhi = kInfinity;
373 G4double distSTheta = kInfinity, distETheta = kInfinity;
374 G4ThreeVector nR, nPs, nPe, nTs, nTe, nZ(0.,0.,1.);
375 G4ThreeVector norm, sumnorm(0.,0.,0.);
376
377 rho2 = p.x()*p.x()+p.y()*p.y();
378 radius = std::sqrt(rho2+p.z()*p.z());
379 rho = std::sqrt(rho2);
380
381 G4double distRMax = std::fabs(radius-fRmax);
382 if (fRmin != 0.0) { distRMin = std::fabs(radius-fRmin); }
383
384 if ( (rho != 0.0) && !fFullSphere )
385 {
386 pPhi = std::atan2(p.y(),p.x());
387
388 if (pPhi < fSPhi-halfAngTolerance) { pPhi += twopi; }
389 else if (pPhi > ePhi+halfAngTolerance) { pPhi -= twopi; }
390 }
391 if ( !fFullPhiSphere )
392 {
393 if ( rho != 0.0 )
394 {
395 distSPhi = std::fabs( pPhi-fSPhi );
396 distEPhi = std::fabs( pPhi-ePhi );
397 }
398 else if( fRmin == 0.0 )
399 {
400 distSPhi = 0.;
401 distEPhi = 0.;
402 }
403 nPs = G4ThreeVector(sinSPhi,-cosSPhi,0);
404 nPe = G4ThreeVector(-sinEPhi,cosEPhi,0);
405 }
406 if ( !fFullThetaSphere )
407 {
408 if ( rho != 0.0 )
409 {
410 pTheta = std::atan2(rho,p.z());
411 distSTheta = std::fabs(pTheta-fSTheta);
412 distETheta = std::fabs(pTheta-eTheta);
413
414 nTs = G4ThreeVector(-cosSTheta*p.x()/rho,
415 -cosSTheta*p.y()/rho,
416 sinSTheta );
417
418 nTe = G4ThreeVector( cosETheta*p.x()/rho,
419 cosETheta*p.y()/rho,
420 -sinETheta );
421 }
422 else if( fRmin == 0.0 )
423 {
424 if ( fSTheta != 0.0 )
425 {
426 distSTheta = 0.;
427 nTs = G4ThreeVector(0.,0.,-1.);
428 }
429 if ( eTheta < pi )
430 {
431 distETheta = 0.;
432 nTe = G4ThreeVector(0.,0.,1.);
433 }
434 }
435 }
436 if( radius != 0.0 ) { nR = G4ThreeVector(p.x()/radius,p.y()/radius,p.z()/radius); }
437
438 if( distRMax <= halfCarTolerance )
439 {
440 ++noSurfaces;
441 sumnorm += nR;
442 }
443 if( (fRmin != 0.0) && (distRMin <= halfCarTolerance) )
444 {
445 ++noSurfaces;
446 sumnorm -= nR;
447 }
448 if( !fFullPhiSphere )
449 {
450 if (distSPhi <= halfAngTolerance)
451 {
452 ++noSurfaces;
453 sumnorm += nPs;
454 }
455 if (distEPhi <= halfAngTolerance)
456 {
457 ++noSurfaces;
458 sumnorm += nPe;
459 }
460 }
461 if ( !fFullThetaSphere )
462 {
463 if ((distSTheta <= halfAngTolerance) && (fSTheta > 0.))
464 {
465 ++noSurfaces;
466 if ((radius <= halfCarTolerance) && fFullPhiSphere) { sumnorm += nZ; }
467 else { sumnorm += nTs; }
468 }
469 if ((distETheta <= halfAngTolerance) && (eTheta < pi))
470 {
471 ++noSurfaces;
472 if ((radius <= halfCarTolerance) && fFullPhiSphere) { sumnorm -= nZ; }
473 else { sumnorm += nTe; }
474 if(sumnorm.z() == 0.) { sumnorm += nZ; }
475 }
476 }
477 if ( noSurfaces == 0 )
478 {
479#ifdef G4CSGDEBUG
480 G4Exception("G4Sphere::SurfaceNormal(p)", "GeomSolids1002",
481 JustWarning, "Point p is not on surface !?" );
482#endif
483 norm = ApproxSurfaceNormal(p);
484 }
485 else if ( noSurfaces == 1 ) { norm = sumnorm; }
486 else { norm = sumnorm.unit(); }
487 return norm;
488}
G4int
int G4int
Definition G4Types.hh:85
CLHEP::Hep3Vector::unit
Hep3Vector unit() const
The documentation for this class was generated from the following files: