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

#include <G4Cons.hh>

+ Inheritance 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 () override
 
G4double GetInnerRadiusMinusZ () const
 
G4double GetOuterRadiusMinusZ () const
 
G4double GetInnerRadiusPlusZ () const
 
G4double GetOuterRadiusPlusZ () const
 
G4double GetZHalfLength () const
 
G4double GetStartPhiAngle () const
 
G4double GetDeltaPhiAngle () const
 
G4double GetSinStartPhi () const
 
G4double GetCosStartPhi () const
 
G4double GetSinEndPhi () const
 
G4double GetCosEndPhi () 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 () override
 
G4double GetSurfaceArea () override
 
void ComputeDimensions (G4VPVParameterisation *p, const G4int n, const G4VPhysicalVolume *pRep) override
 
void BoundingLimits (G4ThreeVector &pMin, G4ThreeVector &pMax) const override
 
G4bool CalculateExtent (const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pMin, G4double &pMax) const override
 
EInside Inside (const G4ThreeVector &p) const override
 
G4ThreeVector SurfaceNormal (const G4ThreeVector &p) const override
 
G4double DistanceToIn (const G4ThreeVector &p, const G4ThreeVector &v) const override
 
G4double DistanceToIn (const G4ThreeVector &p) const override
 
G4double DistanceToOut (const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=false, G4bool *validNorm=nullptr, G4ThreeVector *n=nullptr) const override
 
G4double DistanceToOut (const G4ThreeVector &p) const override
 
G4GeometryType GetEntityType () const override
 
G4ThreeVector GetPointOnSurface () const override
 
G4VSolid * Clone () const override
 
std::ostream & StreamInfo (std::ostream &os) const override
 
void DescribeYourselfTo (G4VGraphicsScene &scene) const override
 
G4Polyhedron * CreatePolyhedron () const override
 
 G4Cons (__void__ &)
 
 G4Cons (const G4Cons &rhs)
 
G4Cons & operator= (const G4Cons &rhs)
 
- Public Member Functions inherited from G4CSGSolid
 G4CSGSolid (const G4String &pName)
 
virtual ~G4CSGSolid ()
 
std::ostream & StreamInfo (std::ostream &os) const override
 
G4Polyhedron * GetPolyhedron () const override
 
 G4CSGSolid (__void__ &)
 
 G4CSGSolid (const G4CSGSolid &rhs)
 
G4CSGSolid & operator= (const G4CSGSolid &rhs)
 
- Public Member Functions inherited from G4VSolid
 G4VSolid (const G4String &name)
 
virtual ~G4VSolid ()
 
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 G4VSolid * GetConstituentSolid (G4int no) const
 
virtual G4VSolid * GetConstituentSolid (G4int no)
 
virtual const G4DisplacedSolid * GetDisplacedSolidPtr () const
 
virtual G4DisplacedSolid * GetDisplacedSolidPtr ()
 
 G4VSolid (__void__ &)
 
 G4VSolid (const G4VSolid &rhs)
 
G4VSolid & operator= (const G4VSolid &rhs)
 
G4double EstimateCubicVolume (G4int nStat, G4double epsilon) const
 
G4double EstimateSurfaceArea (G4int nStat, G4double ell) const
 

Additional Inherited Members

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

Detailed Description

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

82 : G4CSGSolid(pName), fRmin1(pRmin1), fRmin2(pRmin2),
83 fRmax1(pRmax1), fRmax2(pRmax2), fDz(pDz), fSPhi(0.), fDPhi(0.)
84{
85 kRadTolerance = G4GeometryTolerance::GetInstance()->GetRadialTolerance();
86 kAngTolerance = G4GeometryTolerance::GetInstance()->GetAngularTolerance();
87
88 halfCarTolerance=kCarTolerance*0.5;
89 halfRadTolerance=kRadTolerance*0.5;
90 halfAngTolerance=kAngTolerance*0.5;
91
92 // Check z-len
93 //
94 if ( pDz < 0 )
95 {
96 std::ostringstream message;
97 message << "Invalid Z half-length for Solid: " << GetName() << G4endl
98 << " hZ = " << pDz;
99 G4Exception("G4Cons::G4Cons()", "GeomSolids0002",
100 FatalException, message);
101 }
102
103 // Check radii
104 //
105 if (((pRmin1>=pRmax1) || (pRmin2>=pRmax2) || (pRmin1<0)) && (pRmin2<0))
106 {
107 std::ostringstream message;
108 message << "Invalid values of radii for Solid: " << GetName() << G4endl
109 << " pRmin1 = " << pRmin1 << ", pRmin2 = " << pRmin2
110 << ", pRmax1 = " << pRmax1 << ", pRmax2 = " << pRmax2;
111 G4Exception("G4Cons::G4Cons()", "GeomSolids0002",
112 FatalException, message) ;
113 }
114 if( (pRmin1 == 0.0) && (pRmin2 > 0.0) ) { fRmin1 = 1e3*kRadTolerance ; }
115 if( (pRmin2 == 0.0) && (pRmin1 > 0.0) ) { fRmin2 = 1e3*kRadTolerance ; }
116
117 // Check angles
118 //
119 CheckPhiAngles(pSPhi, pDPhi);
120}
FatalException
@ FatalException
Definition G4ExceptionSeverity.hh:67
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition G4Exception.cc:59
G4endl
#define G4endl
Definition G4ios.hh:67
G4CSGSolid::G4CSGSolid
G4CSGSolid(const G4String &pName)
Definition G4CSGSolid.cc:49
G4GeometryTolerance::GetRadialTolerance
G4double GetRadialTolerance() const
Definition G4GeometryTolerance.cc:82
G4GeometryTolerance::GetInstance
static G4GeometryTolerance * GetInstance()
Definition G4GeometryTolerance.cc:58
G4GeometryTolerance::GetAngularTolerance
G4double GetAngularTolerance() const
Definition G4GeometryTolerance.cc:77
G4VSolid::GetName
G4String GetName() const
G4VSolid::kCarTolerance
G4double kCarTolerance
Definition G4VSolid.hh:299

Referenced by Clone().

◆ ~G4Cons()

G4Cons::~G4Cons ( )
overridedefault

◆ G4Cons() [2/3]

G4Cons::G4Cons ( __void__ & a)

Definition at line 127 of file G4Cons.cc.

128 : G4CSGSolid(a)
129{
130}

◆ G4Cons() [3/3]

G4Cons::G4Cons ( const G4Cons & rhs)
default

Member Function Documentation

◆ BoundingLimits()

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

Reimplemented from G4VSolid.

Definition at line 247 of file G4Cons.cc.

248{
249 G4double rmin = std::min(GetInnerRadiusMinusZ(),GetInnerRadiusPlusZ());
250 G4double rmax = std::max(GetOuterRadiusMinusZ(),GetOuterRadiusPlusZ());
251 G4double dz = GetZHalfLength();
252
253 // Find bounding box
254 //
255 if (GetDeltaPhiAngle() < twopi)
256 {
257 G4TwoVector vmin,vmax;
258 G4GeomTools::DiskExtent(rmin,rmax,
259 GetSinStartPhi(),GetCosStartPhi(),
260 GetSinEndPhi(),GetCosEndPhi(),
261 vmin,vmax);
262 pMin.set(vmin.x(),vmin.y(),-dz);
263 pMax.set(vmax.x(),vmax.y(), dz);
264 }
265 else
266 {
267 pMin.set(-rmax,-rmax,-dz);
268 pMax.set( rmax, rmax, dz);
269 }
270
271 // Check correctness of the bounding box
272 //
273 if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
274 {
275 std::ostringstream message;
276 message << "Bad bounding box (min >= max) for solid: "
277 << GetName() << " !"
278 << "\npMin = " << pMin
279 << "\npMax = " << pMax;
280 G4Exception("G4Cons::BoundingLimits()", "GeomMgt0001",
281 JustWarning, message);
282 DumpInfo();
283 }
284}
JustWarning
@ JustWarning
Definition G4ExceptionSeverity.hh:71
G4double
double G4double
Definition G4Types.hh:83
CLHEP::Hep2Vector::x
double x() const
CLHEP::Hep2Vector::y
double y() const
CLHEP::Hep3Vector::z
double z() const
CLHEP::Hep3Vector::x
double x() const
CLHEP::Hep3Vector::y
double y() const
CLHEP::Hep3Vector::set
void set(double x, double y, double z)
G4Cons::GetOuterRadiusPlusZ
G4double GetOuterRadiusPlusZ() const
G4Cons::GetCosStartPhi
G4double GetCosStartPhi() const
G4Cons::GetDeltaPhiAngle
G4double GetDeltaPhiAngle() const
G4Cons::GetInnerRadiusMinusZ
G4double GetInnerRadiusMinusZ() const
G4Cons::GetInnerRadiusPlusZ
G4double GetInnerRadiusPlusZ() const
G4Cons::GetCosEndPhi
G4double GetCosEndPhi() const
G4Cons::GetOuterRadiusMinusZ
G4double GetOuterRadiusMinusZ() const
G4Cons::GetSinEndPhi
G4double GetSinEndPhi() const
G4Cons::GetSinStartPhi
G4double GetSinStartPhi() const
G4Cons::GetZHalfLength
G4double GetZHalfLength() const
G4GeomTools::DiskExtent
static G4bool DiskExtent(G4double rmin, G4double rmax, G4double startPhi, G4double delPhi, G4TwoVector &pmin, G4TwoVector &pmax)
Definition G4GeomTools.cc:390
G4VSolid::DumpInfo
void DumpInfo() const

Referenced by CalculateExtent().

◆ CalculateExtent()

G4bool G4Cons::CalculateExtent ( const EAxis pAxis,
const G4VoxelLimits & pVoxelLimit,
const G4AffineTransform & pTransform,
G4double & pMin,
G4double & pMax ) const
overridevirtual

Implements G4VSolid.

Definition at line 290 of file G4Cons.cc.

295{
296 G4ThreeVector bmin, bmax;
297 G4bool exist;
298
299 // Get bounding box
300 BoundingLimits(bmin,bmax);
301
302 // Check bounding box
303 G4BoundingEnvelope bbox(bmin,bmax);
304#ifdef G4BBOX_EXTENT
305 if (true) return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
306#endif
307 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVoxelLimit,pTransform,pMin,pMax))
308 {
309 return exist = pMin < pMax;
310 }
311
312 // Get parameters of the solid
313 G4double rmin1 = GetInnerRadiusMinusZ();
314 G4double rmax1 = GetOuterRadiusMinusZ();
315 G4double rmin2 = GetInnerRadiusPlusZ();
316 G4double rmax2 = GetOuterRadiusPlusZ();
317 G4double dz = GetZHalfLength();
318 G4double dphi = GetDeltaPhiAngle();
319
320 // Find bounding envelope and calculate extent
321 //
322 const G4int NSTEPS = 24; // number of steps for whole circle
323 G4double astep = twopi/NSTEPS; // max angle for one step
324 G4int ksteps = (dphi <= astep) ? 1 : (G4int)((dphi-deg)/astep) + 1;
325 G4double ang = dphi/ksteps;
326
327 G4double sinHalf = std::sin(0.5*ang);
328 G4double cosHalf = std::cos(0.5*ang);
329 G4double sinStep = 2.*sinHalf*cosHalf;
330 G4double cosStep = 1. - 2.*sinHalf*sinHalf;
331 G4double rext1 = rmax1/cosHalf;
332 G4double rext2 = rmax2/cosHalf;
333
334 // bounding envelope for full cone without hole consists of two polygons,
335 // in other cases it is a sequence of quadrilaterals
336 if (rmin1 == 0 && rmin2 == 0 && dphi == twopi)
337 {
338 G4double sinCur = sinHalf;
339 G4double cosCur = cosHalf;
340
341 G4ThreeVectorList baseA(NSTEPS),baseB(NSTEPS);
342 for (G4int k=0; k<NSTEPS; ++k)
343 {
344 baseA[k].set(rext1*cosCur,rext1*sinCur,-dz);
345 baseB[k].set(rext2*cosCur,rext2*sinCur, dz);
346
347 G4double sinTmp = sinCur;
348 sinCur = sinCur*cosStep + cosCur*sinStep;
349 cosCur = cosCur*cosStep - sinTmp*sinStep;
350 }
351 std::vector<const G4ThreeVectorList *> polygons(2);
352 polygons[0] = &baseA;
353 polygons[1] = &baseB;
354 G4BoundingEnvelope benv(bmin,bmax,polygons);
355 exist = benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
356 }
357 else
358 {
359 G4double sinStart = GetSinStartPhi();
360 G4double cosStart = GetCosStartPhi();
361 G4double sinEnd = GetSinEndPhi();
362 G4double cosEnd = GetCosEndPhi();
363 G4double sinCur = sinStart*cosHalf + cosStart*sinHalf;
364 G4double cosCur = cosStart*cosHalf - sinStart*sinHalf;
365
366 // set quadrilaterals
367 G4ThreeVectorList pols[NSTEPS+2];
368 for (G4int k=0; k<ksteps+2; ++k) pols[k].resize(4);
369 pols[0][0].set(rmin2*cosStart,rmin2*sinStart, dz);
370 pols[0][1].set(rmin1*cosStart,rmin1*sinStart,-dz);
371 pols[0][2].set(rmax1*cosStart,rmax1*sinStart,-dz);
372 pols[0][3].set(rmax2*cosStart,rmax2*sinStart, dz);
373 for (G4int k=1; k<ksteps+1; ++k)
374 {
375 pols[k][0].set(rmin2*cosCur,rmin2*sinCur, dz);
376 pols[k][1].set(rmin1*cosCur,rmin1*sinCur,-dz);
377 pols[k][2].set(rext1*cosCur,rext1*sinCur,-dz);
378 pols[k][3].set(rext2*cosCur,rext2*sinCur, dz);
379
380 G4double sinTmp = sinCur;
381 sinCur = sinCur*cosStep + cosCur*sinStep;
382 cosCur = cosCur*cosStep - sinTmp*sinStep;
383 }
384 pols[ksteps+1][0].set(rmin2*cosEnd,rmin2*sinEnd, dz);
385 pols[ksteps+1][1].set(rmin1*cosEnd,rmin1*sinEnd,-dz);
386 pols[ksteps+1][2].set(rmax1*cosEnd,rmax1*sinEnd,-dz);
387 pols[ksteps+1][3].set(rmax2*cosEnd,rmax2*sinEnd, dz);
388
389 // set envelope and calculate extent
390 std::vector<const G4ThreeVectorList *> polygons;
391 polygons.resize(ksteps+2);
392 for (G4int k=0; k<ksteps+2; ++k) polygons[k] = &pols[k];
393 G4BoundingEnvelope benv(bmin,bmax,polygons);
394 exist = benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
395 }
396 return exist;
397}
G4ThreeVectorList
std::vector< G4ThreeVector > G4ThreeVectorList
Definition G4BoundingEnvelope.hh:54
G4bool
bool G4bool
Definition G4Types.hh:86
G4int
int G4int
Definition G4Types.hh:85
G4Cons::BoundingLimits
void BoundingLimits(G4ThreeVector &pMin, G4ThreeVector &pMax) const override
Definition G4Cons.cc:247

◆ Clone()

G4VSolid * G4Cons::Clone ( ) const
overridevirtual

Reimplemented from G4VSolid.

Definition at line 2089 of file G4Cons.cc.

2090{
2091 return new G4Cons(*this);
2092}
G4Cons::G4Cons
G4Cons(const G4String &pName, G4double pRmin1, G4double pRmax1, G4double pRmin2, G4double pRmax2, G4double pDz, G4double pSPhi, G4double pDPhi)
Definition G4Cons.cc:77

◆ ComputeDimensions()

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

Reimplemented from G4VSolid.

Definition at line 236 of file G4Cons.cc.

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

◆ CreatePolyhedron()

G4Polyhedron * G4Cons::CreatePolyhedron ( ) const
overridevirtual

Reimplemented from G4VSolid.

Definition at line 2210 of file G4Cons.cc.

2211{
2212 return new G4PolyhedronCons(fRmin1,fRmax1,fRmin2,fRmax2,fDz,fSPhi,fDPhi);
2213}

◆ DescribeYourselfTo()

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

Implements G4VSolid.

Definition at line 2205 of file G4Cons.cc.

2206{
2207 scene.AddSolid (*this);
2208}
G4VGraphicsScene::AddSolid
virtual void AddSolid(const G4Box &)=0

◆ DistanceToIn() [1/2]

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

Implements G4VSolid.

Definition at line 1317 of file G4Cons.cc.

1318{
1319 G4double safe=0.0, rho, safeR1, safeR2, safeZ, safePhi, cosPsi ;
1320 G4double tanRMin, secRMin, pRMin ;
1321 G4double tanRMax, secRMax, pRMax ;
1322
1323 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
1324 safeZ = std::fabs(p.z()) - fDz ;
1325
1326 if ( (fRmin1 != 0.0) || (fRmin2 != 0.0) )
1327 {
1328 tanRMin = (fRmin2 - fRmin1)*0.5/fDz ;
1329 secRMin = std::sqrt(1.0 + tanRMin*tanRMin) ;
1330 pRMin = tanRMin*p.z() + (fRmin1 + fRmin2)*0.5 ;
1331 safeR1 = (pRMin - rho)/secRMin ;
1332
1333 tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;
1334 secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
1335 pRMax = tanRMax*p.z() + (fRmax1 + fRmax2)*0.5 ;
1336 safeR2 = (rho - pRMax)/secRMax ;
1337
1338 if ( safeR1 > safeR2) { safe = safeR1; }
1339 else { safe = safeR2; }
1340 }
1341 else
1342 {
1343 tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;
1344 secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
1345 pRMax = tanRMax*p.z() + (fRmax1 + fRmax2)*0.5 ;
1346 safe = (rho - pRMax)/secRMax ;
1347 }
1348 if ( safeZ > safe ) { safe = safeZ; }
1349
1350 if ( !fPhiFullCone && (rho != 0.0) )
1351 {
1352 // Psi=angle from central phi to point
1353
1354 cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/rho ;
1355
1356 if ( cosPsi < cosHDPhi ) // Point lies outside phi range
1357 {
1358 if ( (p.y()*cosCPhi - p.x()*sinCPhi) <= 0.0 )
1359 {
1360 safePhi = std::fabs(p.x()*sinSPhi-p.y()*cosSPhi);
1361 }
1362 else
1363 {
1364 safePhi = std::fabs(p.x()*sinEPhi-p.y()*cosEPhi);
1365 }
1366 if ( safePhi > safe ) { safe = safePhi; }
1367 }
1368 }
1369 if ( safe < 0.0 ) { safe = 0.0; }
1370
1371 return safe ;
1372}

◆ DistanceToIn() [2/2]

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

Implements G4VSolid.

Definition at line 647 of file G4Cons.cc.

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

Referenced by DistanceToIn().

◆ DistanceToOut() [1/2]

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

Implements G4VSolid.

Definition at line 2002 of file G4Cons.cc.

2003{
2004 G4double safe=0.0, rho, safeR1, safeR2, safeZ, safePhi;
2005 G4double tanRMin, secRMin, pRMin;
2006 G4double tanRMax, secRMax, pRMax;
2007
2008#ifdef G4CSGDEBUG
2009 if( Inside(p) == kOutside )
2010 {
2011 G4int oldprc=G4cout.precision(16) ;
2012 G4cout << G4endl ;
2013 DumpInfo();
2014 G4cout << "Position:" << G4endl << G4endl ;
2015 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl ;
2016 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ;
2017 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ;
2018 G4cout << "pho at z = " << std::sqrt( p.x()*p.x()+p.y()*p.y() )/mm
2019 << " mm" << G4endl << G4endl ;
2020 if( (p.x() != 0.) || (p.x() != 0.) )
2021 {
2022 G4cout << "point phi = " << std::atan2(p.y(),p.x())/degree
2023 << " degrees" << G4endl << G4endl ;
2024 }
2025 G4cout.precision(oldprc) ;
2026 G4Exception("G4Cons::DistanceToOut(p)", "GeomSolids1002",
2027 JustWarning, "Point p is outside !?" );
2028 }
2029#endif
2030
2031 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
2032 safeZ = fDz - std::fabs(p.z()) ;
2033
2034 if ((fRmin1 != 0.0) || (fRmin2 != 0.0))
2035 {
2036 tanRMin = (fRmin2 - fRmin1)*0.5/fDz ;
2037 secRMin = std::sqrt(1.0 + tanRMin*tanRMin) ;
2038 pRMin = tanRMin*p.z() + (fRmin1 + fRmin2)*0.5 ;
2039 safeR1 = (rho - pRMin)/secRMin ;
2040 }
2041 else
2042 {
2043 safeR1 = kInfinity ;
2044 }
2045
2046 tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;
2047 secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
2048 pRMax = tanRMax*p.z() + (fRmax1+fRmax2)*0.5 ;
2049 safeR2 = (pRMax - rho)/secRMax ;
2050
2051 if (safeR1 < safeR2) { safe = safeR1; }
2052 else { safe = safeR2; }
2053 if (safeZ < safe) { safe = safeZ ; }
2054
2055 // Check if phi divided, Calc distances closest phi plane
2056
2057 if (!fPhiFullCone)
2058 {
2059 // Above/below central phi of G4Cons?
2060
2061 if ( (p.y()*cosCPhi - p.x()*sinCPhi) <= 0 )
2062 {
2063 safePhi = -(p.x()*sinSPhi - p.y()*cosSPhi) ;
2064 }
2065 else
2066 {
2067 safePhi = (p.x()*sinEPhi - p.y()*cosEPhi) ;
2068 }
2069 if (safePhi < safe) { safe = safePhi; }
2070 }
2071 if ( safe < 0 ) { safe = 0; }
2072
2073 return safe ;
2074}
G4cout
G4GLOB_DLL std::ostream G4cout
G4Cons::Inside
EInside Inside(const G4ThreeVector &p) const override
Definition G4Cons.cc:181
kOutside
@ kOutside
Definition geomdefs.hh:68

◆ DistanceToOut() [2/2]

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

Implements G4VSolid.

Definition at line 1379 of file G4Cons.cc.

1384{
1385 ESide side = kNull, sider = kNull, sidephi = kNull;
1386
1387 G4double snxt,srd,sphi,pdist ;
1388
1389 G4double tanRMax, secRMax, rMaxAv ; // Data for outer cone
1390 G4double tanRMin, secRMin, rMinAv ; // Data for inner cone
1391
1392 G4double t1, t2, t3, rout, rin, nt1, nt2, nt3 ;
1393 G4double b, c, d, sr2, sr3 ;
1394
1395 // Vars for intersection within tolerance
1396
1397 ESide sidetol = kNull ;
1398 G4double slentol = kInfinity ;
1399
1400 // Vars for phi intersection:
1401
1402 G4double pDistS, compS, pDistE, compE, sphi2, xi, yi, risec, vphi ;
1403 G4double zi, ri, deltaRoi2 ;
1404
1405 // Z plane intersection
1406
1407 if ( v.z() > 0.0 )
1408 {
1409 pdist = fDz - p.z() ;
1410
1411 if (pdist > halfCarTolerance)
1412 {
1413 snxt = pdist/v.z() ;
1414 side = kPZ ;
1415 }
1416 else
1417 {
1418 if (calcNorm)
1419 {
1420 *n = G4ThreeVector(0,0,1) ;
1421 *validNorm = true ;
1422 }
1423 return snxt = 0.0;
1424 }
1425 }
1426 else if ( v.z() < 0.0 )
1427 {
1428 pdist = fDz + p.z() ;
1429
1430 if ( pdist > halfCarTolerance)
1431 {
1432 snxt = -pdist/v.z() ;
1433 side = kMZ ;
1434 }
1435 else
1436 {
1437 if ( calcNorm )
1438 {
1439 *n = G4ThreeVector(0,0,-1) ;
1440 *validNorm = true ;
1441 }
1442 return snxt = 0.0 ;
1443 }
1444 }
1445 else // Travel perpendicular to z axis
1446 {
1447 snxt = kInfinity ;
1448 side = kNull ;
1449 }
1450
1451 // Radial Intersections
1452 //
1453 // Intersection with outer cone (possible return) and
1454 // inner cone (must also check phi)
1455 //
1456 // Intersection point (xi,yi,zi) on line x=p.x+t*v.x etc.
1457 //
1458 // Intersects with x^2+y^2=(a*z+b)^2
1459 //
1460 // where a=tanRMax or tanRMin
1461 // b=rMaxAv or rMinAv
1462 //
1463 // (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 ;
1464 // t1 t2 t3
1465 //
1466 // \--------u-------/ \-----------v----------/ \---------w--------/
1467
1468 tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;
1469 secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
1470 rMaxAv = (fRmax1 + fRmax2)*0.5 ;
1471
1472
1473 t1 = 1.0 - v.z()*v.z() ; // since v normalised
1474 t2 = p.x()*v.x() + p.y()*v.y() ;
1475 t3 = p.x()*p.x() + p.y()*p.y() ;
1476 rout = tanRMax*p.z() + rMaxAv ;
1477
1478 nt1 = t1 - (tanRMax*v.z())*(tanRMax*v.z()) ;
1479 nt2 = t2 - tanRMax*v.z()*rout ;
1480 nt3 = t3 - rout*rout ;
1481
1482 if (v.z() > 0.0)
1483 {
1484 deltaRoi2 = snxt*snxt*t1 + 2*snxt*t2 + t3
1485 - fRmax2*(fRmax2 + kRadTolerance*secRMax);
1486 }
1487 else if (v.z() < 0.0)
1488 {
1489 deltaRoi2 = snxt*snxt*t1 + 2*snxt*t2 + t3
1490 - fRmax1*(fRmax1 + kRadTolerance*secRMax);
1491 }
1492 else
1493 {
1494 deltaRoi2 = 1.0;
1495 }
1496
1497 if ( (nt1 != 0.0) && (deltaRoi2 > 0.0) )
1498 {
1499 // Equation quadratic => 2 roots : second root must be leaving
1500
1501 b = nt2/nt1 ;
1502 c = nt3/nt1 ;
1503 d = b*b - c ;
1504
1505 if ( d >= 0 )
1506 {
1507 // Check if on outer cone & heading outwards
1508 // NOTE: Should use rho-rout>-kRadTolerance*0.5
1509
1510 if (nt3 > -halfRadTolerance && nt2 >= 0 )
1511 {
1512 if (calcNorm)
1513 {
1514 risec = std::sqrt(t3)*secRMax ;
1515 *validNorm = true ;
1516 *n = G4ThreeVector(p.x()/risec,p.y()/risec,-tanRMax/secRMax);
1517 }
1518 return snxt=0 ;
1519 }
1520 else
1521 {
1522 sider = kRMax ;
1523 if (b>0) { srd = -b - std::sqrt(d); }
1524 else { srd = c/(-b+std::sqrt(d)) ; }
1525
1526 zi = p.z() + srd*v.z() ;
1527 ri = tanRMax*zi + rMaxAv ;
1528
1529 if ((ri >= 0) && (-halfRadTolerance <= srd) && (srd <= halfRadTolerance))
1530 {
1531 // An intersection within the tolerance
1532 // we will Store it in case it is good -
1533 //
1534 slentol = srd ;
1535 sidetol = kRMax ;
1536 }
1537 if ( (ri < 0) || (srd < halfRadTolerance) )
1538 {
1539 // Safety: if both roots -ve ensure that srd cannot `win'
1540 // distance to out
1541
1542 if (b>0) { sr2 = c/(-b-std::sqrt(d)); }
1543 else { sr2 = -b + std::sqrt(d); }
1544 zi = p.z() + sr2*v.z() ;
1545 ri = tanRMax*zi + rMaxAv ;
1546
1547 if ((ri >= 0) && (sr2 > halfRadTolerance))
1548 {
1549 srd = sr2;
1550 }
1551 else
1552 {
1553 srd = kInfinity ;
1554
1555 if( (-halfRadTolerance <= sr2) && ( sr2 <= halfRadTolerance) )
1556 {
1557 // An intersection within the tolerance.
1558 // Storing it in case it is good.
1559
1560 slentol = sr2 ;
1561 sidetol = kRMax ;
1562 }
1563 }
1564 }
1565 }
1566 }
1567 else
1568 {
1569 // No intersection with outer cone & not parallel
1570 // -> already outside, no intersection
1571
1572 if ( calcNorm )
1573 {
1574 risec = std::sqrt(t3)*secRMax;
1575 *validNorm = true;
1576 *n = G4ThreeVector(p.x()/risec,p.y()/risec,-tanRMax/secRMax);
1577 }
1578 return snxt = 0.0 ;
1579 }
1580 }
1581 else if ( (nt2 != 0.0) && (deltaRoi2 > 0.0) )
1582 {
1583 // Linear case (only one intersection) => point outside outer cone
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.0 ;
1592 }
1593 else
1594 {
1595 // No intersection -> parallel to outer cone
1596 // => Z or inner cone intersection
1597
1598 srd = kInfinity ;
1599 }
1600
1601 // Check possible intersection within tolerance
1602
1603 if ( slentol <= halfCarTolerance )
1604 {
1605 // An intersection within the tolerance was found.
1606 // We must accept it only if the momentum points outwards.
1607 //
1608 // G4ThreeVector ptTol ; // The point of the intersection
1609 // ptTol= p + slentol*v ;
1610 // ri=tanRMax*zi+rMaxAv ;
1611 //
1612 // Calculate a normal vector, as below
1613
1614 xi = p.x() + slentol*v.x();
1615 yi = p.y() + slentol*v.y();
1616 risec = std::sqrt(xi*xi + yi*yi)*secRMax;
1617 G4ThreeVector Normal = G4ThreeVector(xi/risec,yi/risec,-tanRMax/secRMax);
1618
1619 if ( Normal.dot(v) > 0 ) // We will leave the Cone immediatelly
1620 {
1621 if ( calcNorm )
1622 {
1623 *n = Normal.unit() ;
1624 *validNorm = true ;
1625 }
1626 return snxt = 0.0 ;
1627 }
1628 else // On the surface, but not heading out so we ignore this intersection
1629 { // (as it is within tolerance).
1630 slentol = kInfinity ;
1631 }
1632 }
1633
1634 // Inner Cone intersection
1635
1636 if ( (fRmin1 != 0.0) || (fRmin2 != 0.0) )
1637 {
1638 tanRMin = (fRmin2 - fRmin1)*0.5/fDz ;
1639 nt1 = t1 - (tanRMin*v.z())*(tanRMin*v.z()) ;
1640
1641 if ( nt1 != 0.0 )
1642 {
1643 secRMin = std::sqrt(1.0 + tanRMin*tanRMin) ;
1644 rMinAv = (fRmin1 + fRmin2)*0.5 ;
1645 rin = tanRMin*p.z() + rMinAv ;
1646 nt2 = t2 - tanRMin*v.z()*rin ;
1647 nt3 = t3 - rin*rin ;
1648
1649 // Equation quadratic => 2 roots : first root must be leaving
1650
1651 b = nt2/nt1 ;
1652 c = nt3/nt1 ;
1653 d = b*b - c ;
1654
1655 if ( d >= 0.0 )
1656 {
1657 // NOTE: should be rho-rin<kRadTolerance*0.5,
1658 // but using squared versions for efficiency
1659
1660 if (nt3 < kRadTolerance*(rin + kRadTolerance*0.25))
1661 {
1662 if ( nt2 < 0.0 )
1663 {
1664 if (calcNorm) { *validNorm = false; }
1665 return snxt = 0.0;
1666 }
1667 }
1668 else
1669 {
1670 if (b>0) { sr2 = -b - std::sqrt(d); }
1671 else { sr2 = c/(-b+std::sqrt(d)); }
1672 zi = p.z() + sr2*v.z() ;
1673 ri = tanRMin*zi + rMinAv ;
1674
1675 if( (ri>=0.0)&&(-halfRadTolerance<=sr2)&&(sr2<=halfRadTolerance) )
1676 {
1677 // An intersection within the tolerance
1678 // storing it in case it is good.
1679
1680 slentol = sr2 ;
1681 sidetol = kRMax ;
1682 }
1683 if( (ri<0) || (sr2 < halfRadTolerance) )
1684 {
1685 if (b>0) { sr3 = c/(-b-std::sqrt(d)); }
1686 else { sr3 = -b + std::sqrt(d) ; }
1687
1688 // Safety: if both roots -ve ensure that srd cannot `win'
1689 // distancetoout
1690
1691 if ( sr3 > halfRadTolerance )
1692 {
1693 if( sr3 < srd )
1694 {
1695 zi = p.z() + sr3*v.z() ;
1696 ri = tanRMin*zi + rMinAv ;
1697
1698 if ( ri >= 0.0 )
1699 {
1700 srd=sr3 ;
1701 sider=kRMin ;
1702 }
1703 }
1704 }
1705 else if ( sr3 > -halfRadTolerance )
1706 {
1707 // Intersection in tolerance. Store to check if it's good
1708
1709 slentol = sr3 ;
1710 sidetol = kRMin ;
1711 }
1712 }
1713 else if ( (sr2 < srd) && (sr2 > halfCarTolerance) )
1714 {
1715 srd = sr2 ;
1716 sider = kRMin ;
1717 }
1718 else if (sr2 > -halfCarTolerance)
1719 {
1720 // Intersection in tolerance. Store to check if it's good
1721
1722 slentol = sr2 ;
1723 sidetol = kRMin ;
1724 }
1725 if( slentol <= halfCarTolerance )
1726 {
1727 // An intersection within the tolerance was found.
1728 // We must accept it only if the momentum points outwards.
1729
1730 G4ThreeVector Normal ;
1731
1732 // Calculate a normal vector, as below
1733
1734 xi = p.x() + slentol*v.x() ;
1735 yi = p.y() + slentol*v.y() ;
1736 if( sidetol==kRMax )
1737 {
1738 risec = std::sqrt(xi*xi + yi*yi)*secRMax ;
1739 Normal = G4ThreeVector(xi/risec,yi/risec,-tanRMax/secRMax) ;
1740 }
1741 else
1742 {
1743 risec = std::sqrt(xi*xi + yi*yi)*secRMin ;
1744 Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin) ;
1745 }
1746 if( Normal.dot(v) > 0 )
1747 {
1748 // We will leave the cone immediately
1749
1750 if( calcNorm )
1751 {
1752 *n = Normal.unit() ;
1753 *validNorm = true ;
1754 }
1755 return snxt = 0.0 ;
1756 }
1757 else
1758 {
1759 // On the surface, but not heading out so we ignore this
1760 // intersection (as it is within tolerance).
1761
1762 slentol = kInfinity ;
1763 }
1764 }
1765 }
1766 }
1767 }
1768 }
1769
1770 // Linear case => point outside inner cone ---> outer cone intersect
1771 //
1772 // Phi Intersection
1773
1774 if ( !fPhiFullCone )
1775 {
1776 // add angle calculation with correction
1777 // of the difference in domain of atan2 and Sphi
1778
1779 vphi = std::atan2(v.y(),v.x()) ;
1780
1781 if ( vphi < fSPhi - halfAngTolerance ) { vphi += twopi; }
1782 else if ( vphi > fSPhi + fDPhi + halfAngTolerance ) { vphi -= twopi; }
1783
1784 if ( (p.x() != 0.0) || (p.y() != 0.0) ) // Check if on z axis (rho not needed later)
1785 {
1786 // pDist -ve when inside
1787
1788 pDistS = p.x()*sinSPhi - p.y()*cosSPhi ;
1789 pDistE = -p.x()*sinEPhi + p.y()*cosEPhi ;
1790
1791 // Comp -ve when in direction of outwards normal
1792
1793 compS = -sinSPhi*v.x() + cosSPhi*v.y() ;
1794 compE = sinEPhi*v.x() - cosEPhi*v.y() ;
1795
1796 sidephi = kNull ;
1797
1798 if( ( (fDPhi <= pi) && ( (pDistS <= halfCarTolerance)
1799 && (pDistE <= halfCarTolerance) ) )
1800 || ( (fDPhi > pi) && ((pDistS <= halfCarTolerance)
1801 || (pDistE <= halfCarTolerance) ) ) )
1802 {
1803 // Inside both phi *full* planes
1804 if ( compS < 0 )
1805 {
1806 sphi = pDistS/compS ;
1807 if (sphi >= -halfCarTolerance)
1808 {
1809 xi = p.x() + sphi*v.x() ;
1810 yi = p.y() + sphi*v.y() ;
1811
1812 // Check intersecting with correct half-plane
1813 // (if not -> no intersect)
1814 //
1815 if ( (std::fabs(xi)<=kCarTolerance)
1816 && (std::fabs(yi)<=kCarTolerance) )
1817 {
1818 sidephi= kSPhi;
1819 if ( ( fSPhi-halfAngTolerance <= vphi )
1820 && ( fSPhi+fDPhi+halfAngTolerance >=vphi ) )
1821 {
1822 sphi = kInfinity;
1823 }
1824 }
1825 else
1826 if ( (yi*cosCPhi-xi*sinCPhi)>=0 )
1827 {
1828 sphi = kInfinity ;
1829 }
1830 else
1831 {
1832 sidephi = kSPhi ;
1833 if ( pDistS > -halfCarTolerance )
1834 {
1835 sphi = 0.0 ; // Leave by sphi immediately
1836 }
1837 }
1838 }
1839 else
1840 {
1841 sphi = kInfinity ;
1842 }
1843 }
1844 else
1845 {
1846 sphi = kInfinity ;
1847 }
1848
1849 if ( compE < 0 )
1850 {
1851 sphi2 = pDistE/compE ;
1852
1853 // Only check further if < starting phi intersection
1854 //
1855 if ( (sphi2 > -halfCarTolerance) && (sphi2 < sphi) )
1856 {
1857 xi = p.x() + sphi2*v.x() ;
1858 yi = p.y() + sphi2*v.y() ;
1859
1860 // Check intersecting with correct half-plane
1861
1862 if ( (std::fabs(xi)<=kCarTolerance)
1863 && (std::fabs(yi)<=kCarTolerance) )
1864 {
1865 // Leaving via ending phi
1866
1867 if( (fSPhi-halfAngTolerance > vphi)
1868 || (fSPhi+fDPhi+halfAngTolerance < vphi) )
1869 {
1870 sidephi = kEPhi ;
1871 if ( pDistE <= -halfCarTolerance ) { sphi = sphi2; }
1872 else { sphi = 0.0; }
1873 }
1874 }
1875 else // Check intersecting with correct half-plane
1876 if ( yi*cosCPhi-xi*sinCPhi >= 0 )
1877 {
1878 // Leaving via ending phi
1879
1880 sidephi = kEPhi ;
1881 if ( pDistE <= -halfCarTolerance ) { sphi = sphi2; }
1882 else { sphi = 0.0; }
1883 }
1884 }
1885 }
1886 }
1887 else
1888 {
1889 sphi = kInfinity ;
1890 }
1891 }
1892 else
1893 {
1894 // On z axis + travel not || to z axis -> if phi of vector direction
1895 // within phi of shape, Step limited by rmax, else Step =0
1896
1897 if ( (fSPhi-halfAngTolerance <= vphi)
1898 && (vphi <= fSPhi+fDPhi+halfAngTolerance) )
1899 {
1900 sphi = kInfinity ;
1901 }
1902 else
1903 {
1904 sidephi = kSPhi ; // arbitrary
1905 sphi = 0.0 ;
1906 }
1907 }
1908 if ( sphi < snxt ) // Order intersecttions
1909 {
1910 snxt = sphi ;
1911 side = sidephi ;
1912 }
1913 }
1914 if ( srd < snxt ) // Order intersections
1915 {
1916 snxt = srd ;
1917 side = sider ;
1918 }
1919 if (calcNorm)
1920 {
1921 switch(side)
1922 { // Note: returned vector not normalised
1923 case kRMax: // (divide by frmax for unit vector)
1924 xi = p.x() + snxt*v.x() ;
1925 yi = p.y() + snxt*v.y() ;
1926 risec = std::sqrt(xi*xi + yi*yi)*secRMax ;
1927 *n = G4ThreeVector(xi/risec,yi/risec,-tanRMax/secRMax) ;
1928 *validNorm = true ;
1929 break ;
1930 case kRMin:
1931 *validNorm = false ; // Rmin is inconvex
1932 break ;
1933 case kSPhi:
1934 if ( fDPhi <= pi )
1935 {
1936 *n = G4ThreeVector(sinSPhi, -cosSPhi, 0);
1937 *validNorm = true ;
1938 }
1939 else
1940 {
1941 *validNorm = false ;
1942 }
1943 break ;
1944 case kEPhi:
1945 if ( fDPhi <= pi )
1946 {
1947 *n = G4ThreeVector(-sinEPhi, cosEPhi, 0);
1948 *validNorm = true ;
1949 }
1950 else
1951 {
1952 *validNorm = false ;
1953 }
1954 break ;
1955 case kPZ:
1956 *n = G4ThreeVector(0,0,1) ;
1957 *validNorm = true ;
1958 break ;
1959 case kMZ:
1960 *n = G4ThreeVector(0,0,-1) ;
1961 *validNorm = true ;
1962 break ;
1963 default:
1964 G4cout << G4endl ;
1965 DumpInfo();
1966 std::ostringstream message;
1967 G4long oldprc = message.precision(16) ;
1968 message << "Undefined side for valid surface normal to solid."
1969 << G4endl
1970 << "Position:" << G4endl << G4endl
1971 << "p.x() = " << p.x()/mm << " mm" << G4endl
1972 << "p.y() = " << p.y()/mm << " mm" << G4endl
1973 << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl
1974 << "pho at z = " << std::sqrt( p.x()*p.x()+p.y()*p.y() )/mm
1975 << " mm" << G4endl << G4endl ;
1976 if( p.x() != 0. || p.y() != 0.)
1977 {
1978 message << "point phi = " << std::atan2(p.y(),p.x())/degree
1979 << " degrees" << G4endl << G4endl ;
1980 }
1981 message << "Direction:" << G4endl << G4endl
1982 << "v.x() = " << v.x() << G4endl
1983 << "v.y() = " << v.y() << G4endl
1984 << "v.z() = " << v.z() << G4endl<< G4endl
1985 << "Proposed distance :" << G4endl<< G4endl
1986 << "snxt = " << snxt/mm << " mm" << G4endl ;
1987 message.precision(oldprc) ;
1988 G4Exception("G4Cons::DistanceToOut(p,v,..)","GeomSolids1002",
1989 JustWarning, message) ;
1990 break ;
1991 }
1992 }
1993 if (snxt < halfCarTolerance) { snxt = 0.; }
1994
1995 return snxt ;
1996}
ESide
ESide
Definition G4Sphere.cc:61
G4long
long G4long
Definition G4Types.hh:87
CLHEP::Hep3Vector::unit
Hep3Vector unit() const

◆ GetCosEndPhi()

G4double G4Cons::GetCosEndPhi ( ) const
inline

Referenced by BoundingLimits(), and CalculateExtent().

◆ GetCosStartPhi()

G4double G4Cons::GetCosStartPhi ( ) const
inline

Referenced by BoundingLimits(), and CalculateExtent().

◆ GetCubicVolume()

G4double G4Cons::GetCubicVolume ( )
inlineoverridevirtual

Reimplemented from G4VSolid.

◆ GetDeltaPhiAngle()

G4double G4Cons::GetDeltaPhiAngle ( ) const
inline

Referenced by G4HepRepFileSceneHandler::AddSolid(), BoundingLimits(), G4tgbVolume::BuildSolidForDivision(), CalculateExtent(), G4GDMLWriteParamvol::Cone_dimensionsWrite(), G4GDMLWriteSolids::ConeWrite(), and G4tgbGeometryDumper::GetSolidParams().

◆ GetEntityType()

G4GeometryType G4Cons::GetEntityType ( ) const
overridevirtual

Implements G4VSolid.

Definition at line 2080 of file G4Cons.cc.

2081{
2082 return {"G4Cons"};
2083}

◆ GetInnerRadiusMinusZ()

G4double G4Cons::GetInnerRadiusMinusZ ( ) const
inline

Referenced by G4HepRepFileSceneHandler::AddSolid(), BoundingLimits(), G4tgbVolume::BuildSolidForDivision(), CalculateExtent(), G4GDMLWriteParamvol::Cone_dimensionsWrite(), G4GDMLWriteSolids::ConeWrite(), and G4tgbGeometryDumper::GetSolidParams().

◆ GetInnerRadiusPlusZ()

G4double G4Cons::GetInnerRadiusPlusZ ( ) const
inline

Referenced by G4HepRepFileSceneHandler::AddSolid(), BoundingLimits(), G4tgbVolume::BuildSolidForDivision(), CalculateExtent(), G4GDMLWriteParamvol::Cone_dimensionsWrite(), G4GDMLWriteSolids::ConeWrite(), and G4tgbGeometryDumper::GetSolidParams().

◆ GetOuterRadiusMinusZ()

G4double G4Cons::GetOuterRadiusMinusZ ( ) const
inline

Referenced by G4HepRepFileSceneHandler::AddSolid(), BoundingLimits(), G4tgbVolume::BuildSolidForDivision(), CalculateExtent(), G4GDMLWriteParamvol::Cone_dimensionsWrite(), G4GDMLWriteSolids::ConeWrite(), and G4tgbGeometryDumper::GetSolidParams().

◆ GetOuterRadiusPlusZ()

G4double G4Cons::GetOuterRadiusPlusZ ( ) const
inline

Referenced by G4HepRepFileSceneHandler::AddSolid(), BoundingLimits(), G4tgbVolume::BuildSolidForDivision(), CalculateExtent(), G4GDMLWriteParamvol::Cone_dimensionsWrite(), G4GDMLWriteSolids::ConeWrite(), and G4tgbGeometryDumper::GetSolidParams().

◆ GetPointOnSurface()

G4ThreeVector G4Cons::GetPointOnSurface ( ) const
overridevirtual

Reimplemented from G4VSolid.

Definition at line 2125 of file G4Cons.cc.

2126{
2127 // declare working variables
2128 //
2129 G4double rone = (fRmax1-fRmax2)/(2.*fDz);
2130 G4double rtwo = (fRmin1-fRmin2)/(2.*fDz);
2131 G4double qone = (fRmax1 == fRmax2) ? 0. : fDz*(fRmax1+fRmax2)/(fRmax1-fRmax2);
2132 G4double qtwo = (fRmin1 == fRmin2) ? 0. : fDz*(fRmin1+fRmin2)/(fRmin1-fRmin2);
2133
2134 G4double slin = std::hypot(fRmin1-fRmin2, 2.*fDz);
2135 G4double slout = std::hypot(fRmax1-fRmax2, 2.*fDz);
2136 G4double Aone = 0.5*fDPhi*(fRmax2 + fRmax1)*slout; // outer surface
2137 G4double Atwo = 0.5*fDPhi*(fRmin2 + fRmin1)*slin; // inner surface
2138 G4double Athree = 0.5*fDPhi*(fRmax1*fRmax1-fRmin1*fRmin1); // base at -Dz
2139 G4double Afour = 0.5*fDPhi*(fRmax2*fRmax2-fRmin2*fRmin2); // base at +Dz
2140 G4double Afive = fDz*(fRmax1-fRmin1+fRmax2-fRmin2); // phi section
2141
2142 G4double phi = G4RandFlat::shoot(fSPhi,fSPhi+fDPhi);
2143 G4double cosu = std::cos(phi);
2144 G4double sinu = std::sin(phi);
2145 G4double rRand1 = GetRadiusInRing(fRmin1, fRmax1);
2146 G4double rRand2 = GetRadiusInRing(fRmin2, fRmax2);
2147
2148 if ( (fSPhi == 0.) && fPhiFullCone ) { Afive = 0.; }
2149 G4double chose = G4RandFlat::shoot(0.,Aone+Atwo+Athree+Afour+2.*Afive);
2150
2151 if( (chose >= 0.) && (chose < Aone) ) // outer surface
2152 {
2153 if(fRmax1 != fRmax2)
2154 {
2155 G4double zRand = G4RandFlat::shoot(-1.*fDz,fDz);
2156 return { rone*cosu*(qone-zRand), rone*sinu*(qone-zRand), zRand };
2157 }
2158 else
2159 {
2160 return { fRmax1*cosu, fRmax2*sinu, G4RandFlat::shoot(-1.*fDz,fDz) };
2161 }
2162 }
2163 else if( (chose >= Aone) && (chose < Aone + Atwo) ) // inner surface
2164 {
2165 if(fRmin1 != fRmin2)
2166 {
2167 G4double zRand = G4RandFlat::shoot(-1.*fDz,fDz);
2168 return { rtwo*cosu*(qtwo-zRand), rtwo*sinu*(qtwo-zRand), zRand };
2169 }
2170 else
2171 {
2172 return { fRmin1*cosu, fRmin2*sinu, G4RandFlat::shoot(-1.*fDz,fDz) };
2173 }
2174 }
2175 else if( (chose >= Aone + Atwo) && (chose < Aone + Atwo + Athree) ) // base at -Dz
2176 {
2177 return {rRand1*cosu, rRand1*sinu, -1*fDz};
2178 }
2179 else if( (chose >= Aone + Atwo + Athree)
2180 && (chose < Aone + Atwo + Athree + Afour) ) // base at +Dz
2181 {
2182 return { rRand2*cosu, rRand2*sinu, fDz };
2183 }
2184 else if( (chose >= Aone + Atwo + Athree + Afour) // SPhi section
2185 && (chose < Aone + Atwo + Athree + Afour + Afive) )
2186 {
2187 G4double zRand = G4RandFlat::shoot(-1.*fDz,fDz);
2188 rRand1 = G4RandFlat::shoot(fRmin2-((zRand-fDz)/(2.*fDz))*(fRmin1-fRmin2),
2189 fRmax2-((zRand-fDz)/(2.*fDz))*(fRmax1-fRmax2));
2190 return { rRand1*cosSPhi, rRand1*sinSPhi, zRand };
2191 }
2192 else // SPhi+DPhi section
2193 {
2194 G4double zRand = G4RandFlat::shoot(-1.*fDz,fDz);
2195 rRand1 = G4RandFlat::shoot(fRmin2-((zRand-fDz)/(2.*fDz))*(fRmin1-fRmin2),
2196 fRmax2-((zRand-fDz)/(2.*fDz))*(fRmax1-fRmax2));
2197 return { rRand1*cosEPhi, rRand1*sinEPhi, zRand };
2198 }
2199}
G4CSGSolid::GetRadiusInRing
G4double GetRadiusInRing(G4double rmin, G4double rmax) const
Definition G4CSGSolid.cc:109

◆ GetSinEndPhi()

G4double G4Cons::GetSinEndPhi ( ) const
inline

Referenced by BoundingLimits(), and CalculateExtent().

◆ GetSinStartPhi()

G4double G4Cons::GetSinStartPhi ( ) const
inline

Referenced by BoundingLimits(), and CalculateExtent().

◆ GetStartPhiAngle()

G4double G4Cons::GetStartPhiAngle ( ) const
inline

Referenced by G4tgbVolume::BuildSolidForDivision(), G4GDMLWriteParamvol::Cone_dimensionsWrite(), G4GDMLWriteSolids::ConeWrite(), and G4tgbGeometryDumper::GetSolidParams().

◆ GetSurfaceArea()

G4double G4Cons::GetSurfaceArea ( )
inlineoverridevirtual

Reimplemented from G4VSolid.

◆ GetZHalfLength()

G4double G4Cons::GetZHalfLength ( ) const
inline

Referenced by G4HepRepFileSceneHandler::AddSolid(), BoundingLimits(), G4tgbVolume::BuildSolidForDivision(), CalculateExtent(), G4GDMLWriteParamvol::Cone_dimensionsWrite(), G4GDMLWriteSolids::ConeWrite(), and G4tgbGeometryDumper::GetSolidParams().

◆ Inside()

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

Implements G4VSolid.

Definition at line 181 of file G4Cons.cc.

182{
183 G4double r2, rl, rh, pPhi, tolRMin, tolRMax; // rh2, rl2 ;
184 EInside in;
185
186 if (std::fabs(p.z()) > fDz + halfCarTolerance ) { return in = kOutside; }
187 else if(std::fabs(p.z()) >= fDz - halfCarTolerance ) { in = kSurface; }
188 else { in = kInside; }
189
190 r2 = p.x()*p.x() + p.y()*p.y() ;
191 rl = 0.5*(fRmin2*(p.z() + fDz) + fRmin1*(fDz - p.z()))/fDz ;
192 rh = 0.5*(fRmax2*(p.z()+fDz)+fRmax1*(fDz-p.z()))/fDz;
193
194 // rh2 = rh*rh;
195
196 tolRMin = rl - halfRadTolerance;
197 if ( tolRMin < 0 ) { tolRMin = 0; }
198 tolRMax = rh + halfRadTolerance;
199
200 if ( (r2<tolRMin*tolRMin) || (r2>tolRMax*tolRMax) ) { return in = kOutside; }
201
202 if (rl != 0.0) { tolRMin = rl + halfRadTolerance; }
203 else { tolRMin = 0.0; }
204 tolRMax = rh - halfRadTolerance;
205
206 if (in == kInside) // else it's kSurface already
207 {
208 if ( (r2 < tolRMin*tolRMin) || (r2 >= tolRMax*tolRMax) ) { in = kSurface; }
209 }
210 if ( !fPhiFullCone && ((p.x() != 0.0) || (p.y() != 0.0)) )
211 {
212 pPhi = std::atan2(p.y(),p.x()) ;
213
214 if ( pPhi < fSPhi - halfAngTolerance ) { pPhi += twopi; }
215 else if ( pPhi > fSPhi + fDPhi + halfAngTolerance ) { pPhi -= twopi; }
216
217 if ( (pPhi < fSPhi - halfAngTolerance) ||
218 (pPhi > fSPhi + fDPhi + halfAngTolerance) ) { return in = kOutside; }
219
220 else if (in == kInside) // else it's kSurface anyway already
221 {
222 if ( (pPhi < fSPhi + halfAngTolerance) ||
223 (pPhi > fSPhi + fDPhi - halfAngTolerance) ) { in = kSurface; }
224 }
225 }
226 else if ( !fPhiFullCone ) { in = kSurface; }
227
228 return in ;
229}
EInside
EInside
Definition geomdefs.hh:67
kInside
@ kInside
Definition geomdefs.hh:70
kSurface
@ kSurface
Definition geomdefs.hh:69

Referenced by DistanceToOut().

◆ operator=()

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

Definition at line 148 of file G4Cons.cc.

149{
150 // Check assignment to self
151 //
152 if (this == &rhs) { return *this; }
153
154 // Copy base class data
155 //
156 G4CSGSolid::operator=(rhs);
157
158 // Copy data
159 //
160 kRadTolerance = rhs.kRadTolerance;
161 kAngTolerance = rhs.kAngTolerance;
162 fRmin1 = rhs.fRmin1; fRmin2 = rhs.fRmin2;
163 fRmax1 = rhs.fRmax1; fRmax2 = rhs.fRmax2;
164 fDz = rhs.fDz; fSPhi = rhs.fSPhi; fDPhi = rhs.fDPhi;
165 sinCPhi = rhs.sinCPhi; cosCPhi = rhs.cosCPhi; cosHDPhi = rhs.cosHDPhi;
166 cosHDPhiOT = rhs.cosHDPhiOT; cosHDPhiIT = rhs.cosHDPhiIT;
167 sinSPhi = rhs.sinSPhi; cosSPhi = rhs.cosSPhi;
168 sinEPhi = rhs.sinEPhi; cosEPhi = rhs.cosEPhi;
169 fPhiFullCone = rhs.fPhiFullCone;
170 halfCarTolerance = rhs.halfCarTolerance;
171 halfRadTolerance = rhs.halfRadTolerance;
172 halfAngTolerance = rhs.halfAngTolerance;
173
174 return *this;
175}
G4CSGSolid::operator=
G4CSGSolid & operator=(const G4CSGSolid &rhs)
Definition G4CSGSolid.cc:89

◆ SetDeltaPhiAngle()

void G4Cons::SetDeltaPhiAngle ( G4double newDPhi)
inline

Referenced by G4ParameterisationConsPhi::ComputeDimensions(), G4ParameterisationConsRho::ComputeDimensions(), and G4ParameterisationConsZ::ComputeDimensions().

◆ SetInnerRadiusMinusZ()

void G4Cons::SetInnerRadiusMinusZ ( G4double Rmin1)
inline

Referenced by G4ParameterisationConsPhi::ComputeDimensions(), G4ParameterisationConsRho::ComputeDimensions(), and G4ParameterisationConsZ::ComputeDimensions().

◆ SetInnerRadiusPlusZ()

void G4Cons::SetInnerRadiusPlusZ ( G4double Rmin2)
inline

Referenced by G4ParameterisationConsPhi::ComputeDimensions(), G4ParameterisationConsRho::ComputeDimensions(), and G4ParameterisationConsZ::ComputeDimensions().

◆ SetOuterRadiusMinusZ()

void G4Cons::SetOuterRadiusMinusZ ( G4double Rmax1)
inline

Referenced by G4ParameterisationConsPhi::ComputeDimensions(), G4ParameterisationConsRho::ComputeDimensions(), and G4ParameterisationConsZ::ComputeDimensions().

◆ SetOuterRadiusPlusZ()

void G4Cons::SetOuterRadiusPlusZ ( G4double Rmax2)
inline

Referenced by G4ParameterisationConsPhi::ComputeDimensions(), G4ParameterisationConsRho::ComputeDimensions(), and G4ParameterisationConsZ::ComputeDimensions().

◆ SetStartPhiAngle()

void G4Cons::SetStartPhiAngle ( G4double newSPhi,
G4bool trig = true )
inline

Referenced by G4ParameterisationConsPhi::ComputeDimensions(), G4ParameterisationConsRho::ComputeDimensions(), and G4ParameterisationConsZ::ComputeDimensions().

◆ SetZHalfLength()

void G4Cons::SetZHalfLength ( G4double newDz)
inline

Referenced by G4ParameterisationConsPhi::ComputeDimensions(), G4ParameterisationConsRho::ComputeDimensions(), and G4ParameterisationConsZ::ComputeDimensions().

◆ StreamInfo()

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

Implements G4VSolid.

Definition at line 2098 of file G4Cons.cc.

2099{
2100 G4long oldprc = os.precision(16);
2101 os << "-----------------------------------------------------------\n"
2102 << " *** Dump for solid - " << GetName() << " ***\n"
2103 << " ===================================================\n"
2104 << " Solid type: G4Cons\n"
2105 << " Parameters: \n"
2106 << " inside -fDz radius: " << fRmin1/mm << " mm \n"
2107 << " outside -fDz radius: " << fRmax1/mm << " mm \n"
2108 << " inside +fDz radius: " << fRmin2/mm << " mm \n"
2109 << " outside +fDz radius: " << fRmax2/mm << " mm \n"
2110 << " half length in Z : " << fDz/mm << " mm \n"
2111 << " starting angle of segment: " << fSPhi/degree << " degrees \n"
2112 << " delta angle of segment : " << fDPhi/degree << " degrees \n"
2113 << "-----------------------------------------------------------\n";
2114 os.precision(oldprc);
2115
2116 return os;
2117}

◆ SurfaceNormal()

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

Implements G4VSolid.

Definition at line 405 of file G4Cons.cc.

406{
407 G4int noSurfaces = 0;
408 G4double rho, pPhi;
409 G4double distZ, distRMin, distRMax;
410 G4double distSPhi = kInfinity, distEPhi = kInfinity;
411 G4double tanRMin, secRMin, pRMin, widRMin;
412 G4double tanRMax, secRMax, pRMax, widRMax;
413
414 G4ThreeVector norm, sumnorm(0.,0.,0.), nZ = G4ThreeVector(0.,0.,1.);
415 G4ThreeVector nR, nr(0.,0.,0.), nPs, nPe;
416
417 distZ = std::fabs(std::fabs(p.z()) - fDz);
418 rho = std::sqrt(p.x()*p.x() + p.y()*p.y());
419
420 tanRMin = (fRmin2 - fRmin1)*0.5/fDz;
421 secRMin = std::sqrt(1 + tanRMin*tanRMin);
422 pRMin = rho - p.z()*tanRMin;
423 widRMin = fRmin2 - fDz*tanRMin;
424 distRMin = std::fabs(pRMin - widRMin)/secRMin;
425
426 tanRMax = (fRmax2 - fRmax1)*0.5/fDz;
427 secRMax = std::sqrt(1+tanRMax*tanRMax);
428 pRMax = rho - p.z()*tanRMax;
429 widRMax = fRmax2 - fDz*tanRMax;
430 distRMax = std::fabs(pRMax - widRMax)/secRMax;
431
432 if (!fPhiFullCone) // Protected against (0,0,z)
433 {
434 if ( rho != 0.0 )
435 {
436 pPhi = std::atan2(p.y(),p.x());
437
438 if (pPhi < fSPhi-halfCarTolerance) { pPhi += twopi; }
439 else if (pPhi > fSPhi+fDPhi+halfCarTolerance) { pPhi -= twopi; }
440
441 distSPhi = std::fabs( pPhi - fSPhi );
442 distEPhi = std::fabs( pPhi - fSPhi - fDPhi );
443 }
444 else if( ((fRmin1) == 0.0) || ((fRmin2) == 0.0) )
445 {
446 distSPhi = 0.;
447 distEPhi = 0.;
448 }
449 nPs = G4ThreeVector( sinSPhi, -cosSPhi, 0 );
450 nPe = G4ThreeVector( -sinEPhi, cosEPhi, 0 );
451 }
452 if ( rho > halfCarTolerance )
453 {
454 nR = G4ThreeVector(p.x()/rho/secRMax, p.y()/rho/secRMax, -tanRMax/secRMax);
455 if ((fRmin1 != 0.0) || (fRmin2 != 0.0))
456 {
457 nr = G4ThreeVector(-p.x()/rho/secRMin,-p.y()/rho/secRMin,tanRMin/secRMin);
458 }
459 }
460
461 if( distRMax <= halfCarTolerance )
462 {
463 ++noSurfaces;
464 sumnorm += nR;
465 }
466 if( ((fRmin1 != 0.0) || (fRmin2 != 0.0)) && (distRMin <= halfCarTolerance) )
467 {
468 ++noSurfaces;
469 sumnorm += nr;
470 }
471 if( !fPhiFullCone )
472 {
473 if (distSPhi <= halfAngTolerance)
474 {
475 ++noSurfaces;
476 sumnorm += nPs;
477 }
478 if (distEPhi <= halfAngTolerance)
479 {
480 ++noSurfaces;
481 sumnorm += nPe;
482 }
483 }
484 if (distZ <= halfCarTolerance)
485 {
486 ++noSurfaces;
487 if ( p.z() >= 0.) { sumnorm += nZ; }
488 else { sumnorm -= nZ; }
489 }
490 if ( noSurfaces == 0 )
491 {
492#ifdef G4CSGDEBUG
493 G4Exception("G4Cons::SurfaceNormal(p)", "GeomSolids1002",
494 JustWarning, "Point p is not on surface !?" );
495#endif
496 norm = ApproxSurfaceNormal(p);
497 }
498 else if ( noSurfaces == 1 ) { norm = sumnorm; }
499 else { norm = sumnorm.unit(); }
500
501 return norm ;
502}
The documentation for this class was generated from the following files: