Geant4 11.2.2
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
Public Member Functions | List of all members
G4Sphere Class Reference

#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
 
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)
 
G4Sphere & operator= (const G4Sphere &rhs)
 
- Public Member Functions inherited from G4CSGSolid
 G4CSGSolid (const G4String &pName)
 
virtual ~G4CSGSolid ()
 
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 ()
 
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
 

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

Definition at line 80 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 )

Definition at line 72 of file G4Sphere.cc.

76 : G4CSGSolid(pName), fSPhi(0.0), fFullPhiSphere(true), fFullThetaSphere(true)
77{
78 kAngTolerance = G4GeometryTolerance::GetInstance()->GetAngularTolerance();
79 kRadTolerance = G4GeometryTolerance::GetInstance()->GetRadialTolerance();
80
81 halfCarTolerance = 0.5*kCarTolerance;
82 halfAngTolerance = 0.5*kAngTolerance;
83
84 // Check radii and set radial tolerances
85
86 if ( (pRmin >= pRmax) || (pRmax < 1.1*kRadTolerance) || (pRmin < 0) )
87 {
88 std::ostringstream message;
89 message << "Invalid radii for Solid: " << GetName() << G4endl
90 << " pRmin = " << pRmin << ", pRmax = " << pRmax;
91 G4Exception("G4Sphere::G4Sphere()", "GeomSolids0002",
92 FatalException, message);
93 }
94 fRmin=pRmin; fRmax=pRmax;
95 fRminTolerance = (fRmin) != 0.0 ? std::max( kRadTolerance, fEpsilon*fRmin ) : 0;
96 fRmaxTolerance = std::max( kRadTolerance, fEpsilon*fRmax );
97
98 // Check angles
99
100 CheckPhiAngles(pSPhi, pDPhi);
101 CheckThetaAngles(pSTheta, pDTheta);
102}
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:299

Referenced by Clone().

◆ ~G4Sphere()

G4Sphere::~G4Sphere ( )
overridedefault

◆ G4Sphere() [2/3]

G4Sphere::G4Sphere ( __void__ & a)

Definition at line 109 of file G4Sphere.cc.

110 : G4CSGSolid(a)
111{
112}

◆ G4Sphere() [3/3]

G4Sphere::G4Sphere ( const G4Sphere & rhs)
default

Member Function Documentation

◆ BoundingLimits()

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

Reimplemented from G4VSolid.

Definition at line 180 of file G4Sphere.cc.

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

Implements G4VSolid.

Definition at line 237 of file G4Sphere.cc.

241{
242 G4ThreeVector bmin, bmax;
243
244 // Get bounding box
245 BoundingLimits(bmin,bmax);
246
247 // Find extent
248 G4BoundingEnvelope bbox(bmin,bmax);
249 return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
250}
G4Sphere::BoundingLimits
void BoundingLimits(G4ThreeVector &pMin, G4ThreeVector &pMax) const override
Definition G4Sphere.cc:180

◆ Clone()

G4VSolid * G4Sphere::Clone ( ) const
overridevirtual

Reimplemented from G4VSolid.

Definition at line 2712 of file G4Sphere.cc.

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

◆ ComputeDimensions()

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

Reimplemented from G4VSolid.

Definition at line 169 of file G4Sphere.cc.

172{
173 p->ComputeDimensions(*this,n,pRep);
174}
G4VPVParameterisation::ComputeDimensions
virtual void ComputeDimensions(G4Box &, const G4int, const G4VPhysicalVolume *) const
Definition G4VPVParameterisation.hh:89

◆ CreatePolyhedron()

G4Polyhedron * G4Sphere::CreatePolyhedron ( ) const
overridevirtual

Reimplemented from G4VSolid.

Definition at line 2840 of file G4Sphere.cc.

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

◆ DescribeYourselfTo()

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

Implements G4VSolid.

Definition at line 2835 of file G4Sphere.cc.

2836{
2837 scene.AddSolid (*this);
2838}
G4VGraphicsScene::AddSolid
virtual void AddSolid(const G4Box &)=0

◆ DistanceToIn() [1/2]

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

Implements G4VSolid.

Definition at line 1623 of file G4Sphere.cc.

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

◆ DistanceToIn() [2/2]

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

Implements G4VSolid.

Definition at line 674 of file G4Sphere.cc.

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

Referenced by DistanceToIn().

◆ DistanceToOut() [1/2]

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

Implements G4VSolid.

Definition at line 2610 of file G4Sphere.cc.

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

Implements G4VSolid.

Definition at line 1719 of file G4Sphere.cc.

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

◆ 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

Reimplemented from G4VSolid.

Definition at line 2745 of file G4Sphere.cc.

2746{
2747 if (fCubicVolume == 0.)
2748 {
2749 G4double RRR = fRmax*fRmax*fRmax;
2750 G4double rrr = fRmin*fRmin*fRmin;
2751 fCubicVolume = fDPhi*(cosSTheta - cosETheta)*(RRR - rrr)/3.;
2752 }
2753 return fCubicVolume;
2754}
G4CSGSolid::fCubicVolume
G4double fCubicVolume
Definition G4CSGSolid.hh:68

◆ 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

Implements G4VSolid.

Definition at line 2703 of file G4Sphere.cc.

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

◆ GetExtent()

G4VisExtent G4Sphere::GetExtent ( ) const
overridevirtual

Reimplemented from G4VSolid.

Definition at line 2829 of file G4Sphere.cc.

2830{
2831 return { -fRmax, fRmax,-fRmax, fRmax,-fRmax, fRmax };
2832}

◆ GetInnerRadius()

G4double G4Sphere::GetInnerRadius ( ) const
inline

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

Reimplemented from G4VSolid.

Definition at line 2778 of file G4Sphere.cc.

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

Reimplemented from G4VSolid.

Definition at line 2760 of file G4Sphere.cc.

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

◆ Inside()

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

Implements G4VSolid.

Definition at line 258 of file G4Sphere.cc.

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

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

◆ 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

Reimplemented from G4CSGSolid.

Definition at line 2721 of file G4Sphere.cc.

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

◆ SurfaceNormal()

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

Implements G4VSolid.

Definition at line 372 of file G4Sphere.cc.

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