Geant4 9.6.0
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 ()
 
G4double GetInnerRadius () const
 
G4double GetOuterRadius () const
 
G4double GetStartPhiAngle () const
 
G4double GetDeltaPhiAngle () const
 
G4double GetStartThetaAngle () const
 
G4double GetDeltaThetaAngle () const
 
void SetInnerRadius (G4double newRMin)
 
void SetOuterRadius (G4double newRmax)
 
void SetStartPhiAngle (G4double newSphi, G4bool trig=true)
 
void SetDeltaPhiAngle (G4double newDphi)
 
void SetStartThetaAngle (G4double newSTheta)
 
void SetDeltaThetaAngle (G4double newDTheta)
 
G4double GetCubicVolume ()
 
G4double GetSurfaceArea ()
 
void ComputeDimensions (G4VPVParameterisation *p, const G4int n, const G4VPhysicalVolume *pRep)
 
G4bool CalculateExtent (const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pmin, G4double &pmax) const
 
EInside Inside (const G4ThreeVector &p) const
 
G4ThreeVector SurfaceNormal (const G4ThreeVector &p) const
 
G4double DistanceToIn (const G4ThreeVector &p, const G4ThreeVector &v) const
 
G4double DistanceToIn (const G4ThreeVector &p) const
 
G4double DistanceToOut (const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=G4bool(false), G4bool *validNorm=0, G4ThreeVector *n=0) const
 
G4double DistanceToOut (const G4ThreeVector &p) const
 
G4GeometryType GetEntityType () const
 
G4ThreeVector GetPointOnSurface () const
 
G4VSolid * Clone () const
 
std::ostream & StreamInfo (std::ostream &os) const
 
G4VisExtent GetExtent () const
 
void DescribeYourselfTo (G4VGraphicsScene &scene) const
 
G4Polyhedron * CreatePolyhedron () const
 
G4NURBS * CreateNURBS () const
 
 G4Sphere (__void__ &)
 
 G4Sphere (const G4Sphere &rhs)
 
G4Sphere & operator= (const G4Sphere &rhs)
 
G4double GetRmin () const
 
G4double GetRmax () const
 
G4double GetSPhi () const
 
G4double GetDPhi () const
 
G4double GetSTheta () const
 
G4double GetDTheta () const
 
G4double GetInsideRadius () const
 
void SetInsideRadius (G4double newRmin)
 
- Public Member Functions inherited from G4CSGSolid
 G4CSGSolid (const G4String &pName)
 
virtual ~G4CSGSolid ()
 
virtual std::ostream & StreamInfo (std::ostream &os) const
 
virtual G4Polyhedron * GetPolyhedron () const
 
 G4CSGSolid (__void__ &)
 
 G4CSGSolid (const G4CSGSolid &rhs)
 
G4CSGSolid & operator= (const G4CSGSolid &rhs)
 
- Public Member Functions inherited from G4VSolid
 G4VSolid (const G4String &name)
 
virtual ~G4VSolid ()
 
G4bool operator== (const G4VSolid &s) const
 
G4String GetName () const
 
void SetName (const G4String &name)
 
G4double GetTolerance () const
 
virtual G4bool CalculateExtent (const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pMin, G4double &pMax) const =0
 
virtual EInside Inside (const G4ThreeVector &p) const =0
 
virtual G4ThreeVector SurfaceNormal (const G4ThreeVector &p) const =0
 
virtual G4double DistanceToIn (const G4ThreeVector &p, const G4ThreeVector &v) const =0
 
virtual G4double DistanceToIn (const G4ThreeVector &p) const =0
 
virtual G4double DistanceToOut (const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=false, G4bool *validNorm=0, G4ThreeVector *n=0) const =0
 
virtual G4double DistanceToOut (const G4ThreeVector &p) const =0
 
virtual void ComputeDimensions (G4VPVParameterisation *p, const G4int n, const G4VPhysicalVolume *pRep)
 
virtual G4double GetCubicVolume ()
 
virtual G4double GetSurfaceArea ()
 
virtual G4GeometryType GetEntityType () const =0
 
virtual G4ThreeVector GetPointOnSurface () const
 
virtual G4VSolid * Clone () const
 
virtual std::ostream & StreamInfo (std::ostream &os) const =0
 
void DumpInfo () const
 
virtual void DescribeYourselfTo (G4VGraphicsScene &scene) const =0
 
virtual G4VisExtent GetExtent () const
 
virtual G4Polyhedron * CreatePolyhedron () const
 
virtual G4NURBS * CreateNURBS () const
 
virtual G4Polyhedron * GetPolyhedron () const
 
virtual const G4VSolid * GetConstituentSolid (G4int no) const
 
virtual G4VSolid * GetConstituentSolid (G4int no)
 
virtual const G4DisplacedSolid * GetDisplacedSolidPtr () const
 
virtual G4DisplacedSolid * GetDisplacedSolidPtr ()
 
 G4VSolid (__void__ &)
 
 G4VSolid (const G4VSolid &rhs)
 
G4VSolid & operator= (const G4VSolid &rhs)
 

Additional Inherited Members

- Protected Member Functions inherited from G4CSGSolid
G4double GetRadiusInRing (G4double rmin, G4double rmax) const
 
- Protected Member Functions inherited from G4VSolid
void CalculateClippedPolygonExtent (G4ThreeVectorList &pPolygon, const G4VoxelLimits &pVoxelLimit, const EAxis pAxis, G4double &pMin, G4double &pMax) const
 
void ClipCrossSection (G4ThreeVectorList *pVertices, const G4int pSectionIndex, const G4VoxelLimits &pVoxelLimit, const EAxis pAxis, G4double &pMin, G4double &pMax) const
 
void ClipBetweenSections (G4ThreeVectorList *pVertices, const G4int pSectionIndex, const G4VoxelLimits &pVoxelLimit, const EAxis pAxis, G4double &pMin, G4double &pMax) const
 
void ClipPolygon (G4ThreeVectorList &pPolygon, const G4VoxelLimits &pVoxelLimit, const EAxis pAxis) const
 
G4double EstimateCubicVolume (G4int nStat, G4double epsilon) const
 
G4double EstimateSurfaceArea (G4int nStat, G4double ell) const
 
- Protected Attributes inherited from G4CSGSolid
G4double fCubicVolume
 
G4double fSurfaceArea
 
G4Polyhedron * fpPolyhedron
 
- Protected Attributes inherited from G4VSolid
G4double kCarTolerance
 

Detailed Description

Definition at line 77 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 90 of file G4Sphere.cc.

94 : G4CSGSolid(pName), fEpsilon(2.e-11),
95 fFullPhiSphere(true), fFullThetaSphere(true)
96{
97 kAngTolerance = G4GeometryTolerance::GetInstance()->GetAngularTolerance();
98
99 // Check radii and set radial tolerances
100
101 kRadTolerance = G4GeometryTolerance::GetInstance()->GetRadialTolerance();
102 if ( (pRmin >= pRmax) || (pRmax < 1.1*kRadTolerance) || (pRmin < 0) )
103 {
104 std::ostringstream message;
105 message << "Invalid radii for Solid: " << GetName() << G4endl
106 << " pRmin = " << pRmin << ", pRmax = " << pRmax;
107 G4Exception("G4Sphere::G4Sphere()", "GeomSolids0002",
108 FatalException, message);
109 }
110 fRmin=pRmin; fRmax=pRmax;
111 fRminTolerance = (fRmin) ? std::max( kRadTolerance, fEpsilon*fRmin ) : 0;
112 fRmaxTolerance = std::max( kRadTolerance, fEpsilon*fRmax );
113
114 // Check angles
115
116 CheckPhiAngles(pSPhi, pDPhi);
117 CheckThetaAngles(pSTheta, pDTheta);
118}
FatalException
@ FatalException
Definition: G4ExceptionSeverity.hh:60
G4endl
#define G4endl
Definition: G4ios.hh:52
G4GeometryTolerance::GetRadialTolerance
G4double GetRadialTolerance() const
Definition: G4GeometryTolerance.cc:97
G4GeometryTolerance::GetInstance
static G4GeometryTolerance * GetInstance()
Definition: G4GeometryTolerance.cc:73
G4GeometryTolerance::GetAngularTolerance
G4double GetAngularTolerance() const
Definition: G4GeometryTolerance.cc:92
G4VSolid::GetName
G4String GetName() const
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41

◆ ~G4Sphere()

G4Sphere::~G4Sphere ( )

Definition at line 141 of file G4Sphere.cc.

142{
143}

◆ G4Sphere() [2/3]

G4Sphere::G4Sphere ( __void__ &  a)

Definition at line 125 of file G4Sphere.cc.

126 : G4CSGSolid(a), fRminTolerance(0.), fRmaxTolerance(0.),
127 kAngTolerance(0.), kRadTolerance(0.), fEpsilon(0.),
128 fRmin(0.), fRmax(0.), fSPhi(0.), fDPhi(0.), fSTheta(0.),
129 fDTheta(0.), sinCPhi(0.), cosCPhi(0.), cosHDPhiOT(0.), cosHDPhiIT(0.),
130 sinSPhi(0.), cosSPhi(0.), sinEPhi(0.), cosEPhi(0.), hDPhi(0.), cPhi(0.),
131 ePhi(0.), sinSTheta(0.), cosSTheta(0.), sinETheta(0.), cosETheta(0.),
132 tanSTheta(0.), tanSTheta2(0.), tanETheta(0.), tanETheta2(0.), eTheta(0.),
133 fFullPhiSphere(false), fFullThetaSphere(false), fFullSphere(true)
134{
135}

◆ G4Sphere() [3/3]

G4Sphere::G4Sphere ( const G4Sphere &  rhs)

Definition at line 149 of file G4Sphere.cc.

150 : G4CSGSolid(rhs), fRminTolerance(rhs.fRminTolerance),
151 fRmaxTolerance(rhs.fRmaxTolerance), kAngTolerance(rhs.kAngTolerance),
152 kRadTolerance(rhs.kRadTolerance), fEpsilon(rhs.fEpsilon),
153 fRmin(rhs.fRmin), fRmax(rhs.fRmax), fSPhi(rhs.fSPhi), fDPhi(rhs.fDPhi),
154 fSTheta(rhs.fSTheta), fDTheta(rhs.fDTheta),
155 sinCPhi(rhs.sinCPhi), cosCPhi(rhs.cosCPhi),
156 cosHDPhiOT(rhs.cosHDPhiOT), cosHDPhiIT(rhs.cosHDPhiIT),
157 sinSPhi(rhs.sinSPhi), cosSPhi(rhs.cosSPhi),
158 sinEPhi(rhs.sinEPhi), cosEPhi(rhs.cosEPhi),
159 hDPhi(rhs.hDPhi), cPhi(rhs.cPhi), ePhi(rhs.ePhi),
160 sinSTheta(rhs.sinSTheta), cosSTheta(rhs.cosSTheta),
161 sinETheta(rhs.sinETheta), cosETheta(rhs.cosETheta),
162 tanSTheta(rhs.tanSTheta), tanSTheta2(rhs.tanSTheta2),
163 tanETheta(rhs.tanETheta), tanETheta2(rhs.tanETheta2), eTheta(rhs.eTheta),
164 fFullPhiSphere(rhs.fFullPhiSphere), fFullThetaSphere(rhs.fFullThetaSphere),
165 fFullSphere(rhs.fFullSphere)
166{
167}

Member Function Documentation

◆ CalculateExtent()

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

Implements G4VSolid.

Definition at line 220 of file G4Sphere.cc.

224{
225 if ( fFullSphere )
226 {
227 // Special case handling for solid spheres-shells
228 // (rotation doesn't influence).
229 // Compute x/y/z mins and maxs for bounding box respecting limits,
230 // with early returns if outside limits. Then switch() on pAxis,
231 // and compute exact x and y limit for x/y case
232
233 G4double xoffset,xMin,xMax;
234 G4double yoffset,yMin,yMax;
235 G4double zoffset,zMin,zMax;
236
237 G4double diff1,diff2,delta,maxDiff,newMin,newMax;
238 G4double xoff1,xoff2,yoff1,yoff2;
239
240 xoffset=pTransform.NetTranslation().x();
241 xMin=xoffset-fRmax;
242 xMax=xoffset+fRmax;
243 if (pVoxelLimit.IsXLimited())
244 {
245 if ( (xMin>pVoxelLimit.GetMaxXExtent()+kCarTolerance)
246 || (xMax<pVoxelLimit.GetMinXExtent()-kCarTolerance) )
247 {
248 return false;
249 }
250 else
251 {
252 if (xMin<pVoxelLimit.GetMinXExtent())
253 {
254 xMin=pVoxelLimit.GetMinXExtent();
255 }
256 if (xMax>pVoxelLimit.GetMaxXExtent())
257 {
258 xMax=pVoxelLimit.GetMaxXExtent();
259 }
260 }
261 }
262
263 yoffset=pTransform.NetTranslation().y();
264 yMin=yoffset-fRmax;
265 yMax=yoffset+fRmax;
266 if (pVoxelLimit.IsYLimited())
267 {
268 if ( (yMin>pVoxelLimit.GetMaxYExtent()+kCarTolerance)
269 || (yMax<pVoxelLimit.GetMinYExtent()-kCarTolerance) )
270 {
271 return false;
272 }
273 else
274 {
275 if (yMin<pVoxelLimit.GetMinYExtent())
276 {
277 yMin=pVoxelLimit.GetMinYExtent();
278 }
279 if (yMax>pVoxelLimit.GetMaxYExtent())
280 {
281 yMax=pVoxelLimit.GetMaxYExtent();
282 }
283 }
284 }
285
286 zoffset=pTransform.NetTranslation().z();
287 zMin=zoffset-fRmax;
288 zMax=zoffset+fRmax;
289 if (pVoxelLimit.IsZLimited())
290 {
291 if ( (zMin>pVoxelLimit.GetMaxZExtent()+kCarTolerance)
292 || (zMax<pVoxelLimit.GetMinZExtent()-kCarTolerance) )
293 {
294 return false;
295 }
296 else
297 {
298 if (zMin<pVoxelLimit.GetMinZExtent())
299 {
300 zMin=pVoxelLimit.GetMinZExtent();
301 }
302 if (zMax>pVoxelLimit.GetMaxZExtent())
303 {
304 zMax=pVoxelLimit.GetMaxZExtent();
305 }
306 }
307 }
308
309 // Known to cut sphere
310
311 switch (pAxis)
312 {
313 case kXAxis:
314 yoff1=yoffset-yMin;
315 yoff2=yMax-yoffset;
316 if ((yoff1>=0) && (yoff2>=0))
317 {
318 // Y limits cross max/min x => no change
319 //
320 pMin=xMin;
321 pMax=xMax;
322 }
323 else
324 {
325 // Y limits don't cross max/min x => compute max delta x,
326 // hence new mins/maxs
327 //
328 delta=fRmax*fRmax-yoff1*yoff1;
329 diff1=(delta>0.) ? std::sqrt(delta) : 0.;
330 delta=fRmax*fRmax-yoff2*yoff2;
331 diff2=(delta>0.) ? std::sqrt(delta) : 0.;
332 maxDiff=(diff1>diff2) ? diff1:diff2;
333 newMin=xoffset-maxDiff;
334 newMax=xoffset+maxDiff;
335 pMin=(newMin<xMin) ? xMin : newMin;
336 pMax=(newMax>xMax) ? xMax : newMax;
337 }
338 break;
339 case kYAxis:
340 xoff1=xoffset-xMin;
341 xoff2=xMax-xoffset;
342 if ((xoff1>=0) && (xoff2>=0))
343 {
344 // X limits cross max/min y => no change
345 //
346 pMin=yMin;
347 pMax=yMax;
348 }
349 else
350 {
351 // X limits don't cross max/min y => compute max delta y,
352 // hence new mins/maxs
353 //
354 delta=fRmax*fRmax-xoff1*xoff1;
355 diff1=(delta>0.) ? std::sqrt(delta) : 0.;
356 delta=fRmax*fRmax-xoff2*xoff2;
357 diff2=(delta>0.) ? std::sqrt(delta) : 0.;
358 maxDiff=(diff1>diff2) ? diff1:diff2;
359 newMin=yoffset-maxDiff;
360 newMax=yoffset+maxDiff;
361 pMin=(newMin<yMin) ? yMin : newMin;
362 pMax=(newMax>yMax) ? yMax : newMax;
363 }
364 break;
365 case kZAxis:
366 pMin=zMin;
367 pMax=zMax;
368 break;
369 default:
370 break;
371 }
372 pMin-=kCarTolerance;
373 pMax+=kCarTolerance;
374
375 return true;
376 }
377 else // Transformed cutted sphere
378 {
379 G4int i,j,noEntries,noBetweenSections;
380 G4bool existsAfterClip=false;
381
382 // Calculate rotated vertex coordinates
383
384 G4ThreeVectorList* vertices;
385 G4int noPolygonVertices ;
386 vertices=CreateRotatedVertices(pTransform,noPolygonVertices);
387
388 pMin=+kInfinity;
389 pMax=-kInfinity;
390
391 noEntries=vertices->size(); // noPolygonVertices*noPhiCrossSections
392 noBetweenSections=noEntries-noPolygonVertices;
393
394 G4ThreeVectorList ThetaPolygon ;
395 for (i=0;i<noEntries;i+=noPolygonVertices)
396 {
397 for(j=0;j<(noPolygonVertices/2)-1;j++)
398 {
399 ThetaPolygon.push_back((*vertices)[i+j]) ;
400 ThetaPolygon.push_back((*vertices)[i+j+1]) ;
401 ThetaPolygon.push_back((*vertices)[i+noPolygonVertices-2-j]) ;
402 ThetaPolygon.push_back((*vertices)[i+noPolygonVertices-1-j]) ;
403 CalculateClippedPolygonExtent(ThetaPolygon,pVoxelLimit,pAxis,pMin,pMax);
404 ThetaPolygon.clear() ;
405 }
406 }
407 for (i=0;i<noBetweenSections;i+=noPolygonVertices)
408 {
409 for(j=0;j<noPolygonVertices-1;j++)
410 {
411 ThetaPolygon.push_back((*vertices)[i+j]) ;
412 ThetaPolygon.push_back((*vertices)[i+j+1]) ;
413 ThetaPolygon.push_back((*vertices)[i+noPolygonVertices+j+1]) ;
414 ThetaPolygon.push_back((*vertices)[i+noPolygonVertices+j]) ;
415 CalculateClippedPolygonExtent(ThetaPolygon,pVoxelLimit,pAxis,pMin,pMax);
416 ThetaPolygon.clear() ;
417 }
418 ThetaPolygon.push_back((*vertices)[i+noPolygonVertices-1]) ;
419 ThetaPolygon.push_back((*vertices)[i]) ;
420 ThetaPolygon.push_back((*vertices)[i+noPolygonVertices]) ;
421 ThetaPolygon.push_back((*vertices)[i+2*noPolygonVertices-1]) ;
422 CalculateClippedPolygonExtent(ThetaPolygon,pVoxelLimit,pAxis,pMin,pMax);
423 ThetaPolygon.clear() ;
424 }
425
426 if ((pMin!=kInfinity) || (pMax!=-kInfinity))
427 {
428 existsAfterClip=true;
429
430 // Add 2*tolerance to avoid precision troubles
431 //
432 pMin-=kCarTolerance;
433 pMax+=kCarTolerance;
434 }
435 else
436 {
437 // Check for case where completely enveloping clipping volume
438 // If point inside then we are confident that the solid completely
439 // envelopes the clipping volume. Hence set min/max extents according
440 // to clipping volume extents along the specified axis.
441
442 G4ThreeVector clipCentre(
443 (pVoxelLimit.GetMinXExtent()+pVoxelLimit.GetMaxXExtent())*0.5,
444 (pVoxelLimit.GetMinYExtent()+pVoxelLimit.GetMaxYExtent())*0.5,
445 (pVoxelLimit.GetMinZExtent()+pVoxelLimit.GetMaxZExtent())*0.5);
446
447 if (Inside(pTransform.Inverse().TransformPoint(clipCentre))!=kOutside)
448 {
449 existsAfterClip=true;
450 pMin=pVoxelLimit.GetMinExtent(pAxis);
451 pMax=pVoxelLimit.GetMaxExtent(pAxis);
452 }
453 }
454 delete vertices;
455 return existsAfterClip;
456 }
457}
G4double
double G4double
Definition: G4Types.hh:64
G4int
int G4int
Definition: G4Types.hh:66
G4bool
bool G4bool
Definition: G4Types.hh:67
G4ThreeVectorList
std::vector< G4ThreeVector > G4ThreeVectorList
Definition: G4VSolid.hh:85
CLHEP::Hep3Vector::z
double z() const
CLHEP::Hep3Vector::x
double x() const
CLHEP::Hep3Vector::y
double y() const
G4AffineTransform::Inverse
G4AffineTransform Inverse() const
G4AffineTransform::NetTranslation
G4ThreeVector NetTranslation() const
G4AffineTransform::TransformPoint
G4ThreeVector TransformPoint(const G4ThreeVector &vec) const
G4Sphere::Inside
EInside Inside(const G4ThreeVector &p) const
Definition: G4Sphere.cc:465
G4VSolid::kCarTolerance
G4double kCarTolerance
Definition: G4VSolid.hh:307
G4VSolid::CalculateClippedPolygonExtent
void CalculateClippedPolygonExtent(G4ThreeVectorList &pPolygon, const G4VoxelLimits &pVoxelLimit, const EAxis pAxis, G4double &pMin, G4double &pMax) const
Definition: G4VSolid.cc:425
G4VoxelLimits::GetMinExtent
G4double GetMinExtent(const EAxis pAxis) const
G4VoxelLimits::IsYLimited
G4bool IsYLimited() const
G4VoxelLimits::GetMinZExtent
G4double GetMinZExtent() const
G4VoxelLimits::IsXLimited
G4bool IsXLimited() const
G4VoxelLimits::GetMaxExtent
G4double GetMaxExtent(const EAxis pAxis) const
G4VoxelLimits::GetMaxYExtent
G4double GetMaxYExtent() const
G4VoxelLimits::GetMaxZExtent
G4double GetMaxZExtent() const
G4VoxelLimits::GetMinYExtent
G4double GetMinYExtent() const
G4VoxelLimits::GetMinXExtent
G4double GetMinXExtent() const
G4VoxelLimits::IsZLimited
G4bool IsZLimited() const
G4VoxelLimits::GetMaxXExtent
G4double GetMaxXExtent() const
kYAxis
@ kYAxis
Definition: geomdefs.hh:54
kXAxis
@ kXAxis
Definition: geomdefs.hh:54
kZAxis
@ kZAxis
Definition: geomdefs.hh:54
kOutside
@ kOutside
Definition: geomdefs.hh:58

◆ Clone()

G4VSolid * G4Sphere::Clone ( ) const
virtual

Reimplemented from G4VSolid.

Definition at line 3009 of file G4Sphere.cc.

3010{
3011 return new G4Sphere(*this);
3012}

◆ ComputeDimensions()

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

Reimplemented from G4VSolid.

Definition at line 209 of file G4Sphere.cc.

212{
213 p->ComputeDimensions(*this,n,pRep);
214}
G4VPVParameterisation::ComputeDimensions
virtual void ComputeDimensions(G4Box &, const G4int, const G4VPhysicalVolume *) const
Definition: G4VPVParameterisation.hh:94

◆ CreateNURBS()

G4NURBS * G4Sphere::CreateNURBS ( ) const
virtual

Reimplemented from G4VSolid.

Definition at line 3187 of file G4Sphere.cc.

3188{
3189 return new G4NURBSbox (fRmax, fRmax, fRmax); // Box for now!!!
3190}

◆ CreatePolyhedron()

G4Polyhedron * G4Sphere::CreatePolyhedron ( ) const
virtual

Reimplemented from G4VSolid.

Definition at line 3182 of file G4Sphere.cc.

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

◆ DescribeYourselfTo()

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

Implements G4VSolid.

Definition at line 3177 of file G4Sphere.cc.

3178{
3179 scene.AddSolid (*this);
3180}
G4VGraphicsScene::AddSolid
virtual void AddSolid(const G4Box &)=0

◆ DistanceToIn() [1/2]

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

Implements G4VSolid.

Definition at line 1801 of file G4Sphere.cc.

1802{
1803 G4double safe=0.0,safeRMin,safeRMax,safePhi,safeTheta;
1804 G4double rho2,rds,rho;
1805 G4double cosPsi;
1806 G4double pTheta,dTheta1,dTheta2;
1807 rho2=p.x()*p.x()+p.y()*p.y();
1808 rds=std::sqrt(rho2+p.z()*p.z());
1809 rho=std::sqrt(rho2);
1810
1811 //
1812 // Distance to r shells
1813 //
1814 if (fRmin)
1815 {
1816 safeRMin=fRmin-rds;
1817 safeRMax=rds-fRmax;
1818 if (safeRMin>safeRMax)
1819 {
1820 safe=safeRMin;
1821 }
1822 else
1823 {
1824 safe=safeRMax;
1825 }
1826 }
1827 else
1828 {
1829 safe=rds-fRmax;
1830 }
1831
1832 //
1833 // Distance to phi extent
1834 //
1835 if (!fFullPhiSphere && rho)
1836 {
1837 // Psi=angle from central phi to point
1838 //
1839 cosPsi=(p.x()*cosCPhi+p.y()*sinCPhi)/rho;
1840 if (cosPsi<std::cos(hDPhi))
1841 {
1842 // Point lies outside phi range
1843 //
1844 if ((p.y()*cosCPhi-p.x()*sinCPhi)<=0)
1845 {
1846 safePhi=std::fabs(p.x()*sinSPhi-p.y()*cosSPhi);
1847 }
1848 else
1849 {
1850 safePhi=std::fabs(p.x()*sinEPhi-p.y()*cosEPhi);
1851 }
1852 if (safePhi>safe) { safe=safePhi; }
1853 }
1854 }
1855 //
1856 // Distance to Theta extent
1857 //
1858 if ((rds!=0.0) && (!fFullThetaSphere))
1859 {
1860 pTheta=std::acos(p.z()/rds);
1861 if (pTheta<0) { pTheta+=pi; }
1862 dTheta1=fSTheta-pTheta;
1863 dTheta2=pTheta-eTheta;
1864 if (dTheta1>dTheta2)
1865 {
1866 if (dTheta1>=0) // WHY ???????????
1867 {
1868 safeTheta=rds*std::sin(dTheta1);
1869 if (safe<=safeTheta)
1870 {
1871 safe=safeTheta;
1872 }
1873 }
1874 }
1875 else
1876 {
1877 if (dTheta2>=0)
1878 {
1879 safeTheta=rds*std::sin(dTheta2);
1880 if (safe<=safeTheta)
1881 {
1882 safe=safeTheta;
1883 }
1884 }
1885 }
1886 }
1887
1888 if (safe<0) { safe=0; }
1889 return safe;
1890}
G4INCL::Math::pi
const G4double pi
Definition: G4INCLGlobals.hh:70

◆ DistanceToIn() [2/2]

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

Implements G4VSolid.

Definition at line 867 of file G4Sphere.cc.

869{
870 G4double snxt = kInfinity ; // snxt = default return value
871 G4double rho2, rad2, pDotV2d, pDotV3d, pTheta ;
872 G4double tolSTheta=0., tolETheta=0. ;
873 const G4double dRmax = 100.*fRmax;
874
875 static const G4double halfCarTolerance = kCarTolerance*0.5;
876 static const G4double halfAngTolerance = kAngTolerance*0.5;
877 const G4double halfRmaxTolerance = fRmaxTolerance*0.5;
878 const G4double halfRminTolerance = fRminTolerance*0.5;
879 const G4double tolORMin2 = (fRmin>halfRminTolerance)
880 ? (fRmin-halfRminTolerance)*(fRmin-halfRminTolerance) : 0;
881 const G4double tolIRMin2 =
882 (fRmin+halfRminTolerance)*(fRmin+halfRminTolerance);
883 const G4double tolORMax2 =
884 (fRmax+halfRmaxTolerance)*(fRmax+halfRmaxTolerance);
885 const G4double tolIRMax2 =
886 (fRmax-halfRmaxTolerance)*(fRmax-halfRmaxTolerance);
887
888 // Intersection point
889 //
890 G4double xi, yi, zi, rhoi, rhoi2, radi2, iTheta ;
891
892 // Phi intersection
893 //
894 G4double Comp ;
895
896 // Phi precalcs
897 //
898 G4double Dist, cosPsi ;
899
900 // Theta precalcs
901 //
902 G4double dist2STheta, dist2ETheta ;
903 G4double t1, t2, b, c, d2, d, sd = kInfinity ;
904
905 // General Precalcs
906 //
907 rho2 = p.x()*p.x() + p.y()*p.y() ;
908 rad2 = rho2 + p.z()*p.z() ;
909 pTheta = std::atan2(std::sqrt(rho2),p.z()) ;
910
911 pDotV2d = p.x()*v.x() + p.y()*v.y() ;
912 pDotV3d = pDotV2d + p.z()*v.z() ;
913
914 // Theta precalcs
915 //
916 if (!fFullThetaSphere)
917 {
918 tolSTheta = fSTheta - halfAngTolerance ;
919 tolETheta = eTheta + halfAngTolerance ;
920 }
921
922 // Outer spherical shell intersection
923 // - Only if outside tolerant fRmax
924 // - Check for if inside and outer G4Sphere heading through solid (-> 0)
925 // - No intersect -> no intersection with G4Sphere
926 //
927 // Shell eqn: x^2+y^2+z^2=RSPH^2
928 //
929 // => (px+svx)^2+(py+svy)^2+(pz+svz)^2=R^2
930 //
931 // => (px^2+py^2+pz^2) +2sd(pxvx+pyvy+pzvz)+sd^2(vx^2+vy^2+vz^2)=R^2
932 // => rad2 +2sd(pDotV3d) +sd^2 =R^2
933 //
934 // => sd=-pDotV3d+-std::sqrt(pDotV3d^2-(rad2-R^2))
935
936 c = rad2 - fRmax*fRmax ;
937
938 if (c > fRmaxTolerance*fRmax)
939 {
940 // If outside tolerant boundary of outer G4Sphere
941 // [should be std::sqrt(rad2)-fRmax > halfRmaxTolerance]
942
943 d2 = pDotV3d*pDotV3d - c ;
944
945 if ( d2 >= 0 )
946 {
947 sd = -pDotV3d - std::sqrt(d2) ;
948
949 if (sd >= 0 )
950 {
951 if ( sd>dRmax ) // Avoid rounding errors due to precision issues seen on
952 { // 64 bits systems. Split long distances and recompute
953 G4double fTerm = sd-std::fmod(sd,dRmax);
954 sd = fTerm + DistanceToIn(p+fTerm*v,v);
955 }
956 xi = p.x() + sd*v.x() ;
957 yi = p.y() + sd*v.y() ;
958 rhoi = std::sqrt(xi*xi + yi*yi) ;
959
960 if (!fFullPhiSphere && rhoi) // Check phi intersection
961 {
962 cosPsi = (xi*cosCPhi + yi*sinCPhi)/rhoi ;
963
964 if (cosPsi >= cosHDPhiOT)
965 {
966 if (!fFullThetaSphere) // Check theta intersection
967 {
968 zi = p.z() + sd*v.z() ;
969
970 // rhoi & zi can never both be 0
971 // (=>intersect at origin =>fRmax=0)
972 //
973 iTheta = std::atan2(rhoi,zi) ;
974 if ( (iTheta >= tolSTheta) && (iTheta <= tolETheta) )
975 {
976 return snxt = sd ;
977 }
978 }
979 else
980 {
981 return snxt=sd;
982 }
983 }
984 }
985 else
986 {
987 if (!fFullThetaSphere) // Check theta intersection
988 {
989 zi = p.z() + sd*v.z() ;
990
991 // rhoi & zi can never both be 0
992 // (=>intersect at origin => fRmax=0 !)
993 //
994 iTheta = std::atan2(rhoi,zi) ;
995 if ( (iTheta >= tolSTheta) && (iTheta <= tolETheta) )
996 {
997 return snxt=sd;
998 }
999 }
1000 else
1001 {
1002 return snxt = sd;
1003 }
1004 }
1005 }
1006 }
1007 else // No intersection with G4Sphere
1008 {
1009 return snxt=kInfinity;
1010 }
1011 }
1012 else
1013 {
1014 // Inside outer radius
1015 // check not inside, and heading through G4Sphere (-> 0 to in)
1016
1017 d2 = pDotV3d*pDotV3d - c ;
1018
1019 if ( (rad2 > tolIRMax2)
1020 && ( (d2 >= fRmaxTolerance*fRmax) && (pDotV3d < 0) ) )
1021 {
1022 if (!fFullPhiSphere)
1023 {
1024 // Use inner phi tolerant boundary -> if on tolerant
1025 // phi boundaries, phi intersect code handles leaving/entering checks
1026
1027 cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/std::sqrt(rho2) ;
1028
1029 if (cosPsi>=cosHDPhiIT)
1030 {
1031 // inside radii, delta r -ve, inside phi
1032
1033 if ( !fFullThetaSphere )
1034 {
1035 if ( (pTheta >= tolSTheta + kAngTolerance)
1036 && (pTheta <= tolETheta - kAngTolerance) )
1037 {
1038 return snxt=0;
1039 }
1040 }
1041 else // strictly inside Theta in both cases
1042 {
1043 return snxt=0;
1044 }
1045 }
1046 }
1047 else
1048 {
1049 if ( !fFullThetaSphere )
1050 {
1051 if ( (pTheta >= tolSTheta + kAngTolerance)
1052 && (pTheta <= tolETheta - kAngTolerance) )
1053 {
1054 return snxt=0;
1055 }
1056 }
1057 else // strictly inside Theta in both cases
1058 {
1059 return snxt=0;
1060 }
1061 }
1062 }
1063 }
1064
1065 // Inner spherical shell intersection
1066 // - Always farthest root, because would have passed through outer
1067 // surface first.
1068 // - Tolerant check if travelling through solid
1069
1070 if (fRmin)
1071 {
1072 c = rad2 - fRmin*fRmin ;
1073 d2 = pDotV3d*pDotV3d - c ;
1074
1075 // Within tolerance inner radius of inner G4Sphere
1076 // Check for immediate entry/already inside and travelling outwards
1077
1078 if ( (c > -halfRminTolerance) && (rad2 < tolIRMin2)
1079 && ( (d2 < fRmin*kCarTolerance) || (pDotV3d >= 0) ) )
1080 {
1081 if ( !fFullPhiSphere )
1082 {
1083 // Use inner phi tolerant boundary -> if on tolerant
1084 // phi boundaries, phi intersect code handles leaving/entering checks
1085
1086 cosPsi = (p.x()*cosCPhi+p.y()*sinCPhi)/std::sqrt(rho2) ;
1087 if (cosPsi >= cosHDPhiIT)
1088 {
1089 // inside radii, delta r -ve, inside phi
1090 //
1091 if ( !fFullThetaSphere )
1092 {
1093 if ( (pTheta >= tolSTheta + kAngTolerance)
1094 && (pTheta <= tolETheta - kAngTolerance) )
1095 {
1096 return snxt=0;
1097 }
1098 }
1099 else
1100 {
1101 return snxt = 0 ;
1102 }
1103 }
1104 }
1105 else
1106 {
1107 if ( !fFullThetaSphere )
1108 {
1109 if ( (pTheta >= tolSTheta + kAngTolerance)
1110 && (pTheta <= tolETheta - kAngTolerance) )
1111 {
1112 return snxt = 0 ;
1113 }
1114 }
1115 else
1116 {
1117 return snxt=0;
1118 }
1119 }
1120 }
1121 else // Not special tolerant case
1122 {
1123 if (d2 >= 0)
1124 {
1125 sd = -pDotV3d + std::sqrt(d2) ;
1126 if ( sd >= halfRminTolerance ) // It was >= 0 ??
1127 {
1128 xi = p.x() + sd*v.x() ;
1129 yi = p.y() + sd*v.y() ;
1130 rhoi = std::sqrt(xi*xi+yi*yi) ;
1131
1132 if ( !fFullPhiSphere && rhoi ) // Check phi intersection
1133 {
1134 cosPsi = (xi*cosCPhi + yi*sinCPhi)/rhoi ;
1135
1136 if (cosPsi >= cosHDPhiOT)
1137 {
1138 if ( !fFullThetaSphere ) // Check theta intersection
1139 {
1140 zi = p.z() + sd*v.z() ;
1141
1142 // rhoi & zi can never both be 0
1143 // (=>intersect at origin =>fRmax=0)
1144 //
1145 iTheta = std::atan2(rhoi,zi) ;
1146 if ( (iTheta >= tolSTheta) && (iTheta<=tolETheta) )
1147 {
1148 snxt = sd;
1149 }
1150 }
1151 else
1152 {
1153 snxt=sd;
1154 }
1155 }
1156 }
1157 else
1158 {
1159 if ( !fFullThetaSphere ) // Check theta intersection
1160 {
1161 zi = p.z() + sd*v.z() ;
1162
1163 // rhoi & zi can never both be 0
1164 // (=>intersect at origin => fRmax=0 !)
1165 //
1166 iTheta = std::atan2(rhoi,zi) ;
1167 if ( (iTheta >= tolSTheta) && (iTheta <= tolETheta) )
1168 {
1169 snxt = sd;
1170 }
1171 }
1172 else
1173 {
1174 snxt = sd;
1175 }
1176 }
1177 }
1178 }
1179 }
1180 }
1181
1182 // Phi segment intersection
1183 //
1184 // o Tolerant of points inside phi planes by up to kCarTolerance*0.5
1185 //
1186 // o NOTE: Large duplication of code between sphi & ephi checks
1187 // -> only diffs: sphi -> ephi, Comp -> -Comp and half-plane
1188 // intersection check <=0 -> >=0
1189 // -> Should use some form of loop Construct
1190 //
1191 if ( !fFullPhiSphere )
1192 {
1193 // First phi surface ('S'tarting phi)
1194 // Comp = Component in outwards normal dirn
1195 //
1196 Comp = v.x()*sinSPhi - v.y()*cosSPhi ;
1197
1198 if ( Comp < 0 )
1199 {
1200 Dist = p.y()*cosSPhi - p.x()*sinSPhi ;
1201
1202 if (Dist < halfCarTolerance)
1203 {
1204 sd = Dist/Comp ;
1205
1206 if (sd < snxt)
1207 {
1208 if ( sd > 0 )
1209 {
1210 xi = p.x() + sd*v.x() ;
1211 yi = p.y() + sd*v.y() ;
1212 zi = p.z() + sd*v.z() ;
1213 rhoi2 = xi*xi + yi*yi ;
1214 radi2 = rhoi2 + zi*zi ;
1215 }
1216 else
1217 {
1218 sd = 0 ;
1219 xi = p.x() ;
1220 yi = p.y() ;
1221 zi = p.z() ;
1222 rhoi2 = rho2 ;
1223 radi2 = rad2 ;
1224 }
1225 if ( (radi2 <= tolORMax2)
1226 && (radi2 >= tolORMin2)
1227 && ((yi*cosCPhi-xi*sinCPhi) <= 0) )
1228 {
1229 // Check theta intersection
1230 // rhoi & zi can never both be 0
1231 // (=>intersect at origin =>fRmax=0)
1232 //
1233 if ( !fFullThetaSphere )
1234 {
1235 iTheta = std::atan2(std::sqrt(rhoi2),zi) ;
1236 if ( (iTheta >= tolSTheta) && (iTheta <= tolETheta) )
1237 {
1238 // r and theta intersections good
1239 // - check intersecting with correct half-plane
1240
1241 if ((yi*cosCPhi-xi*sinCPhi) <= 0)
1242 {
1243 snxt = sd;
1244 }
1245 }
1246 }
1247 else
1248 {
1249 snxt = sd;
1250 }
1251 }
1252 }
1253 }
1254 }
1255
1256 // Second phi surface ('E'nding phi)
1257 // Component in outwards normal dirn
1258
1259 Comp = -( v.x()*sinEPhi-v.y()*cosEPhi ) ;
1260
1261 if (Comp < 0)
1262 {
1263 Dist = -(p.y()*cosEPhi-p.x()*sinEPhi) ;
1264 if ( Dist < halfCarTolerance )
1265 {
1266 sd = Dist/Comp ;
1267
1268 if ( sd < snxt )
1269 {
1270 if (sd > 0)
1271 {
1272 xi = p.x() + sd*v.x() ;
1273 yi = p.y() + sd*v.y() ;
1274 zi = p.z() + sd*v.z() ;
1275 rhoi2 = xi*xi + yi*yi ;
1276 radi2 = rhoi2 + zi*zi ;
1277 }
1278 else
1279 {
1280 sd = 0 ;
1281 xi = p.x() ;
1282 yi = p.y() ;
1283 zi = p.z() ;
1284 rhoi2 = rho2 ;
1285 radi2 = rad2 ;
1286 }
1287 if ( (radi2 <= tolORMax2)
1288 && (radi2 >= tolORMin2)
1289 && ((yi*cosCPhi-xi*sinCPhi) >= 0) )
1290 {
1291 // Check theta intersection
1292 // rhoi & zi can never both be 0
1293 // (=>intersect at origin =>fRmax=0)
1294 //
1295 if ( !fFullThetaSphere )
1296 {
1297 iTheta = std::atan2(std::sqrt(rhoi2),zi) ;
1298 if ( (iTheta >= tolSTheta) && (iTheta <= tolETheta) )
1299 {
1300 // r and theta intersections good
1301 // - check intersecting with correct half-plane
1302
1303 if ((yi*cosCPhi-xi*sinCPhi) >= 0)
1304 {
1305 snxt = sd;
1306 }
1307 }
1308 }
1309 else
1310 {
1311 snxt = sd;
1312 }
1313 }
1314 }
1315 }
1316 }
1317 }
1318
1319 // Theta segment intersection
1320
1321 if ( !fFullThetaSphere )
1322 {
1323
1324 // Intersection with theta surfaces
1325 // Known failure cases:
1326 // o Inside tolerance of stheta surface, skim
1327 // ~parallel to cone and Hit & enter etheta surface [& visa versa]
1328 //
1329 // To solve: Check 2nd root of etheta surface in addition to stheta
1330 //
1331 // o start/end theta is exactly pi/2
1332 // Intersections with cones
1333 //
1334 // Cone equation: x^2+y^2=z^2tan^2(t)
1335 //
1336 // => (px+svx)^2+(py+svy)^2=(pz+svz)^2tan^2(t)
1337 //
1338 // => (px^2+py^2-pz^2tan^2(t))+2sd(pxvx+pyvy-pzvztan^2(t))
1339 // + sd^2(vx^2+vy^2-vz^2tan^2(t)) = 0
1340 //
1341 // => sd^2(1-vz^2(1+tan^2(t))+2sd(pdotv2d-pzvztan^2(t))+(rho2-pz^2tan^2(t))=0
1342
1343 if (fSTheta)
1344 {
1345 dist2STheta = rho2 - p.z()*p.z()*tanSTheta2 ;
1346 }
1347 else
1348 {
1349 dist2STheta = kInfinity ;
1350 }
1351 if ( eTheta < pi )
1352 {
1353 dist2ETheta=rho2-p.z()*p.z()*tanETheta2;
1354 }
1355 else
1356 {
1357 dist2ETheta=kInfinity;
1358 }
1359 if ( pTheta < tolSTheta )
1360 {
1361 // Inside (theta<stheta-tol) stheta cone
1362 // First root of stheta cone, second if first root -ve
1363
1364 t1 = 1 - v.z()*v.z()*(1 + tanSTheta2) ;
1365 t2 = pDotV2d - p.z()*v.z()*tanSTheta2 ;
1366 if (t1)
1367 {
1368 b = t2/t1 ;
1369 c = dist2STheta/t1 ;
1370 d2 = b*b - c ;
1371
1372 if ( d2 >= 0 )
1373 {
1374 d = std::sqrt(d2) ;
1375 sd = -b - d ; // First root
1376 zi = p.z() + sd*v.z();
1377
1378 if ( (sd < 0) || (zi*(fSTheta - halfpi) > 0) )
1379 {
1380 sd = -b+d; // Second root
1381 }
1382 if ((sd >= 0) && (sd < snxt))
1383 {
1384 xi = p.x() + sd*v.x();
1385 yi = p.y() + sd*v.y();
1386 zi = p.z() + sd*v.z();
1387 rhoi2 = xi*xi + yi*yi;
1388 radi2 = rhoi2 + zi*zi;
1389 if ( (radi2 <= tolORMax2)
1390 && (radi2 >= tolORMin2)
1391 && (zi*(fSTheta - halfpi) <= 0) )
1392 {
1393 if ( !fFullPhiSphere && rhoi2 ) // Check phi intersection
1394 {
1395 cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
1396 if (cosPsi >= cosHDPhiOT)
1397 {
1398 snxt = sd;
1399 }
1400 }
1401 else
1402 {
1403 snxt = sd;
1404 }
1405 }
1406 }
1407 }
1408 }
1409
1410 // Possible intersection with ETheta cone.
1411 // Second >= 0 root should be considered
1412
1413 if ( eTheta < pi )
1414 {
1415 t1 = 1 - v.z()*v.z()*(1 + tanETheta2) ;
1416 t2 = pDotV2d - p.z()*v.z()*tanETheta2 ;
1417 if (t1)
1418 {
1419 b = t2/t1 ;
1420 c = dist2ETheta/t1 ;
1421 d2 = b*b - c ;
1422
1423 if (d2 >= 0)
1424 {
1425 d = std::sqrt(d2) ;
1426 sd = -b + d ; // Second root
1427
1428 if ( (sd >= 0) && (sd < snxt) )
1429 {
1430 xi = p.x() + sd*v.x() ;
1431 yi = p.y() + sd*v.y() ;
1432 zi = p.z() + sd*v.z() ;
1433 rhoi2 = xi*xi + yi*yi ;
1434 radi2 = rhoi2 + zi*zi ;
1435
1436 if ( (radi2 <= tolORMax2)
1437 && (radi2 >= tolORMin2)
1438 && (zi*(eTheta - halfpi) <= 0) )
1439 {
1440 if (!fFullPhiSphere && rhoi2) // Check phi intersection
1441 {
1442 cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
1443 if (cosPsi >= cosHDPhiOT)
1444 {
1445 snxt = sd;
1446 }
1447 }
1448 else
1449 {
1450 snxt = sd;
1451 }
1452 }
1453 }
1454 }
1455 }
1456 }
1457 }
1458 else if ( pTheta > tolETheta )
1459 {
1460 // dist2ETheta<-kRadTolerance*0.5 && dist2STheta>0)
1461 // Inside (theta > etheta+tol) e-theta cone
1462 // First root of etheta cone, second if first root 'imaginary'
1463
1464 t1 = 1 - v.z()*v.z()*(1 + tanETheta2) ;
1465 t2 = pDotV2d - p.z()*v.z()*tanETheta2 ;
1466 if (t1)
1467 {
1468 b = t2/t1 ;
1469 c = dist2ETheta/t1 ;
1470 d2 = b*b - c ;
1471
1472 if (d2 >= 0)
1473 {
1474 d = std::sqrt(d2) ;
1475 sd = -b - d ; // First root
1476 zi = p.z() + sd*v.z();
1477
1478 if ( (sd < 0) || (zi*(eTheta - halfpi) > 0) )
1479 {
1480 sd = -b + d ; // second root
1481 }
1482 if ( (sd >= 0) && (sd < snxt) )
1483 {
1484 xi = p.x() + sd*v.x() ;
1485 yi = p.y() + sd*v.y() ;
1486 zi = p.z() + sd*v.z() ;
1487 rhoi2 = xi*xi + yi*yi ;
1488 radi2 = rhoi2 + zi*zi ;
1489
1490 if ( (radi2 <= tolORMax2)
1491 && (radi2 >= tolORMin2)
1492 && (zi*(eTheta - halfpi) <= 0) )
1493 {
1494 if (!fFullPhiSphere && rhoi2) // Check phi intersection
1495 {
1496 cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
1497 if (cosPsi >= cosHDPhiOT)
1498 {
1499 snxt = sd;
1500 }
1501 }
1502 else
1503 {
1504 snxt = sd;
1505 }
1506 }
1507 }
1508 }
1509 }
1510
1511 // Possible intersection with STheta cone.
1512 // Second >= 0 root should be considered
1513
1514 if ( fSTheta )
1515 {
1516 t1 = 1 - v.z()*v.z()*(1 + tanSTheta2) ;
1517 t2 = pDotV2d - p.z()*v.z()*tanSTheta2 ;
1518 if (t1)
1519 {
1520 b = t2/t1 ;
1521 c = dist2STheta/t1 ;
1522 d2 = b*b - c ;
1523
1524 if (d2 >= 0)
1525 {
1526 d = std::sqrt(d2) ;
1527 sd = -b + d ; // Second root
1528
1529 if ( (sd >= 0) && (sd < snxt) )
1530 {
1531 xi = p.x() + sd*v.x() ;
1532 yi = p.y() + sd*v.y() ;
1533 zi = p.z() + sd*v.z() ;
1534 rhoi2 = xi*xi + yi*yi ;
1535 radi2 = rhoi2 + zi*zi ;
1536
1537 if ( (radi2 <= tolORMax2)
1538 && (radi2 >= tolORMin2)
1539 && (zi*(fSTheta - halfpi) <= 0) )
1540 {
1541 if (!fFullPhiSphere && rhoi2) // Check phi intersection
1542 {
1543 cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
1544 if (cosPsi >= cosHDPhiOT)
1545 {
1546 snxt = sd;
1547 }
1548 }
1549 else
1550 {
1551 snxt = sd;
1552 }
1553 }
1554 }
1555 }
1556 }
1557 }
1558 }
1559 else if ( (pTheta < tolSTheta + kAngTolerance)
1560 && (fSTheta > halfAngTolerance) )
1561 {
1562 // In tolerance of stheta
1563 // If entering through solid [r,phi] => 0 to in
1564 // else try 2nd root
1565
1566 t2 = pDotV2d - p.z()*v.z()*tanSTheta2 ;
1567 if ( (t2>=0 && tolIRMin2<rad2 && rad2<tolIRMax2 && fSTheta<halfpi)
1568 || (t2<0 && tolIRMin2<rad2 && rad2<tolIRMax2 && fSTheta>halfpi)
1569 || (v.z()<0 && tolIRMin2<rad2 && rad2<tolIRMax2 && fSTheta==halfpi) )
1570 {
1571 if (!fFullPhiSphere && rho2) // Check phi intersection
1572 {
1573 cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/std::sqrt(rho2) ;
1574 if (cosPsi >= cosHDPhiIT)
1575 {
1576 return 0 ;
1577 }
1578 }
1579 else
1580 {
1581 return 0 ;
1582 }
1583 }
1584
1585 // Not entering immediately/travelling through
1586
1587 t1 = 1 - v.z()*v.z()*(1 + tanSTheta2) ;
1588 if (t1)
1589 {
1590 b = t2/t1 ;
1591 c = dist2STheta/t1 ;
1592 d2 = b*b - c ;
1593
1594 if (d2 >= 0)
1595 {
1596 d = std::sqrt(d2) ;
1597 sd = -b + d ;
1598 if ( (sd >= halfCarTolerance) && (sd < snxt) && (fSTheta < halfpi) )
1599 { // ^^^^^^^^^^^^^^^^^^^^^ shouldn't it be >=0 instead ?
1600 xi = p.x() + sd*v.x() ;
1601 yi = p.y() + sd*v.y() ;
1602 zi = p.z() + sd*v.z() ;
1603 rhoi2 = xi*xi + yi*yi ;
1604 radi2 = rhoi2 + zi*zi ;
1605
1606 if ( (radi2 <= tolORMax2)
1607 && (radi2 >= tolORMin2)
1608 && (zi*(fSTheta - halfpi) <= 0) )
1609 {
1610 if ( !fFullPhiSphere && rhoi2 ) // Check phi intersection
1611 {
1612 cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
1613 if ( cosPsi >= cosHDPhiOT )
1614 {
1615 snxt = sd;
1616 }
1617 }
1618 else
1619 {
1620 snxt = sd;
1621 }
1622 }
1623 }
1624 }
1625 }
1626 }
1627 else if ((pTheta > tolETheta-kAngTolerance) && (eTheta < pi-kAngTolerance))
1628 {
1629
1630 // In tolerance of etheta
1631 // If entering through solid [r,phi] => 0 to in
1632 // else try 2nd root
1633
1634 t2 = pDotV2d - p.z()*v.z()*tanETheta2 ;
1635
1636 if ( ((t2<0) && (eTheta < halfpi)
1637 && (tolIRMin2 < rad2) && (rad2 < tolIRMax2))
1638 || ((t2>=0) && (eTheta > halfpi)
1639 && (tolIRMin2 < rad2) && (rad2 < tolIRMax2))
1640 || ((v.z()>0) && (eTheta == halfpi)
1641 && (tolIRMin2 < rad2) && (rad2 < tolIRMax2)) )
1642 {
1643 if (!fFullPhiSphere && rho2) // Check phi intersection
1644 {
1645 cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/std::sqrt(rho2) ;
1646 if (cosPsi >= cosHDPhiIT)
1647 {
1648 return 0 ;
1649 }
1650 }
1651 else
1652 {
1653 return 0 ;
1654 }
1655 }
1656
1657 // Not entering immediately/travelling through
1658
1659 t1 = 1 - v.z()*v.z()*(1 + tanETheta2) ;
1660 if (t1)
1661 {
1662 b = t2/t1 ;
1663 c = dist2ETheta/t1 ;
1664 d2 = b*b - c ;
1665
1666 if (d2 >= 0)
1667 {
1668 d = std::sqrt(d2) ;
1669 sd = -b + d ;
1670
1671 if ( (sd >= halfCarTolerance)
1672 && (sd < snxt) && (eTheta > halfpi) )
1673 {
1674 xi = p.x() + sd*v.x() ;
1675 yi = p.y() + sd*v.y() ;
1676 zi = p.z() + sd*v.z() ;
1677 rhoi2 = xi*xi + yi*yi ;
1678 radi2 = rhoi2 + zi*zi ;
1679
1680 if ( (radi2 <= tolORMax2)
1681 && (radi2 >= tolORMin2)
1682 && (zi*(eTheta - halfpi) <= 0) )
1683 {
1684 if (!fFullPhiSphere && rhoi2) // Check phi intersection
1685 {
1686 cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
1687 if (cosPsi >= cosHDPhiOT)
1688 {
1689 snxt = sd;
1690 }
1691 }
1692 else
1693 {
1694 snxt = sd;
1695 }
1696 }
1697 }
1698 }
1699 }
1700 }
1701 else
1702 {
1703 // stheta+tol<theta<etheta-tol
1704 // For BOTH stheta & etheta check 2nd root for validity [r,phi]
1705
1706 t1 = 1 - v.z()*v.z()*(1 + tanSTheta2) ;
1707 t2 = pDotV2d - p.z()*v.z()*tanSTheta2 ;
1708 if (t1)
1709 {
1710 b = t2/t1;
1711 c = dist2STheta/t1 ;
1712 d2 = b*b - c ;
1713
1714 if (d2 >= 0)
1715 {
1716 d = std::sqrt(d2) ;
1717 sd = -b + d ; // second root
1718
1719 if ((sd >= 0) && (sd < snxt))
1720 {
1721 xi = p.x() + sd*v.x() ;
1722 yi = p.y() + sd*v.y() ;
1723 zi = p.z() + sd*v.z() ;
1724 rhoi2 = xi*xi + yi*yi ;
1725 radi2 = rhoi2 + zi*zi ;
1726
1727 if ( (radi2 <= tolORMax2)
1728 && (radi2 >= tolORMin2)
1729 && (zi*(fSTheta - halfpi) <= 0) )
1730 {
1731 if (!fFullPhiSphere && rhoi2) // Check phi intersection
1732 {
1733 cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
1734 if (cosPsi >= cosHDPhiOT)
1735 {
1736 snxt = sd;
1737 }
1738 }
1739 else
1740 {
1741 snxt = sd;
1742 }
1743 }
1744 }
1745 }
1746 }
1747 t1 = 1 - v.z()*v.z()*(1 + tanETheta2) ;
1748 t2 = pDotV2d - p.z()*v.z()*tanETheta2 ;
1749 if (t1)
1750 {
1751 b = t2/t1 ;
1752 c = dist2ETheta/t1 ;
1753 d2 = b*b - c ;
1754
1755 if (d2 >= 0)
1756 {
1757 d = std::sqrt(d2) ;
1758 sd = -b + d; // second root
1759
1760 if ((sd >= 0) && (sd < snxt))
1761 {
1762 xi = p.x() + sd*v.x() ;
1763 yi = p.y() + sd*v.y() ;
1764 zi = p.z() + sd*v.z() ;
1765 rhoi2 = xi*xi + yi*yi ;
1766 radi2 = rhoi2 + zi*zi ;
1767
1768 if ( (radi2 <= tolORMax2)
1769 && (radi2 >= tolORMin2)
1770 && (zi*(eTheta - halfpi) <= 0) )
1771 {
1772 if (!fFullPhiSphere && rhoi2) // Check phi intersection
1773 {
1774 cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
1775 if ( cosPsi >= cosHDPhiOT )
1776 {
1777 snxt = sd;
1778 }
1779 }
1780 else
1781 {
1782 snxt = sd;
1783 }
1784 }
1785 }
1786 }
1787 }
1788 }
1789 }
1790 return snxt;
1791}
G4Sphere::DistanceToIn
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
Definition: G4Sphere.cc:867

Referenced by DistanceToIn().

◆ DistanceToOut() [1/2]

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

Implements G4VSolid.

Definition at line 2781 of file G4Sphere.cc.

2782{
2783 G4double safe=0.0,safeRMin,safeRMax,safePhi,safeTheta;
2784 G4double rho2,rds,rho;
2785 G4double pTheta,dTheta1,dTheta2;
2786 rho2=p.x()*p.x()+p.y()*p.y();
2787 rds=std::sqrt(rho2+p.z()*p.z());
2788 rho=std::sqrt(rho2);
2789
2790#ifdef G4CSGDEBUG
2791 if( Inside(p) == kOutside )
2792 {
2793 G4int old_prc = G4cout.precision(16);
2794 G4cout << G4endl;
2795 DumpInfo();
2796 G4cout << "Position:" << G4endl << G4endl ;
2797 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl ;
2798 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ;
2799 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ;
2800 G4cout.precision(old_prc) ;
2801 G4Exception("G4Sphere::DistanceToOut(p)",
2802 "GeomSolids1002", JustWarning, "Point p is outside !?" );
2803 }
2804#endif
2805
2806 //
2807 // Distance to r shells
2808 //
2809 if (fRmin)
2810 {
2811 safeRMin=rds-fRmin;
2812 safeRMax=fRmax-rds;
2813 if (safeRMin<safeRMax)
2814 {
2815 safe=safeRMin;
2816 }
2817 else
2818 {
2819 safe=safeRMax;
2820 }
2821 }
2822 else
2823 {
2824 safe=fRmax-rds;
2825 }
2826
2827 //
2828 // Distance to phi extent
2829 //
2830 if (!fFullPhiSphere && rho)
2831 {
2832 if ((p.y()*cosCPhi-p.x()*sinCPhi)<=0)
2833 {
2834 safePhi=-(p.x()*sinSPhi-p.y()*cosSPhi);
2835 }
2836 else
2837 {
2838 safePhi=(p.x()*sinEPhi-p.y()*cosEPhi);
2839 }
2840 if (safePhi<safe) { safe=safePhi; }
2841 }
2842
2843 //
2844 // Distance to Theta extent
2845 //
2846 if (rds)
2847 {
2848 pTheta=std::acos(p.z()/rds);
2849 if (pTheta<0) { pTheta+=pi; }
2850 dTheta1=pTheta-fSTheta;
2851 dTheta2=eTheta-pTheta;
2852 if (dTheta1<dTheta2) { safeTheta=rds*std::sin(dTheta1); }
2853 else { safeTheta=rds*std::sin(dTheta2); }
2854 if (safe>safeTheta) { safe=safeTheta; }
2855 }
2856
2857 if (safe<0) { safe=0; }
2858 return safe;
2859}
JustWarning
@ JustWarning
Definition: G4ExceptionSeverity.hh:64
G4cout
G4DLLIMPORT std::ostream G4cout
G4VSolid::DumpInfo
void DumpInfo() const

◆ DistanceToOut() [2/2]

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

Implements G4VSolid.

Definition at line 1897 of file G4Sphere.cc.

1902{
1903 G4double snxt = kInfinity; // snxt is default return value
1904 G4double sphi= kInfinity,stheta= kInfinity;
1905 ESide side=kNull,sidephi=kNull,sidetheta=kNull;
1906
1907 static const G4double halfCarTolerance = kCarTolerance*0.5;
1908 static const G4double halfAngTolerance = kAngTolerance*0.5;
1909 const G4double halfRmaxTolerance = fRmaxTolerance*0.5;
1910 const G4double halfRminTolerance = fRminTolerance*0.5;
1911 const G4double Rmax_plus = fRmax + halfRmaxTolerance;
1912 const G4double Rmin_minus = (fRmin) ? fRmin-halfRminTolerance : 0;
1913 G4double t1,t2;
1914 G4double b,c,d;
1915
1916 // Variables for phi intersection:
1917
1918 G4double pDistS,compS,pDistE,compE,sphi2,vphi;
1919
1920 G4double rho2,rad2,pDotV2d,pDotV3d;
1921
1922 G4double xi,yi,zi; // Intersection point
1923
1924 // Theta precals
1925 //
1926 G4double rhoSecTheta;
1927 G4double dist2STheta, dist2ETheta, distTheta;
1928 G4double d2,sd;
1929
1930 // General Precalcs
1931 //
1932 rho2 = p.x()*p.x()+p.y()*p.y();
1933 rad2 = rho2+p.z()*p.z();
1934
1935 pDotV2d = p.x()*v.x()+p.y()*v.y();
1936 pDotV3d = pDotV2d+p.z()*v.z();
1937
1938 // Radial Intersections from G4Sphere::DistanceToIn
1939 //
1940 // Outer spherical shell intersection
1941 // - Only if outside tolerant fRmax
1942 // - Check for if inside and outer G4Sphere heading through solid (-> 0)
1943 // - No intersect -> no intersection with G4Sphere
1944 //
1945 // Shell eqn: x^2+y^2+z^2=RSPH^2
1946 //
1947 // => (px+svx)^2+(py+svy)^2+(pz+svz)^2=R^2
1948 //
1949 // => (px^2+py^2+pz^2) +2sd(pxvx+pyvy+pzvz)+sd^2(vx^2+vy^2+vz^2)=R^2
1950 // => rad2 +2sd(pDotV3d) +sd^2 =R^2
1951 //
1952 // => sd=-pDotV3d+-std::sqrt(pDotV3d^2-(rad2-R^2))
1953
1954 if( (rad2 <= Rmax_plus*Rmax_plus) && (rad2 >= Rmin_minus*Rmin_minus) )
1955 {
1956 c = rad2 - fRmax*fRmax;
1957
1958 if (c < fRmaxTolerance*fRmax)
1959 {
1960 // Within tolerant Outer radius
1961 //
1962 // The test is
1963 // rad - fRmax < 0.5*kRadTolerance
1964 // => rad < fRmax + 0.5*kRadTol
1965 // => rad2 < (fRmax + 0.5*kRadTol)^2
1966 // => rad2 < fRmax^2 + 2.*0.5*fRmax*kRadTol + 0.25*kRadTol*kRadTol
1967 // => rad2 - fRmax^2 <~ fRmax*kRadTol
1968
1969 d2 = pDotV3d*pDotV3d - c;
1970
1971 if( (c >- fRmaxTolerance*fRmax) // on tolerant surface
1972 && ((pDotV3d >=0) || (d2 < 0)) ) // leaving outside from Rmax
1973 // not re-entering
1974 {
1975 if(calcNorm)
1976 {
1977 *validNorm = true ;
1978 *n = G4ThreeVector(p.x()/fRmax,p.y()/fRmax,p.z()/fRmax) ;
1979 }
1980 return snxt = 0;
1981 }
1982 else
1983 {
1984 snxt = -pDotV3d+std::sqrt(d2); // second root since inside Rmax
1985 side = kRMax ;
1986 }
1987 }
1988
1989 // Inner spherical shell intersection:
1990 // Always first >=0 root, because would have passed
1991 // from outside of Rmin surface .
1992
1993 if (fRmin)
1994 {
1995 c = rad2 - fRmin*fRmin;
1996 d2 = pDotV3d*pDotV3d - c;
1997
1998 if (c >- fRminTolerance*fRmin) // 2.0 * (0.5*kRadTolerance) * fRmin
1999 {
2000 if ( (c < fRminTolerance*fRmin) // leaving from Rmin
2001 && (d2 >= fRminTolerance*fRmin) && (pDotV3d < 0) )
2002 {
2003 if(calcNorm) { *validNorm = false; } // Rmin surface is concave
2004 return snxt = 0 ;
2005 }
2006 else
2007 {
2008 if ( d2 >= 0. )
2009 {
2010 sd = -pDotV3d-std::sqrt(d2);
2011
2012 if ( sd >= 0. ) // Always intersect Rmin first
2013 {
2014 snxt = sd ;
2015 side = kRMin ;
2016 }
2017 }
2018 }
2019 }
2020 }
2021 }
2022
2023 // Theta segment intersection
2024
2025 if ( !fFullThetaSphere )
2026 {
2027 // Intersection with theta surfaces
2028 //
2029 // Known failure cases:
2030 // o Inside tolerance of stheta surface, skim
2031 // ~parallel to cone and Hit & enter etheta surface [& visa versa]
2032 //
2033 // To solve: Check 2nd root of etheta surface in addition to stheta
2034 //
2035 // o start/end theta is exactly pi/2
2036 //
2037 // Intersections with cones
2038 //
2039 // Cone equation: x^2+y^2=z^2tan^2(t)
2040 //
2041 // => (px+svx)^2+(py+svy)^2=(pz+svz)^2tan^2(t)
2042 //
2043 // => (px^2+py^2-pz^2tan^2(t))+2sd(pxvx+pyvy-pzvztan^2(t))
2044 // + sd^2(vx^2+vy^2-vz^2tan^2(t)) = 0
2045 //
2046 // => sd^2(1-vz^2(1+tan^2(t))+2sd(pdotv2d-pzvztan^2(t))+(rho2-pz^2tan^2(t))=0
2047 //
2048
2049 if(fSTheta) // intersection with first cons
2050 {
2051 if( std::fabs(tanSTheta) > 5./kAngTolerance ) // kons is plane z=0
2052 {
2053 if( v.z() > 0. )
2054 {
2055 if ( std::fabs( p.z() ) <= halfRmaxTolerance )
2056 {
2057 if(calcNorm)
2058 {
2059 *validNorm = true;
2060 *n = G4ThreeVector(0.,0.,1.);
2061 }
2062 return snxt = 0 ;
2063 }
2064 stheta = -p.z()/v.z();
2065 sidetheta = kSTheta;
2066 }
2067 }
2068 else // kons is not plane
2069 {
2070 t1 = 1-v.z()*v.z()*(1+tanSTheta2);
2071 t2 = pDotV2d-p.z()*v.z()*tanSTheta2; // ~vDotN if p on cons
2072 dist2STheta = rho2-p.z()*p.z()*tanSTheta2; // t3
2073
2074 distTheta = std::sqrt(rho2)-p.z()*tanSTheta;
2075
2076 if( std::fabs(t1) < halfAngTolerance ) // 1st order equation,
2077 { // v parallel to kons
2078 if( v.z() > 0. )
2079 {
2080 if(std::fabs(distTheta) < halfRmaxTolerance) // p on surface
2081 {
2082 if( (fSTheta < halfpi) && (p.z() > 0.) )
2083 {
2084 if( calcNorm ) { *validNorm = false; }
2085 return snxt = 0.;
2086 }
2087 else if( (fSTheta > halfpi) && (p.z() <= 0) )
2088 {
2089 if( calcNorm )
2090 {
2091 *validNorm = true;
2092 if (rho2)
2093 {
2094 rhoSecTheta = std::sqrt(rho2*(1+tanSTheta2));
2095
2096 *n = G4ThreeVector( p.x()/rhoSecTheta,
2097 p.y()/rhoSecTheta,
2098 std::sin(fSTheta) );
2099 }
2100 else *n = G4ThreeVector(0.,0.,1.);
2101 }
2102 return snxt = 0.;
2103 }
2104 }
2105 stheta = -0.5*dist2STheta/t2;
2106 sidetheta = kSTheta;
2107 }
2108 } // 2nd order equation, 1st root of fSTheta cone,
2109 else // 2nd if 1st root -ve
2110 {
2111 if( std::fabs(distTheta) < halfRmaxTolerance )
2112 {
2113 if( (fSTheta > halfpi) && (t2 >= 0.) ) // leave
2114 {
2115 if( calcNorm )
2116 {
2117 *validNorm = true;
2118 if (rho2)
2119 {
2120 rhoSecTheta = std::sqrt(rho2*(1+tanSTheta2));
2121
2122 *n = G4ThreeVector( p.x()/rhoSecTheta,
2123 p.y()/rhoSecTheta,
2124 std::sin(fSTheta) );
2125 }
2126 else { *n = G4ThreeVector(0.,0.,1.); }
2127 }
2128 return snxt = 0.;
2129 }
2130 else if( (fSTheta < halfpi) && (t2 < 0.) && (p.z() >=0.) ) // leave
2131 {
2132 if( calcNorm ) { *validNorm = false; }
2133 return snxt = 0.;
2134 }
2135 }
2136 b = t2/t1;
2137 c = dist2STheta/t1;
2138 d2 = b*b - c ;
2139
2140 if ( d2 >= 0. )
2141 {
2142 d = std::sqrt(d2);
2143
2144 if( fSTheta > halfpi )
2145 {
2146 sd = -b - d; // First root
2147
2148 if ( ((std::fabs(s) < halfRmaxTolerance) && (t2 < 0.))
2149 || (sd < 0.) || ( (sd > 0.) && (p.z() + sd*v.z() > 0.) ) )
2150 {
2151 sd = -b + d ; // 2nd root
2152 }
2153 if( (sd > halfRmaxTolerance) && (p.z() + sd*v.z() <= 0.) )
2154 {
2155 stheta = sd;
2156 sidetheta = kSTheta;
2157 }
2158 }
2159 else // sTheta < pi/2, concave surface, no normal
2160 {
2161 sd = -b - d; // First root
2162
2163 if ( ( (std::fabs(sd) < halfRmaxTolerance) && (t2 >= 0.) )
2164 || (sd < 0.) || ( (sd > 0.) && (p.z() + sd*v.z() < 0.) ) )
2165 {
2166 sd = -b + d ; // 2nd root
2167 }
2168 if( (sd > halfRmaxTolerance) && (p.z() + sd*v.z() >= 0.) )
2169 {
2170 stheta = sd;
2171 sidetheta = kSTheta;
2172 }
2173 }
2174 }
2175 }
2176 }
2177 }
2178 if (eTheta < pi) // intersection with second cons
2179 {
2180 if( std::fabs(tanETheta) > 5./kAngTolerance ) // kons is plane z=0
2181 {
2182 if( v.z() < 0. )
2183 {
2184 if ( std::fabs( p.z() ) <= halfRmaxTolerance )
2185 {
2186 if(calcNorm)
2187 {
2188 *validNorm = true;
2189 *n = G4ThreeVector(0.,0.,-1.);
2190 }
2191 return snxt = 0 ;
2192 }
2193 sd = -p.z()/v.z();
2194
2195 if( sd < stheta )
2196 {
2197 stheta = sd;
2198 sidetheta = kETheta;
2199 }
2200 }
2201 }
2202 else // kons is not plane
2203 {
2204 t1 = 1-v.z()*v.z()*(1+tanETheta2);
2205 t2 = pDotV2d-p.z()*v.z()*tanETheta2; // ~vDotN if p on cons
2206 dist2ETheta = rho2-p.z()*p.z()*tanETheta2; // t3
2207
2208 distTheta = std::sqrt(rho2)-p.z()*tanETheta;
2209
2210 if( std::fabs(t1) < halfAngTolerance ) // 1st order equation,
2211 { // v parallel to kons
2212 if( v.z() < 0. )
2213 {
2214 if(std::fabs(distTheta) < halfRmaxTolerance) // p on surface
2215 {
2216 if( (eTheta > halfpi) && (p.z() < 0.) )
2217 {
2218 if( calcNorm ) { *validNorm = false; }
2219 return snxt = 0.;
2220 }
2221 else if ( (eTheta < halfpi) && (p.z() >= 0) )
2222 {
2223 if( calcNorm )
2224 {
2225 *validNorm = true;
2226 if (rho2)
2227 {
2228 rhoSecTheta = std::sqrt(rho2*(1+tanETheta2));
2229 *n = G4ThreeVector( p.x()/rhoSecTheta,
2230 p.y()/rhoSecTheta,
2231 -sinETheta );
2232 }
2233 else { *n = G4ThreeVector(0.,0.,-1.); }
2234 }
2235 return snxt = 0.;
2236 }
2237 }
2238 sd = -0.5*dist2ETheta/t2;
2239
2240 if( sd < stheta )
2241 {
2242 stheta = sd;
2243 sidetheta = kETheta;
2244 }
2245 }
2246 } // 2nd order equation, 1st root of fSTheta cone
2247 else // 2nd if 1st root -ve
2248 {
2249 if ( std::fabs(distTheta) < halfRmaxTolerance )
2250 {
2251 if( (eTheta < halfpi) && (t2 >= 0.) ) // leave
2252 {
2253 if( calcNorm )
2254 {
2255 *validNorm = true;
2256 if (rho2)
2257 {
2258 rhoSecTheta = std::sqrt(rho2*(1+tanETheta2));
2259 *n = G4ThreeVector( p.x()/rhoSecTheta,
2260 p.y()/rhoSecTheta,
2261 -sinETheta );
2262 }
2263 else *n = G4ThreeVector(0.,0.,-1.);
2264 }
2265 return snxt = 0.;
2266 }
2267 else if ( (eTheta > halfpi)
2268 && (t2 < 0.) && (p.z() <=0.) ) // leave
2269 {
2270 if( calcNorm ) { *validNorm = false; }
2271 return snxt = 0.;
2272 }
2273 }
2274 b = t2/t1;
2275 c = dist2ETheta/t1;
2276 d2 = b*b - c ;
2277
2278 if ( d2 >= 0. )
2279 {
2280 d = std::sqrt(d2);
2281
2282 if( eTheta < halfpi )
2283 {
2284 sd = -b - d; // First root
2285
2286 if( ((std::fabs(sd) < halfRmaxTolerance) && (t2 < 0.))
2287 || (sd < 0.) )
2288 {
2289 sd = -b + d ; // 2nd root
2290 }
2291 if( sd > halfRmaxTolerance )
2292 {
2293 if( sd < stheta )
2294 {
2295 stheta = sd;
2296 sidetheta = kETheta;
2297 }
2298 }
2299 }
2300 else // sTheta+fDTheta > pi/2, concave surface, no normal
2301 {
2302 sd = -b - d; // First root
2303
2304 if ( ((std::fabs(sd) < halfRmaxTolerance) && (t2 >= 0.))
2305 || (sd < 0.) || ( (sd > 0.) && (p.z() + sd*v.z() > 0.) ) )
2306 {
2307 sd = -b + d ; // 2nd root
2308 }
2309 if( (sd > halfRmaxTolerance) && (p.z() + sd*v.z() <= 0.) )
2310 {
2311 if( sd < stheta )
2312 {
2313 stheta = sd;
2314 sidetheta = kETheta;
2315 }
2316 }
2317 }
2318 }
2319 }
2320 }
2321 }
2322
2323 } // end theta intersections
2324
2325 // Phi Intersection
2326
2327 if ( !fFullPhiSphere )
2328 {
2329 if ( p.x() || p.y() ) // Check if on z axis (rho not needed later)
2330 {
2331 // pDist -ve when inside
2332
2333 pDistS=p.x()*sinSPhi-p.y()*cosSPhi;
2334 pDistE=-p.x()*sinEPhi+p.y()*cosEPhi;
2335
2336 // Comp -ve when in direction of outwards normal
2337
2338 compS = -sinSPhi*v.x()+cosSPhi*v.y() ;
2339 compE = sinEPhi*v.x()-cosEPhi*v.y() ;
2340 sidephi = kNull ;
2341
2342 if ( (pDistS <= 0) && (pDistE <= 0) )
2343 {
2344 // Inside both phi *full* planes
2345
2346 if ( compS < 0 )
2347 {
2348 sphi = pDistS/compS ;
2349 xi = p.x()+sphi*v.x() ;
2350 yi = p.y()+sphi*v.y() ;
2351
2352 // Check intersection with correct half-plane (if not -> no intersect)
2353 //
2354 if( (std::fabs(xi)<=kCarTolerance) && (std::fabs(yi)<=kCarTolerance) )
2355 {
2356 vphi = std::atan2(v.y(),v.x());
2357 sidephi = kSPhi;
2358 if ( ( (fSPhi-halfAngTolerance) <= vphi)
2359 && ( (ePhi+halfAngTolerance) >= vphi) )
2360 {
2361 sphi = kInfinity;
2362 }
2363 }
2364 else if ( ( yi*cosCPhi - xi*sinCPhi ) >= 0 )
2365 {
2366 sphi=kInfinity;
2367 }
2368 else
2369 {
2370 sidephi = kSPhi ;
2371 if ( pDistS > -halfCarTolerance) { sphi = 0; } // Leave by sphi
2372 }
2373 }
2374 else { sphi = kInfinity; }
2375
2376 if ( compE < 0 )
2377 {
2378 sphi2=pDistE/compE ;
2379 if (sphi2 < sphi) // Only check further if < starting phi intersection
2380 {
2381 xi = p.x()+sphi2*v.x() ;
2382 yi = p.y()+sphi2*v.y() ;
2383
2384 // Check intersection with correct half-plane
2385 //
2386 if ((std::fabs(xi)<=kCarTolerance) && (std::fabs(yi)<=kCarTolerance))
2387 {
2388 // Leaving via ending phi
2389 //
2390 vphi = std::atan2(v.y(),v.x()) ;
2391
2392 if( !((fSPhi-halfAngTolerance <= vphi)
2393 &&(fSPhi+fDPhi+halfAngTolerance >= vphi)) )
2394 {
2395 sidephi = kEPhi;
2396 if ( pDistE <= -halfCarTolerance ) { sphi = sphi2; }
2397 else { sphi = 0.0; }
2398 }
2399 }
2400 else if ((yi*cosCPhi-xi*sinCPhi)>=0) // Leaving via ending phi
2401 {
2402 sidephi = kEPhi ;
2403 if ( pDistE <= -halfCarTolerance )
2404 {
2405 sphi=sphi2;
2406 }
2407 else
2408 {
2409 sphi = 0 ;
2410 }
2411 }
2412 }
2413 }
2414 }
2415 else if ((pDistS >= 0) && (pDistE >= 0)) // Outside both *full* phi planes
2416 {
2417 if ( pDistS <= pDistE )
2418 {
2419 sidephi = kSPhi ;
2420 }
2421 else
2422 {
2423 sidephi = kEPhi ;
2424 }
2425 if ( fDPhi > pi )
2426 {
2427 if ( (compS < 0) && (compE < 0) ) { sphi = 0; }
2428 else { sphi = kInfinity; }
2429 }
2430 else
2431 {
2432 // if towards both >=0 then once inside (after error)
2433 // will remain inside
2434
2435 if ( (compS >= 0) && (compE >= 0) ) { sphi = kInfinity; }
2436 else { sphi = 0; }
2437 }
2438 }
2439 else if ( (pDistS > 0) && (pDistE < 0) )
2440 {
2441 // Outside full starting plane, inside full ending plane
2442
2443 if ( fDPhi > pi )
2444 {
2445 if ( compE < 0 )
2446 {
2447 sphi = pDistE/compE ;
2448 xi = p.x() + sphi*v.x() ;
2449 yi = p.y() + sphi*v.y() ;
2450
2451 // Check intersection in correct half-plane
2452 // (if not -> not leaving phi extent)
2453 //
2454 if( (std::fabs(xi)<=kCarTolerance)&&(std::fabs(yi)<=kCarTolerance) )
2455 {
2456 vphi = std::atan2(v.y(),v.x());
2457 sidephi = kSPhi;
2458 if ( ( (fSPhi-halfAngTolerance) <= vphi)
2459 && ( (ePhi+halfAngTolerance) >= vphi) )
2460 {
2461 sphi = kInfinity;
2462 }
2463 }
2464 else if ( ( yi*cosCPhi - xi*sinCPhi ) <= 0 )
2465 {
2466 sphi = kInfinity ;
2467 }
2468 else // Leaving via Ending phi
2469 {
2470 sidephi = kEPhi ;
2471 if ( pDistE > -halfCarTolerance ) { sphi = 0.; }
2472 }
2473 }
2474 else
2475 {
2476 sphi = kInfinity ;
2477 }
2478 }
2479 else
2480 {
2481 if ( compS >= 0 )
2482 {
2483 if ( compE < 0 )
2484 {
2485 sphi = pDistE/compE ;
2486 xi = p.x() + sphi*v.x() ;
2487 yi = p.y() + sphi*v.y() ;
2488
2489 // Check intersection in correct half-plane
2490 // (if not -> remain in extent)
2491 //
2492 if( (std::fabs(xi)<=kCarTolerance)&&(std::fabs(yi)<=kCarTolerance) )
2493 {
2494 vphi = std::atan2(v.y(),v.x());
2495 sidephi = kSPhi;
2496 if ( ( (fSPhi-halfAngTolerance) <= vphi)
2497 && ( (ePhi+halfAngTolerance) >= vphi) )
2498 {
2499 sphi = kInfinity;
2500 }
2501 }
2502 else if ( ( yi*cosCPhi - xi*sinCPhi) <= 0 )
2503 {
2504 sphi=kInfinity;
2505 }
2506 else // otherwise leaving via Ending phi
2507 {
2508 sidephi = kEPhi ;
2509 }
2510 }
2511 else sphi=kInfinity;
2512 }
2513 else // leaving immediately by starting phi
2514 {
2515 sidephi = kSPhi ;
2516 sphi = 0 ;
2517 }
2518 }
2519 }
2520 else
2521 {
2522 // Must be pDistS < 0 && pDistE > 0
2523 // Inside full starting plane, outside full ending plane
2524
2525 if ( fDPhi > pi )
2526 {
2527 if ( compS < 0 )
2528 {
2529 sphi=pDistS/compS;
2530 xi=p.x()+sphi*v.x();
2531 yi=p.y()+sphi*v.y();
2532
2533 // Check intersection in correct half-plane
2534 // (if not -> not leaving phi extent)
2535 //
2536 if( (std::fabs(xi)<=kCarTolerance)&&(std::fabs(yi)<=kCarTolerance) )
2537 {
2538 vphi = std::atan2(v.y(),v.x()) ;
2539 sidephi = kSPhi;
2540 if ( ( (fSPhi-halfAngTolerance) <= vphi)
2541 && ( (ePhi+halfAngTolerance) >= vphi) )
2542 {
2543 sphi = kInfinity;
2544 }
2545 }
2546 else if ( ( yi*cosCPhi - xi*sinCPhi ) >= 0 )
2547 {
2548 sphi = kInfinity ;
2549 }
2550 else // Leaving via Starting phi
2551 {
2552 sidephi = kSPhi ;
2553 if ( pDistS > -halfCarTolerance ) { sphi = 0; }
2554 }
2555 }
2556 else
2557 {
2558 sphi = kInfinity ;
2559 }
2560 }
2561 else
2562 {
2563 if ( compE >= 0 )
2564 {
2565 if ( compS < 0 )
2566 {
2567 sphi = pDistS/compS ;
2568 xi = p.x()+sphi*v.x() ;
2569 yi = p.y()+sphi*v.y() ;
2570
2571 // Check intersection in correct half-plane
2572 // (if not -> remain in extent)
2573 //
2574 if((std::fabs(xi)<=kCarTolerance) && (std::fabs(yi)<=kCarTolerance))
2575 {
2576 vphi = std::atan2(v.y(),v.x()) ;
2577 sidephi = kSPhi;
2578 if ( ( (fSPhi-halfAngTolerance) <= vphi)
2579 && ( (ePhi+halfAngTolerance) >= vphi) )
2580 {
2581 sphi = kInfinity;
2582 }
2583 }
2584 else if ( ( yi*cosCPhi - xi*sinCPhi ) >= 0 )
2585 {
2586 sphi = kInfinity ;
2587 }
2588 else // otherwise leaving via Starting phi
2589 {
2590 sidephi = kSPhi ;
2591 }
2592 }
2593 else
2594 {
2595 sphi = kInfinity ;
2596 }
2597 }
2598 else // leaving immediately by ending
2599 {
2600 sidephi = kEPhi ;
2601 sphi = 0 ;
2602 }
2603 }
2604 }
2605 }
2606 else
2607 {
2608 // On z axis + travel not || to z axis -> if phi of vector direction
2609 // within phi of shape, Step limited by rmax, else Step =0
2610
2611 if ( v.x() || v.y() )
2612 {
2613 vphi = std::atan2(v.y(),v.x()) ;
2614 if ((fSPhi-halfAngTolerance < vphi) && (vphi < ePhi+halfAngTolerance))
2615 {
2616 sphi = kInfinity;
2617 }
2618 else
2619 {
2620 sidephi = kSPhi ; // arbitrary
2621 sphi = 0 ;
2622 }
2623 }
2624 else // travel along z - no phi intersection
2625 {
2626 sphi = kInfinity ;
2627 }
2628 }
2629 if ( sphi < snxt ) // Order intersecttions
2630 {
2631 snxt = sphi ;
2632 side = sidephi ;
2633 }
2634 }
2635 if (stheta < snxt ) // Order intersections
2636 {
2637 snxt = stheta ;
2638 side = sidetheta ;
2639 }
2640
2641 if (calcNorm) // Output switch operator
2642 {
2643 switch( side )
2644 {
2645 case kRMax:
2646 xi=p.x()+snxt*v.x();
2647 yi=p.y()+snxt*v.y();
2648 zi=p.z()+snxt*v.z();
2649 *n=G4ThreeVector(xi/fRmax,yi/fRmax,zi/fRmax);
2650 *validNorm=true;
2651 break;
2652
2653 case kRMin:
2654 *validNorm=false; // Rmin is concave
2655 break;
2656
2657 case kSPhi:
2658 if ( fDPhi <= pi ) // Normal to Phi-
2659 {
2660 *n=G4ThreeVector(sinSPhi,-cosSPhi,0);
2661 *validNorm=true;
2662 }
2663 else { *validNorm=false; }
2664 break ;
2665
2666 case kEPhi:
2667 if ( fDPhi <= pi ) // Normal to Phi+
2668 {
2669 *n=G4ThreeVector(-sinEPhi,cosEPhi,0);
2670 *validNorm=true;
2671 }
2672 else { *validNorm=false; }
2673 break;
2674
2675 case kSTheta:
2676 if( fSTheta == halfpi )
2677 {
2678 *n=G4ThreeVector(0.,0.,1.);
2679 *validNorm=true;
2680 }
2681 else if ( fSTheta > halfpi )
2682 {
2683 xi = p.x() + snxt*v.x();
2684 yi = p.y() + snxt*v.y();
2685 rho2=xi*xi+yi*yi;
2686 if (rho2)
2687 {
2688 rhoSecTheta = std::sqrt(rho2*(1+tanSTheta2));
2689 *n = G4ThreeVector( xi/rhoSecTheta, yi/rhoSecTheta,
2690 -tanSTheta/std::sqrt(1+tanSTheta2));
2691 }
2692 else
2693 {
2694 *n = G4ThreeVector(0.,0.,1.);
2695 }
2696 *validNorm=true;
2697 }
2698 else { *validNorm=false; } // Concave STheta cone
2699 break;
2700
2701 case kETheta:
2702 if( eTheta == halfpi )
2703 {
2704 *n = G4ThreeVector(0.,0.,-1.);
2705 *validNorm = true;
2706 }
2707 else if ( eTheta < halfpi )
2708 {
2709 xi=p.x()+snxt*v.x();
2710 yi=p.y()+snxt*v.y();
2711 rho2=xi*xi+yi*yi;
2712 if (rho2)
2713 {
2714 rhoSecTheta = std::sqrt(rho2*(1+tanETheta2));
2715 *n = G4ThreeVector( xi/rhoSecTheta, yi/rhoSecTheta,
2716 -tanETheta/std::sqrt(1+tanETheta2) );
2717 }
2718 else
2719 {
2720 *n = G4ThreeVector(0.,0.,-1.);
2721 }
2722 *validNorm=true;
2723 }
2724 else { *validNorm=false; } // Concave ETheta cone
2725 break;
2726
2727 default:
2728 G4cout << G4endl;
2729 DumpInfo();
2730 std::ostringstream message;
2731 G4int oldprc = message.precision(16);
2732 message << "Undefined side for valid surface normal to solid."
2733 << G4endl
2734 << "Position:" << G4endl << G4endl
2735 << "p.x() = " << p.x()/mm << " mm" << G4endl
2736 << "p.y() = " << p.y()/mm << " mm" << G4endl
2737 << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl
2738 << "Direction:" << G4endl << G4endl
2739 << "v.x() = " << v.x() << G4endl
2740 << "v.y() = " << v.y() << G4endl
2741 << "v.z() = " << v.z() << G4endl << G4endl
2742 << "Proposed distance :" << G4endl << G4endl
2743 << "snxt = " << snxt/mm << " mm" << G4endl;
2744 message.precision(oldprc);
2745 G4Exception("G4Sphere::DistanceToOut(p,v,..)",
2746 "GeomSolids1002", JustWarning, message);
2747 break;
2748 }
2749 }
2750 if (snxt == kInfinity)
2751 {
2752 G4cout << G4endl;
2753 DumpInfo();
2754 std::ostringstream message;
2755 G4int oldprc = message.precision(16);
2756 message << "Logic error: snxt = kInfinity ???" << G4endl
2757 << "Position:" << G4endl << G4endl
2758 << "p.x() = " << p.x()/mm << " mm" << G4endl
2759 << "p.y() = " << p.y()/mm << " mm" << G4endl
2760 << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl
2761 << "Rp = "<< std::sqrt( p.x()*p.x()+p.y()*p.y()+p.z()*p.z() )/mm
2762 << " mm" << G4endl << G4endl
2763 << "Direction:" << G4endl << G4endl
2764 << "v.x() = " << v.x() << G4endl
2765 << "v.y() = " << v.y() << G4endl
2766 << "v.z() = " << v.z() << G4endl << G4endl
2767 << "Proposed distance :" << G4endl << G4endl
2768 << "snxt = " << snxt/mm << " mm" << G4endl;
2769 message.precision(oldprc);
2770 G4Exception("G4Sphere::DistanceToOut(p,v,..)",
2771 "GeomSolids1002", JustWarning, message);
2772 }
2773
2774 return snxt;
2775}
ESide
ESide
Definition: G4Cons.cc:68
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:42

◆ GetCubicVolume()

G4double G4Sphere::GetCubicVolume ( )
inlinevirtual

Reimplemented from G4VSolid.

◆ GetDeltaPhiAngle()

G4double G4Sphere::GetDeltaPhiAngle ( ) const
inline

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

◆ GetDeltaThetaAngle()

G4double G4Sphere::GetDeltaThetaAngle ( ) const
inline

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

◆ GetDPhi()

G4double G4Sphere::GetDPhi ( ) const
inline

◆ GetDTheta()

G4double G4Sphere::GetDTheta ( ) const
inline

◆ GetEntityType()

G4GeometryType G4Sphere::GetEntityType ( ) const
virtual

Implements G4VSolid.

Definition at line 3000 of file G4Sphere.cc.

3001{
3002 return G4String("G4Sphere");
3003}

◆ GetExtent()

G4VisExtent G4Sphere::GetExtent ( ) const
virtual

Reimplemented from G4VSolid.

Definition at line 3171 of file G4Sphere.cc.

3172{
3173 return G4VisExtent(-fRmax, fRmax,-fRmax, fRmax,-fRmax, fRmax );
3174}

◆ GetInnerRadius()

G4double G4Sphere::GetInnerRadius ( ) const
inline

Referenced by G4tgbGeometryDumper::GetSolidParams().

◆ GetInsideRadius()

G4double G4Sphere::GetInsideRadius ( ) const
inline

Referenced by G4PSSphereSurfaceCurrent::IsSelectedSurface(), G4PSSphereSurfaceFlux::IsSelectedSurface(), G4PSSphereSurfaceCurrent::ProcessHits(), G4PSSphereSurfaceFlux::ProcessHits(), G4GDMLWriteParamvol::Sphere_dimensionsWrite(), and G4GDMLWriteSolids::SphereWrite().

◆ GetOuterRadius()

G4double G4Sphere::GetOuterRadius ( ) const
inline

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

◆ GetPointOnSurface()

G4ThreeVector G4Sphere::GetPointOnSurface ( ) const
virtual

Reimplemented from G4VSolid.

Definition at line 3042 of file G4Sphere.cc.

3043{
3044 G4double zRand, aOne, aTwo, aThr, aFou, aFiv, chose, phi, sinphi, cosphi;
3045 G4double height1, height2, slant1, slant2, costheta, sintheta, rRand;
3046
3047 height1 = (fRmax-fRmin)*cosSTheta;
3048 height2 = (fRmax-fRmin)*cosETheta;
3049 slant1 = std::sqrt(sqr((fRmax - fRmin)*sinSTheta) + height1*height1);
3050 slant2 = std::sqrt(sqr((fRmax - fRmin)*sinETheta) + height2*height2);
3051 rRand = GetRadiusInRing(fRmin,fRmax);
3052
3053 aOne = fRmax*fRmax*fDPhi*(cosSTheta-cosETheta);
3054 aTwo = fRmin*fRmin*fDPhi*(cosSTheta-cosETheta);
3055 aThr = fDPhi*((fRmax + fRmin)*sinSTheta)*slant1;
3056 aFou = fDPhi*((fRmax + fRmin)*sinETheta)*slant2;
3057 aFiv = 0.5*fDTheta*(fRmax*fRmax-fRmin*fRmin);
3058
3059 phi = RandFlat::shoot(fSPhi, ePhi);
3060 cosphi = std::cos(phi);
3061 sinphi = std::sin(phi);
3062 costheta = RandFlat::shoot(cosETheta,cosSTheta);
3063 sintheta = std::sqrt(1.-sqr(costheta));
3064
3065 if(fFullPhiSphere) { aFiv = 0; }
3066 if(fSTheta == 0) { aThr=0; }
3067 if(eTheta == pi) { aFou = 0; }
3068 if(fSTheta == halfpi) { aThr = pi*(fRmax*fRmax-fRmin*fRmin); }
3069 if(eTheta == halfpi) { aFou = pi*(fRmax*fRmax-fRmin*fRmin); }
3070
3071 chose = RandFlat::shoot(0.,aOne+aTwo+aThr+aFou+2.*aFiv);
3072 if( (chose>=0.) && (chose<aOne) )
3073 {
3074 return G4ThreeVector(fRmax*sintheta*cosphi,
3075 fRmax*sintheta*sinphi, fRmax*costheta);
3076 }
3077 else if( (chose>=aOne) && (chose<aOne+aTwo) )
3078 {
3079 return G4ThreeVector(fRmin*sintheta*cosphi,
3080 fRmin*sintheta*sinphi, fRmin*costheta);
3081 }
3082 else if( (chose>=aOne+aTwo) && (chose<aOne+aTwo+aThr) )
3083 {
3084 if (fSTheta != halfpi)
3085 {
3086 zRand = RandFlat::shoot(fRmin*cosSTheta,fRmax*cosSTheta);
3087 return G4ThreeVector(tanSTheta*zRand*cosphi,
3088 tanSTheta*zRand*sinphi,zRand);
3089 }
3090 else
3091 {
3092 return G4ThreeVector(rRand*cosphi, rRand*sinphi, 0.);
3093 }
3094 }
3095 else if( (chose>=aOne+aTwo+aThr) && (chose<aOne+aTwo+aThr+aFou) )
3096 {
3097 if(eTheta != halfpi)
3098 {
3099 zRand = RandFlat::shoot(fRmin*cosETheta, fRmax*cosETheta);
3100 return G4ThreeVector (tanETheta*zRand*cosphi,
3101 tanETheta*zRand*sinphi,zRand);
3102 }
3103 else
3104 {
3105 return G4ThreeVector(rRand*cosphi, rRand*sinphi, 0.);
3106 }
3107 }
3108 else if( (chose>=aOne+aTwo+aThr+aFou) && (chose<aOne+aTwo+aThr+aFou+aFiv) )
3109 {
3110 return G4ThreeVector(rRand*sintheta*cosSPhi,
3111 rRand*sintheta*sinSPhi,rRand*costheta);
3112 }
3113 else
3114 {
3115 return G4ThreeVector(rRand*sintheta*cosEPhi,
3116 rRand*sintheta*sinEPhi,rRand*costheta);
3117 }
3118}
CLHEP::RandFlat::shoot
static double shoot()
Definition: RandFlat.cc:59
G4CSGSolid::GetRadiusInRing
G4double GetRadiusInRing(G4double rmin, G4double rmax) const
Definition: G4CSGSolid.cc:101
sqr
T sqr(const T &x)
Definition: templates.hh:145

◆ GetRmax()

G4double G4Sphere::GetRmax ( ) const
inline

◆ GetRmin()

G4double G4Sphere::GetRmin ( ) const
inline

◆ GetSPhi()

G4double G4Sphere::GetSPhi ( ) const
inline

◆ GetStartPhiAngle()

G4double G4Sphere::GetStartPhiAngle ( ) const
inline

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

◆ GetStartThetaAngle()

G4double G4Sphere::GetStartThetaAngle ( ) const
inline

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

◆ GetSTheta()

G4double G4Sphere::GetSTheta ( ) const
inline

◆ GetSurfaceArea()

G4double G4Sphere::GetSurfaceArea ( )
virtual

Reimplemented from G4VSolid.

Definition at line 3124 of file G4Sphere.cc.

3125{
3126 if(fSurfaceArea != 0.) {;}
3127 else
3128 {
3129 G4double Rsq=fRmax*fRmax;
3130 G4double rsq=fRmin*fRmin;
3131
3132 fSurfaceArea = fDPhi*(rsq+Rsq)*(cosSTheta - cosETheta);
3133 if(!fFullPhiSphere)
3134 {
3135 fSurfaceArea = fSurfaceArea + fDTheta*(Rsq-rsq);
3136 }
3137 if(fSTheta >0)
3138 {
3139 G4double acos1=std::acos( std::pow(sinSTheta,2) * std::cos(fDPhi)
3140 + std::pow(cosSTheta,2));
3141 if(fDPhi>pi)
3142 {
3143 fSurfaceArea = fSurfaceArea + 0.5*(Rsq-rsq)*(twopi-acos1);
3144 }
3145 else
3146 {
3147 fSurfaceArea = fSurfaceArea + 0.5*(Rsq-rsq)*acos1;
3148 }
3149 }
3150 if(eTheta < pi)
3151 {
3152 G4double acos2=std::acos( std::pow(sinETheta,2) * std::cos(fDPhi)
3153 + std::pow(cosETheta,2));
3154 if(fDPhi>pi)
3155 {
3156 fSurfaceArea = fSurfaceArea + 0.5*(Rsq-rsq)*(twopi-acos2);
3157 }
3158 else
3159 {
3160 fSurfaceArea = fSurfaceArea + 0.5*(Rsq-rsq)*acos2;
3161 }
3162 }
3163 }
3164 return fSurfaceArea;
3165}
G4CSGSolid::fSurfaceArea
G4double fSurfaceArea
Definition: G4CSGSolid.hh:79

◆ Inside()

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

Implements G4VSolid.

Definition at line 465 of file G4Sphere.cc.

466{
467 G4double rho,rho2,rad2,tolRMin,tolRMax;
468 G4double pPhi,pTheta;
469 EInside in = kOutside;
470 static const G4double halfAngTolerance = kAngTolerance*0.5;
471 const G4double halfRmaxTolerance = fRmaxTolerance*0.5;
472 const G4double halfRminTolerance = fRminTolerance*0.5;
473 const G4double Rmax_minus = fRmax - halfRmaxTolerance;
474 const G4double Rmin_plus = (fRmin > 0) ? fRmin+halfRminTolerance : 0;
475
476 rho2 = p.x()*p.x() + p.y()*p.y() ;
477 rad2 = rho2 + p.z()*p.z() ;
478
479 // Check radial surfaces. Sets 'in'
480
481 tolRMin = Rmin_plus;
482 tolRMax = Rmax_minus;
483
484 if ( (rad2 <= Rmax_minus*Rmax_minus) && (rad2 >= Rmin_plus*Rmin_plus) )
485 {
486 in = kInside;
487 }
488 else
489 {
490 tolRMax = fRmax + halfRmaxTolerance; // outside case
491 tolRMin = std::max(fRmin-halfRminTolerance, 0.); // outside case
492 if ( (rad2 <= tolRMax*tolRMax) && (rad2 >= tolRMin*tolRMin) )
493 {
494 in = kSurface;
495 }
496 else
497 {
498 return in = kOutside;
499 }
500 }
501
502 // Phi boundaries : Do not check if it has no phi boundary!
503
504 if ( !fFullPhiSphere && rho2 ) // [fDPhi < twopi] and [p.x or p.y]
505 {
506 pPhi = std::atan2(p.y(),p.x()) ;
507
508 if ( pPhi < fSPhi - halfAngTolerance ) { pPhi += twopi; }
509 else if ( pPhi > ePhi + halfAngTolerance ) { pPhi -= twopi; }
510
511 if ( (pPhi < fSPhi - halfAngTolerance)
512 || (pPhi > ePhi + halfAngTolerance) ) { return in = kOutside; }
513
514 else if (in == kInside) // else it's kSurface anyway already
515 {
516 if ( (pPhi < fSPhi + halfAngTolerance)
517 || (pPhi > ePhi - halfAngTolerance) ) { in = kSurface; }
518 }
519 }
520
521 // Theta bondaries
522
523 if ( (rho2 || p.z()) && (!fFullThetaSphere) )
524 {
525 rho = std::sqrt(rho2);
526 pTheta = std::atan2(rho,p.z());
527
528 if ( in == kInside )
529 {
530 if ( (pTheta < fSTheta + halfAngTolerance)
531 || (pTheta > eTheta - halfAngTolerance) )
532 {
533 if ( (pTheta >= fSTheta - halfAngTolerance)
534 && (pTheta <= eTheta + halfAngTolerance) )
535 {
536 in = kSurface;
537 }
538 else
539 {
540 in = kOutside;
541 }
542 }
543 }
544 else
545 {
546 if ( (pTheta < fSTheta - halfAngTolerance)
547 || (pTheta > eTheta + halfAngTolerance) )
548 {
549 in = kOutside;
550 }
551 }
552 }
553 return in;
554}
EInside
EInside
Definition: geomdefs.hh:58
kInside
@ kInside
Definition: geomdefs.hh:58
kSurface
@ kSurface
Definition: geomdefs.hh:58

Referenced by CalculateExtent(), and DistanceToOut().

◆ operator=()

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

Definition at line 173 of file G4Sphere.cc.

174{
175 // Check assignment to self
176 //
177 if (this == &rhs) { return *this; }
178
179 // Copy base class data
180 //
181 G4CSGSolid::operator=(rhs);
182
183 // Copy data
184 //
185 fRminTolerance = rhs.fRminTolerance; fRmaxTolerance = rhs.fRmaxTolerance;
186 kAngTolerance = rhs.kAngTolerance; kRadTolerance = rhs.kRadTolerance;
187 fEpsilon = rhs.fEpsilon; fRmin = rhs.fRmin; fRmax = rhs.fRmax;
188 fSPhi = rhs.fSPhi; fDPhi = rhs.fDPhi; fSTheta = rhs.fSTheta;
189 fDTheta = rhs.fDTheta; sinCPhi = rhs.sinCPhi; cosCPhi = rhs.cosCPhi;
190 cosHDPhiOT = rhs.cosHDPhiOT; cosHDPhiIT = rhs.cosHDPhiIT;
191 sinSPhi = rhs.sinSPhi; cosSPhi = rhs.cosSPhi;
192 sinEPhi = rhs.sinEPhi; cosEPhi = rhs.cosEPhi;
193 hDPhi = rhs.hDPhi; cPhi = rhs.cPhi; ePhi = rhs.ePhi;
194 sinSTheta = rhs.sinSTheta; cosSTheta = rhs.cosSTheta;
195 sinETheta = rhs.sinETheta; cosETheta = rhs.cosETheta;
196 tanSTheta = rhs.tanSTheta; tanSTheta2 = rhs.tanSTheta2;
197 tanETheta = rhs.tanETheta; tanETheta2 = rhs.tanETheta2;
198 eTheta = rhs.eTheta; fFullPhiSphere = rhs.fFullPhiSphere;
199 fFullThetaSphere = rhs.fFullThetaSphere; fFullSphere = rhs.fFullSphere;
200
201 return *this;
202}
G4CSGSolid::operator=
G4CSGSolid & operator=(const G4CSGSolid &rhs)
Definition: G4CSGSolid.cc:82

◆ SetDeltaPhiAngle()

void G4Sphere::SetDeltaPhiAngle ( G4double  newDphi)
inline

◆ SetDeltaThetaAngle()

void G4Sphere::SetDeltaThetaAngle ( G4double  newDTheta)
inline

◆ SetInnerRadius()

void G4Sphere::SetInnerRadius ( G4double  newRMin)
inline

◆ SetInsideRadius()

void G4Sphere::SetInsideRadius ( 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
virtual

Reimplemented from G4CSGSolid.

Definition at line 3018 of file G4Sphere.cc.

3019{
3020 G4int oldprc = os.precision(16);
3021 os << "-----------------------------------------------------------\n"
3022 << " *** Dump for solid - " << GetName() << " ***\n"
3023 << " ===================================================\n"
3024 << " Solid type: G4Sphere\n"
3025 << " Parameters: \n"
3026 << " inner radius: " << fRmin/mm << " mm \n"
3027 << " outer radius: " << fRmax/mm << " mm \n"
3028 << " starting phi of segment : " << fSPhi/degree << " degrees \n"
3029 << " delta phi of segment : " << fDPhi/degree << " degrees \n"
3030 << " starting theta of segment: " << fSTheta/degree << " degrees \n"
3031 << " delta theta of segment : " << fDTheta/degree << " degrees \n"
3032 << "-----------------------------------------------------------\n";
3033 os.precision(oldprc);
3034
3035 return os;
3036}

◆ SurfaceNormal()

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

Implements G4VSolid.

Definition at line 562 of file G4Sphere.cc.

563{
564 G4int noSurfaces = 0;
565 G4double rho, rho2, radius, pTheta, pPhi=0.;
566 G4double distRMin = kInfinity;
567 G4double distSPhi = kInfinity, distEPhi = kInfinity;
568 G4double distSTheta = kInfinity, distETheta = kInfinity;
569 G4ThreeVector nR, nPs, nPe, nTs, nTe, nZ(0.,0.,1.);
570 G4ThreeVector norm, sumnorm(0.,0.,0.);
571
572 static const G4double halfCarTolerance = 0.5*kCarTolerance;
573 static const G4double halfAngTolerance = 0.5*kAngTolerance;
574
575 rho2 = p.x()*p.x()+p.y()*p.y();
576 radius = std::sqrt(rho2+p.z()*p.z());
577 rho = std::sqrt(rho2);
578
579 G4double distRMax = std::fabs(radius-fRmax);
580 if (fRmin) distRMin = std::fabs(radius-fRmin);
581
582 if ( rho && !fFullSphere )
583 {
584 pPhi = std::atan2(p.y(),p.x());
585
586 if (pPhi < fSPhi-halfAngTolerance) { pPhi += twopi; }
587 else if (pPhi > ePhi+halfAngTolerance) { pPhi -= twopi; }
588 }
589 if ( !fFullPhiSphere )
590 {
591 if ( rho )
592 {
593 distSPhi = std::fabs( pPhi-fSPhi );
594 distEPhi = std::fabs( pPhi-ePhi );
595 }
596 else if( !fRmin )
597 {
598 distSPhi = 0.;
599 distEPhi = 0.;
600 }
601 nPs = G4ThreeVector(sinSPhi,-cosSPhi,0);
602 nPe = G4ThreeVector(-sinEPhi,cosEPhi,0);
603 }
604 if ( !fFullThetaSphere )
605 {
606 if ( rho )
607 {
608 pTheta = std::atan2(rho,p.z());
609 distSTheta = std::fabs(pTheta-fSTheta);
610 distETheta = std::fabs(pTheta-eTheta);
611
612 nTs = G4ThreeVector(-cosSTheta*p.x()/rho,
613 -cosSTheta*p.y()/rho,
614 sinSTheta );
615
616 nTe = G4ThreeVector( cosETheta*p.x()/rho,
617 cosETheta*p.y()/rho,
618 -sinETheta );
619 }
620 else if( !fRmin )
621 {
622 if ( fSTheta )
623 {
624 distSTheta = 0.;
625 nTs = G4ThreeVector(0.,0.,-1.);
626 }
627 if ( eTheta < pi )
628 {
629 distETheta = 0.;
630 nTe = G4ThreeVector(0.,0.,1.);
631 }
632 }
633 }
634 if( radius ) { nR = G4ThreeVector(p.x()/radius,p.y()/radius,p.z()/radius); }
635
636 if( distRMax <= halfCarTolerance )
637 {
638 noSurfaces ++;
639 sumnorm += nR;
640 }
641 if( fRmin && (distRMin <= halfCarTolerance) )
642 {
643 noSurfaces ++;
644 sumnorm -= nR;
645 }
646 if( !fFullPhiSphere )
647 {
648 if (distSPhi <= halfAngTolerance)
649 {
650 noSurfaces ++;
651 sumnorm += nPs;
652 }
653 if (distEPhi <= halfAngTolerance)
654 {
655 noSurfaces ++;
656 sumnorm += nPe;
657 }
658 }
659 if ( !fFullThetaSphere )
660 {
661 if ((distSTheta <= halfAngTolerance) && (fSTheta > 0.))
662 {
663 noSurfaces ++;
664 if ((radius <= halfCarTolerance) && fFullPhiSphere) { sumnorm += nZ; }
665 else { sumnorm += nTs; }
666 }
667 if ((distETheta <= halfAngTolerance) && (eTheta < pi))
668 {
669 noSurfaces ++;
670 if ((radius <= halfCarTolerance) && fFullPhiSphere) { sumnorm -= nZ; }
671 else { sumnorm += nTe; }
672 if(sumnorm.z() == 0.) { sumnorm += nZ; }
673 }
674 }
675 if ( noSurfaces == 0 )
676 {
677#ifdef G4CSGDEBUG
678 G4Exception("G4Sphere::SurfaceNormal(p)", "GeomSolids1002",
679 JustWarning, "Point p is not on surface !?" );
680#endif
681 norm = ApproxSurfaceNormal(p);
682 }
683 else if ( noSurfaces == 1 ) { norm = sumnorm; }
684 else { norm = sumnorm.unit(); }
685 return norm;
686}
CLHEP::Hep3Vector::unit
Hep3Vector unit() const
The documentation for this class was generated from the following files: