Geant4 11.4.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4Cons Class Reference

G4Cons is, in the general case, a Phi segment of a cone, with half-length fDz, inner and outer radii specified at -fDz and +fDz. The Phi segment is described by a starting fSPhi angle, and the +fDPhi delta angle for the shape. If the delta angle is >=2*pi, the shape is treated as continuous in Phi. More...

#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=default
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)=default
G4Cons & operator= (const G4Cons &rhs)
Public Member Functions inherited from G4CSGSolid
 G4CSGSolid (const G4String &pName)
 ~G4CSGSolid () override
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 ()
 G4VSolid (const G4VSolid &rhs)
G4VSolid & operator= (const G4VSolid &rhs)
G4bool operator== (const G4VSolid &s) const
G4String GetName () const
void SetName (const G4String &name)
G4double GetTolerance () const
virtual G4int GetNumOfConstituents () const
virtual G4bool IsFaceted () 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__ &)
G4double EstimateCubicVolume (G4int nStat, G4double epsilon) const
G4double EstimateSurfaceArea (G4int nStat, G4double epsilon) 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

G4Cons is, in the general case, a Phi segment of a cone, with half-length fDz, inner and outer radii specified at -fDz and +fDz. The Phi segment is described by a starting fSPhi angle, and the +fDPhi delta angle for the shape. If the delta angle is >=2*pi, the shape is treated as continuous in Phi.

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

Constructs a cone with the given name and dimensions.

Parameters
[in]pNameThe name of the solid.
[in]pRmin1Inside radius at -fDz.
[in]pRmax1Outside radius at -fDz
[in]pRmin2Inside radius at +fDz.
[in]pRmax2Outside radius at +fDz
[in]pDZHalf length in Z.
[in]pSPhiStarting angle of the segment in radians.
[in]pDPhiDelta angle of the segment in radians.

Definition at line 83 of file G4Cons.cc.

88 : G4CSGSolid(pName), fRmin1(pRmin1), fRmin2(pRmin2),
89 fRmax1(pRmax1), fRmax2(pRmax2), fDz(pDz), fSPhi(0.), fDPhi(0.)
90{
91 kRadTolerance = G4GeometryTolerance::GetInstance()->GetRadialTolerance();
92 kAngTolerance = G4GeometryTolerance::GetInstance()->GetAngularTolerance();
93
94 halfCarTolerance=kCarTolerance*0.5;
95 halfRadTolerance=kRadTolerance*0.5;
96 halfAngTolerance=kAngTolerance*0.5;
97
98 // Check z-len
99 //
100 if ( pDz < 0 )
101 {
102 std::ostringstream message;
103 message << "Invalid Z half-length for Solid: " << GetName() << G4endl
104 << " hZ = " << pDz;
105 G4Exception("G4Cons::G4Cons()", "GeomSolids0002",
106 FatalException, message);
107 }
108
109 // Check radii
110 //
111 if (((pRmin1>=pRmax1) || (pRmin2>=pRmax2) || (pRmin1<0)) && (pRmin2<0))
112 {
113 std::ostringstream message;
114 message << "Invalid values of radii for Solid: " << GetName() << G4endl
115 << " pRmin1 = " << pRmin1 << ", pRmin2 = " << pRmin2
116 << ", pRmax1 = " << pRmax1 << ", pRmax2 = " << pRmax2;
117 G4Exception("G4Cons::G4Cons()", "GeomSolids0002",
118 FatalException, message) ;
119 }
120 if( (pRmin1 == 0.0) && (pRmin2 > 0.0) ) { fRmin1 = 1e3*kRadTolerance ; }
121 if( (pRmin2 == 0.0) && (pRmin1 > 0.0) ) { fRmin2 = 1e3*kRadTolerance ; }
122
123 // Check angles
124 //
125 CheckPhiAngles(pSPhi, pDPhi);
126}
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:418

Referenced by Clone(), G4Cons(), and operator=().

◆ ~G4Cons()

G4Cons::~G4Cons ( )
overridedefault

Default destructor.

◆ G4Cons() [2/3]

G4Cons::G4Cons ( __void__ & a)

Fake default constructor for usage restricted to direct object persistency for clients requiring preallocation of memory for persistifiable objects.

Definition at line 133 of file G4Cons.cc.

134 : G4CSGSolid(a)
135{
136}

◆ G4Cons() [3/3]

G4Cons::G4Cons ( const G4Cons & rhs)
default

Copy constructor and assignment operator.

Member Function Documentation

◆ BoundingLimits()

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

Computes the bounding limits of the solid.

Parameters
[out]pMinThe minimum bounding limit point.
[out]pMaxThe maximum bounding limit point.

Reimplemented from G4VSolid.

Definition at line 241 of file G4Cons.cc.

242{
243 G4double rmin = std::min(GetInnerRadiusMinusZ(),GetInnerRadiusPlusZ());
244 G4double rmax = std::max(GetOuterRadiusMinusZ(),GetOuterRadiusPlusZ());
245 G4double dz = GetZHalfLength();
246
247 // Find bounding box
248 //
249 if (GetDeltaPhiAngle() < twopi)
250 {
251 G4TwoVector vmin,vmax;
252 G4GeomTools::DiskExtent(rmin,rmax,
253 GetSinStartPhi(),GetCosStartPhi(),
254 GetSinEndPhi(),GetCosEndPhi(),
255 vmin,vmax);
256 pMin.set(vmin.x(),vmin.y(),-dz);
257 pMax.set(vmax.x(),vmax.y(), dz);
258 }
259 else
260 {
261 pMin.set(-rmax,-rmax,-dz);
262 pMax.set( rmax, rmax, dz);
263 }
264
265 // Check correctness of the bounding box
266 //
267 if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
268 {
269 std::ostringstream message;
270 message << "Bad bounding box (min >= max) for solid: "
271 << GetName() << " !"
272 << "\npMin = " << pMin
273 << "\npMax = " << pMax;
274 G4Exception("G4Cons::BoundingLimits()", "GeomMgt0001",
275 JustWarning, message);
276 DumpInfo();
277 }
278}
JustWarning
@ JustWarning
Definition G4ExceptionSeverity.hh:71
G4TwoVector
CLHEP::Hep2Vector G4TwoVector
Definition G4TwoVector.hh:36
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:399
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

Calculates the minimum and maximum extent of the solid, when under the specified transform, and within the specified limits.

Parameters
[in]pAxisThe axis along which compute the extent.
[in]pVoxelLimitThe limiting space dictated by voxels.
[in]pTransformThe internal transformation applied to the solid.
[out]pMinThe minimum extent value.
[out]pMaxThe maximum extent value.
Returns
True if the solid is intersected by the extent region.

Implements G4VSolid.

Definition at line 284 of file G4Cons.cc.

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

◆ Clone()

G4VSolid * G4Cons::Clone ( ) const
overridevirtual

Makes a clone of the object for use in multi-treading.

Returns
A pointer to the new cloned allocated solid.

Reimplemented from G4VSolid.

Definition at line 2073 of file G4Cons.cc.

2074{
2075 return new G4Cons(*this);
2076}
G4Cons::G4Cons
G4Cons(const G4String &pName, G4double pRmin1, G4double pRmax1, G4double pRmin2, G4double pRmax2, G4double pDz, G4double pSPhi, G4double pDPhi)
Definition G4Cons.cc:83

◆ ComputeDimensions()

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

Dispatch method for parameterisation replication mechanism and dimension computation.

Reimplemented from G4VSolid.

Definition at line 230 of file G4Cons.cc.

233{
234 p->ComputeDimensions(*this,n,pRep) ;
235}
G4VPVParameterisation::ComputeDimensions
virtual void ComputeDimensions(G4Box &, const G4int, const G4VPhysicalVolume *) const
Definition G4VPVParameterisation.hh:123

◆ CreatePolyhedron()

G4Polyhedron * G4Cons::CreatePolyhedron ( ) const
overridevirtual

Creates a Polyhedron used for Visualisation. It is the caller's responsibility to delete it. A null pointer means "not created".

Reimplemented from G4VSolid.

Definition at line 2239 of file G4Cons.cc.

2240{
2241 return new G4PolyhedronCons(fRmin1,fRmax1,fRmin2,fRmax2,fDz,fSPhi,fDPhi);
2242}

◆ DescribeYourselfTo()

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

Methods for creating graphical representations (i.e. for visualisation).

Implements G4VSolid.

Definition at line 2234 of file G4Cons.cc.

2235{
2236 scene.AddSolid (*this);
2237}
G4VGraphicsScene::AddSolid
virtual void AddSolid(const G4Box &)=0

◆ DistanceToIn() [1/2]

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

Calculates the distance to the nearest surface of a shape from an outside point. The distance can be an underestimate.

Parameters
[in]pThe point at offset p.
Returns
The safety distance to enter the shape.

Implements G4VSolid.

Definition at line 1306 of file G4Cons.cc.

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

◆ DistanceToIn() [2/2]

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

Returns the distance along the normalised vector 'v' to the shape, from the point at offset 'p'. If there is no intersection, returns kInfinity. The first intersection resulting from 'leaving' a surface/volume is discarded. Hence, it is tolerant of points on the surface of the shape.

Parameters
[in]pThe point at offset p.
[in]vThe normalised direction vector.
Returns
The distance to enter the shape.

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
847 xi = p.x() + sd*v.x() ;
848 yi = p.y() + sd*v.y() ;
849 ri = rMaxAv + zi*tanRMax ;
850 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
851
852 if ( cosPsi >= cosHDPhiIT ) { return sd; }
853 }
854 } // end if (sd>0)
855 }
856 }
857 else
858 {
859 // Inside outer cone
860 // check not inside, and heading through G4Cons (-> 0 to in)
861
862 if ( ( t3 > (rin + halfRadTolerance*secRMin)*
863 (rin + halfRadTolerance*secRMin) )
864 && (nt2 < 0) && (d >= 0) && (std::fabs(p.z()) <= tolIDz) )
865 {
866 // Inside cones, delta r -ve, inside z extent
867 // Point is on the Surface => check Direction using Normal.dot(v)
868
869 xi = p.x() ;
870 yi = p.y() ;
871 risec = std::sqrt(xi*xi + yi*yi)*secRMax ;
872 Normal = G4ThreeVector(xi/risec,yi/risec,-tanRMax/secRMax) ;
873 if ( !fPhiFullCone )
874 {
875 cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/std::sqrt(t3) ;
876 if ( cosPsi >= cosHDPhiIT )
877 {
878 if ( Normal.dot(v) <= 0 ) { return 0.0; }
879 }
880 }
881 else
882 {
883 if ( Normal.dot(v) <= 0 ) { return 0.0; }
884 }
885 }
886 }
887 }
888 else // Single root case
889 {
890 if ( std::fabs(nt2) > kRadTolerance )
891 {
892 sd = -0.5*nt3/nt2 ;
893
894 if ( sd < 0 ) { return kInfinity; } // travel away
895
896 // sd >= 0, If 'forwards'. Check z intersection
897 zi = p.z() + sd*v.z() ;
898
899 if ((std::fabs(zi) <= tolODz) && (nt2 < 0))
900 {
901 // Z ok. Check phi intersection if reqd
902
903 if ( fPhiFullCone ) { return sd; }
904
905 xi = p.x() + sd*v.x() ;
906 yi = p.y() + sd*v.y() ;
907 ri = rMaxAv + zi*tanRMax ;
908 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
909
910 if (cosPsi >= cosHDPhiIT) { return sd; }
911 }
912 }
913 else // travel || cone surface from its origin
914 {
915 sd = kInfinity ;
916 }
917 }
918
919 // Inner Cone Intersection
920 // o Space is divided into 3 areas:
921 // 1) Radius greater than real inner cone & imaginary cone & outside
922 // tolerance
923 // 2) Radius less than inner or imaginary cone & outside tolarance
924 // 3) Within tolerance of real or imaginary cones
925 // - Extra checks needed for 3's intersections
926 // => lots of duplicated code
927
928 if (rMinAv != 0.0)
929 {
930 nt1 = t1 - (tanRMin*v.z())*(tanRMin*v.z()) ;
931 nt2 = t2 - tanRMin*v.z()*rin ;
932 nt3 = t3 - rin*rin ;
933
934 if ( nt1 != 0.0 )
935 {
936 if ( nt3 > rin*kRadTolerance*secRMin )
937 {
938 // At radius greater than real & imaginary cones
939 // -> 2nd root, with zi check
940
941 b = nt2/nt1 ;
942 c = nt3/nt1 ;
943 d = b*b-c ;
944 if (d >= 0) // > 0
945 {
946 if(b>0){sd = c/( -b-std::sqrt(d));}
947 else {sd = -b + std::sqrt(d) ;}
948
949 if ( sd >= 0 ) // > 0
950 {
951 if ( sd>dRmax ) // Avoid rounding errors due to precision issues on
952 { // 64 bits systems. Split long distance and recompute
953 G4double fTerm = sd-std::fmod(sd,dRmax);
954 sd = fTerm + DistanceToIn(p+fTerm*v,v);
955 }
956 zi = p.z() + sd*v.z() ;
957
958 if ( std::fabs(zi) <= tolODz )
959 {
960 if ( !fPhiFullCone )
961 {
962 xi = p.x() + sd*v.x() ;
963 yi = p.y() + sd*v.y() ;
964 ri = rMinAv + zi*tanRMin ;
965 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
966
967 if (cosPsi >= cosHDPhiIT)
968 {
969 if ( sd > halfRadTolerance ) { snxt=sd; }
970 else
971 {
972 // Calculate a normal vector in order to check Direction
973
974 risec = std::sqrt(xi*xi + yi*yi)*secRMin ;
975 Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin);
976 if ( Normal.dot(v) <= 0 ) { snxt = sd; }
977 }
978 }
979 }
980 else
981 {
982 if ( sd > halfRadTolerance ) { return sd; }
983
984 // Calculate a normal vector in order to check Direction
985
986 xi = p.x() + sd*v.x() ;
987 yi = p.y() + sd*v.y() ;
988 risec = std::sqrt(xi*xi + yi*yi)*secRMin ;
989 Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin) ;
990 if ( Normal.dot(v) <= 0 ) { return sd; }
991 }
992 }
993 }
994 }
995 }
996 else if ( nt3 < -rin*kRadTolerance*secRMin )
997 {
998 // Within radius of inner cone (real or imaginary)
999 // -> Try 2nd root, with checking intersection is with real cone
1000 // -> If check fails, try 1st root, also checking intersection is
1001 // on real cone
1002
1003 b = nt2/nt1 ;
1004 c = nt3/nt1 ;
1005 d = b*b - c ;
1006
1007 if ( d >= 0 ) // > 0
1008 {
1009 if (b>0) { sd = c/(-b-std::sqrt(d)); }
1010 else { sd = -b + std::sqrt(d); }
1011 zi = p.z() + sd*v.z() ;
1012 ri = rMinAv + zi*tanRMin ;
1013
1014 if ( ri > 0 )
1015 {
1016 if ( (sd >= 0) && (std::fabs(zi) <= tolODz) ) // sd > 0
1017 {
1018 if ( sd>dRmax ) // Avoid rounding errors due to precision issues
1019 { // seen on 64 bits systems. Split and recompute
1020 G4double fTerm = sd-std::fmod(sd,dRmax);
1021 sd = fTerm + DistanceToIn(p+fTerm*v,v);
1022 }
1023 if ( !fPhiFullCone )
1024 {
1025 xi = p.x() + sd*v.x() ;
1026 yi = p.y() + sd*v.y() ;
1027 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
1028
1029 if (cosPsi >= cosHDPhiOT)
1030 {
1031 if ( sd > halfRadTolerance ) { snxt=sd; }
1032 else
1033 {
1034 // Calculate a normal vector in order to check Direction
1035
1036 risec = std::sqrt(xi*xi + yi*yi)*secRMin ;
1037 Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin);
1038 if ( Normal.dot(v) <= 0 ) { snxt = sd; }
1039 }
1040 }
1041 }
1042 else
1043 {
1044 if( sd > halfRadTolerance ) { return sd; }
1045
1046 // Calculate a normal vector in order to check Direction
1047
1048 xi = p.x() + sd*v.x() ;
1049 yi = p.y() + sd*v.y() ;
1050 risec = std::sqrt(xi*xi + yi*yi)*secRMin ;
1051 Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin) ;
1052 if ( Normal.dot(v) <= 0 ) { return sd; }
1053 }
1054 }
1055 }
1056 else
1057 {
1058 if (b>0) { sd = -b - std::sqrt(d); }
1059 else { sd = c/(-b+std::sqrt(d)); }
1060 zi = p.z() + sd*v.z() ;
1061 ri = rMinAv + zi*tanRMin ;
1062
1063 if ( (sd >= 0) && (ri > 0) && (std::fabs(zi) <= tolODz) ) // sd>0
1064 {
1065 if ( sd>dRmax ) // Avoid rounding errors due to precision issues
1066 { // seen on 64 bits systems. Split and recompute
1067 G4double fTerm = sd-std::fmod(sd,dRmax);
1068 sd = fTerm + DistanceToIn(p+fTerm*v,v);
1069 }
1070 if ( !fPhiFullCone )
1071 {
1072 xi = p.x() + sd*v.x() ;
1073 yi = p.y() + sd*v.y() ;
1074 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
1075
1076 if (cosPsi >= cosHDPhiIT)
1077 {
1078 if ( sd > halfRadTolerance ) { snxt=sd; }
1079 else
1080 {
1081 // Calculate a normal vector in order to check Direction
1082
1083 risec = std::sqrt(xi*xi + yi*yi)*secRMin ;
1084 Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin);
1085 if ( Normal.dot(v) <= 0 ) { snxt = sd; }
1086 }
1087 }
1088 }
1089 else
1090 {
1091 if ( sd > halfRadTolerance ) { return sd; }
1092
1093 // Calculate a normal vector in order to check Direction
1094
1095 xi = p.x() + sd*v.x() ;
1096 yi = p.y() + sd*v.y() ;
1097 risec = std::sqrt(xi*xi + yi*yi)*secRMin ;
1098 Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin) ;
1099 if ( Normal.dot(v) <= 0 ) { return sd; }
1100 }
1101 }
1102 }
1103 }
1104 }
1105 else
1106 {
1107 // Within kRadTol*0.5 of inner cone (real OR imaginary)
1108 // ----> Check not travelling through (=>0 to in)
1109 // ----> if not:
1110 // -2nd root with validity check
1111
1112 if ( std::fabs(p.z()) <= tolODz )
1113 {
1114 if ( nt2 > 0 )
1115 {
1116 // Inside inner real cone, heading outwards, inside z range
1117
1118 if ( !fPhiFullCone )
1119 {
1120 cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/std::sqrt(t3) ;
1121
1122 if (cosPsi >= cosHDPhiIT) { return 0.0; }
1123 }
1124 else { return 0.0; }
1125 }
1126 else
1127 {
1128 // Within z extent, but not travelling through
1129 // -> 2nd root or kInfinity if 1st root on imaginary cone
1130
1131 b = nt2/nt1 ;
1132 c = nt3/nt1 ;
1133 d = b*b - c ;
1134
1135 if ( d >= 0 ) // > 0
1136 {
1137 if (b>0) { sd = -b - std::sqrt(d); }
1138 else { sd = c/(-b+std::sqrt(d)); }
1139 zi = p.z() + sd*v.z() ;
1140 ri = rMinAv + zi*tanRMin ;
1141
1142 if ( ri > 0 ) // 2nd root
1143 {
1144 if (b>0) { sd = c/(-b-std::sqrt(d)); }
1145 else { sd = -b + std::sqrt(d); }
1146
1147 zi = p.z() + sd*v.z() ;
1148
1149 if ( (sd >= 0) && (std::fabs(zi) <= tolODz) ) // sd>0
1150 {
1151 if ( sd>dRmax ) // Avoid rounding errors due to precision issue
1152 { // seen on 64 bits systems. Split and recompute
1153 G4double fTerm = sd-std::fmod(sd,dRmax);
1154 sd = fTerm + DistanceToIn(p+fTerm*v,v);
1155 }
1156 if ( !fPhiFullCone )
1157 {
1158 xi = p.x() + sd*v.x() ;
1159 yi = p.y() + sd*v.y() ;
1160 ri = rMinAv + zi*tanRMin ;
1161 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
1162
1163 if ( cosPsi >= cosHDPhiIT ) { snxt = sd; }
1164 }
1165 else { return sd; }
1166 }
1167 }
1168 else { return kInfinity; }
1169 }
1170 }
1171 }
1172 else // 2nd root
1173 {
1174 b = nt2/nt1 ;
1175 c = nt3/nt1 ;
1176 d = b*b - c ;
1177
1178 if ( d > 0 )
1179 {
1180 if (b>0) { sd = c/(-b-std::sqrt(d)); }
1181 else { sd = -b + std::sqrt(d) ; }
1182 zi = p.z() + sd*v.z() ;
1183
1184 if ( (sd >= 0) && (std::fabs(zi) <= tolODz) ) // sd>0
1185 {
1186 if ( sd>dRmax ) // Avoid rounding errors due to precision issues
1187 { // seen on 64 bits systems. Split and recompute
1188 G4double fTerm = sd-std::fmod(sd,dRmax);
1189 sd = fTerm + DistanceToIn(p+fTerm*v,v);
1190 }
1191 if ( !fPhiFullCone )
1192 {
1193 xi = p.x() + sd*v.x();
1194 yi = p.y() + sd*v.y();
1195 ri = rMinAv + zi*tanRMin ;
1196 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri;
1197
1198 if (cosPsi >= cosHDPhiIT) { snxt = sd; }
1199 }
1200 else { return sd; }
1201 }
1202 }
1203 }
1204 }
1205 }
1206 }
1207
1208 // Phi segment intersection
1209 //
1210 // o Tolerant of points inside phi planes by up to kCarTolerance*0.5
1211 //
1212 // o NOTE: Large duplication of code between sphi & ephi checks
1213 // -> only diffs: sphi -> ephi, Comp -> -Comp and half-plane
1214 // intersection check <=0 -> >=0
1215 // -> Should use some form of loop Construct
1216
1217 if ( !fPhiFullCone )
1218 {
1219 // First phi surface (starting phi)
1220
1221 Comp = v.x()*sinSPhi - v.y()*cosSPhi ;
1222
1223 if ( Comp < 0 ) // Component in outwards normal dirn
1224 {
1225 Dist = (p.y()*cosSPhi - p.x()*sinSPhi) ;
1226
1227 if (Dist < halfCarTolerance)
1228 {
1229 sd = Dist/Comp ;
1230
1231 if ( sd < snxt )
1232 {
1233 if ( sd < 0 ) { sd = 0.0; }
1234
1235 zi = p.z() + sd*v.z() ;
1236
1237 if ( std::fabs(zi) <= tolODz )
1238 {
1239 xi = p.x() + sd*v.x() ;
1240 yi = p.y() + sd*v.y() ;
1241 rhoi2 = xi*xi + yi*yi ;
1242 tolORMin2 = (rMinOAv + zi*tanRMin)*(rMinOAv + zi*tanRMin) ;
1243 tolORMax2 = (rMaxOAv + zi*tanRMax)*(rMaxOAv + zi*tanRMax) ;
1244
1245 if ( (rhoi2 >= tolORMin2) && (rhoi2 <= tolORMax2) )
1246 {
1247 // z and r intersections good - check intersecting with
1248 // correct half-plane
1249
1250 if ((yi*cosCPhi - xi*sinCPhi) <= 0 ) { snxt = sd; }
1251 }
1252 }
1253 }
1254 }
1255 }
1256
1257 // Second phi surface (Ending phi)
1258
1259 Comp = -(v.x()*sinEPhi - v.y()*cosEPhi) ;
1260
1261 if ( Comp < 0 ) // Component in outwards normal dirn
1262 {
1263 Dist = -(p.y()*cosEPhi - p.x()*sinEPhi) ;
1264 if (Dist < halfCarTolerance)
1265 {
1266 sd = Dist/Comp ;
1267
1268 if ( sd < snxt )
1269 {
1270 if ( sd < 0 ) { sd = 0.0; }
1271
1272 zi = p.z() + sd*v.z() ;
1273
1274 if (std::fabs(zi) <= tolODz)
1275 {
1276 xi = p.x() + sd*v.x() ;
1277 yi = p.y() + sd*v.y() ;
1278 rhoi2 = xi*xi + yi*yi ;
1279 tolORMin2 = (rMinOAv + zi*tanRMin)*(rMinOAv + zi*tanRMin) ;
1280 tolORMax2 = (rMaxOAv + zi*tanRMax)*(rMaxOAv + zi*tanRMax) ;
1281
1282 if ( (rhoi2 >= tolORMin2) && (rhoi2 <= tolORMax2) )
1283 {
1284 // z and r intersections good - check intersecting with
1285 // correct half-plane
1286
1287 if ( (yi*cosCPhi - xi*sinCPhi) >= 0.0 ) { snxt = sd; }
1288 }
1289 }
1290 }
1291 }
1292 }
1293 }
1294 if (snxt < halfCarTolerance) { snxt = 0.; }
1295
1296 return snxt ;
1297}
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

Calculates the distance to the nearest surface of a shape from an inside point 'p'. The distance can be an underestimate.

Parameters
[in]pThe point at offset p.
Returns
The safety distance to exit the shape.

Implements G4VSolid.

Definition at line 1986 of file G4Cons.cc.

1987{
1988 G4double safe=0.0, rho, safeR1, safeR2, safeZ, safePhi;
1989 G4double tanRMin, secRMin, pRMin;
1990 G4double tanRMax, secRMax, pRMax;
1991
1992#ifdef G4CSGDEBUG
1993 if( Inside(p) == kOutside )
1994 {
1995 G4int oldprc=G4cout.precision(16) ;
1996 G4cout << G4endl ;
1997 DumpInfo();
1998 G4cout << "Position:" << G4endl << G4endl ;
1999 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl ;
2000 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ;
2001 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ;
2002 G4cout << "pho at z = " << std::sqrt( p.x()*p.x()+p.y()*p.y() )/mm
2003 << " mm" << G4endl << G4endl ;
2004 if( (p.x() != 0.) || (p.x() != 0.) )
2005 {
2006 G4cout << "point phi = " << std::atan2(p.y(),p.x())/degree
2007 << " degrees" << G4endl << G4endl ;
2008 }
2009 G4cout.precision(oldprc) ;
2010 G4Exception("G4Cons::DistanceToOut(p)", "GeomSolids1002",
2011 JustWarning, "Point p is outside !?" );
2012 }
2013#endif
2014
2015 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
2016 safeZ = fDz - std::fabs(p.z()) ;
2017
2018 if ((fRmin1 != 0.0) || (fRmin2 != 0.0))
2019 {
2020 tanRMin = (fRmin2 - fRmin1)*0.5/fDz ;
2021 secRMin = std::sqrt(1.0 + tanRMin*tanRMin) ;
2022 pRMin = tanRMin*p.z() + (fRmin1 + fRmin2)*0.5 ;
2023 safeR1 = (rho - pRMin)/secRMin ;
2024 }
2025 else
2026 {
2027 safeR1 = kInfinity ;
2028 }
2029
2030 tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;
2031 secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
2032 pRMax = tanRMax*p.z() + (fRmax1+fRmax2)*0.5 ;
2033 safeR2 = (pRMax - rho)/secRMax ;
2034
2035 if (safeR1 < safeR2) { safe = safeR1; }
2036 else { safe = safeR2; }
2037 if (safeZ < safe) { safe = safeZ ; }
2038
2039 // Check if phi divided, Calc distances closest phi plane
2040
2041 if (!fPhiFullCone)
2042 {
2043 // Above/below central phi of G4Cons?
2044
2045 if ( (p.y()*cosCPhi - p.x()*sinCPhi) <= 0 )
2046 {
2047 safePhi = -(p.x()*sinSPhi - p.y()*cosSPhi) ;
2048 }
2049 else
2050 {
2051 safePhi = (p.x()*sinEPhi - p.y()*cosEPhi) ;
2052 }
2053 if (safePhi < safe) { safe = safePhi; }
2054 }
2055 if ( safe < 0 ) { safe = 0; }
2056
2057 return safe ;
2058}
G4cout
G4GLOB_DLL std::ostream G4cout
G4Cons::Inside
EInside Inside(const G4ThreeVector &p) const override
Definition G4Cons.cc:175
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

Returns the distance along the normalised vector 'v' to the shape, from a point at an offset 'p' inside or on the surface of the shape. Intersections with surfaces, when the point is less than Tolerance/2 from a surface must be ignored.

Parameters
[in]pThe point at offset p.
[in]vThe normalised direction vector.
[in]calcNormFlag to indicate if to calculate the normal or not.
[out]validNormFlag set to true if the solid lies entirely behind or on the exiting surface. It is set false if the solid does not lie entirely behind or on the exiting surface. 'calcNorm' must be true, otherwise it is unused.
[out]nThe exiting outwards normal vector (undefined Magnitude). 'calcNorm' must be true, otherwise it is unused.
Returns
The distance to exit the shape.

Implements G4VSolid.

Definition at line 1368 of file G4Cons.cc.

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

Returning an estimation of the solid volume (capacity) and surface area, in internal units.

Reimplemented from G4VSolid.

Definition at line 2182 of file G4Cons.cc.

2183{
2184 if (fCubicVolume == 0)
2185 {
2186 G4AutoLock l(&consMutex);
2187 G4double Rmean, rMean, deltaR, deltar;
2188
2189 Rmean = 0.5*(fRmax1+fRmax2);
2190 deltaR = fRmax1-fRmax2;
2191
2192 rMean = 0.5*(fRmin1+fRmin2);
2193 deltar = fRmin1-fRmin2;
2194 fCubicVolume = fDPhi*fDz*(Rmean*Rmean-rMean*rMean
2195 +(deltaR*deltaR-deltar*deltar)/12);
2196 l.unlock();
2197 }
2198 return fCubicVolume;
2199}
G4AutoLock
G4TemplateAutoLock< G4Mutex > G4AutoLock
Definition G4AutoLock.hh:582
G4CSGSolid::fCubicVolume
G4double fCubicVolume
Definition G4CSGSolid.hh:92

◆ GetDeltaPhiAngle()

G4double G4Cons::GetDeltaPhiAngle ( ) const
inline

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

◆ GetEntityType()

G4GeometryType G4Cons::GetEntityType ( ) const
overridevirtual

Returns the type ID, "G4Cons" of the solid.

Implements G4VSolid.

Definition at line 2064 of file G4Cons.cc.

2065{
2066 return {"G4Cons"};
2067}

◆ GetInnerRadiusMinusZ()

G4double G4Cons::GetInnerRadiusMinusZ ( ) const
inline

Accessors.

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

◆ GetInnerRadiusPlusZ()

G4double G4Cons::GetInnerRadiusPlusZ ( ) const
inline

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

◆ GetOuterRadiusMinusZ()

G4double G4Cons::GetOuterRadiusMinusZ ( ) const
inline

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

◆ GetOuterRadiusPlusZ()

G4double G4Cons::GetOuterRadiusPlusZ ( ) const
inline

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

◆ GetPointOnSurface()

G4ThreeVector G4Cons::GetPointOnSurface ( ) const
overridevirtual

Returns a random point located and uniformly distributed on the surface of the solid.

Reimplemented from G4VSolid.

Definition at line 2107 of file G4Cons.cc.

2108{
2109 // declare working variables
2110 //
2111 G4double rone = (fRmax1-fRmax2)/(2.*fDz);
2112 G4double rtwo = (fRmin1-fRmin2)/(2.*fDz);
2113 G4double qone = (fRmax1 == fRmax2) ? 0. : fDz*(fRmax1+fRmax2)/(fRmax1-fRmax2);
2114 G4double qtwo = (fRmin1 == fRmin2) ? 0. : fDz*(fRmin1+fRmin2)/(fRmin1-fRmin2);
2115
2116 G4double slin = std::hypot(fRmin1-fRmin2, 2.*fDz);
2117 G4double slout = std::hypot(fRmax1-fRmax2, 2.*fDz);
2118 G4double Aone = 0.5*fDPhi*(fRmax2 + fRmax1)*slout; // outer surface
2119 G4double Atwo = 0.5*fDPhi*(fRmin2 + fRmin1)*slin; // inner surface
2120 G4double Athree = 0.5*fDPhi*(fRmax1*fRmax1-fRmin1*fRmin1); // base at -Dz
2121 G4double Afour = 0.5*fDPhi*(fRmax2*fRmax2-fRmin2*fRmin2); // base at +Dz
2122 G4double Afive = fDz*(fRmax1-fRmin1+fRmax2-fRmin2); // phi section
2123
2124 G4double phi = G4RandFlat::shoot(fSPhi,fSPhi+fDPhi);
2125 G4double cosu = std::cos(phi);
2126 G4double sinu = std::sin(phi);
2127 G4double rRand1 = GetRadiusInRing(fRmin1, fRmax1);
2128 G4double rRand2 = GetRadiusInRing(fRmin2, fRmax2);
2129
2130 if ( (fSPhi == 0.) && fPhiFullCone ) { Afive = 0.; }
2131 G4double chose = G4RandFlat::shoot(0.,Aone+Atwo+Athree+Afour+2.*Afive);
2132
2133 if( (chose >= 0.) && (chose < Aone) ) // outer surface
2134 {
2135 if(fRmax1 != fRmax2)
2136 {
2137 G4double zRand = G4RandFlat::shoot(-1.*fDz,fDz);
2138 return { rone*cosu*(qone-zRand), rone*sinu*(qone-zRand), zRand };
2139 }
2140
2141 return { fRmax1*cosu, fRmax2*sinu, G4RandFlat::shoot(-1.*fDz,fDz) };
2142 }
2143 if( (chose >= Aone) && (chose < Aone + Atwo) ) // inner surface
2144 {
2145 if(fRmin1 != fRmin2)
2146 {
2147 G4double zRand = G4RandFlat::shoot(-1.*fDz,fDz);
2148 return { rtwo*cosu*(qtwo-zRand), rtwo*sinu*(qtwo-zRand), zRand };
2149 }
2150
2151 return { fRmin1*cosu, fRmin2*sinu, G4RandFlat::shoot(-1.*fDz,fDz) };
2152 }
2153 if( (chose >= Aone + Atwo) && (chose < Aone + Atwo + Athree) ) // base at -Dz
2154 {
2155 return {rRand1*cosu, rRand1*sinu, -1*fDz};
2156 }
2157 if( (chose >= Aone + Atwo + Athree)
2158 && (chose < Aone + Atwo + Athree + Afour) ) // base at +Dz
2159 {
2160 return { rRand2*cosu, rRand2*sinu, fDz };
2161 }
2162 if( (chose >= Aone + Atwo + Athree + Afour) // SPhi section
2163 && (chose < Aone + Atwo + Athree + Afour + Afive) )
2164 {
2165 G4double zRand = G4RandFlat::shoot(-1.*fDz,fDz);
2166 rRand1 = G4RandFlat::shoot(fRmin2-((zRand-fDz)/(2.*fDz))*(fRmin1-fRmin2),
2167 fRmax2-((zRand-fDz)/(2.*fDz))*(fRmax1-fRmax2));
2168 return { rRand1*cosSPhi, rRand1*sinSPhi, zRand };
2169 }
2170
2171 // SPhi+DPhi section
2172 G4double zRand = G4RandFlat::shoot(-1.*fDz,fDz);
2173 rRand1 = G4RandFlat::shoot(fRmin2-((zRand-fDz)/(2.*fDz))*(fRmin1-fRmin2),
2174 fRmax2-((zRand-fDz)/(2.*fDz))*(fRmax1-fRmax2));
2175 return { rRand1*cosEPhi, rRand1*sinEPhi, zRand };
2176}
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 ( )
overridevirtual

Returns an estimation of the solid surface area in internal units. This method may be overloaded by derived classes to compute the exact geometrical quantity for solids where this is possible, or anyway to cache the computed value. Note: the computed value is NOT cached.

Reimplemented from G4VSolid.

Definition at line 2205 of file G4Cons.cc.

2206{
2207 if (fSurfaceArea == 0)
2208 {
2209 G4AutoLock l(&consMutex);
2210 G4double mmin, mmax, dmin, dmax;
2211
2212 mmin= (fRmin1+fRmin2)*0.5;
2213 mmax= (fRmax1+fRmax2)*0.5;
2214 dmin= (fRmin2-fRmin1);
2215 dmax= (fRmax2-fRmax1);
2216
2217 fSurfaceArea = fDPhi*( mmin * std::sqrt(dmin*dmin+4*fDz*fDz)
2218 + mmax * std::sqrt(dmax*dmax+4*fDz*fDz)
2219 + 0.5*(fRmax1*fRmax1-fRmin1*fRmin1
2220 +fRmax2*fRmax2-fRmin2*fRmin2 ));
2221 if(!fPhiFullCone)
2222 {
2223 fSurfaceArea = fSurfaceArea+4*fDz*(mmax-mmin);
2224 }
2225 l.unlock();
2226 }
2227 return fSurfaceArea;
2228}
G4CSGSolid::fSurfaceArea
G4double fSurfaceArea
Definition G4CSGSolid.hh:93

◆ GetZHalfLength()

G4double G4Cons::GetZHalfLength ( ) const
inline

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

◆ Inside()

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

Concrete implementations of the expected query interfaces for solids, as defined in the base class G4VSolid.

Implements G4VSolid.

Definition at line 175 of file G4Cons.cc.

176{
177 G4double r2, rl, rh, pPhi, tolRMin, tolRMax; // rh2, rl2 ;
178 EInside in;
179
180 if (std::fabs(p.z()) > fDz + halfCarTolerance ) { return in = kOutside; }
181 if(std::fabs(p.z()) >= fDz - halfCarTolerance ) { in = kSurface; }
182 else { in = kInside; }
183
184 r2 = p.x()*p.x() + p.y()*p.y() ;
185 rl = 0.5*(fRmin2*(p.z() + fDz) + fRmin1*(fDz - p.z()))/fDz ;
186 rh = 0.5*(fRmax2*(p.z()+fDz)+fRmax1*(fDz-p.z()))/fDz;
187
188 // rh2 = rh*rh;
189
190 tolRMin = rl - halfRadTolerance;
191 if ( tolRMin < 0 ) { tolRMin = 0; }
192 tolRMax = rh + halfRadTolerance;
193
194 if ( (r2<tolRMin*tolRMin) || (r2>tolRMax*tolRMax) ) { return in = kOutside; }
195
196 if (rl != 0.0) { tolRMin = rl + halfRadTolerance; }
197 else { tolRMin = 0.0; }
198 tolRMax = rh - halfRadTolerance;
199
200 if (in == kInside) // else it's kSurface already
201 {
202 if ( (r2 < tolRMin*tolRMin) || (r2 >= tolRMax*tolRMax) ) { in = kSurface; }
203 }
204 if ( !fPhiFullCone && ((p.x() != 0.0) || (p.y() != 0.0)) )
205 {
206 pPhi = std::atan2(p.y(),p.x()) ;
207
208 if ( pPhi < fSPhi - halfAngTolerance ) { pPhi += twopi; }
209 else if ( pPhi > fSPhi + fDPhi + halfAngTolerance ) { pPhi -= twopi; }
210
211 if ( (pPhi < fSPhi - halfAngTolerance) ||
212 (pPhi > fSPhi + fDPhi + halfAngTolerance) ) { return in = kOutside; }
213
214 if (in == kInside) // else it's kSurface anyway already
215 {
216 if ( (pPhi < fSPhi + halfAngTolerance) ||
217 (pPhi > fSPhi + fDPhi - halfAngTolerance) ) { in = kSurface; }
218 }
219 }
220 else if ( !fPhiFullCone ) { in = kSurface; }
221
222 return in ;
223}
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 142 of file G4Cons.cc.

143{
144 // Check assignment to self
145 //
146 if (this == &rhs) { return *this; }
147
148 // Copy base class data
149 //
150 G4CSGSolid::operator=(rhs);
151
152 // Copy data
153 //
154 kRadTolerance = rhs.kRadTolerance;
155 kAngTolerance = rhs.kAngTolerance;
156 fRmin1 = rhs.fRmin1; fRmin2 = rhs.fRmin2;
157 fRmax1 = rhs.fRmax1; fRmax2 = rhs.fRmax2;
158 fDz = rhs.fDz; fSPhi = rhs.fSPhi; fDPhi = rhs.fDPhi;
159 sinCPhi = rhs.sinCPhi; cosCPhi = rhs.cosCPhi; cosHDPhi = rhs.cosHDPhi;
160 cosHDPhiOT = rhs.cosHDPhiOT; cosHDPhiIT = rhs.cosHDPhiIT;
161 sinSPhi = rhs.sinSPhi; cosSPhi = rhs.cosSPhi;
162 sinEPhi = rhs.sinEPhi; cosEPhi = rhs.cosEPhi;
163 fPhiFullCone = rhs.fPhiFullCone;
164 halfCarTolerance = rhs.halfCarTolerance;
165 halfRadTolerance = rhs.halfRadTolerance;
166 halfAngTolerance = rhs.halfAngTolerance;
167
168 return *this;
169}
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

Modifiers.

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

Streams the object contents to an output stream.

Implements G4VSolid.

Definition at line 2082 of file G4Cons.cc.

2083{
2084 G4long oldprc = os.precision(16);
2085 os << "-----------------------------------------------------------\n"
2086 << " *** Dump for solid - " << GetName() << " ***\n"
2087 << " ===================================================\n"
2088 << " Solid type: G4Cons\n"
2089 << " Parameters: \n"
2090 << " inside -fDz radius: " << fRmin1/mm << " mm \n"
2091 << " outside -fDz radius: " << fRmax1/mm << " mm \n"
2092 << " inside +fDz radius: " << fRmin2/mm << " mm \n"
2093 << " outside +fDz radius: " << fRmax2/mm << " mm \n"
2094 << " half length in Z : " << fDz/mm << " mm \n"
2095 << " starting angle of segment: " << fSPhi/degree << " degrees \n"
2096 << " delta angle of segment : " << fDPhi/degree << " degrees \n"
2097 << "-----------------------------------------------------------\n";
2098 os.precision(oldprc);
2099
2100 return os;
2101}

◆ SurfaceNormal()

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

Returns the outwards pointing unit normal of the shape for the surface closest to the point at offset 'p'.

Parameters
[in]pThe point at offset p.
Returns
The outwards pointing unit normal.

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: