Geant4.10
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Public Member Functions
G4Sphere Class Reference

#include <G4Sphere.hh>

Inheritance diagram for G4Sphere:
G4CSGSolid G4VSolid

Public Member Functions

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

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 86 of file G4Sphere.hh.

Constructor & Destructor Documentation

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.

References FatalException, G4endl, G4Exception(), G4GeometryTolerance::GetAngularTolerance(), G4GeometryTolerance::GetInstance(), G4VSolid::GetName(), G4GeometryTolerance::GetRadialTolerance(), G4VSolid::kCarTolerance, and G4INCL::Math::max().

Referenced by Clone().

94  : G4CSGSolid(pName), fEpsilon(2.e-11), fSPhi(0.0),
95  fFullPhiSphere(true), fFullThetaSphere(true)
96 {
97  kAngTolerance = G4GeometryTolerance::GetInstance()->GetAngularTolerance();
98  kRadTolerance = G4GeometryTolerance::GetInstance()->GetRadialTolerance();
99 
100  halfCarTolerance = 0.5*kCarTolerance;
101  halfAngTolerance = 0.5*kAngTolerance;
102 
103  // Check radii and set radial tolerances
104 
105  if ( (pRmin >= pRmax) || (pRmax < 1.1*kRadTolerance) || (pRmin < 0) )
106  {
107  std::ostringstream message;
108  message << "Invalid radii for Solid: " << GetName() << G4endl
109  << " pRmin = " << pRmin << ", pRmax = " << pRmax;
110  G4Exception("G4Sphere::G4Sphere()", "GeomSolids0002",
111  FatalException, message);
112  }
113  fRmin=pRmin; fRmax=pRmax;
114  fRminTolerance = (fRmin) ? std::max( kRadTolerance, fEpsilon*fRmin ) : 0;
115  fRmaxTolerance = std::max( kRadTolerance, fEpsilon*fRmax );
116 
117  // Check angles
118 
119  CheckPhiAngles(pSPhi, pDPhi);
120  CheckThetaAngles(pSTheta, pDTheta);
121 }
G4VSolid::GetName
G4String GetName() const
G4GeometryTolerance::GetRadialTolerance
G4double GetRadialTolerance() const
Definition: G4GeometryTolerance.cc:97
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4CSGSolid::G4CSGSolid
G4CSGSolid(const G4String &pName)
Definition: G4CSGSolid.cc:42
G4INCL::Math::max
T max(const T t1, const T t2)
brief Return the largest of the two arguments
Definition: G4INCLGlobals.hh:116
FatalException
Definition: G4ExceptionSeverity.hh:60
G4endl
#define G4endl
Definition: G4ios.hh:61
G4VSolid::kCarTolerance
G4double kCarTolerance
Definition: G4VSolid.hh:305
G4GeometryTolerance::GetAngularTolerance
G4double GetAngularTolerance() const
Definition: G4GeometryTolerance.cc:92
G4GeometryTolerance::GetInstance
static G4GeometryTolerance * GetInstance()
Definition: G4GeometryTolerance.cc:74
G4Sphere::~G4Sphere ( )

Definition at line 145 of file G4Sphere.cc.

146 {
147 }
G4Sphere::G4Sphere ( __void__ &  a)

Definition at line 128 of file G4Sphere.cc.

129  : G4CSGSolid(a), fRminTolerance(0.), fRmaxTolerance(0.),
130  kAngTolerance(0.), kRadTolerance(0.), fEpsilon(0.),
131  fRmin(0.), fRmax(0.), fSPhi(0.), fDPhi(0.), fSTheta(0.),
132  fDTheta(0.), sinCPhi(0.), cosCPhi(0.), cosHDPhiOT(0.), cosHDPhiIT(0.),
133  sinSPhi(0.), cosSPhi(0.), sinEPhi(0.), cosEPhi(0.), hDPhi(0.), cPhi(0.),
134  ePhi(0.), sinSTheta(0.), cosSTheta(0.), sinETheta(0.), cosETheta(0.),
135  tanSTheta(0.), tanSTheta2(0.), tanETheta(0.), tanETheta2(0.), eTheta(0.),
136  fFullPhiSphere(false), fFullThetaSphere(false), fFullSphere(true),
137  halfCarTolerance(0.), halfAngTolerance(0.)
138 {
139 }
G4CSGSolid::G4CSGSolid
G4CSGSolid(const G4String &pName)
Definition: G4CSGSolid.cc:42
G4Sphere::G4Sphere ( const G4Sphere &  rhs)

Definition at line 153 of file G4Sphere.cc.

154  : G4CSGSolid(rhs), fRminTolerance(rhs.fRminTolerance),
155  fRmaxTolerance(rhs.fRmaxTolerance), kAngTolerance(rhs.kAngTolerance),
156  kRadTolerance(rhs.kRadTolerance), fEpsilon(rhs.fEpsilon),
157  fRmin(rhs.fRmin), fRmax(rhs.fRmax), fSPhi(rhs.fSPhi), fDPhi(rhs.fDPhi),
158  fSTheta(rhs.fSTheta), fDTheta(rhs.fDTheta),
159  sinCPhi(rhs.sinCPhi), cosCPhi(rhs.cosCPhi),
160  cosHDPhiOT(rhs.cosHDPhiOT), cosHDPhiIT(rhs.cosHDPhiIT),
161  sinSPhi(rhs.sinSPhi), cosSPhi(rhs.cosSPhi),
162  sinEPhi(rhs.sinEPhi), cosEPhi(rhs.cosEPhi),
163  hDPhi(rhs.hDPhi), cPhi(rhs.cPhi), ePhi(rhs.ePhi),
164  sinSTheta(rhs.sinSTheta), cosSTheta(rhs.cosSTheta),
165  sinETheta(rhs.sinETheta), cosETheta(rhs.cosETheta),
166  tanSTheta(rhs.tanSTheta), tanSTheta2(rhs.tanSTheta2),
167  tanETheta(rhs.tanETheta), tanETheta2(rhs.tanETheta2), eTheta(rhs.eTheta),
168  fFullPhiSphere(rhs.fFullPhiSphere), fFullThetaSphere(rhs.fFullThetaSphere),
169  fFullSphere(rhs.fFullSphere),
170  halfCarTolerance(rhs.halfCarTolerance),
171  halfAngTolerance(rhs.halfAngTolerance)
172 {
173 }
G4CSGSolid::G4CSGSolid
G4CSGSolid(const G4String &pName)
Definition: G4CSGSolid.cc:42

Member Function Documentation

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

Implements G4VSolid.

Definition at line 228 of file G4Sphere.cc.

References G4VSolid::CalculateClippedPolygonExtent(), G4VoxelLimits::GetMaxExtent(), G4VoxelLimits::GetMaxXExtent(), G4VoxelLimits::GetMaxYExtent(), G4VoxelLimits::GetMaxZExtent(), G4VoxelLimits::GetMinExtent(), G4VoxelLimits::GetMinXExtent(), G4VoxelLimits::GetMinYExtent(), G4VoxelLimits::GetMinZExtent(), Inside(), G4AffineTransform::Inverse(), G4VoxelLimits::IsXLimited(), G4VoxelLimits::IsYLimited(), G4VoxelLimits::IsZLimited(), G4VSolid::kCarTolerance, kOutside, kXAxis, kYAxis, kZAxis, G4AffineTransform::NetTranslation(), G4AffineTransform::TransformPoint(), CLHEP::Hep3Vector::x(), CLHEP::Hep3Vector::y(), and CLHEP::Hep3Vector::z().

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

Reimplemented from G4VSolid.

Definition at line 3010 of file G4Sphere.cc.

References G4Sphere().

3011 {
3012  return new G4Sphere(*this);
3013 }
G4Sphere::G4Sphere
G4Sphere(const G4String &pName, G4double pRmin, G4double pRmax, G4double pSPhi, G4double pDPhi, G4double pSTheta, G4double pDTheta)
Definition: G4Sphere.cc:90
void G4Sphere::ComputeDimensions ( G4VPVParameterisation *  p,
const G4int  n,
const G4VPhysicalVolume *  pRep 
)
virtual

Reimplemented from G4VSolid.

Definition at line 217 of file G4Sphere.cc.

References G4VPVParameterisation::ComputeDimensions().

220 {
221  p->ComputeDimensions(*this,n,pRep);
222 }
n
const G4int n
Definition: G4UrQMD1_3Interface.hh:144
G4VPVParameterisation::ComputeDimensions
virtual void ComputeDimensions(G4Box &, const G4int, const G4VPhysicalVolume *) const
Definition: G4VPVParameterisation.hh:95
G4Polyhedron * G4Sphere::CreatePolyhedron ( ) const
virtual

Reimplemented from G4VSolid.

Definition at line 3183 of file G4Sphere.cc.

3184 {
3185  return new G4PolyhedronSphere (fRmin, fRmax, fSPhi, fDPhi, fSTheta, fDTheta);
3186 }
void G4Sphere::DescribeYourselfTo ( G4VGraphicsScene &  scene) const
virtual

Implements G4VSolid.

Definition at line 3178 of file G4Sphere.cc.

References G4VGraphicsScene::AddSolid().

3179 {
3180  scene.AddSolid (*this);
3181 }
G4VGraphicsScene::AddSolid
virtual void AddSolid(const G4Box &)=0
G4double G4Sphere::DistanceToIn ( const G4ThreeVector &  p,
const G4ThreeVector &  v 
) const
virtual

Implements G4VSolid.

Definition at line 872 of file G4Sphere.cc.

References test::b, test::c, python.hepunit::halfpi, G4VSolid::kCarTolerance, python.hepunit::pi, plottest35::t1, CLHEP::Hep3Vector::x(), CLHEP::Hep3Vector::y(), and CLHEP::Hep3Vector::z().

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

Implements G4VSolid.

Definition at line 1804 of file G4Sphere.cc.

References python.hepunit::pi, CLHEP::Hep3Vector::x(), CLHEP::Hep3Vector::y(), and CLHEP::Hep3Vector::z().

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

References test::b, test::c, G4VSolid::DumpInfo(), G4cout, G4endl, G4Exception(), python.hepunit::halfpi, JustWarning, G4VSolid::kCarTolerance, python.hepunit::mm, python.hepunit::pi, plottest35::t1, CLHEP::Hep3Vector::x(), CLHEP::Hep3Vector::y(), and CLHEP::Hep3Vector::z().

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

Implements G4VSolid.

Definition at line 2782 of file G4Sphere.cc.

References G4VSolid::DumpInfo(), G4cout, G4endl, G4Exception(), Inside(), JustWarning, kOutside, python.hepunit::mm, python.hepunit::pi, CLHEP::Hep3Vector::x(), CLHEP::Hep3Vector::y(), and CLHEP::Hep3Vector::z().

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

Reimplemented from G4VSolid.

G4double G4Sphere::GetDeltaPhiAngle ( ) const
inline

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

G4double G4Sphere::GetDeltaThetaAngle ( ) const
inline

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

G4double G4Sphere::GetDPhi ( ) const
inline
G4double G4Sphere::GetDTheta ( ) const
inline
G4GeometryType G4Sphere::GetEntityType ( ) const
virtual

Implements G4VSolid.

Definition at line 3001 of file G4Sphere.cc.

3002 {
3003  return G4String("G4Sphere");
3004 }
G4VisExtent G4Sphere::GetExtent ( ) const
virtual

Reimplemented from G4VSolid.

Definition at line 3172 of file G4Sphere.cc.

3173 {
3174  return G4VisExtent(-fRmax, fRmax,-fRmax, fRmax,-fRmax, fRmax );
3175 }
G4double G4Sphere::GetInnerRadius ( ) const
inline

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

G4double G4Sphere::GetInsideRadius ( ) const
inline

Referenced by export_G4Sphere().

G4double G4Sphere::GetOuterRadius ( ) const
inline

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

G4ThreeVector G4Sphere::GetPointOnSurface ( ) const
virtual

Reimplemented from G4VSolid.

Definition at line 3043 of file G4Sphere.cc.

References G4CSGSolid::GetRadiusInRing(), python.hepunit::halfpi, python.hepunit::pi, G4INCL::DeJongSpin::shoot(), and sqr().

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

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

G4double G4Sphere::GetStartThetaAngle ( ) const
inline

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

G4double G4Sphere::GetSTheta ( ) const
inline
G4double G4Sphere::GetSurfaceArea ( )
virtual

Reimplemented from G4VSolid.

Definition at line 3125 of file G4Sphere.cc.

References G4CSGSolid::fSurfaceArea, python.hepunit::pi, and python.hepunit::twopi.

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

Implements G4VSolid.

Definition at line 473 of file G4Sphere.cc.

References kInside, kOutside, kSurface, G4INCL::Math::max(), python.hepunit::twopi, CLHEP::Hep3Vector::x(), CLHEP::Hep3Vector::y(), and CLHEP::Hep3Vector::z().

Referenced by CalculateExtent(), and DistanceToOut().

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

Definition at line 179 of file G4Sphere.cc.

References G4CSGSolid::operator=().

180 {
181  // Check assignment to self
182  //
183  if (this == &rhs) { return *this; }
184 
185  // Copy base class data
186  //
187  G4CSGSolid::operator=(rhs);
188 
189  // Copy data
190  //
191  fRminTolerance = rhs.fRminTolerance; fRmaxTolerance = rhs.fRmaxTolerance;
192  kAngTolerance = rhs.kAngTolerance; kRadTolerance = rhs.kRadTolerance;
193  fEpsilon = rhs.fEpsilon; fRmin = rhs.fRmin; fRmax = rhs.fRmax;
194  fSPhi = rhs.fSPhi; fDPhi = rhs.fDPhi; fSTheta = rhs.fSTheta;
195  fDTheta = rhs.fDTheta; sinCPhi = rhs.sinCPhi; cosCPhi = rhs.cosCPhi;
196  cosHDPhiOT = rhs.cosHDPhiOT; cosHDPhiIT = rhs.cosHDPhiIT;
197  sinSPhi = rhs.sinSPhi; cosSPhi = rhs.cosSPhi;
198  sinEPhi = rhs.sinEPhi; cosEPhi = rhs.cosEPhi;
199  hDPhi = rhs.hDPhi; cPhi = rhs.cPhi; ePhi = rhs.ePhi;
200  sinSTheta = rhs.sinSTheta; cosSTheta = rhs.cosSTheta;
201  sinETheta = rhs.sinETheta; cosETheta = rhs.cosETheta;
202  tanSTheta = rhs.tanSTheta; tanSTheta2 = rhs.tanSTheta2;
203  tanETheta = rhs.tanETheta; tanETheta2 = rhs.tanETheta2;
204  eTheta = rhs.eTheta; fFullPhiSphere = rhs.fFullPhiSphere;
205  fFullThetaSphere = rhs.fFullThetaSphere; fFullSphere = rhs.fFullSphere;
206  halfCarTolerance = rhs.halfCarTolerance;
207  halfAngTolerance = rhs.halfAngTolerance;
208 
209  return *this;
210 }
G4CSGSolid::operator=
G4CSGSolid & operator=(const G4CSGSolid &rhs)
Definition: G4CSGSolid.cc:82
void G4Sphere::SetDeltaPhiAngle ( G4double  newDphi)
inline

Referenced by export_G4Sphere().

void G4Sphere::SetDeltaThetaAngle ( G4double  newDTheta)
inline

Referenced by export_G4Sphere().

void G4Sphere::SetInnerRadius ( G4double  newRMin)
inline
void G4Sphere::SetInsideRadius ( G4double  newRmin)
inline

Referenced by export_G4Sphere().

void G4Sphere::SetOuterRadius ( G4double  newRmax)
inline

Referenced by export_G4Sphere().

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

Referenced by export_G4Sphere().

void G4Sphere::SetStartThetaAngle ( G4double  newSTheta)
inline

Referenced by export_G4Sphere().

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

Reimplemented from G4CSGSolid.

Definition at line 3019 of file G4Sphere.cc.

References python.hepunit::degree, G4VSolid::GetName(), and python.hepunit::mm.

3020 {
3021  G4int oldprc = os.precision(16);
3022  os << "-----------------------------------------------------------\n"
3023  << " *** Dump for solid - " << GetName() << " ***\n"
3024  << " ===================================================\n"
3025  << " Solid type: G4Sphere\n"
3026  << " Parameters: \n"
3027  << " inner radius: " << fRmin/mm << " mm \n"
3028  << " outer radius: " << fRmax/mm << " mm \n"
3029  << " starting phi of segment : " << fSPhi/degree << " degrees \n"
3030  << " delta phi of segment : " << fDPhi/degree << " degrees \n"
3031  << " starting theta of segment: " << fSTheta/degree << " degrees \n"
3032  << " delta theta of segment : " << fDTheta/degree << " degrees \n"
3033  << "-----------------------------------------------------------\n";
3034  os.precision(oldprc);
3035 
3036  return os;
3037 }
G4VSolid::GetName
G4String GetName() const
G4int
int G4int
Definition: G4Types.hh:78
python.hepunit.degree
tuple degree
Definition: hepunit.py:69
G4ThreeVector G4Sphere::SurfaceNormal ( const G4ThreeVector &  p) const
virtual

Implements G4VSolid.

Definition at line 570 of file G4Sphere.cc.

References G4Exception(), JustWarning, python.hepunit::pi, python.hepunit::twopi, CLHEP::Hep3Vector::unit(), CLHEP::Hep3Vector::x(), CLHEP::Hep3Vector::y(), and CLHEP::Hep3Vector::z().

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