Geant4  10.02.p03
G4Cons Class Reference

#include <G4Cons.hh>

Inheritance diagram for G4Cons:
Collaboration diagram for G4Cons:

Public Member Functions

 G4Cons (const G4String &pName, G4double pRmin1, G4double pRmax1, G4double pRmin2, G4double pRmax2, G4double pDz, G4double pSPhi, G4double pDPhi)
 
 ~G4Cons ()
 
G4double GetInnerRadiusMinusZ () const
 
G4double GetOuterRadiusMinusZ () const
 
G4double GetInnerRadiusPlusZ () const
 
G4double GetOuterRadiusPlusZ () const
 
G4double GetZHalfLength () const
 
G4double GetStartPhiAngle () const
 
G4double GetDeltaPhiAngle () const
 
void SetInnerRadiusMinusZ (G4double Rmin1)
 
void SetOuterRadiusMinusZ (G4double Rmax1)
 
void SetInnerRadiusPlusZ (G4double Rmin2)
 
void SetOuterRadiusPlusZ (G4double Rmax2)
 
void SetZHalfLength (G4double newDz)
 
void SetStartPhiAngle (G4double newSPhi, G4bool trig=true)
 
void SetDeltaPhiAngle (G4double newDPhi)
 
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
 
G4VSolidClone () const
 
std::ostream & StreamInfo (std::ostream &os) const
 
void DescribeYourselfTo (G4VGraphicsScene &scene) const
 
G4PolyhedronCreatePolyhedron () const
 
 G4Cons (__void__ &)
 
 G4Cons (const G4Cons &rhs)
 
G4Consoperator= (const G4Cons &rhs)
 
G4double GetRmin1 () const
 
G4double GetRmax1 () const
 
G4double GetRmin2 () const
 
G4double GetRmax2 () const
 
G4double GetDz () const
 
G4double GetSPhi () const
 
G4double GetDPhi () const
 
- Public Member Functions inherited from G4CSGSolid
 G4CSGSolid (const G4String &pName)
 
virtual ~G4CSGSolid ()
 
virtual G4PolyhedronGetPolyhedron () const
 
 G4CSGSolid (__void__ &)
 
 G4CSGSolid (const G4CSGSolid &rhs)
 
G4CSGSolidoperator= (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 G4VisExtent GetExtent () const
 
virtual const G4VSolidGetConstituentSolid (G4int no) const
 
virtual G4VSolidGetConstituentSolid (G4int no)
 
virtual const G4DisplacedSolidGetDisplacedSolidPtr () const
 
virtual G4DisplacedSolidGetDisplacedSolidPtr ()
 
 G4VSolid (__void__ &)
 
 G4VSolid (const G4VSolid &rhs)
 
G4VSolidoperator= (const G4VSolid &rhs)
 
G4double EstimateCubicVolume (G4int nStat, G4double epsilon) const
 
G4double EstimateSurfaceArea (G4int nStat, G4double ell) const
 

Private Types

enum  ESide {
  kNull, kRMin, kRMax, kSPhi,
  kEPhi, kPZ, kMZ
}
 
enum  ENorm {
  kNRMin, kNRMax, kNSPhi, kNEPhi,
  kNZ
}
 

Private Member Functions

G4ThreeVectorListCreateRotatedVertices (const G4AffineTransform &pTransform) const
 
void Initialize ()
 
void CheckSPhiAngle (G4double sPhi)
 
void CheckDPhiAngle (G4double dPhi)
 
void CheckPhiAngles (G4double sPhi, G4double dPhi)
 
void InitializeTrigonometry ()
 
G4ThreeVector ApproxSurfaceNormal (const G4ThreeVector &p) const
 

Private Attributes

G4double kRadTolerance
 
G4double kAngTolerance
 
G4double fRmin1
 
G4double fRmin2
 
G4double fRmax1
 
G4double fRmax2
 
G4double fDz
 
G4double fSPhi
 
G4double fDPhi
 
G4double sinCPhi
 
G4double cosCPhi
 
G4double cosHDPhiOT
 
G4double cosHDPhiIT
 
G4double sinSPhi
 
G4double cosSPhi
 
G4double sinEPhi
 
G4double cosEPhi
 
G4bool fPhiFullCone
 
G4double halfCarTolerance
 
G4double halfRadTolerance
 
G4double halfAngTolerance
 

Additional Inherited Members

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

Detailed Description

Definition at line 83 of file G4Cons.hh.

Member Enumeration Documentation

◆ ENorm

enum G4Cons::ENorm
private
Enumerator
kNRMin 
kNRMax 
kNSPhi 
kNEPhi 
kNZ 

Definition at line 214 of file G4Cons.hh.

◆ ESide

enum G4Cons::ESide
private
Enumerator
kNull 
kRMin 
kRMax 
kSPhi 
kEPhi 
kPZ 
kMZ 

Definition at line 210 of file G4Cons.hh.

Constructor & Destructor Documentation

◆ G4Cons() [1/3]

G4Cons::G4Cons ( const G4String pName,
G4double  pRmin1,
G4double  pRmax1,
G4double  pRmin2,
G4double  pRmax2,
G4double  pDz,
G4double  pSPhi,
G4double  pDPhi 
)

Definition at line 82 of file G4Cons.cc.

87  : G4CSGSolid(pName), fRmin1(pRmin1), fRmin2(pRmin2),
88  fRmax1(pRmax1), fRmax2(pRmax2), fDz(pDz), fSPhi(0.), fDPhi(0.)
89 {
92 
96 
97  // Check z-len
98  //
99  if ( pDz < 0 )
100  {
101  std::ostringstream message;
102  message << "Invalid Z half-length for Solid: " << GetName() << G4endl
103  << " hZ = " << pDz;
104  G4Exception("G4Cons::G4Cons()", "GeomSolids0002",
105  FatalException, message);
106  }
107 
108  // Check radii
109  //
110  if (((pRmin1>=pRmax1) || (pRmin2>=pRmax2) || (pRmin1<0)) && (pRmin2<0))
111  {
112  std::ostringstream message;
113  message << "Invalid values of radii for Solid: " << GetName() << G4endl
114  << " pRmin1 = " << pRmin1 << ", pRmin2 = " << pRmin2
115  << ", pRmax1 = " << pRmax1 << ", pRmax2 = " << pRmax2;
116  G4Exception("G4Cons::G4Cons()", "GeomSolids0002",
117  FatalException, message) ;
118  }
119  if( (pRmin1 == 0.0) && (pRmin2 > 0.0) ) { fRmin1 = 1e3*kRadTolerance ; }
120  if( (pRmin2 == 0.0) && (pRmin1 > 0.0) ) { fRmin2 = 1e3*kRadTolerance ; }
121 
122  // Check angles
123  //
124  CheckPhiAngles(pSPhi, pDPhi);
125 }
G4double halfRadTolerance
Definition: G4Cons.hh:233
G4double fDz
Definition: G4Cons.hh:220
G4double fRmin1
Definition: G4Cons.hh:220
G4double fRmax2
Definition: G4Cons.hh:220
G4String GetName() const
void CheckPhiAngles(G4double sPhi, G4double dPhi)
G4double kRadTolerance
Definition: G4Cons.hh:216
G4double fRmin2
Definition: G4Cons.hh:220
G4double GetRadialTolerance() const
G4double fSPhi
Definition: G4Cons.hh:220
G4double halfCarTolerance
Definition: G4Cons.hh:233
G4double GetAngularTolerance() const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4CSGSolid(const G4String &pName)
Definition: G4CSGSolid.cc:49
G4double fDPhi
Definition: G4Cons.hh:220
G4double halfAngTolerance
Definition: G4Cons.hh:233
G4double kAngTolerance
Definition: G4Cons.hh:216
#define G4endl
Definition: G4ios.hh:61
G4double kCarTolerance
Definition: G4VSolid.hh:304
static const G4double e3
static G4GeometryTolerance * GetInstance()
G4double fRmax1
Definition: G4Cons.hh:220
Here is the call graph for this function:
Here is the caller graph for this function:

◆ ~G4Cons()

G4Cons::~G4Cons ( )

Definition at line 146 of file G4Cons.cc.

147 {
148 }

◆ G4Cons() [2/3]

G4Cons::G4Cons ( __void__ &  a)

Definition at line 132 of file G4Cons.cc.

134  fRmin1(0.), fRmin2(0.), fRmax1(0.), fRmax2(0.), fDz(0.),
135  fSPhi(0.), fDPhi(0.), sinCPhi(0.), cosCPhi(0.), cosHDPhiOT(0.),
136  cosHDPhiIT(0.), sinSPhi(0.), cosSPhi(0.), sinEPhi(0.), cosEPhi(0.),
138  halfAngTolerance(0.)
139 {
140 }
G4double halfRadTolerance
Definition: G4Cons.hh:233
G4double fDz
Definition: G4Cons.hh:220
G4bool fPhiFullCone
Definition: G4Cons.hh:229
G4double sinEPhi
Definition: G4Cons.hh:224
G4double fRmin1
Definition: G4Cons.hh:220
G4double fRmax2
Definition: G4Cons.hh:220
G4double cosCPhi
Definition: G4Cons.hh:224
G4double kRadTolerance
Definition: G4Cons.hh:216
G4double fRmin2
Definition: G4Cons.hh:220
G4double sinSPhi
Definition: G4Cons.hh:224
G4double fSPhi
Definition: G4Cons.hh:220
G4double halfCarTolerance
Definition: G4Cons.hh:233
G4double cosSPhi
Definition: G4Cons.hh:224
G4CSGSolid(const G4String &pName)
Definition: G4CSGSolid.cc:49
G4double fDPhi
Definition: G4Cons.hh:220
G4double halfAngTolerance
Definition: G4Cons.hh:233
G4double sinCPhi
Definition: G4Cons.hh:224
G4double kAngTolerance
Definition: G4Cons.hh:216
G4double cosEPhi
Definition: G4Cons.hh:224
G4double cosHDPhiOT
Definition: G4Cons.hh:224
G4double cosHDPhiIT
Definition: G4Cons.hh:224
G4double fRmax1
Definition: G4Cons.hh:220

◆ G4Cons() [3/3]

G4Cons::G4Cons ( const G4Cons rhs)

Definition at line 154 of file G4Cons.cc.

157  fRmax1(rhs.fRmax1), fRmax2(rhs.fRmax2), fDz(rhs.fDz), fSPhi(rhs.fSPhi),
158  fDPhi(rhs.fDPhi), sinCPhi(rhs.sinCPhi), cosCPhi(rhs.cosCPhi),
160  sinSPhi(rhs.sinSPhi), cosSPhi(rhs.cosSPhi), sinEPhi(rhs.sinEPhi),
165 {
166 }
G4double halfRadTolerance
Definition: G4Cons.hh:233
G4double fDz
Definition: G4Cons.hh:220
G4bool fPhiFullCone
Definition: G4Cons.hh:229
G4double sinEPhi
Definition: G4Cons.hh:224
G4double fRmin1
Definition: G4Cons.hh:220
G4double fRmax2
Definition: G4Cons.hh:220
G4double cosCPhi
Definition: G4Cons.hh:224
G4double kRadTolerance
Definition: G4Cons.hh:216
G4double fRmin2
Definition: G4Cons.hh:220
G4double sinSPhi
Definition: G4Cons.hh:224
G4double fSPhi
Definition: G4Cons.hh:220
G4double halfCarTolerance
Definition: G4Cons.hh:233
G4double cosSPhi
Definition: G4Cons.hh:224
G4CSGSolid(const G4String &pName)
Definition: G4CSGSolid.cc:49
G4double fDPhi
Definition: G4Cons.hh:220
G4double halfAngTolerance
Definition: G4Cons.hh:233
G4double sinCPhi
Definition: G4Cons.hh:224
G4double kAngTolerance
Definition: G4Cons.hh:216
G4double cosEPhi
Definition: G4Cons.hh:224
G4double cosHDPhiOT
Definition: G4Cons.hh:224
G4double cosHDPhiIT
Definition: G4Cons.hh:224
G4double fRmax1
Definition: G4Cons.hh:220

Member Function Documentation

◆ ApproxSurfaceNormal()

G4ThreeVector G4Cons::ApproxSurfaceNormal ( const G4ThreeVector p) const
private

Definition at line 594 of file G4Cons.cc.

595 {
596  ENorm side ;
598  G4double rho, phi ;
599  G4double distZ, distRMin, distRMax, distSPhi, distEPhi, distMin ;
600  G4double tanRMin, secRMin, pRMin, widRMin ;
601  G4double tanRMax, secRMax, pRMax, widRMax ;
602 
603  distZ = std::fabs(std::fabs(p.z()) - fDz) ;
604  rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
605 
606  tanRMin = (fRmin2 - fRmin1)*0.5/fDz ;
607  secRMin = std::sqrt(1 + tanRMin*tanRMin) ;
608  pRMin = rho - p.z()*tanRMin ;
609  widRMin = fRmin2 - fDz*tanRMin ;
610  distRMin = std::fabs(pRMin - widRMin)/secRMin ;
611 
612  tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;
613  secRMax = std::sqrt(1+tanRMax*tanRMax) ;
614  pRMax = rho - p.z()*tanRMax ;
615  widRMax = fRmax2 - fDz*tanRMax ;
616  distRMax = std::fabs(pRMax - widRMax)/secRMax ;
617 
618  if (distRMin < distRMax) // First minimum
619  {
620  if (distZ < distRMin)
621  {
622  distMin = distZ ;
623  side = kNZ ;
624  }
625  else
626  {
627  distMin = distRMin ;
628  side = kNRMin ;
629  }
630  }
631  else
632  {
633  if (distZ < distRMax)
634  {
635  distMin = distZ ;
636  side = kNZ ;
637  }
638  else
639  {
640  distMin = distRMax ;
641  side = kNRMax ;
642  }
643  }
644  if ( !fPhiFullCone && rho ) // Protected against (0,0,z)
645  {
646  phi = std::atan2(p.y(),p.x()) ;
647 
648  if (phi < 0) { phi += twopi; }
649 
650  if (fSPhi < 0) { distSPhi = std::fabs(phi - (fSPhi + twopi))*rho; }
651  else { distSPhi = std::fabs(phi - fSPhi)*rho; }
652 
653  distEPhi = std::fabs(phi - fSPhi - fDPhi)*rho ;
654 
655  // Find new minimum
656 
657  if (distSPhi < distEPhi)
658  {
659  if (distSPhi < distMin) { side = kNSPhi; }
660  }
661  else
662  {
663  if (distEPhi < distMin) { side = kNEPhi; }
664  }
665  }
666  switch (side)
667  {
668  case kNRMin: // Inner radius
669  rho *= secRMin ;
670  norm = G4ThreeVector(-p.x()/rho, -p.y()/rho, tanRMin/secRMin) ;
671  break ;
672  case kNRMax: // Outer radius
673  rho *= secRMax ;
674  norm = G4ThreeVector(p.x()/rho, p.y()/rho, -tanRMax/secRMax) ;
675  break ;
676  case kNZ: // +/- dz
677  if (p.z() > 0) { norm = G4ThreeVector(0,0,1); }
678  else { norm = G4ThreeVector(0,0,-1); }
679  break ;
680  case kNSPhi:
681  norm = G4ThreeVector(std::sin(fSPhi), -std::cos(fSPhi), 0) ;
682  break ;
683  case kNEPhi:
684  norm=G4ThreeVector(-std::sin(fSPhi+fDPhi), std::cos(fSPhi+fDPhi), 0) ;
685  break ;
686  default: // Should never reach this case...
687  DumpInfo();
688  G4Exception("G4Cons::ApproxSurfaceNormal()",
689  "GeomSolids1002", JustWarning,
690  "Undefined side for valid surface normal to solid.");
691  break ;
692  }
693  return norm ;
694 }
G4double fDz
Definition: G4Cons.hh:220
G4bool fPhiFullCone
Definition: G4Cons.hh:229
CLHEP::Hep3Vector G4ThreeVector
Float_t norm
G4double fRmin1
Definition: G4Cons.hh:220
G4double fRmax2
Definition: G4Cons.hh:220
void DumpInfo() const
ENorm
Definition: G4Cons.cc:75
G4double fRmin2
Definition: G4Cons.hh:220
G4double fSPhi
Definition: G4Cons.hh:220
static const double twopi
Definition: G4SIunits.hh:75
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
double x() const
G4double fDPhi
Definition: G4Cons.hh:220
double y() const
double z() const
double G4double
Definition: G4Types.hh:76
G4double fRmax1
Definition: G4Cons.hh:220
Here is the call graph for this function:
Here is the caller graph for this function:

◆ CalculateExtent()

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

Implements G4VSolid.

Definition at line 272 of file G4Cons.cc.

277 {
278  if ( !pTransform.IsRotated() && (fDPhi == twopi)
279  && (fRmin1 == 0) && (fRmin2 == 0) )
280  {
281  // Special case handling for unrotated solid cones
282  // Compute z/x/y mins and maxs for bounding box respecting limits,
283  // with early returns if outside limits. Then switch() on pAxis,
284  // and compute exact x and y limit for x/y case
285 
286  G4double xoffset, xMin, xMax ;
287  G4double yoffset, yMin, yMax ;
288  G4double zoffset, zMin, zMax ;
289 
290  G4double diff1, diff2, delta, maxDiff, newMin, newMax, RMax ;
291  G4double xoff1, xoff2, yoff1, yoff2 ;
292 
293  zoffset = pTransform.NetTranslation().z();
294  zMin = zoffset - fDz ;
295  zMax = zoffset + fDz ;
296 
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  xoffset = pTransform.NetTranslation().x() ;
317  RMax = (fRmax2 >= fRmax1) ? zMax : zMin ;
318  xMax = xoffset + (fRmax1 + fRmax2)*0.5 +
319  (RMax - zoffset)*(fRmax2 - fRmax1)/(2*fDz) ;
320  xMin = 2*xoffset-xMax ;
321 
322  if (pVoxelLimit.IsXLimited())
323  {
324  if ( (xMin > pVoxelLimit.GetMaxXExtent() + kCarTolerance) ||
325  (xMax < pVoxelLimit.GetMinXExtent() - kCarTolerance) )
326  {
327  return false ;
328  }
329  else
330  {
331  if ( xMin < pVoxelLimit.GetMinXExtent() )
332  {
333  xMin = pVoxelLimit.GetMinXExtent() ;
334  }
335  if ( xMax > pVoxelLimit.GetMaxXExtent() )
336  {
337  xMax=pVoxelLimit.GetMaxXExtent() ;
338  }
339  }
340  }
341  yoffset = pTransform.NetTranslation().y() ;
342  yMax = yoffset + (fRmax1 + fRmax2)*0.5 +
343  (RMax - zoffset)*(fRmax2 - fRmax1)/(2*fDz) ;
344  yMin = 2*yoffset-yMax ;
345  RMax = yMax - yoffset ; // = max radius due to Zmax/Zmin cuttings
346 
347  if (pVoxelLimit.IsYLimited())
348  {
349  if ( (yMin > pVoxelLimit.GetMaxYExtent() + kCarTolerance) ||
350  (yMax < pVoxelLimit.GetMinYExtent() - kCarTolerance) )
351  {
352  return false ;
353  }
354  else
355  {
356  if ( yMin < pVoxelLimit.GetMinYExtent() )
357  {
358  yMin = pVoxelLimit.GetMinYExtent() ;
359  }
360  if ( yMax > pVoxelLimit.GetMaxYExtent() )
361  {
362  yMax = pVoxelLimit.GetMaxYExtent() ;
363  }
364  }
365  }
366  switch (pAxis) // Known to cut cones
367  {
368  case kXAxis:
369  yoff1 = yoffset - yMin ;
370  yoff2 = yMax - yoffset ;
371 
372  if ((yoff1 >= 0) && (yoff2 >= 0)) // Y limits cross max/min x
373  { // => no change
374  pMin = xMin ;
375  pMax = xMax ;
376  }
377  else
378  {
379  // Y limits don't cross max/min x => compute max delta x,
380  // hence new mins/maxs
381  delta=RMax*RMax-yoff1*yoff1;
382  diff1=(delta>0.) ? std::sqrt(delta) : 0.;
383  delta=RMax*RMax-yoff2*yoff2;
384  diff2=(delta>0.) ? std::sqrt(delta) : 0.;
385  maxDiff = (diff1>diff2) ? diff1:diff2 ;
386  newMin = xoffset - maxDiff ;
387  newMax = xoffset + maxDiff ;
388  pMin = ( newMin < xMin ) ? xMin : newMin ;
389  pMax = ( newMax > xMax) ? xMax : newMax ;
390  }
391  break ;
392 
393  case kYAxis:
394  xoff1 = xoffset - xMin ;
395  xoff2 = xMax - xoffset ;
396 
397  if ((xoff1 >= 0) && (xoff2 >= 0) ) // X limits cross max/min y
398  { // => no change
399  pMin = yMin ;
400  pMax = yMax ;
401  }
402  else
403  {
404  // X limits don't cross max/min y => compute max delta y,
405  // hence new mins/maxs
406  delta=RMax*RMax-xoff1*xoff1;
407  diff1=(delta>0.) ? std::sqrt(delta) : 0.;
408  delta=RMax*RMax-xoff2*xoff2;
409  diff2=(delta>0.) ? std::sqrt(delta) : 0.;
410  maxDiff = (diff1 > diff2) ? diff1:diff2 ;
411  newMin = yoffset - maxDiff ;
412  newMax = yoffset + maxDiff ;
413  pMin = (newMin < yMin) ? yMin : newMin ;
414  pMax = (newMax > yMax) ? yMax : newMax ;
415  }
416  break ;
417 
418  case kZAxis:
419  pMin = zMin ;
420  pMax = zMax ;
421  break ;
422 
423  default:
424  break ;
425  }
426  pMin -= kCarTolerance ;
427  pMax += kCarTolerance ;
428 
429  return true ;
430  }
431  else // Calculate rotated vertex coordinates
432  {
433  G4int i, noEntries, noBetweenSections4 ;
434  G4bool existsAfterClip = false ;
435  G4ThreeVectorList* vertices = CreateRotatedVertices(pTransform) ;
436 
437  pMin = +kInfinity ;
438  pMax = -kInfinity ;
439 
440  noEntries = vertices->size() ;
441  noBetweenSections4 = noEntries-4 ;
442 
443  for ( i = 0 ; i < noEntries ; i += 4 )
444  {
445  ClipCrossSection(vertices, i, pVoxelLimit, pAxis, pMin, pMax) ;
446  }
447  for ( i = 0 ; i < noBetweenSections4 ; i += 4 )
448  {
449  ClipBetweenSections(vertices, i, pVoxelLimit, pAxis, pMin, pMax) ;
450  }
451  if ( (pMin != kInfinity) || (pMax != -kInfinity) )
452  {
453  existsAfterClip = true ;
454 
455  // Add 2*tolerance to avoid precision troubles
456 
457  pMin -= kCarTolerance ;
458  pMax += kCarTolerance ;
459  }
460  else
461  {
462  // Check for case where completely enveloping clipping volume
463  // If point inside then we are confident that the solid completely
464  // envelopes the clipping volume. Hence set min/max extents according
465  // to clipping volume extents along the specified axis.
466 
467  G4ThreeVector clipCentre(
468  (pVoxelLimit.GetMinXExtent() + pVoxelLimit.GetMaxXExtent())*0.5,
469  (pVoxelLimit.GetMinYExtent() + pVoxelLimit.GetMaxYExtent())*0.5,
470  (pVoxelLimit.GetMinZExtent() + pVoxelLimit.GetMaxZExtent())*0.5 ) ;
471 
472  if (Inside(pTransform.Inverse().TransformPoint(clipCentre)) != kOutside)
473  {
474  existsAfterClip = true ;
475  pMin = pVoxelLimit.GetMinExtent(pAxis) ;
476  pMax = pVoxelLimit.GetMaxExtent(pAxis) ;
477  }
478  }
479  delete vertices ;
480  return existsAfterClip ;
481  }
482 }
void ClipBetweenSections(G4ThreeVectorList *pVertices, const G4int pSectionIndex, const G4VoxelLimits &pVoxelLimit, const EAxis pAxis, G4double &pMin, G4double &pMax) const
Definition: G4VSolid.cc:378
G4double GetMinXExtent() const
G4double fDz
Definition: G4Cons.hh:220
static const G4double kInfinity
Definition: geomdefs.hh:42
G4double GetMinYExtent() const
G4ThreeVectorList * CreateRotatedVertices(const G4AffineTransform &pTransform) const
Definition: G4Cons.cc:2161
G4double GetMinZExtent() const
G4bool IsYLimited() const
G4ThreeVector TransformPoint(const G4ThreeVector &vec) const
G4double fRmin1
Definition: G4Cons.hh:220
void ClipCrossSection(G4ThreeVectorList *pVertices, const G4int pSectionIndex, const G4VoxelLimits &pVoxelLimit, const EAxis pAxis, G4double &pMin, G4double &pMax) const
Definition: G4VSolid.cc:347
G4bool IsXLimited() const
G4bool IsRotated() const
int G4int
Definition: G4Types.hh:78
G4double fRmax2
Definition: G4Cons.hh:220
G4double GetMaxZExtent() const
G4double GetMaxXExtent() const
G4AffineTransform Inverse() const
G4ThreeVector NetTranslation() const
G4double fRmin2
Definition: G4Cons.hh:220
bool G4bool
Definition: G4Types.hh:79
static const double twopi
Definition: G4SIunits.hh:75
std::vector< G4ThreeVector > G4ThreeVectorList
Definition: G4VSolid.hh:79
G4double GetMinExtent(const EAxis pAxis) const
EInside Inside(const G4ThreeVector &p) const
Definition: G4Cons.cc:205
double x() const
G4bool IsZLimited() const
G4double fDPhi
Definition: G4Cons.hh:220
double y() const
double z() const
G4double GetMaxExtent(const EAxis pAxis) const
G4double kCarTolerance
Definition: G4VSolid.hh:304
double G4double
Definition: G4Types.hh:76
G4double GetMaxYExtent() const
G4double fRmax1
Definition: G4Cons.hh:220
Here is the call graph for this function:

◆ CheckDPhiAngle()

void G4Cons::CheckDPhiAngle ( G4double  dPhi)
inlineprivate

◆ CheckPhiAngles()

void G4Cons::CheckPhiAngles ( G4double  sPhi,
G4double  dPhi 
)
inlineprivate
Here is the caller graph for this function:

◆ CheckSPhiAngle()

void G4Cons::CheckSPhiAngle ( G4double  sPhi)
inlineprivate

◆ Clone()

G4VSolid * G4Cons::Clone ( ) const
virtual

Reimplemented from G4VSolid.

Definition at line 2256 of file G4Cons.cc.

2257 {
2258  return new G4Cons(*this);
2259 }
G4Cons(const G4String &pName, G4double pRmin1, G4double pRmax1, G4double pRmin2, G4double pRmax2, G4double pDz, G4double pSPhi, G4double pDPhi)
Definition: G4Cons.cc:82
Here is the call graph for this function:

◆ ComputeDimensions()

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

Reimplemented from G4VSolid.

Definition at line 260 of file G4Cons.cc.

263 {
264  p->ComputeDimensions(*this,n,pRep) ;
265 }
Char_t n[5]
virtual void ComputeDimensions(G4Box &, const G4int, const G4VPhysicalVolume *) const
Here is the call graph for this function:

◆ CreatePolyhedron()

G4Polyhedron * G4Cons::CreatePolyhedron ( ) const
virtual

Reimplemented from G4VSolid.

Definition at line 2384 of file G4Cons.cc.

2385 {
2387 }
G4double fDz
Definition: G4Cons.hh:220
G4double fRmin1
Definition: G4Cons.hh:220
G4double fRmax2
Definition: G4Cons.hh:220
G4double fRmin2
Definition: G4Cons.hh:220
G4double fSPhi
Definition: G4Cons.hh:220
G4double fDPhi
Definition: G4Cons.hh:220
G4double fRmax1
Definition: G4Cons.hh:220

◆ CreateRotatedVertices()

G4ThreeVectorList * G4Cons::CreateRotatedVertices ( const G4AffineTransform pTransform) const
private

Definition at line 2161 of file G4Cons.cc.

2162 {
2163  G4ThreeVectorList* vertices ;
2164  G4ThreeVector vertex0, vertex1, vertex2, vertex3 ;
2165  G4double meshAngle, meshRMax1, meshRMax2, crossAngle;
2166  G4double cosCrossAngle, sinCrossAngle, sAngle ;
2167  G4double rMaxX1, rMaxX2, rMaxY1, rMaxY2, rMinX1, rMinX2, rMinY1, rMinY2 ;
2168  G4int crossSection, noCrossSections ;
2169 
2170  // Compute no of cross-sections necessary to mesh cone
2171 
2172  noCrossSections = G4int(fDPhi/kMeshAngleDefault) + 1 ;
2173 
2174  if (noCrossSections < kMinMeshSections)
2175  {
2176  noCrossSections = kMinMeshSections ;
2177  }
2178  else if (noCrossSections > kMaxMeshSections)
2179  {
2180  noCrossSections = kMaxMeshSections ;
2181  }
2182  meshAngle = fDPhi/(noCrossSections - 1) ;
2183 
2184  meshRMax1 = fRmax1/std::cos(meshAngle*0.5) ;
2185  meshRMax2 = fRmax2/std::cos(meshAngle*0.5) ;
2186 
2187  // If complete in phi, set start angle such that mesh will be at RMax
2188  // on the x axis. Will give better extent calculations when not rotated.
2189 
2190  if ( fPhiFullCone && (fSPhi == 0.0) )
2191  {
2192  sAngle = -meshAngle*0.5 ;
2193  }
2194  else
2195  {
2196  sAngle = fSPhi ;
2197  }
2198  vertices = new G4ThreeVectorList();
2199 
2200  if (vertices)
2201  {
2202  vertices->reserve(noCrossSections*4) ;
2203  for (crossSection = 0 ; crossSection < noCrossSections ; crossSection++)
2204  {
2205  // Compute coordinates of cross section at section crossSection
2206 
2207  crossAngle = sAngle + crossSection*meshAngle ;
2208  cosCrossAngle = std::cos(crossAngle) ;
2209  sinCrossAngle = std::sin(crossAngle) ;
2210 
2211  rMaxX1 = meshRMax1*cosCrossAngle ;
2212  rMaxY1 = meshRMax1*sinCrossAngle ;
2213  rMaxX2 = meshRMax2*cosCrossAngle ;
2214  rMaxY2 = meshRMax2*sinCrossAngle ;
2215 
2216  rMinX1 = fRmin1*cosCrossAngle ;
2217  rMinY1 = fRmin1*sinCrossAngle ;
2218  rMinX2 = fRmin2*cosCrossAngle ;
2219  rMinY2 = fRmin2*sinCrossAngle ;
2220 
2221  vertex0 = G4ThreeVector(rMinX1,rMinY1,-fDz) ;
2222  vertex1 = G4ThreeVector(rMaxX1,rMaxY1,-fDz) ;
2223  vertex2 = G4ThreeVector(rMaxX2,rMaxY2,+fDz) ;
2224  vertex3 = G4ThreeVector(rMinX2,rMinY2,+fDz) ;
2225 
2226  vertices->push_back(pTransform.TransformPoint(vertex0)) ;
2227  vertices->push_back(pTransform.TransformPoint(vertex1)) ;
2228  vertices->push_back(pTransform.TransformPoint(vertex2)) ;
2229  vertices->push_back(pTransform.TransformPoint(vertex3)) ;
2230  }
2231  }
2232  else
2233  {
2234  DumpInfo();
2235  G4Exception("G4Cons::CreateRotatedVertices()",
2236  "GeomSolids0003", FatalException,
2237  "Error in allocation of vertices. Out of memory !");
2238  }
2239 
2240  return vertices ;
2241 }
G4double fDz
Definition: G4Cons.hh:220
G4bool fPhiFullCone
Definition: G4Cons.hh:229
CLHEP::Hep3Vector G4ThreeVector
G4ThreeVector TransformPoint(const G4ThreeVector &vec) const
G4double fRmin1
Definition: G4Cons.hh:220
const G4int kMinMeshSections
Definition: meshdefs.hh:45
int G4int
Definition: G4Types.hh:78
G4double fRmax2
Definition: G4Cons.hh:220
void DumpInfo() const
G4double fRmin2
Definition: G4Cons.hh:220
G4double fSPhi
Definition: G4Cons.hh:220
std::vector< G4ThreeVector > G4ThreeVectorList
Definition: G4VSolid.hh:79
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4double fDPhi
Definition: G4Cons.hh:220
double G4double
Definition: G4Types.hh:76
const G4double kMeshAngleDefault
Definition: meshdefs.hh:42
G4double fRmax1
Definition: G4Cons.hh:220
const G4int kMaxMeshSections
Definition: meshdefs.hh:46
Here is the call graph for this function:
Here is the caller graph for this function:

◆ DescribeYourselfTo()

void G4Cons::DescribeYourselfTo ( G4VGraphicsScene scene) const
virtual

Implements G4VSolid.

Definition at line 2379 of file G4Cons.cc.

2380 {
2381  scene.AddSolid (*this);
2382 }
virtual void AddSolid(const G4Box &)=0
Here is the call graph for this function:

◆ DistanceToIn() [1/2]

G4double G4Cons::DistanceToIn ( const G4ThreeVector p,
const G4ThreeVector v 
) const
virtual

Implements G4VSolid.

Definition at line 720 of file G4Cons.cc.

722 {
723  G4double snxt = kInfinity ; // snxt = default return value
724  const G4double dRmax = 50*(fRmax1+fRmax2);// 100*(Rmax1+Rmax2)/2.
725 
726  G4double tanRMax,secRMax,rMaxAv,rMaxOAv ; // Data for cones
727  G4double tanRMin,secRMin,rMinAv,rMinOAv ;
728  G4double rout,rin ;
729 
730  G4double tolORMin,tolORMin2,tolIRMin,tolIRMin2 ; // `generous' radii squared
731  G4double tolORMax2,tolIRMax,tolIRMax2 ;
732  G4double tolODz,tolIDz ;
733 
734  G4double Dist,sd,xi,yi,zi,ri=0.,risec,rhoi2,cosPsi ; // Intersection point vars
735 
736  G4double t1,t2,t3,b,c,d ; // Quadratic solver variables
737  G4double nt1,nt2,nt3 ;
738  G4double Comp ;
739 
740  G4ThreeVector Normal;
741 
742  // Cone Precalcs
743 
744  tanRMin = (fRmin2 - fRmin1)*0.5/fDz ;
745  secRMin = std::sqrt(1.0 + tanRMin*tanRMin) ;
746  rMinAv = (fRmin1 + fRmin2)*0.5 ;
747 
748  if (rMinAv > halfRadTolerance)
749  {
750  rMinOAv = rMinAv - halfRadTolerance ;
751  }
752  else
753  {
754  rMinOAv = 0.0 ;
755  }
756  tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;
757  secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
758  rMaxAv = (fRmax1 + fRmax2)*0.5 ;
759  rMaxOAv = rMaxAv + halfRadTolerance ;
760 
761  // Intersection with z-surfaces
762 
763  tolIDz = fDz - halfCarTolerance ;
764  tolODz = fDz + halfCarTolerance ;
765 
766  if (std::fabs(p.z()) >= tolIDz)
767  {
768  if ( p.z()*v.z() < 0 ) // at +Z going in -Z or visa versa
769  {
770  sd = (std::fabs(p.z()) - fDz)/std::fabs(v.z()) ; // Z intersect distance
771 
772  if( sd < 0.0 ) { sd = 0.0; } // negative dist -> zero
773 
774  xi = p.x() + sd*v.x() ; // Intersection coords
775  yi = p.y() + sd*v.y() ;
776  rhoi2 = xi*xi + yi*yi ;
777 
778  // Check validity of intersection
779  // Calculate (outer) tolerant radi^2 at intersecion
780 
781  if (v.z() > 0)
782  {
783  tolORMin = fRmin1 - halfRadTolerance*secRMin ;
784  tolIRMin = fRmin1 + halfRadTolerance*secRMin ;
785  tolIRMax = fRmax1 - halfRadTolerance*secRMin ;
786  // tolORMax2 = (fRmax1 + halfRadTolerance*secRMax)*
787  // (fRmax1 + halfRadTolerance*secRMax) ;
788  }
789  else
790  {
791  tolORMin = fRmin2 - halfRadTolerance*secRMin ;
792  tolIRMin = fRmin2 + halfRadTolerance*secRMin ;
793  tolIRMax = fRmax2 - halfRadTolerance*secRMin ;
794  // tolORMax2 = (fRmax2 + halfRadTolerance*secRMax)*
795  // (fRmax2 + halfRadTolerance*secRMax) ;
796  }
797  if ( tolORMin > 0 )
798  {
799  // tolORMin2 = tolORMin*tolORMin ;
800  tolIRMin2 = tolIRMin*tolIRMin ;
801  }
802  else
803  {
804  // tolORMin2 = 0.0 ;
805  tolIRMin2 = 0.0 ;
806  }
807  if ( tolIRMax > 0 ) { tolIRMax2 = tolIRMax*tolIRMax; }
808  else { tolIRMax2 = 0.0; }
809 
810  if ( (tolIRMin2 <= rhoi2) && (rhoi2 <= tolIRMax2) )
811  {
812  if ( !fPhiFullCone && rhoi2 )
813  {
814  // Psi = angle made with central (average) phi of shape
815 
816  cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
817 
818  if (cosPsi >= cosHDPhiIT) { return sd; }
819  }
820  else
821  {
822  return sd;
823  }
824  }
825  }
826  else // On/outside extent, and heading away -> cannot intersect
827  {
828  return snxt ;
829  }
830  }
831 
832 // ----> Can not intersect z surfaces
833 
834 
835 // Intersection with outer cone (possible return) and
836 // inner cone (must also check phi)
837 //
838 // Intersection point (xi,yi,zi) on line x=p.x+t*v.x etc.
839 //
840 // Intersects with x^2+y^2=(a*z+b)^2
841 //
842 // where a=tanRMax or tanRMin
843 // b=rMaxAv or rMinAv
844 //
845 // (vx^2+vy^2-(a*vz)^2)t^2+2t(pxvx+pyvy-a*vz(a*pz+b))+px^2+py^2-(a*pz+b)^2=0 ;
846 // t1 t2 t3
847 //
848 // \--------u-------/ \-----------v----------/ \---------w--------/
849 //
850 
851  t1 = 1.0 - v.z()*v.z() ;
852  t2 = p.x()*v.x() + p.y()*v.y() ;
853  t3 = p.x()*p.x() + p.y()*p.y() ;
854  rin = tanRMin*p.z() + rMinAv ;
855  rout = tanRMax*p.z() + rMaxAv ;
856 
857  // Outer Cone Intersection
858  // Must be outside/on outer cone for valid intersection
859 
860  nt1 = t1 - (tanRMax*v.z())*(tanRMax*v.z()) ;
861  nt2 = t2 - tanRMax*v.z()*rout ;
862  nt3 = t3 - rout*rout ;
863 
864  if (std::fabs(nt1) > kRadTolerance) // Equation quadratic => 2 roots
865  {
866  b = nt2/nt1;
867  c = nt3/nt1;
868  d = b*b-c ;
869  if ( (nt3 > rout*rout*kRadTolerance*kRadTolerance*secRMax*secRMax)
870  || (rout < 0) )
871  {
872  // If outside real cone (should be rho-rout>kRadTolerance*0.5
873  // NOT rho^2 etc) saves a std::sqrt() at expense of accuracy
874 
875  if (d >= 0)
876  {
877 
878  if ((rout < 0) && (nt3 <= 0))
879  {
880  // Inside `shadow cone' with -ve radius
881  // -> 2nd root could be on real cone
882 
883  if (b>0) { sd = c/(-b-std::sqrt(d)); }
884  else { sd = -b + std::sqrt(d); }
885  }
886  else
887  {
888  if ((b <= 0) && (c >= 0)) // both >=0, try smaller root
889  {
890  sd=c/(-b+std::sqrt(d));
891  }
892  else
893  {
894  if ( c <= 0 ) // second >=0
895  {
896  sd = -b + std::sqrt(d) ;
897  if((sd<0) & (sd>-halfRadTolerance)) sd=0;
898  }
899  else // both negative, travel away
900  {
901  return kInfinity ;
902  }
903  }
904  }
905  if ( sd >= 0 ) // If 'forwards'. Check z intersection
906  {
907  if ( sd>dRmax ) // Avoid rounding errors due to precision issues on
908  { // 64 bits systems. Split long distances and recompute
909  G4double fTerm = sd-std::fmod(sd,dRmax);
910  sd = fTerm + DistanceToIn(p+fTerm*v,v);
911  }
912  zi = p.z() + sd*v.z() ;
913 
914  if (std::fabs(zi) <= tolODz)
915  {
916  // Z ok. Check phi intersection if reqd
917 
918  if ( fPhiFullCone ) { return sd; }
919  else
920  {
921  xi = p.x() + sd*v.x() ;
922  yi = p.y() + sd*v.y() ;
923  ri = rMaxAv + zi*tanRMax ;
924  cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
925 
926  if ( cosPsi >= cosHDPhiIT ) { return sd; }
927  }
928  }
929  } // end if (sd>0)
930  }
931  }
932  else
933  {
934  // Inside outer cone
935  // check not inside, and heading through G4Cons (-> 0 to in)
936 
937  if ( ( t3 > (rin + halfRadTolerance*secRMin)*
938  (rin + halfRadTolerance*secRMin) )
939  && (nt2 < 0) && (d >= 0) && (std::fabs(p.z()) <= tolIDz) )
940  {
941  // Inside cones, delta r -ve, inside z extent
942  // Point is on the Surface => check Direction using Normal.dot(v)
943 
944  xi = p.x() ;
945  yi = p.y() ;
946  risec = std::sqrt(xi*xi + yi*yi)*secRMax ;
947  Normal = G4ThreeVector(xi/risec,yi/risec,-tanRMax/secRMax) ;
948  if ( !fPhiFullCone )
949  {
950  cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/std::sqrt(t3) ;
951  if ( cosPsi >= cosHDPhiIT )
952  {
953  if ( Normal.dot(v) <= 0 ) { return 0.0; }
954  }
955  }
956  else
957  {
958  if ( Normal.dot(v) <= 0 ) { return 0.0; }
959  }
960  }
961  }
962  }
963  else // Single root case
964  {
965  if ( std::fabs(nt2) > kRadTolerance )
966  {
967  sd = -0.5*nt3/nt2 ;
968 
969  if ( sd < 0 ) { return kInfinity; } // travel away
970  else // sd >= 0, If 'forwards'. Check z intersection
971  {
972  zi = p.z() + sd*v.z() ;
973 
974  if ((std::fabs(zi) <= tolODz) && (nt2 < 0))
975  {
976  // Z ok. Check phi intersection if reqd
977 
978  if ( fPhiFullCone ) { return sd; }
979  else
980  {
981  xi = p.x() + sd*v.x() ;
982  yi = p.y() + sd*v.y() ;
983  ri = rMaxAv + zi*tanRMax ;
984  cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
985 
986  if (cosPsi >= cosHDPhiIT) { return sd; }
987  }
988  }
989  }
990  }
991  else // travel || cone surface from its origin
992  {
993  sd = kInfinity ;
994  }
995  }
996 
997  // Inner Cone Intersection
998  // o Space is divided into 3 areas:
999  // 1) Radius greater than real inner cone & imaginary cone & outside
1000  // tolerance
1001  // 2) Radius less than inner or imaginary cone & outside tolarance
1002  // 3) Within tolerance of real or imaginary cones
1003  // - Extra checks needed for 3's intersections
1004  // => lots of duplicated code
1005 
1006  if (rMinAv)
1007  {
1008  nt1 = t1 - (tanRMin*v.z())*(tanRMin*v.z()) ;
1009  nt2 = t2 - tanRMin*v.z()*rin ;
1010  nt3 = t3 - rin*rin ;
1011 
1012  if ( nt1 )
1013  {
1014  if ( nt3 > rin*kRadTolerance*secRMin )
1015  {
1016  // At radius greater than real & imaginary cones
1017  // -> 2nd root, with zi check
1018 
1019  b = nt2/nt1 ;
1020  c = nt3/nt1 ;
1021  d = b*b-c ;
1022  if (d >= 0) // > 0
1023  {
1024  if(b>0){sd = c/( -b-std::sqrt(d));}
1025  else {sd = -b + std::sqrt(d) ;}
1026 
1027  if ( sd >= 0 ) // > 0
1028  {
1029  if ( sd>dRmax ) // Avoid rounding errors due to precision issues on
1030  { // 64 bits systems. Split long distance and recompute
1031  G4double fTerm = sd-std::fmod(sd,dRmax);
1032  sd = fTerm + DistanceToIn(p+fTerm*v,v);
1033  }
1034  zi = p.z() + sd*v.z() ;
1035 
1036  if ( std::fabs(zi) <= tolODz )
1037  {
1038  if ( !fPhiFullCone )
1039  {
1040  xi = p.x() + sd*v.x() ;
1041  yi = p.y() + sd*v.y() ;
1042  ri = rMinAv + zi*tanRMin ;
1043  cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
1044 
1045  if (cosPsi >= cosHDPhiIT)
1046  {
1047  if ( sd > halfRadTolerance ) { snxt=sd; }
1048  else
1049  {
1050  // Calculate a normal vector in order to check Direction
1051 
1052  risec = std::sqrt(xi*xi + yi*yi)*secRMin ;
1053  Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin);
1054  if ( Normal.dot(v) <= 0 ) { snxt = sd; }
1055  }
1056  }
1057  }
1058  else
1059  {
1060  if ( sd > halfRadTolerance ) { return sd; }
1061  else
1062  {
1063  // Calculate a normal vector in order to check Direction
1064 
1065  xi = p.x() + sd*v.x() ;
1066  yi = p.y() + sd*v.y() ;
1067  risec = std::sqrt(xi*xi + yi*yi)*secRMin ;
1068  Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin) ;
1069  if ( Normal.dot(v) <= 0 ) { return sd; }
1070  }
1071  }
1072  }
1073  }
1074  }
1075  }
1076  else if ( nt3 < -rin*kRadTolerance*secRMin )
1077  {
1078  // Within radius of inner cone (real or imaginary)
1079  // -> Try 2nd root, with checking intersection is with real cone
1080  // -> If check fails, try 1st root, also checking intersection is
1081  // on real cone
1082 
1083  b = nt2/nt1 ;
1084  c = nt3/nt1 ;
1085  d = b*b - c ;
1086 
1087  if ( d >= 0 ) // > 0
1088  {
1089  if (b>0) { sd = c/(-b-std::sqrt(d)); }
1090  else { sd = -b + std::sqrt(d); }
1091  zi = p.z() + sd*v.z() ;
1092  ri = rMinAv + zi*tanRMin ;
1093 
1094  if ( ri > 0 )
1095  {
1096  if ( (sd >= 0) && (std::fabs(zi) <= tolODz) ) // sd > 0
1097  {
1098  if ( sd>dRmax ) // Avoid rounding errors due to precision issues
1099  { // seen on 64 bits systems. Split and recompute
1100  G4double fTerm = sd-std::fmod(sd,dRmax);
1101  sd = fTerm + DistanceToIn(p+fTerm*v,v);
1102  }
1103  if ( !fPhiFullCone )
1104  {
1105  xi = p.x() + sd*v.x() ;
1106  yi = p.y() + sd*v.y() ;
1107  cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
1108 
1109  if (cosPsi >= cosHDPhiOT)
1110  {
1111  if ( sd > halfRadTolerance ) { snxt=sd; }
1112  else
1113  {
1114  // Calculate a normal vector in order to check Direction
1115 
1116  risec = std::sqrt(xi*xi + yi*yi)*secRMin ;
1117  Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin);
1118  if ( Normal.dot(v) <= 0 ) { snxt = sd; }
1119  }
1120  }
1121  }
1122  else
1123  {
1124  if( sd > halfRadTolerance ) { return sd; }
1125  else
1126  {
1127  // Calculate a normal vector in order to check Direction
1128 
1129  xi = p.x() + sd*v.x() ;
1130  yi = p.y() + sd*v.y() ;
1131  risec = std::sqrt(xi*xi + yi*yi)*secRMin ;
1132  Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin) ;
1133  if ( Normal.dot(v) <= 0 ) { return sd; }
1134  }
1135  }
1136  }
1137  }
1138  else
1139  {
1140  if (b>0) { sd = -b - std::sqrt(d); }
1141  else { sd = c/(-b+std::sqrt(d)); }
1142  zi = p.z() + sd*v.z() ;
1143  ri = rMinAv + zi*tanRMin ;
1144 
1145  if ( (sd >= 0) && (ri > 0) && (std::fabs(zi) <= tolODz) ) // sd>0
1146  {
1147  if ( sd>dRmax ) // Avoid rounding errors due to precision issues
1148  { // seen on 64 bits systems. Split and recompute
1149  G4double fTerm = sd-std::fmod(sd,dRmax);
1150  sd = fTerm + DistanceToIn(p+fTerm*v,v);
1151  }
1152  if ( !fPhiFullCone )
1153  {
1154  xi = p.x() + sd*v.x() ;
1155  yi = p.y() + sd*v.y() ;
1156  cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
1157 
1158  if (cosPsi >= cosHDPhiIT)
1159  {
1160  if ( sd > halfRadTolerance ) { snxt=sd; }
1161  else
1162  {
1163  // Calculate a normal vector in order to check Direction
1164 
1165  risec = std::sqrt(xi*xi + yi*yi)*secRMin ;
1166  Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin);
1167  if ( Normal.dot(v) <= 0 ) { snxt = sd; }
1168  }
1169  }
1170  }
1171  else
1172  {
1173  if ( sd > halfRadTolerance ) { return sd; }
1174  else
1175  {
1176  // Calculate a normal vector in order to check Direction
1177 
1178  xi = p.x() + sd*v.x() ;
1179  yi = p.y() + sd*v.y() ;
1180  risec = std::sqrt(xi*xi + yi*yi)*secRMin ;
1181  Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin) ;
1182  if ( Normal.dot(v) <= 0 ) { return sd; }
1183  }
1184  }
1185  }
1186  }
1187  }
1188  }
1189  else
1190  {
1191  // Within kRadTol*0.5 of inner cone (real OR imaginary)
1192  // ----> Check not travelling through (=>0 to in)
1193  // ----> if not:
1194  // -2nd root with validity check
1195 
1196  if ( std::fabs(p.z()) <= tolODz )
1197  {
1198  if ( nt2 > 0 )
1199  {
1200  // Inside inner real cone, heading outwards, inside z range
1201 
1202  if ( !fPhiFullCone )
1203  {
1204  cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/std::sqrt(t3) ;
1205 
1206  if (cosPsi >= cosHDPhiIT) { return 0.0; }
1207  }
1208  else { return 0.0; }
1209  }
1210  else
1211  {
1212  // Within z extent, but not travelling through
1213  // -> 2nd root or kInfinity if 1st root on imaginary cone
1214 
1215  b = nt2/nt1 ;
1216  c = nt3/nt1 ;
1217  d = b*b - c ;
1218 
1219  if ( d >= 0 ) // > 0
1220  {
1221  if (b>0) { sd = -b - std::sqrt(d); }
1222  else { sd = c/(-b+std::sqrt(d)); }
1223  zi = p.z() + sd*v.z() ;
1224  ri = rMinAv + zi*tanRMin ;
1225 
1226  if ( ri > 0 ) // 2nd root
1227  {
1228  if (b>0) { sd = c/(-b-std::sqrt(d)); }
1229  else { sd = -b + std::sqrt(d); }
1230 
1231  zi = p.z() + sd*v.z() ;
1232 
1233  if ( (sd >= 0) && (std::fabs(zi) <= tolODz) ) // sd>0
1234  {
1235  if ( sd>dRmax ) // Avoid rounding errors due to precision issue
1236  { // seen on 64 bits systems. Split and recompute
1237  G4double fTerm = sd-std::fmod(sd,dRmax);
1238  sd = fTerm + DistanceToIn(p+fTerm*v,v);
1239  }
1240  if ( !fPhiFullCone )
1241  {
1242  xi = p.x() + sd*v.x() ;
1243  yi = p.y() + sd*v.y() ;
1244  ri = rMinAv + zi*tanRMin ;
1245  cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
1246 
1247  if ( cosPsi >= cosHDPhiIT ) { snxt = sd; }
1248  }
1249  else { return sd; }
1250  }
1251  }
1252  else { return kInfinity; }
1253  }
1254  }
1255  }
1256  else // 2nd root
1257  {
1258  b = nt2/nt1 ;
1259  c = nt3/nt1 ;
1260  d = b*b - c ;
1261 
1262  if ( d > 0 )
1263  {
1264  if (b>0) { sd = c/(-b-std::sqrt(d)); }
1265  else { sd = -b + std::sqrt(d) ; }
1266  zi = p.z() + sd*v.z() ;
1267 
1268  if ( (sd >= 0) && (std::fabs(zi) <= tolODz) ) // sd>0
1269  {
1270  if ( sd>dRmax ) // Avoid rounding errors due to precision issues
1271  { // seen on 64 bits systems. Split and recompute
1272  G4double fTerm = sd-std::fmod(sd,dRmax);
1273  sd = fTerm + DistanceToIn(p+fTerm*v,v);
1274  }
1275  if ( !fPhiFullCone )
1276  {
1277  xi = p.x() + sd*v.x();
1278  yi = p.y() + sd*v.y();
1279  ri = rMinAv + zi*tanRMin ;
1280  cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri;
1281 
1282  if (cosPsi >= cosHDPhiIT) { snxt = sd; }
1283  }
1284  else { return sd; }
1285  }
1286  }
1287  }
1288  }
1289  }
1290  }
1291 
1292  // Phi segment intersection
1293  //
1294  // o Tolerant of points inside phi planes by up to kCarTolerance*0.5
1295  //
1296  // o NOTE: Large duplication of code between sphi & ephi checks
1297  // -> only diffs: sphi -> ephi, Comp -> -Comp and half-plane
1298  // intersection check <=0 -> >=0
1299  // -> Should use some form of loop Construct
1300 
1301  if ( !fPhiFullCone )
1302  {
1303  // First phi surface (starting phi)
1304 
1305  Comp = v.x()*sinSPhi - v.y()*cosSPhi ;
1306 
1307  if ( Comp < 0 ) // Component in outwards normal dirn
1308  {
1309  Dist = (p.y()*cosSPhi - p.x()*sinSPhi) ;
1310 
1311  if (Dist < halfCarTolerance)
1312  {
1313  sd = Dist/Comp ;
1314 
1315  if ( sd < snxt )
1316  {
1317  if ( sd < 0 ) { sd = 0.0; }
1318 
1319  zi = p.z() + sd*v.z() ;
1320 
1321  if ( std::fabs(zi) <= tolODz )
1322  {
1323  xi = p.x() + sd*v.x() ;
1324  yi = p.y() + sd*v.y() ;
1325  rhoi2 = xi*xi + yi*yi ;
1326  tolORMin2 = (rMinOAv + zi*tanRMin)*(rMinOAv + zi*tanRMin) ;
1327  tolORMax2 = (rMaxOAv + zi*tanRMax)*(rMaxOAv + zi*tanRMax) ;
1328 
1329  if ( (rhoi2 >= tolORMin2) && (rhoi2 <= tolORMax2) )
1330  {
1331  // z and r intersections good - check intersecting with
1332  // correct half-plane
1333 
1334  if ((yi*cosCPhi - xi*sinCPhi) <= 0 ) { snxt = sd; }
1335  }
1336  }
1337  }
1338  }
1339  }
1340 
1341  // Second phi surface (Ending phi)
1342 
1343  Comp = -(v.x()*sinEPhi - v.y()*cosEPhi) ;
1344 
1345  if ( Comp < 0 ) // Component in outwards normal dirn
1346  {
1347  Dist = -(p.y()*cosEPhi - p.x()*sinEPhi) ;
1348  if (Dist < halfCarTolerance)
1349  {
1350  sd = Dist/Comp ;
1351 
1352  if ( sd < snxt )
1353  {
1354  if ( sd < 0 ) { sd = 0.0; }
1355 
1356  zi = p.z() + sd*v.z() ;
1357 
1358  if (std::fabs(zi) <= tolODz)
1359  {
1360  xi = p.x() + sd*v.x() ;
1361  yi = p.y() + sd*v.y() ;
1362  rhoi2 = xi*xi + yi*yi ;
1363  tolORMin2 = (rMinOAv + zi*tanRMin)*(rMinOAv + zi*tanRMin) ;
1364  tolORMax2 = (rMaxOAv + zi*tanRMax)*(rMaxOAv + zi*tanRMax) ;
1365 
1366  if ( (rhoi2 >= tolORMin2) && (rhoi2 <= tolORMax2) )
1367  {
1368  // z and r intersections good - check intersecting with
1369  // correct half-plane
1370 
1371  if ( (yi*cosCPhi - xi*sinCPhi) >= 0.0 ) { snxt = sd; }
1372  }
1373  }
1374  }
1375  }
1376  }
1377  }
1378  if (snxt < halfCarTolerance) { snxt = 0.; }
1379 
1380  return snxt ;
1381 }
TTree * t1
Definition: plottest35.C:26
G4double halfRadTolerance
Definition: G4Cons.hh:233
G4double fDz
Definition: G4Cons.hh:220
static const G4double kInfinity
Definition: geomdefs.hh:42
G4bool fPhiFullCone
Definition: G4Cons.hh:229
Float_t d
CLHEP::Hep3Vector G4ThreeVector
G4double sinEPhi
Definition: G4Cons.hh:224
G4double fRmin1
Definition: G4Cons.hh:220
G4double fRmax2
Definition: G4Cons.hh:220
G4double cosCPhi
Definition: G4Cons.hh:224
G4double kRadTolerance
Definition: G4Cons.hh:216
G4double fRmin2
Definition: G4Cons.hh:220
G4double sinSPhi
Definition: G4Cons.hh:224
G4double halfCarTolerance
Definition: G4Cons.hh:233
G4double cosSPhi
Definition: G4Cons.hh:224
double x() const
double dot(const Hep3Vector &) const
double y() const
double z() const
G4double sinCPhi
Definition: G4Cons.hh:224
G4double cosEPhi
Definition: G4Cons.hh:224
TTree * t2
Definition: plottest35.C:36
G4double cosHDPhiOT
Definition: G4Cons.hh:224
G4double cosHDPhiIT
Definition: G4Cons.hh:224
double G4double
Definition: G4Types.hh:76
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
Definition: G4Cons.cc:720
G4double fRmax1
Definition: G4Cons.hh:220
Here is the call graph for this function:

◆ DistanceToIn() [2/2]

G4double G4Cons::DistanceToIn ( const G4ThreeVector p) const
virtual

Implements G4VSolid.

Definition at line 1390 of file G4Cons.cc.

1391 {
1392  G4double safe=0.0, rho, safeR1, safeR2, safeZ, safePhi, cosPsi ;
1393  G4double tanRMin, secRMin, pRMin ;
1394  G4double tanRMax, secRMax, pRMax ;
1395 
1396  rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
1397  safeZ = std::fabs(p.z()) - fDz ;
1398 
1399  if ( fRmin1 || fRmin2 )
1400  {
1401  tanRMin = (fRmin2 - fRmin1)*0.5/fDz ;
1402  secRMin = std::sqrt(1.0 + tanRMin*tanRMin) ;
1403  pRMin = tanRMin*p.z() + (fRmin1 + fRmin2)*0.5 ;
1404  safeR1 = (pRMin - rho)/secRMin ;
1405 
1406  tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;
1407  secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
1408  pRMax = tanRMax*p.z() + (fRmax1 + fRmax2)*0.5 ;
1409  safeR2 = (rho - pRMax)/secRMax ;
1410 
1411  if ( safeR1 > safeR2) { safe = safeR1; }
1412  else { safe = safeR2; }
1413  }
1414  else
1415  {
1416  tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;
1417  secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
1418  pRMax = tanRMax*p.z() + (fRmax1 + fRmax2)*0.5 ;
1419  safe = (rho - pRMax)/secRMax ;
1420  }
1421  if ( safeZ > safe ) { safe = safeZ; }
1422 
1423  if ( !fPhiFullCone && rho )
1424  {
1425  // Psi=angle from central phi to point
1426 
1427  cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/rho ;
1428 
1429  if ( cosPsi < std::cos(fDPhi*0.5) ) // Point lies outside phi range
1430  {
1431  if ( (p.y()*cosCPhi - p.x()*sinCPhi) <= 0.0 )
1432  {
1433  safePhi = std::fabs(p.x()*std::sin(fSPhi)-p.y()*std::cos(fSPhi));
1434  }
1435  else
1436  {
1437  safePhi = std::fabs(p.x()*sinEPhi-p.y()*cosEPhi);
1438  }
1439  if ( safePhi > safe ) { safe = safePhi; }
1440  }
1441  }
1442  if ( safe < 0.0 ) { safe = 0.0; }
1443 
1444  return safe ;
1445 }
G4double fDz
Definition: G4Cons.hh:220
G4bool fPhiFullCone
Definition: G4Cons.hh:229
G4double sinEPhi
Definition: G4Cons.hh:224
G4double fRmin1
Definition: G4Cons.hh:220
G4double fRmax2
Definition: G4Cons.hh:220
G4double cosCPhi
Definition: G4Cons.hh:224
G4double fRmin2
Definition: G4Cons.hh:220
G4double fSPhi
Definition: G4Cons.hh:220
double x() const
G4double fDPhi
Definition: G4Cons.hh:220
double y() const
double z() const
G4double sinCPhi
Definition: G4Cons.hh:224
G4double cosEPhi
Definition: G4Cons.hh:224
double G4double
Definition: G4Types.hh:76
G4double fRmax1
Definition: G4Cons.hh:220
Here is the call graph for this function:

◆ DistanceToOut() [1/2]

G4double G4Cons::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 1452 of file G4Cons.cc.

1457 {
1458  ESide side = kNull, sider = kNull, sidephi = kNull;
1459 
1460  G4double snxt,srd,sphi,pdist ;
1461 
1462  G4double tanRMax, secRMax, rMaxAv ; // Data for outer cone
1463  G4double tanRMin, secRMin, rMinAv ; // Data for inner cone
1464 
1465  G4double t1, t2, t3, rout, rin, nt1, nt2, nt3 ;
1466  G4double b, c, d, sr2, sr3 ;
1467 
1468  // Vars for intersection within tolerance
1469 
1470  ESide sidetol = kNull ;
1471  G4double slentol = kInfinity ;
1472 
1473  // Vars for phi intersection:
1474 
1475  G4double pDistS, compS, pDistE, compE, sphi2, xi, yi, risec, vphi ;
1476  G4double zi, ri, deltaRoi2 ;
1477 
1478  // Z plane intersection
1479 
1480  if ( v.z() > 0.0 )
1481  {
1482  pdist = fDz - p.z() ;
1483 
1484  if (pdist > halfCarTolerance)
1485  {
1486  snxt = pdist/v.z() ;
1487  side = kPZ ;
1488  }
1489  else
1490  {
1491  if (calcNorm)
1492  {
1493  *n = G4ThreeVector(0,0,1) ;
1494  *validNorm = true ;
1495  }
1496  return snxt = 0.0;
1497  }
1498  }
1499  else if ( v.z() < 0.0 )
1500  {
1501  pdist = fDz + p.z() ;
1502 
1503  if ( pdist > halfCarTolerance)
1504  {
1505  snxt = -pdist/v.z() ;
1506  side = kMZ ;
1507  }
1508  else
1509  {
1510  if ( calcNorm )
1511  {
1512  *n = G4ThreeVector(0,0,-1) ;
1513  *validNorm = true ;
1514  }
1515  return snxt = 0.0 ;
1516  }
1517  }
1518  else // Travel perpendicular to z axis
1519  {
1520  snxt = kInfinity ;
1521  side = kNull ;
1522  }
1523 
1524  // Radial Intersections
1525  //
1526  // Intersection with outer cone (possible return) and
1527  // inner cone (must also check phi)
1528  //
1529  // Intersection point (xi,yi,zi) on line x=p.x+t*v.x etc.
1530  //
1531  // Intersects with x^2+y^2=(a*z+b)^2
1532  //
1533  // where a=tanRMax or tanRMin
1534  // b=rMaxAv or rMinAv
1535  //
1536  // (vx^2+vy^2-(a*vz)^2)t^2+2t(pxvx+pyvy-a*vz(a*pz+b))+px^2+py^2-(a*pz+b)^2=0 ;
1537  // t1 t2 t3
1538  //
1539  // \--------u-------/ \-----------v----------/ \---------w--------/
1540 
1541  tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;
1542  secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
1543  rMaxAv = (fRmax1 + fRmax2)*0.5 ;
1544 
1545 
1546  t1 = 1.0 - v.z()*v.z() ; // since v normalised
1547  t2 = p.x()*v.x() + p.y()*v.y() ;
1548  t3 = p.x()*p.x() + p.y()*p.y() ;
1549  rout = tanRMax*p.z() + rMaxAv ;
1550 
1551  nt1 = t1 - (tanRMax*v.z())*(tanRMax*v.z()) ;
1552  nt2 = t2 - tanRMax*v.z()*rout ;
1553  nt3 = t3 - rout*rout ;
1554 
1555  if (v.z() > 0.0)
1556  {
1557  deltaRoi2 = snxt*snxt*t1 + 2*snxt*t2 + t3
1558  - fRmax2*(fRmax2 + kRadTolerance*secRMax);
1559  }
1560  else if ( v.z() < 0.0 )
1561  {
1562  deltaRoi2 = snxt*snxt*t1 + 2*snxt*t2 + t3
1563  - fRmax1*(fRmax1 + kRadTolerance*secRMax);
1564  }
1565  else
1566  {
1567  deltaRoi2 = 1.0;
1568  }
1569 
1570  if ( nt1 && (deltaRoi2 > 0.0) )
1571  {
1572  // Equation quadratic => 2 roots : second root must be leaving
1573 
1574  b = nt2/nt1 ;
1575  c = nt3/nt1 ;
1576  d = b*b - c ;
1577 
1578  if ( d >= 0 )
1579  {
1580  // Check if on outer cone & heading outwards
1581  // NOTE: Should use rho-rout>-kRadTolerance*0.5
1582 
1583  if (nt3 > -halfRadTolerance && nt2 >= 0 )
1584  {
1585  if (calcNorm)
1586  {
1587  risec = std::sqrt(t3)*secRMax ;
1588  *validNorm = true ;
1589  *n = G4ThreeVector(p.x()/risec,p.y()/risec,-tanRMax/secRMax);
1590  }
1591  return snxt=0 ;
1592  }
1593  else
1594  {
1595  sider = kRMax ;
1596  if (b>0) { srd = -b - std::sqrt(d); }
1597  else { srd = c/(-b+std::sqrt(d)) ; }
1598 
1599  zi = p.z() + srd*v.z() ;
1600  ri = tanRMax*zi + rMaxAv ;
1601 
1602  if ((ri >= 0) && (-halfRadTolerance <= srd) && (srd <= halfRadTolerance))
1603  {
1604  // An intersection within the tolerance
1605  // we will Store it in case it is good -
1606  //
1607  slentol = srd ;
1608  sidetol = kRMax ;
1609  }
1610  if ( (ri < 0) || (srd < halfRadTolerance) )
1611  {
1612  // Safety: if both roots -ve ensure that srd cannot `win'
1613  // distance to out
1614 
1615  if (b>0) { sr2 = c/(-b-std::sqrt(d)); }
1616  else { sr2 = -b + std::sqrt(d); }
1617  zi = p.z() + sr2*v.z() ;
1618  ri = tanRMax*zi + rMaxAv ;
1619 
1620  if ((ri >= 0) && (sr2 > halfRadTolerance))
1621  {
1622  srd = sr2;
1623  }
1624  else
1625  {
1626  srd = kInfinity ;
1627 
1628  if( (-halfRadTolerance <= sr2) && ( sr2 <= halfRadTolerance) )
1629  {
1630  // An intersection within the tolerance.
1631  // Storing it in case it is good.
1632 
1633  slentol = sr2 ;
1634  sidetol = kRMax ;
1635  }
1636  }
1637  }
1638  }
1639  }
1640  else
1641  {
1642  // No intersection with outer cone & not parallel
1643  // -> already outside, no intersection
1644 
1645  if ( calcNorm )
1646  {
1647  risec = std::sqrt(t3)*secRMax;
1648  *validNorm = true;
1649  *n = G4ThreeVector(p.x()/risec,p.y()/risec,-tanRMax/secRMax);
1650  }
1651  return snxt = 0.0 ;
1652  }
1653  }
1654  else if ( nt2 && (deltaRoi2 > 0.0) )
1655  {
1656  // Linear case (only one intersection) => point outside outer cone
1657 
1658  if ( calcNorm )
1659  {
1660  risec = std::sqrt(t3)*secRMax;
1661  *validNorm = true;
1662  *n = G4ThreeVector(p.x()/risec,p.y()/risec,-tanRMax/secRMax);
1663  }
1664  return snxt = 0.0 ;
1665  }
1666  else
1667  {
1668  // No intersection -> parallel to outer cone
1669  // => Z or inner cone intersection
1670 
1671  srd = kInfinity ;
1672  }
1673 
1674  // Check possible intersection within tolerance
1675 
1676  if ( slentol <= halfCarTolerance )
1677  {
1678  // An intersection within the tolerance was found.
1679  // We must accept it only if the momentum points outwards.
1680  //
1681  // G4ThreeVector ptTol ; // The point of the intersection
1682  // ptTol= p + slentol*v ;
1683  // ri=tanRMax*zi+rMaxAv ;
1684  //
1685  // Calculate a normal vector, as below
1686 
1687  xi = p.x() + slentol*v.x();
1688  yi = p.y() + slentol*v.y();
1689  risec = std::sqrt(xi*xi + yi*yi)*secRMax;
1690  G4ThreeVector Normal = G4ThreeVector(xi/risec,yi/risec,-tanRMax/secRMax);
1691 
1692  if ( Normal.dot(v) > 0 ) // We will leave the Cone immediatelly
1693  {
1694  if ( calcNorm )
1695  {
1696  *n = Normal.unit() ;
1697  *validNorm = true ;
1698  }
1699  return snxt = 0.0 ;
1700  }
1701  else // On the surface, but not heading out so we ignore this intersection
1702  { // (as it is within tolerance).
1703  slentol = kInfinity ;
1704  }
1705  }
1706 
1707  // Inner Cone intersection
1708 
1709  if ( fRmin1 || fRmin2 )
1710  {
1711  tanRMin = (fRmin2 - fRmin1)*0.5/fDz ;
1712  nt1 = t1 - (tanRMin*v.z())*(tanRMin*v.z()) ;
1713 
1714  if ( nt1 )
1715  {
1716  secRMin = std::sqrt(1.0 + tanRMin*tanRMin) ;
1717  rMinAv = (fRmin1 + fRmin2)*0.5 ;
1718  rin = tanRMin*p.z() + rMinAv ;
1719  nt2 = t2 - tanRMin*v.z()*rin ;
1720  nt3 = t3 - rin*rin ;
1721 
1722  // Equation quadratic => 2 roots : first root must be leaving
1723 
1724  b = nt2/nt1 ;
1725  c = nt3/nt1 ;
1726  d = b*b - c ;
1727 
1728  if ( d >= 0.0 )
1729  {
1730  // NOTE: should be rho-rin<kRadTolerance*0.5,
1731  // but using squared versions for efficiency
1732 
1733  if (nt3 < kRadTolerance*(rin + kRadTolerance*0.25))
1734  {
1735  if ( nt2 < 0.0 )
1736  {
1737  if (calcNorm) { *validNorm = false; }
1738  return snxt = 0.0;
1739  }
1740  }
1741  else
1742  {
1743  if (b>0) { sr2 = -b - std::sqrt(d); }
1744  else { sr2 = c/(-b+std::sqrt(d)); }
1745  zi = p.z() + sr2*v.z() ;
1746  ri = tanRMin*zi + rMinAv ;
1747 
1748  if( (ri>=0.0)&&(-halfRadTolerance<=sr2)&&(sr2<=halfRadTolerance) )
1749  {
1750  // An intersection within the tolerance
1751  // storing it in case it is good.
1752 
1753  slentol = sr2 ;
1754  sidetol = kRMax ;
1755  }
1756  if( (ri<0) || (sr2 < halfRadTolerance) )
1757  {
1758  if (b>0) { sr3 = c/(-b-std::sqrt(d)); }
1759  else { sr3 = -b + std::sqrt(d) ; }
1760 
1761  // Safety: if both roots -ve ensure that srd cannot `win'
1762  // distancetoout
1763 
1764  if ( sr3 > halfRadTolerance )
1765  {
1766  if( sr3 < srd )
1767  {
1768  zi = p.z() + sr3*v.z() ;
1769  ri = tanRMin*zi + rMinAv ;
1770 
1771  if ( ri >= 0.0 )
1772  {
1773  srd=sr3 ;
1774  sider=kRMin ;
1775  }
1776  }
1777  }
1778  else if ( sr3 > -halfRadTolerance )
1779  {
1780  // Intersection in tolerance. Store to check if it's good
1781 
1782  slentol = sr3 ;
1783  sidetol = kRMin ;
1784  }
1785  }
1786  else if ( (sr2 < srd) && (sr2 > halfCarTolerance) )
1787  {
1788  srd = sr2 ;
1789  sider = kRMin ;
1790  }
1791  else if (sr2 > -halfCarTolerance)
1792  {
1793  // Intersection in tolerance. Store to check if it's good
1794 
1795  slentol = sr2 ;
1796  sidetol = kRMin ;
1797  }
1798  if( slentol <= halfCarTolerance )
1799  {
1800  // An intersection within the tolerance was found.
1801  // We must accept it only if the momentum points outwards.
1802 
1803  G4ThreeVector Normal ;
1804 
1805  // Calculate a normal vector, as below
1806 
1807  xi = p.x() + slentol*v.x() ;
1808  yi = p.y() + slentol*v.y() ;
1809  if( sidetol==kRMax )
1810  {
1811  risec = std::sqrt(xi*xi + yi*yi)*secRMax ;
1812  Normal = G4ThreeVector(xi/risec,yi/risec,-tanRMax/secRMax) ;
1813  }
1814  else
1815  {
1816  risec = std::sqrt(xi*xi + yi*yi)*secRMin ;
1817  Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin) ;
1818  }
1819  if( Normal.dot(v) > 0 )
1820  {
1821  // We will leave the cone immediately
1822 
1823  if( calcNorm )
1824  {
1825  *n = Normal.unit() ;
1826  *validNorm = true ;
1827  }
1828  return snxt = 0.0 ;
1829  }
1830  else
1831  {
1832  // On the surface, but not heading out so we ignore this
1833  // intersection (as it is within tolerance).
1834 
1835  slentol = kInfinity ;
1836  }
1837  }
1838  }
1839  }
1840  }
1841  }
1842 
1843  // Linear case => point outside inner cone ---> outer cone intersect
1844  //
1845  // Phi Intersection
1846 
1847  if ( !fPhiFullCone )
1848  {
1849  // add angle calculation with correction
1850  // of the difference in domain of atan2 and Sphi
1851 
1852  vphi = std::atan2(v.y(),v.x()) ;
1853 
1854  if ( vphi < fSPhi - halfAngTolerance ) { vphi += twopi; }
1855  else if ( vphi > fSPhi + fDPhi + halfAngTolerance ) { vphi -= twopi; }
1856 
1857  if ( p.x() || p.y() ) // Check if on z axis (rho not needed later)
1858  {
1859  // pDist -ve when inside
1860 
1861  pDistS = p.x()*sinSPhi - p.y()*cosSPhi ;
1862  pDistE = -p.x()*sinEPhi + p.y()*cosEPhi ;
1863 
1864  // Comp -ve when in direction of outwards normal
1865 
1866  compS = -sinSPhi*v.x() + cosSPhi*v.y() ;
1867  compE = sinEPhi*v.x() - cosEPhi*v.y() ;
1868 
1869  sidephi = kNull ;
1870 
1871  if( ( (fDPhi <= pi) && ( (pDistS <= halfCarTolerance)
1872  && (pDistE <= halfCarTolerance) ) )
1873  || ( (fDPhi > pi) && !((pDistS > halfCarTolerance)
1874  && (pDistE > halfCarTolerance) ) ) )
1875  {
1876  // Inside both phi *full* planes
1877  if ( compS < 0 )
1878  {
1879  sphi = pDistS/compS ;
1880  if (sphi >= -halfCarTolerance)
1881  {
1882  xi = p.x() + sphi*v.x() ;
1883  yi = p.y() + sphi*v.y() ;
1884 
1885  // Check intersecting with correct half-plane
1886  // (if not -> no intersect)
1887  //
1888  if ( (std::fabs(xi)<=kCarTolerance)
1889  && (std::fabs(yi)<=kCarTolerance) )
1890  {
1891  sidephi= kSPhi;
1892  if ( ( fSPhi-halfAngTolerance <= vphi )
1893  && ( fSPhi+fDPhi+halfAngTolerance >=vphi ) )
1894  {
1895  sphi = kInfinity;
1896  }
1897  }
1898  else
1899  if ( (yi*cosCPhi-xi*sinCPhi)>=0 )
1900  {
1901  sphi = kInfinity ;
1902  }
1903  else
1904  {
1905  sidephi = kSPhi ;
1906  if ( pDistS > -halfCarTolerance )
1907  {
1908  sphi = 0.0 ; // Leave by sphi immediately
1909  }
1910  }
1911  }
1912  else
1913  {
1914  sphi = kInfinity ;
1915  }
1916  }
1917  else
1918  {
1919  sphi = kInfinity ;
1920  }
1921 
1922  if ( compE < 0 )
1923  {
1924  sphi2 = pDistE/compE ;
1925 
1926  // Only check further if < starting phi intersection
1927  //
1928  if ( (sphi2 > -halfCarTolerance) && (sphi2 < sphi) )
1929  {
1930  xi = p.x() + sphi2*v.x() ;
1931  yi = p.y() + sphi2*v.y() ;
1932 
1933  // Check intersecting with correct half-plane
1934 
1935  if ( (std::fabs(xi)<=kCarTolerance)
1936  && (std::fabs(yi)<=kCarTolerance) )
1937  {
1938  // Leaving via ending phi
1939 
1940  if(!( (fSPhi-halfAngTolerance <= vphi)
1941  && (fSPhi+fDPhi+halfAngTolerance >= vphi) ) )
1942  {
1943  sidephi = kEPhi ;
1944  if ( pDistE <= -halfCarTolerance ) { sphi = sphi2; }
1945  else { sphi = 0.0; }
1946  }
1947  }
1948  else // Check intersecting with correct half-plane
1949  if ( yi*cosCPhi-xi*sinCPhi >= 0 )
1950  {
1951  // Leaving via ending phi
1952 
1953  sidephi = kEPhi ;
1954  if ( pDistE <= -halfCarTolerance ) { sphi = sphi2; }
1955  else { sphi = 0.0; }
1956  }
1957  }
1958  }
1959  }
1960  else
1961  {
1962  sphi = kInfinity ;
1963  }
1964  }
1965  else
1966  {
1967  // On z axis + travel not || to z axis -> if phi of vector direction
1968  // within phi of shape, Step limited by rmax, else Step =0
1969 
1970  if ( (fSPhi-halfAngTolerance <= vphi)
1971  && (vphi <= fSPhi+fDPhi+halfAngTolerance) )
1972  {
1973  sphi = kInfinity ;
1974  }
1975  else
1976  {
1977  sidephi = kSPhi ; // arbitrary
1978  sphi = 0.0 ;
1979  }
1980  }
1981  if ( sphi < snxt ) // Order intersecttions
1982  {
1983  snxt=sphi ;
1984  side=sidephi ;
1985  }
1986  }
1987  if ( srd < snxt ) // Order intersections
1988  {
1989  snxt = srd ;
1990  side = sider ;
1991  }
1992  if (calcNorm)
1993  {
1994  switch(side)
1995  { // Note: returned vector not normalised
1996  case kRMax: // (divide by frmax for unit vector)
1997  xi = p.x() + snxt*v.x() ;
1998  yi = p.y() + snxt*v.y() ;
1999  risec = std::sqrt(xi*xi + yi*yi)*secRMax ;
2000  *n = G4ThreeVector(xi/risec,yi/risec,-tanRMax/secRMax) ;
2001  *validNorm = true ;
2002  break ;
2003  case kRMin:
2004  *validNorm = false ; // Rmin is inconvex
2005  break ;
2006  case kSPhi:
2007  if ( fDPhi <= pi )
2008  {
2009  *n = G4ThreeVector(sinSPhi, -cosSPhi, 0);
2010  *validNorm = true ;
2011  }
2012  else
2013  {
2014  *validNorm = false ;
2015  }
2016  break ;
2017  case kEPhi:
2018  if ( fDPhi <= pi )
2019  {
2020  *n = G4ThreeVector(-sinEPhi, cosEPhi, 0);
2021  *validNorm = true ;
2022  }
2023  else
2024  {
2025  *validNorm = false ;
2026  }
2027  break ;
2028  case kPZ:
2029  *n = G4ThreeVector(0,0,1) ;
2030  *validNorm = true ;
2031  break ;
2032  case kMZ:
2033  *n = G4ThreeVector(0,0,-1) ;
2034  *validNorm = true ;
2035  break ;
2036  default:
2037  G4cout << G4endl ;
2038  DumpInfo();
2039  std::ostringstream message;
2040  G4int oldprc = message.precision(16) ;
2041  message << "Undefined side for valid surface normal to solid."
2042  << G4endl
2043  << "Position:" << G4endl << G4endl
2044  << "p.x() = " << p.x()/mm << " mm" << G4endl
2045  << "p.y() = " << p.y()/mm << " mm" << G4endl
2046  << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl
2047  << "pho at z = " << std::sqrt( p.x()*p.x()+p.y()*p.y() )/mm
2048  << " mm" << G4endl << G4endl ;
2049  if( p.x() != 0. || p.y() != 0.)
2050  {
2051  message << "point phi = " << std::atan2(p.y(),p.x())/degree
2052  << " degree" << G4endl << G4endl ;
2053  }
2054  message << "Direction:" << G4endl << G4endl
2055  << "v.x() = " << v.x() << G4endl
2056  << "v.y() = " << v.y() << G4endl
2057  << "v.z() = " << v.z() << G4endl<< G4endl
2058  << "Proposed distance :" << G4endl<< G4endl
2059  << "snxt = " << snxt/mm << " mm" << G4endl ;
2060  message.precision(oldprc) ;
2061  G4Exception("G4Cons::DistanceToOut(p,v,..)","GeomSolids1002",
2062  JustWarning, message) ;
2063  break ;
2064  }
2065  }
2066  if (snxt < halfCarTolerance) { snxt = 0.; }
2067 
2068  return snxt ;
2069 }
TTree * t1
Definition: plottest35.C:26
G4double halfRadTolerance
Definition: G4Cons.hh:233
G4double fDz
Definition: G4Cons.hh:220
static const G4double kInfinity
Definition: geomdefs.hh:42
G4bool fPhiFullCone
Definition: G4Cons.hh:229
Float_t d
CLHEP::Hep3Vector G4ThreeVector
G4double sinEPhi
Definition: G4Cons.hh:224
G4double fRmin1
Definition: G4Cons.hh:220
int G4int
Definition: G4Types.hh:78
G4double fRmax2
Definition: G4Cons.hh:220
void DumpInfo() const
G4double cosCPhi
Definition: G4Cons.hh:224
G4double kRadTolerance
Definition: G4Cons.hh:216
G4GLOB_DLL std::ostream G4cout
G4double fRmin2
Definition: G4Cons.hh:220
G4double sinSPhi
Definition: G4Cons.hh:224
G4double fSPhi
Definition: G4Cons.hh:220
G4double halfCarTolerance
Definition: G4Cons.hh:233
Hep3Vector unit() const
static const double twopi
Definition: G4SIunits.hh:75
G4double cosSPhi
Definition: G4Cons.hh:224
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
double x() const
double dot(const Hep3Vector &) const
G4double fDPhi
Definition: G4Cons.hh:220
static const double pi
Definition: G4SIunits.hh:74
double y() const
G4double halfAngTolerance
Definition: G4Cons.hh:233
double z() const
G4double sinCPhi
Definition: G4Cons.hh:224
static const double degree
Definition: G4SIunits.hh:143
G4double cosEPhi
Definition: G4Cons.hh:224
TTree * t2
Definition: plottest35.C:36
#define G4endl
Definition: G4ios.hh:61
G4double kCarTolerance
Definition: G4VSolid.hh:304
ESide
Definition: G4Cons.cc:71
double G4double
Definition: G4Types.hh:76
static const double mm
Definition: G4SIunits.hh:114
G4double fRmax1
Definition: G4Cons.hh:220
Here is the call graph for this function:

◆ DistanceToOut() [2/2]

G4double G4Cons::DistanceToOut ( const G4ThreeVector p) const
virtual

Implements G4VSolid.

Definition at line 2075 of file G4Cons.cc.

2076 {
2077  G4double safe=0.0, rho, safeR1, safeR2, safeZ, safePhi;
2078  G4double tanRMin, secRMin, pRMin;
2079  G4double tanRMax, secRMax, pRMax;
2080 
2081 #ifdef G4CSGDEBUG
2082  if( Inside(p) == kOutside )
2083  {
2084  G4int oldprc=G4cout.precision(16) ;
2085  G4cout << G4endl ;
2086  DumpInfo();
2087  G4cout << "Position:" << G4endl << G4endl ;
2088  G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl ;
2089  G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ;
2090  G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ;
2091  G4cout << "pho at z = " << std::sqrt( p.x()*p.x()+p.y()*p.y() )/mm
2092  << " mm" << G4endl << G4endl ;
2093  if( (p.x() != 0.) || (p.x() != 0.) )
2094  {
2095  G4cout << "point phi = " << std::atan2(p.y(),p.x())/degree
2096  << " degree" << G4endl << G4endl ;
2097  }
2098  G4cout.precision(oldprc) ;
2099  G4Exception("G4Cons::DistanceToOut(p)", "GeomSolids1002",
2100  JustWarning, "Point p is outside !?" );
2101  }
2102 #endif
2103 
2104  rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
2105  safeZ = fDz - std::fabs(p.z()) ;
2106 
2107  if (fRmin1 || fRmin2)
2108  {
2109  tanRMin = (fRmin2 - fRmin1)*0.5/fDz ;
2110  secRMin = std::sqrt(1.0 + tanRMin*tanRMin) ;
2111  pRMin = tanRMin*p.z() + (fRmin1 + fRmin2)*0.5 ;
2112  safeR1 = (rho - pRMin)/secRMin ;
2113  }
2114  else
2115  {
2116  safeR1 = kInfinity ;
2117  }
2118 
2119  tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;
2120  secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
2121  pRMax = tanRMax*p.z() + (fRmax1+fRmax2)*0.5 ;
2122  safeR2 = (pRMax - rho)/secRMax ;
2123 
2124  if (safeR1 < safeR2) { safe = safeR1; }
2125  else { safe = safeR2; }
2126  if (safeZ < safe) { safe = safeZ ; }
2127 
2128  // Check if phi divided, Calc distances closest phi plane
2129 
2130  if (!fPhiFullCone)
2131  {
2132  // Above/below central phi of G4Cons?
2133 
2134  if ( (p.y()*cosCPhi - p.x()*sinCPhi) <= 0 )
2135  {
2136  safePhi = -(p.x()*sinSPhi - p.y()*cosSPhi) ;
2137  }
2138  else
2139  {
2140  safePhi = (p.x()*sinEPhi - p.y()*cosEPhi) ;
2141  }
2142  if (safePhi < safe) { safe = safePhi; }
2143  }
2144  if ( safe < 0 ) { safe = 0; }
2145 
2146  return safe ;
2147 }
G4double fDz
Definition: G4Cons.hh:220
static const G4double kInfinity
Definition: geomdefs.hh:42
G4bool fPhiFullCone
Definition: G4Cons.hh:229
G4double sinEPhi
Definition: G4Cons.hh:224
G4double fRmin1
Definition: G4Cons.hh:220
int G4int
Definition: G4Types.hh:78
G4double fRmax2
Definition: G4Cons.hh:220
void DumpInfo() const
G4double cosCPhi
Definition: G4Cons.hh:224
G4GLOB_DLL std::ostream G4cout
G4double fRmin2
Definition: G4Cons.hh:220
G4double sinSPhi
Definition: G4Cons.hh:224
EInside Inside(const G4ThreeVector &p) const
Definition: G4Cons.cc:205
G4double cosSPhi
Definition: G4Cons.hh:224
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
double x() const
double y() const
double z() const
G4double sinCPhi
Definition: G4Cons.hh:224
static const double degree
Definition: G4SIunits.hh:143
G4double cosEPhi
Definition: G4Cons.hh:224
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
static const double mm
Definition: G4SIunits.hh:114
G4double fRmax1
Definition: G4Cons.hh:220
Here is the call graph for this function:

◆ GetCubicVolume()

G4double G4Cons::GetCubicVolume ( )
inlinevirtual

Reimplemented from G4VSolid.

◆ GetDeltaPhiAngle()

G4double G4Cons::GetDeltaPhiAngle ( ) const
inline
Here is the caller graph for this function:

◆ GetDPhi()

G4double G4Cons::GetDPhi ( ) const
inline

◆ GetDz()

G4double G4Cons::GetDz ( ) const
inline

◆ GetEntityType()

G4GeometryType G4Cons::GetEntityType ( ) const
virtual

Implements G4VSolid.

Definition at line 2247 of file G4Cons.cc.

◆ GetInnerRadiusMinusZ()

G4double G4Cons::GetInnerRadiusMinusZ ( ) const
inline
Here is the caller graph for this function:

◆ GetInnerRadiusPlusZ()

G4double G4Cons::GetInnerRadiusPlusZ ( ) const
inline
Here is the caller graph for this function:

◆ GetOuterRadiusMinusZ()

G4double G4Cons::GetOuterRadiusMinusZ ( ) const
inline
Here is the caller graph for this function:

◆ GetOuterRadiusPlusZ()

G4double G4Cons::GetOuterRadiusPlusZ ( ) const
inline
Here is the caller graph for this function:

◆ GetPointOnSurface()

G4ThreeVector G4Cons::GetPointOnSurface ( ) const
virtual

Reimplemented from G4VSolid.

Definition at line 2292 of file G4Cons.cc.

2293 {
2294  // declare working variables
2295  //
2296  G4double Aone, Atwo, Athree, Afour, Afive, slin, slout, phi;
2297  G4double zRand, cosu, sinu, rRand1, rRand2, chose, rone, rtwo, qone, qtwo;
2298  rone = (fRmax1-fRmax2)/(2.*fDz);
2299  rtwo = (fRmin1-fRmin2)/(2.*fDz);
2300  qone=0.; qtwo=0.;
2301  if(fRmax1!=fRmax2) { qone = fDz*(fRmax1+fRmax2)/(fRmax1-fRmax2); }
2302  if(fRmin1!=fRmin2) { qtwo = fDz*(fRmin1+fRmin2)/(fRmin1-fRmin2); }
2303  slin = std::sqrt(sqr(fRmin1-fRmin2)+sqr(2.*fDz));
2304  slout = std::sqrt(sqr(fRmax1-fRmax2)+sqr(2.*fDz));
2305  Aone = 0.5*fDPhi*(fRmax2 + fRmax1)*slout;
2306  Atwo = 0.5*fDPhi*(fRmin2 + fRmin1)*slin;
2307  Athree = 0.5*fDPhi*(fRmax1*fRmax1-fRmin1*fRmin1);
2308  Afour = 0.5*fDPhi*(fRmax2*fRmax2-fRmin2*fRmin2);
2309  Afive = fDz*(fRmax1-fRmin1+fRmax2-fRmin2);
2310 
2312  cosu = std::cos(phi); sinu = std::sin(phi);
2313  rRand1 = GetRadiusInRing(fRmin1, fRmax1);
2314  rRand2 = GetRadiusInRing(fRmin2, fRmax2);
2315 
2316  if ( (fSPhi == 0.) && fPhiFullCone ) { Afive = 0.; }
2317  chose = G4RandFlat::shoot(0.,Aone+Atwo+Athree+Afour+2.*Afive);
2318 
2319  if( (chose >= 0.) && (chose < Aone) )
2320  {
2321  if(fRmin1 != fRmin2)
2322  {
2323  zRand = G4RandFlat::shoot(-1.*fDz,fDz);
2324  return G4ThreeVector (rtwo*cosu*(qtwo-zRand),
2325  rtwo*sinu*(qtwo-zRand), zRand);
2326  }
2327  else
2328  {
2329  return G4ThreeVector(fRmin1*cosu, fRmin2*sinu,
2330  G4RandFlat::shoot(-1.*fDz,fDz));
2331  }
2332  }
2333  else if( (chose >= Aone) && (chose <= Aone + Atwo) )
2334  {
2335  if(fRmax1 != fRmax2)
2336  {
2337  zRand = G4RandFlat::shoot(-1.*fDz,fDz);
2338  return G4ThreeVector (rone*cosu*(qone-zRand),
2339  rone*sinu*(qone-zRand), zRand);
2340  }
2341  else
2342  {
2343  return G4ThreeVector(fRmax1*cosu, fRmax2*sinu,
2344  G4RandFlat::shoot(-1.*fDz,fDz));
2345  }
2346  }
2347  else if( (chose >= Aone + Atwo) && (chose < Aone + Atwo + Athree) )
2348  {
2349  return G4ThreeVector (rRand1*cosu, rRand1*sinu, -1*fDz);
2350  }
2351  else if( (chose >= Aone + Atwo + Athree)
2352  && (chose < Aone + Atwo + Athree + Afour) )
2353  {
2354  return G4ThreeVector (rRand2*cosu,rRand2*sinu,fDz);
2355  }
2356  else if( (chose >= Aone + Atwo + Athree + Afour)
2357  && (chose < Aone + Atwo + Athree + Afour + Afive) )
2358  {
2359  zRand = G4RandFlat::shoot(-1.*fDz,fDz);
2360  rRand1 = G4RandFlat::shoot(fRmin2-((zRand-fDz)/(2.*fDz))*(fRmin1-fRmin2),
2361  fRmax2-((zRand-fDz)/(2.*fDz))*(fRmax1-fRmax2));
2362  return G4ThreeVector (rRand1*std::cos(fSPhi),
2363  rRand1*std::sin(fSPhi), zRand);
2364  }
2365  else
2366  {
2367  zRand = G4RandFlat::shoot(-1.*fDz,fDz);
2368  rRand1 = G4RandFlat::shoot(fRmin2-((zRand-fDz)/(2.*fDz))*(fRmin1-fRmin2),
2369  fRmax2-((zRand-fDz)/(2.*fDz))*(fRmax1-fRmax2));
2370  return G4ThreeVector (rRand1*std::cos(fSPhi+fDPhi),
2371  rRand1*std::sin(fSPhi+fDPhi), zRand);
2372  }
2373 }
ThreeVector shoot(const G4int Ap, const G4int Af)
G4double fDz
Definition: G4Cons.hh:220
G4bool fPhiFullCone
Definition: G4Cons.hh:229
CLHEP::Hep3Vector G4ThreeVector
G4double GetRadiusInRing(G4double rmin, G4double rmax) const
Definition: G4CSGSolid.cc:111
G4double fRmin1
Definition: G4Cons.hh:220
G4double fRmax2
Definition: G4Cons.hh:220
G4double fRmin2
Definition: G4Cons.hh:220
G4double fSPhi
Definition: G4Cons.hh:220
G4double fDPhi
Definition: G4Cons.hh:220
T sqr(const T &x)
Definition: templates.hh:145
double G4double
Definition: G4Types.hh:76
G4double fRmax1
Definition: G4Cons.hh:220
Here is the call graph for this function:

◆ GetRmax1()

G4double G4Cons::GetRmax1 ( ) const
inline

◆ GetRmax2()

G4double G4Cons::GetRmax2 ( ) const
inline

◆ GetRmin1()

G4double G4Cons::GetRmin1 ( ) const
inline

◆ GetRmin2()

G4double G4Cons::GetRmin2 ( ) const
inline

◆ GetSPhi()

G4double G4Cons::GetSPhi ( ) const
inline

◆ GetStartPhiAngle()

G4double G4Cons::GetStartPhiAngle ( ) const
inline
Here is the caller graph for this function:

◆ GetSurfaceArea()

G4double G4Cons::GetSurfaceArea ( )
inlinevirtual

Reimplemented from G4VSolid.

◆ GetZHalfLength()

G4double G4Cons::GetZHalfLength ( ) const
inline
Here is the caller graph for this function:

◆ Initialize()

void G4Cons::Initialize ( )
inlineprivate

◆ InitializeTrigonometry()

void G4Cons::InitializeTrigonometry ( )
inlineprivate

◆ Inside()

EInside G4Cons::Inside ( const G4ThreeVector p) const
virtual

Implements G4VSolid.

Definition at line 205 of file G4Cons.cc.

206 {
207  G4double r2, rl, rh, pPhi, tolRMin, tolRMax; // rh2, rl2 ;
208  EInside in;
209 
210  if (std::fabs(p.z()) > fDz + halfCarTolerance ) { return in = kOutside; }
211  else if(std::fabs(p.z()) >= fDz - halfCarTolerance ) { in = kSurface; }
212  else { in = kInside; }
213 
214  r2 = p.x()*p.x() + p.y()*p.y() ;
215  rl = 0.5*(fRmin2*(p.z() + fDz) + fRmin1*(fDz - p.z()))/fDz ;
216  rh = 0.5*(fRmax2*(p.z()+fDz)+fRmax1*(fDz-p.z()))/fDz;
217 
218  // rh2 = rh*rh;
219 
220  tolRMin = rl - halfRadTolerance;
221  if ( tolRMin < 0 ) { tolRMin = 0; }
222  tolRMax = rh + halfRadTolerance;
223 
224  if ( (r2<tolRMin*tolRMin) || (r2>tolRMax*tolRMax) ) { return in = kOutside; }
225 
226  if (rl) { tolRMin = rl + halfRadTolerance; }
227  else { tolRMin = 0.0; }
228  tolRMax = rh - halfRadTolerance;
229 
230  if (in == kInside) // else it's kSurface already
231  {
232  if ( (r2 < tolRMin*tolRMin) || (r2 >= tolRMax*tolRMax) ) { in = kSurface; }
233  }
234  if ( !fPhiFullCone && ((p.x() != 0.0) || (p.y() != 0.0)) )
235  {
236  pPhi = std::atan2(p.y(),p.x()) ;
237 
238  if ( pPhi < fSPhi - halfAngTolerance ) { pPhi += twopi; }
239  else if ( pPhi > fSPhi + fDPhi + halfAngTolerance ) { pPhi -= twopi; }
240 
241  if ( (pPhi < fSPhi - halfAngTolerance) ||
242  (pPhi > fSPhi + fDPhi + halfAngTolerance) ) { return in = kOutside; }
243 
244  else if (in == kInside) // else it's kSurface anyway already
245  {
246  if ( (pPhi < fSPhi + halfAngTolerance) ||
247  (pPhi > fSPhi + fDPhi - halfAngTolerance) ) { in = kSurface; }
248  }
249  }
250  else if ( !fPhiFullCone ) { in = kSurface; }
251 
252  return in ;
253 }
G4double halfRadTolerance
Definition: G4Cons.hh:233
G4double fDz
Definition: G4Cons.hh:220
G4bool fPhiFullCone
Definition: G4Cons.hh:229
ifstream in
Definition: comparison.C:7
G4double fRmin1
Definition: G4Cons.hh:220
G4double fRmax2
Definition: G4Cons.hh:220
G4double fRmin2
Definition: G4Cons.hh:220
G4double fSPhi
Definition: G4Cons.hh:220
G4double halfCarTolerance
Definition: G4Cons.hh:233
static const double twopi
Definition: G4SIunits.hh:75
double x() const
G4double fDPhi
Definition: G4Cons.hh:220
double y() const
G4double halfAngTolerance
Definition: G4Cons.hh:233
EInside
Definition: geomdefs.hh:58
double z() const
double G4double
Definition: G4Types.hh:76
G4double fRmax1
Definition: G4Cons.hh:220
Here is the call graph for this function:
Here is the caller graph for this function:

◆ operator=()

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

Definition at line 172 of file G4Cons.cc.

173 {
174  // Check assignment to self
175  //
176  if (this == &rhs) { return *this; }
177 
178  // Copy base class data
179  //
181 
182  // Copy data
183  //
186  fRmin1 = rhs.fRmin1; fRmin2 = rhs.fRmin2;
187  fRmax1 = rhs.fRmax1; fRmax2 = rhs.fRmax2;
188  fDz = rhs.fDz; fSPhi = rhs.fSPhi; fDPhi = rhs.fDPhi;
189  sinCPhi = rhs.sinCPhi; cosCPhi = rhs.cosCPhi;
191  sinSPhi = rhs.sinSPhi; cosSPhi = rhs.cosSPhi;
192  sinEPhi = rhs.sinEPhi; cosEPhi = rhs.cosEPhi;
197 
198  return *this;
199 }
G4double halfRadTolerance
Definition: G4Cons.hh:233
G4double fDz
Definition: G4Cons.hh:220
G4bool fPhiFullCone
Definition: G4Cons.hh:229
G4double sinEPhi
Definition: G4Cons.hh:224
G4double fRmin1
Definition: G4Cons.hh:220
G4double fRmax2
Definition: G4Cons.hh:220
G4double cosCPhi
Definition: G4Cons.hh:224
G4double kRadTolerance
Definition: G4Cons.hh:216
G4double fRmin2
Definition: G4Cons.hh:220
G4double sinSPhi
Definition: G4Cons.hh:224
G4double fSPhi
Definition: G4Cons.hh:220
G4double halfCarTolerance
Definition: G4Cons.hh:233
G4double cosSPhi
Definition: G4Cons.hh:224
G4double fDPhi
Definition: G4Cons.hh:220
G4double halfAngTolerance
Definition: G4Cons.hh:233
G4double sinCPhi
Definition: G4Cons.hh:224
G4double kAngTolerance
Definition: G4Cons.hh:216
G4double cosEPhi
Definition: G4Cons.hh:224
G4double cosHDPhiOT
Definition: G4Cons.hh:224
G4double cosHDPhiIT
Definition: G4Cons.hh:224
G4CSGSolid & operator=(const G4CSGSolid &rhs)
Definition: G4CSGSolid.cc:91
G4double fRmax1
Definition: G4Cons.hh:220
Here is the call graph for this function:

◆ SetDeltaPhiAngle()

void G4Cons::SetDeltaPhiAngle ( G4double  newDPhi)
inline
Here is the caller graph for this function:

◆ SetInnerRadiusMinusZ()

void G4Cons::SetInnerRadiusMinusZ ( G4double  Rmin1)
inline
Here is the caller graph for this function:

◆ SetInnerRadiusPlusZ()

void G4Cons::SetInnerRadiusPlusZ ( G4double  Rmin2)
inline
Here is the caller graph for this function:

◆ SetOuterRadiusMinusZ()

void G4Cons::SetOuterRadiusMinusZ ( G4double  Rmax1)
inline
Here is the caller graph for this function:

◆ SetOuterRadiusPlusZ()

void G4Cons::SetOuterRadiusPlusZ ( G4double  Rmax2)
inline
Here is the caller graph for this function:

◆ SetStartPhiAngle()

void G4Cons::SetStartPhiAngle ( G4double  newSPhi,
G4bool  trig = true 
)
inline
Here is the caller graph for this function:

◆ SetZHalfLength()

void G4Cons::SetZHalfLength ( G4double  newDz)
inline
Here is the caller graph for this function:

◆ StreamInfo()

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

Reimplemented from G4CSGSolid.

Definition at line 2265 of file G4Cons.cc.

2266 {
2267  G4int oldprc = os.precision(16);
2268  os << "-----------------------------------------------------------\n"
2269  << " *** Dump for solid - " << GetName() << " ***\n"
2270  << " ===================================================\n"
2271  << " Solid type: G4Cons\n"
2272  << " Parameters: \n"
2273  << " inside -fDz radius: " << fRmin1/mm << " mm \n"
2274  << " outside -fDz radius: " << fRmax1/mm << " mm \n"
2275  << " inside +fDz radius: " << fRmin2/mm << " mm \n"
2276  << " outside +fDz radius: " << fRmax2/mm << " mm \n"
2277  << " half length in Z : " << fDz/mm << " mm \n"
2278  << " starting angle of segment: " << fSPhi/degree << " degrees \n"
2279  << " delta angle of segment : " << fDPhi/degree << " degrees \n"
2280  << "-----------------------------------------------------------\n";
2281  os.precision(oldprc);
2282 
2283  return os;
2284 }
G4double fDz
Definition: G4Cons.hh:220
G4double fRmin1
Definition: G4Cons.hh:220
int G4int
Definition: G4Types.hh:78
G4double fRmax2
Definition: G4Cons.hh:220
G4String GetName() const
G4double fRmin2
Definition: G4Cons.hh:220
G4double fSPhi
Definition: G4Cons.hh:220
G4double fDPhi
Definition: G4Cons.hh:220
static const double degree
Definition: G4SIunits.hh:143
static const double mm
Definition: G4SIunits.hh:114
G4double fRmax1
Definition: G4Cons.hh:220
Here is the call graph for this function:

◆ SurfaceNormal()

G4ThreeVector G4Cons::SurfaceNormal ( const G4ThreeVector p) const
virtual

Implements G4VSolid.

Definition at line 490 of file G4Cons.cc.

491 {
492  G4int noSurfaces = 0;
493  G4double rho, pPhi;
494  G4double distZ, distRMin, distRMax;
495  G4double distSPhi = kInfinity, distEPhi = kInfinity;
496  G4double tanRMin, secRMin, pRMin, widRMin;
497  G4double tanRMax, secRMax, pRMax, widRMax;
498 
499  G4ThreeVector norm, sumnorm(0.,0.,0.), nZ = G4ThreeVector(0.,0.,1.);
500  G4ThreeVector nR, nr(0.,0.,0.), nPs, nPe;
501 
502  distZ = std::fabs(std::fabs(p.z()) - fDz);
503  rho = std::sqrt(p.x()*p.x() + p.y()*p.y());
504 
505  tanRMin = (fRmin2 - fRmin1)*0.5/fDz;
506  secRMin = std::sqrt(1 + tanRMin*tanRMin);
507  pRMin = rho - p.z()*tanRMin;
508  widRMin = fRmin2 - fDz*tanRMin;
509  distRMin = std::fabs(pRMin - widRMin)/secRMin;
510 
511  tanRMax = (fRmax2 - fRmax1)*0.5/fDz;
512  secRMax = std::sqrt(1+tanRMax*tanRMax);
513  pRMax = rho - p.z()*tanRMax;
514  widRMax = fRmax2 - fDz*tanRMax;
515  distRMax = std::fabs(pRMax - widRMax)/secRMax;
516 
517  if (!fPhiFullCone) // Protected against (0,0,z)
518  {
519  if ( rho )
520  {
521  pPhi = std::atan2(p.y(),p.x());
522 
523  if (pPhi < fSPhi-halfCarTolerance) { pPhi += twopi; }
524  else if (pPhi > fSPhi+fDPhi+halfCarTolerance) { pPhi -= twopi; }
525 
526  distSPhi = std::fabs( pPhi - fSPhi );
527  distEPhi = std::fabs( pPhi - fSPhi - fDPhi );
528  }
529  else if( !(fRmin1) || !(fRmin2) )
530  {
531  distSPhi = 0.;
532  distEPhi = 0.;
533  }
534  nPs = G4ThreeVector(std::sin(fSPhi), -std::cos(fSPhi), 0);
535  nPe = G4ThreeVector(-std::sin(fSPhi+fDPhi), std::cos(fSPhi+fDPhi), 0);
536  }
537  if ( rho > halfCarTolerance )
538  {
539  nR = G4ThreeVector(p.x()/rho/secRMax, p.y()/rho/secRMax, -tanRMax/secRMax);
540  if (fRmin1 || fRmin2)
541  {
542  nr = G4ThreeVector(-p.x()/rho/secRMin,-p.y()/rho/secRMin,tanRMin/secRMin);
543  }
544  }
545 
546  if( distRMax <= halfCarTolerance )
547  {
548  noSurfaces ++;
549  sumnorm += nR;
550  }
551  if( (fRmin1 || fRmin2) && (distRMin <= halfCarTolerance) )
552  {
553  noSurfaces ++;
554  sumnorm += nr;
555  }
556  if( !fPhiFullCone )
557  {
558  if (distSPhi <= halfAngTolerance)
559  {
560  noSurfaces ++;
561  sumnorm += nPs;
562  }
563  if (distEPhi <= halfAngTolerance)
564  {
565  noSurfaces ++;
566  sumnorm += nPe;
567  }
568  }
569  if (distZ <= halfCarTolerance)
570  {
571  noSurfaces ++;
572  if ( p.z() >= 0.) { sumnorm += nZ; }
573  else { sumnorm -= nZ; }
574  }
575  if ( noSurfaces == 0 )
576  {
577 #ifdef G4CSGDEBUG
578  G4Exception("G4Cons::SurfaceNormal(p)", "GeomSolids1002",
579  JustWarning, "Point p is not on surface !?" );
580 #endif
581  norm = ApproxSurfaceNormal(p);
582  }
583  else if ( noSurfaces == 1 ) { norm = sumnorm; }
584  else { norm = sumnorm.unit(); }
585 
586  return norm ;
587 }
G4ThreeVector ApproxSurfaceNormal(const G4ThreeVector &p) const
Definition: G4Cons.cc:594
G4double fDz
Definition: G4Cons.hh:220
static const G4double kInfinity
Definition: geomdefs.hh:42
G4bool fPhiFullCone
Definition: G4Cons.hh:229
CLHEP::Hep3Vector G4ThreeVector
Float_t norm
G4double fRmin1
Definition: G4Cons.hh:220
int G4int
Definition: G4Types.hh:78
G4double fRmax2
Definition: G4Cons.hh:220
G4double fRmin2
Definition: G4Cons.hh:220
G4double fSPhi
Definition: G4Cons.hh:220
G4double halfCarTolerance
Definition: G4Cons.hh:233
Hep3Vector unit() const
static const double twopi
Definition: G4SIunits.hh:75
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
double x() const
G4double fDPhi
Definition: G4Cons.hh:220
double y() const
G4double halfAngTolerance
Definition: G4Cons.hh:233
double z() const
double G4double
Definition: G4Types.hh:76
G4double fRmax1
Definition: G4Cons.hh:220
Here is the call graph for this function:

Member Data Documentation

◆ cosCPhi

G4double G4Cons::cosCPhi
private

Definition at line 224 of file G4Cons.hh.

◆ cosEPhi

G4double G4Cons::cosEPhi
private

Definition at line 224 of file G4Cons.hh.

◆ cosHDPhiIT

G4double G4Cons::cosHDPhiIT
private

Definition at line 224 of file G4Cons.hh.

◆ cosHDPhiOT

G4double G4Cons::cosHDPhiOT
private

Definition at line 224 of file G4Cons.hh.

◆ cosSPhi

G4double G4Cons::cosSPhi
private

Definition at line 224 of file G4Cons.hh.

◆ fDPhi

G4double G4Cons::fDPhi
private

Definition at line 220 of file G4Cons.hh.

◆ fDz

G4double G4Cons::fDz
private

Definition at line 220 of file G4Cons.hh.

◆ fPhiFullCone

G4bool G4Cons::fPhiFullCone
private

Definition at line 229 of file G4Cons.hh.

◆ fRmax1

G4double G4Cons::fRmax1
private

Definition at line 220 of file G4Cons.hh.

◆ fRmax2

G4double G4Cons::fRmax2
private

Definition at line 220 of file G4Cons.hh.

◆ fRmin1

G4double G4Cons::fRmin1
private

Definition at line 220 of file G4Cons.hh.

◆ fRmin2

G4double G4Cons::fRmin2
private

Definition at line 220 of file G4Cons.hh.

◆ fSPhi

G4double G4Cons::fSPhi
private

Definition at line 220 of file G4Cons.hh.

◆ halfAngTolerance

G4double G4Cons::halfAngTolerance
private

Definition at line 233 of file G4Cons.hh.

◆ halfCarTolerance

G4double G4Cons::halfCarTolerance
private

Definition at line 233 of file G4Cons.hh.

◆ halfRadTolerance

G4double G4Cons::halfRadTolerance
private

Definition at line 233 of file G4Cons.hh.

◆ kAngTolerance

G4double G4Cons::kAngTolerance
private

Definition at line 216 of file G4Cons.hh.

◆ kRadTolerance

G4double G4Cons::kRadTolerance
private

Definition at line 216 of file G4Cons.hh.

◆ sinCPhi

G4double G4Cons::sinCPhi
private

Definition at line 224 of file G4Cons.hh.

◆ sinEPhi

G4double G4Cons::sinEPhi
private

Definition at line 224 of file G4Cons.hh.

◆ sinSPhi

G4double G4Cons::sinSPhi
private

Definition at line 224 of file G4Cons.hh.


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