Geant4 10.7.0
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 ()
 
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 ()
 
G4double GetSurfaceArea ()
 
void ComputeDimensions (G4VPVParameterisation *p, const G4int n, const G4VPhysicalVolume *pRep)
 
void BoundingLimits (G4ThreeVector &pMin, G4ThreeVector &pMax) const
 
G4bool CalculateExtent (const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pMin, G4double &pMax) const
 
EInside Inside (const G4ThreeVector &p) const
 
G4ThreeVector SurfaceNormal (const G4ThreeVector &p) const
 
G4double DistanceToIn (const G4ThreeVector &p, const G4ThreeVector &v) const
 
G4double DistanceToIn (const G4ThreeVector &p) const
 
G4double DistanceToOut (const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=false, G4bool *validNorm=nullptr, G4ThreeVector *n=nullptr) const
 
G4double DistanceToOut (const G4ThreeVector &p) const
 
G4GeometryType GetEntityType () const
 
G4ThreeVector GetPointOnSurface () const
 
G4VSolid * Clone () const
 
std::ostream & StreamInfo (std::ostream &os) const
 
void DescribeYourselfTo (G4VGraphicsScene &scene) const
 
G4Polyhedron * CreatePolyhedron () const
 
 G4Cons (__void__ &)
 
 G4Cons (const G4Cons &rhs)
 
G4Cons & operator= (const G4Cons &rhs)
 
G4double GetRmin1 () const
 
G4double GetRmax1 () const
 
G4double GetRmin2 () const
 
G4double GetRmax2 () const
 
G4double GetDz () const
 
G4double GetSPhi () const
 
G4double GetDPhi () const
 
- Public Member Functions inherited from G4CSGSolid
 G4CSGSolid (const G4String &pName)
 
virtual ~G4CSGSolid ()
 
virtual std::ostream & StreamInfo (std::ostream &os) const
 
virtual G4Polyhedron * GetPolyhedron () const
 
 G4CSGSolid (__void__ &)
 
 G4CSGSolid (const G4CSGSolid &rhs)
 
G4CSGSolid & operator= (const G4CSGSolid &rhs)
 
- Public Member Functions inherited from G4VSolid
 G4VSolid (const G4String &name)
 
virtual ~G4VSolid ()
 
G4bool operator== (const G4VSolid &s) const
 
G4String GetName () const
 
void SetName (const G4String &name)
 
G4double GetTolerance () const
 
virtual void BoundingLimits (G4ThreeVector &pMin, G4ThreeVector &pMax) const
 
virtual G4bool CalculateExtent (const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pMin, G4double &pMax) const =0
 
virtual EInside Inside (const G4ThreeVector &p) const =0
 
virtual G4ThreeVector SurfaceNormal (const G4ThreeVector &p) const =0
 
virtual G4double DistanceToIn (const G4ThreeVector &p, const G4ThreeVector &v) const =0
 
virtual G4double DistanceToIn (const G4ThreeVector &p) const =0
 
virtual G4double DistanceToOut (const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=false, G4bool *validNorm=nullptr, G4ThreeVector *n=nullptr) const =0
 
virtual G4double DistanceToOut (const G4ThreeVector &p) const =0
 
virtual void ComputeDimensions (G4VPVParameterisation *p, const G4int n, const G4VPhysicalVolume *pRep)
 
virtual G4double GetCubicVolume ()
 
virtual G4double GetSurfaceArea ()
 
virtual G4GeometryType GetEntityType () const =0
 
virtual G4ThreeVector GetPointOnSurface () const
 
virtual G4VSolid * Clone () const
 
virtual std::ostream & StreamInfo (std::ostream &os) const =0
 
void DumpInfo () const
 
virtual void DescribeYourselfTo (G4VGraphicsScene &scene) const =0
 
virtual G4VisExtent GetExtent () const
 
virtual G4Polyhedron * CreatePolyhedron () const
 
virtual G4Polyhedron * GetPolyhedron () const
 
virtual const G4VSolid * GetConstituentSolid (G4int no) const
 
virtual G4VSolid * GetConstituentSolid (G4int no)
 
virtual const G4DisplacedSolid * GetDisplacedSolidPtr () const
 
virtual G4DisplacedSolid * GetDisplacedSolidPtr ()
 
 G4VSolid (__void__ &)
 
 G4VSolid (const G4VSolid &rhs)
 
G4VSolid & operator= (const G4VSolid &rhs)
 
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 72 of file G4Cons.cc.

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

◆ ~G4Cons()

G4Cons::~G4Cons ( )

Definition at line 136 of file G4Cons.cc.

137{
138}

◆ G4Cons() [2/3]

G4Cons::G4Cons ( __void__ &  a)

Definition at line 122 of file G4Cons.cc.

123 : G4CSGSolid(a), kRadTolerance(0.), kAngTolerance(0.),
124 fRmin1(0.), fRmin2(0.), fRmax1(0.), fRmax2(0.), fDz(0.),
125 fSPhi(0.), fDPhi(0.), sinCPhi(0.), cosCPhi(0.), cosHDPhi(0.),
126 cosHDPhiOT(0.), cosHDPhiIT(0.), sinSPhi(0.), cosSPhi(0.),
127 sinEPhi(0.), cosEPhi(0.),
128 halfCarTolerance(0.), halfRadTolerance(0.), halfAngTolerance(0.)
129{
130}

◆ G4Cons() [3/3]

G4Cons::G4Cons ( const G4Cons &  rhs)

Definition at line 144 of file G4Cons.cc.

145 : G4CSGSolid(rhs), kRadTolerance(rhs.kRadTolerance),
146 kAngTolerance(rhs.kAngTolerance), fRmin1(rhs.fRmin1), fRmin2(rhs.fRmin2),
147 fRmax1(rhs.fRmax1), fRmax2(rhs.fRmax2), fDz(rhs.fDz), fSPhi(rhs.fSPhi),
148 fDPhi(rhs.fDPhi), sinCPhi(rhs.sinCPhi), cosCPhi(rhs.cosCPhi),
149 cosHDPhi(rhs.cosHDPhi), cosHDPhiOT(rhs.cosHDPhiOT),
150 cosHDPhiIT(rhs.cosHDPhiIT), sinSPhi(rhs.sinSPhi), cosSPhi(rhs.cosSPhi),
151 sinEPhi(rhs.sinEPhi), cosEPhi(rhs.cosEPhi), fPhiFullCone(rhs.fPhiFullCone),
152 halfCarTolerance(rhs.halfCarTolerance),
153 halfRadTolerance(rhs.halfRadTolerance),
154 halfAngTolerance(rhs.halfAngTolerance)
155{
156}

Member Function Documentation

◆ BoundingLimits()

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

Reimplemented from G4VSolid.

Definition at line 261 of file G4Cons.cc.

262{
263 G4double rmin = std::min(GetInnerRadiusMinusZ(),GetInnerRadiusPlusZ());
264 G4double rmax = std::max(GetOuterRadiusMinusZ(),GetOuterRadiusPlusZ());
265 G4double dz = GetZHalfLength();
266
267 // Find bounding box
268 //
269 if (GetDeltaPhiAngle() < twopi)
270 {
271 G4TwoVector vmin,vmax;
272 G4GeomTools::DiskExtent(rmin,rmax,
273 GetSinStartPhi(),GetCosStartPhi(),
274 GetSinEndPhi(),GetCosEndPhi(),
275 vmin,vmax);
276 pMin.set(vmin.x(),vmin.y(),-dz);
277 pMax.set(vmax.x(),vmax.y(), dz);
278 }
279 else
280 {
281 pMin.set(-rmax,-rmax,-dz);
282 pMax.set( rmax, rmax, dz);
283 }
284
285 // Check correctness of the bounding box
286 //
287 if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
288 {
289 std::ostringstream message;
290 message << "Bad bounding box (min >= max) for solid: "
291 << GetName() << " !"
292 << "\npMin = " << pMin
293 << "\npMax = " << pMax;
294 G4Exception("G4Cons::BoundingLimits()", "GeomMgt0001",
295 JustWarning, message);
296 DumpInfo();
297 }
298}
JustWarning
@ JustWarning
Definition: G4ExceptionSeverity.hh:65
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
virtual

Implements G4VSolid.

Definition at line 304 of file G4Cons.cc.

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

◆ Clone()

G4VSolid * G4Cons::Clone ( ) const
virtual

Reimplemented from G4VSolid.

Definition at line 2103 of file G4Cons.cc.

2104{
2105 return new G4Cons(*this);
2106}
G4Cons
Definition: G4Cons.hh:78

◆ ComputeDimensions()

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

Reimplemented from G4VSolid.

Definition at line 250 of file G4Cons.cc.

253{
254 p->ComputeDimensions(*this,n,pRep) ;
255}
G4VPVParameterisation::ComputeDimensions
virtual void ComputeDimensions(G4Box &, const G4int, const G4VPhysicalVolume *) const
Definition: G4VPVParameterisation.hh:91

◆ CreatePolyhedron()

G4Polyhedron * G4Cons::CreatePolyhedron ( ) const
virtual

Reimplemented from G4VSolid.

Definition at line 2230 of file G4Cons.cc.

2231{
2232 return new G4PolyhedronCons(fRmin1,fRmax1,fRmin2,fRmax2,fDz,fSPhi,fDPhi);
2233}

◆ DescribeYourselfTo()

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

Implements G4VSolid.

Definition at line 2225 of file G4Cons.cc.

2226{
2227 scene.AddSolid (*this);
2228}
G4VGraphicsScene::AddSolid
virtual void AddSolid(const G4Box &)=0

◆ DistanceToIn() [1/2]

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

Implements G4VSolid.

Definition at line 1331 of file G4Cons.cc.

1332{
1333 G4double safe=0.0, rho, safeR1, safeR2, safeZ, safePhi, cosPsi ;
1334 G4double tanRMin, secRMin, pRMin ;
1335 G4double tanRMax, secRMax, pRMax ;
1336
1337 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
1338 safeZ = std::fabs(p.z()) - fDz ;
1339
1340 if ( fRmin1 || fRmin2 )
1341 {
1342 tanRMin = (fRmin2 - fRmin1)*0.5/fDz ;
1343 secRMin = std::sqrt(1.0 + tanRMin*tanRMin) ;
1344 pRMin = tanRMin*p.z() + (fRmin1 + fRmin2)*0.5 ;
1345 safeR1 = (pRMin - rho)/secRMin ;
1346
1347 tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;
1348 secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
1349 pRMax = tanRMax*p.z() + (fRmax1 + fRmax2)*0.5 ;
1350 safeR2 = (rho - pRMax)/secRMax ;
1351
1352 if ( safeR1 > safeR2) { safe = safeR1; }
1353 else { safe = safeR2; }
1354 }
1355 else
1356 {
1357 tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;
1358 secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
1359 pRMax = tanRMax*p.z() + (fRmax1 + fRmax2)*0.5 ;
1360 safe = (rho - pRMax)/secRMax ;
1361 }
1362 if ( safeZ > safe ) { safe = safeZ; }
1363
1364 if ( !fPhiFullCone && rho )
1365 {
1366 // Psi=angle from central phi to point
1367
1368 cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/rho ;
1369
1370 if ( cosPsi < cosHDPhi ) // Point lies outside phi range
1371 {
1372 if ( (p.y()*cosCPhi - p.x()*sinCPhi) <= 0.0 )
1373 {
1374 safePhi = std::fabs(p.x()*sinSPhi-p.y()*cosSPhi);
1375 }
1376 else
1377 {
1378 safePhi = std::fabs(p.x()*sinEPhi-p.y()*cosEPhi);
1379 }
1380 if ( safePhi > safe ) { safe = safePhi; }
1381 }
1382 }
1383 if ( safe < 0.0 ) { safe = 0.0; }
1384
1385 return safe ;
1386}

◆ DistanceToIn() [2/2]

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

Implements G4VSolid.

Definition at line 661 of file G4Cons.cc.

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

Referenced by DistanceToIn().

◆ DistanceToOut() [1/2]

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

Implements G4VSolid.

Definition at line 2016 of file G4Cons.cc.

2017{
2018 G4double safe=0.0, rho, safeR1, safeR2, safeZ, safePhi;
2019 G4double tanRMin, secRMin, pRMin;
2020 G4double tanRMax, secRMax, pRMax;
2021
2022#ifdef G4CSGDEBUG
2023 if( Inside(p) == kOutside )
2024 {
2025 G4int oldprc=G4cout.precision(16) ;
2026 G4cout << G4endl ;
2027 DumpInfo();
2028 G4cout << "Position:" << G4endl << G4endl ;
2029 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl ;
2030 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ;
2031 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ;
2032 G4cout << "pho at z = " << std::sqrt( p.x()*p.x()+p.y()*p.y() )/mm
2033 << " mm" << G4endl << G4endl ;
2034 if( (p.x() != 0.) || (p.x() != 0.) )
2035 {
2036 G4cout << "point phi = " << std::atan2(p.y(),p.x())/degree
2037 << " degree" << G4endl << G4endl ;
2038 }
2039 G4cout.precision(oldprc) ;
2040 G4Exception("G4Cons::DistanceToOut(p)", "GeomSolids1002",
2041 JustWarning, "Point p is outside !?" );
2042 }
2043#endif
2044
2045 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
2046 safeZ = fDz - std::fabs(p.z()) ;
2047
2048 if (fRmin1 || fRmin2)
2049 {
2050 tanRMin = (fRmin2 - fRmin1)*0.5/fDz ;
2051 secRMin = std::sqrt(1.0 + tanRMin*tanRMin) ;
2052 pRMin = tanRMin*p.z() + (fRmin1 + fRmin2)*0.5 ;
2053 safeR1 = (rho - pRMin)/secRMin ;
2054 }
2055 else
2056 {
2057 safeR1 = kInfinity ;
2058 }
2059
2060 tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;
2061 secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
2062 pRMax = tanRMax*p.z() + (fRmax1+fRmax2)*0.5 ;
2063 safeR2 = (pRMax - rho)/secRMax ;
2064
2065 if (safeR1 < safeR2) { safe = safeR1; }
2066 else { safe = safeR2; }
2067 if (safeZ < safe) { safe = safeZ ; }
2068
2069 // Check if phi divided, Calc distances closest phi plane
2070
2071 if (!fPhiFullCone)
2072 {
2073 // Above/below central phi of G4Cons?
2074
2075 if ( (p.y()*cosCPhi - p.x()*sinCPhi) <= 0 )
2076 {
2077 safePhi = -(p.x()*sinSPhi - p.y()*cosSPhi) ;
2078 }
2079 else
2080 {
2081 safePhi = (p.x()*sinEPhi - p.y()*cosEPhi) ;
2082 }
2083 if (safePhi < safe) { safe = safePhi; }
2084 }
2085 if ( safe < 0 ) { safe = 0; }
2086
2087 return safe ;
2088}
G4cout
G4GLOB_DLL std::ostream G4cout
G4Cons::Inside
EInside Inside(const G4ThreeVector &p) const
Definition: G4Cons.cc:195
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
virtual

Implements G4VSolid.

Definition at line 1393 of file G4Cons.cc.

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

Reimplemented from G4VSolid.

◆ GetDeltaPhiAngle()

G4double G4Cons::GetDeltaPhiAngle ( ) const
inline

Referenced by G4HepRepFileSceneHandler::AddSolid(), G4HepRepSceneHandler::AddSolid(), BoundingLimits(), G4tgbVolume::BuildSolidForDivision(), CalculateExtent(), G4ParameterisationConsRho::ComputeDimensions(), G4ParameterisationConsZ::ComputeDimensions(), G4GDMLWriteParamvol::Cone_dimensionsWrite(), G4GDMLWriteSolids::ConeWrite(), G4ParameterisationConsPhi::G4ParameterisationConsPhi(), G4VParameterisationCons::G4VParameterisationCons(), G4ParameterisationConsPhi::GetMaxParameter(), and G4tgbGeometryDumper::GetSolidParams().

◆ GetDPhi()

G4double G4Cons::GetDPhi ( ) const
inline

◆ GetDz()

G4double G4Cons::GetDz ( ) const
inline

◆ GetEntityType()

G4GeometryType G4Cons::GetEntityType ( ) const
virtual

Implements G4VSolid.

Definition at line 2094 of file G4Cons.cc.

2095{
2096 return G4String("G4Cons");
2097}

◆ GetInnerRadiusMinusZ()

G4double G4Cons::GetInnerRadiusMinusZ ( ) const
inline

Referenced by G4HepRepFileSceneHandler::AddSolid(), G4HepRepSceneHandler::AddSolid(), BoundingLimits(), G4tgbVolume::BuildSolidForDivision(), CalculateExtent(), G4ParameterisationConsRho::ComputeDimensions(), G4ParameterisationConsPhi::ComputeDimensions(), G4ParameterisationConsZ::ComputeDimensions(), G4GDMLWriteParamvol::Cone_dimensionsWrite(), G4GDMLWriteSolids::ConeWrite(), G4ParameterisationConsRho::G4ParameterisationConsRho(), G4VParameterisationCons::G4VParameterisationCons(), G4ParameterisationConsRho::GetMaxParameter(), and G4tgbGeometryDumper::GetSolidParams().

◆ GetInnerRadiusPlusZ()

G4double G4Cons::GetInnerRadiusPlusZ ( ) const
inline

Referenced by G4HepRepFileSceneHandler::AddSolid(), G4HepRepSceneHandler::AddSolid(), BoundingLimits(), G4tgbVolume::BuildSolidForDivision(), CalculateExtent(), G4ParameterisationConsRho::ComputeDimensions(), G4ParameterisationConsPhi::ComputeDimensions(), G4ParameterisationConsZ::ComputeDimensions(), G4GDMLWriteParamvol::Cone_dimensionsWrite(), G4GDMLWriteSolids::ConeWrite(), G4ParameterisationConsRho::G4ParameterisationConsRho(), G4VParameterisationCons::G4VParameterisationCons(), and G4tgbGeometryDumper::GetSolidParams().

◆ GetOuterRadiusMinusZ()

G4double G4Cons::GetOuterRadiusMinusZ ( ) const
inline

Referenced by G4HepRepFileSceneHandler::AddSolid(), G4HepRepSceneHandler::AddSolid(), BoundingLimits(), G4tgbVolume::BuildSolidForDivision(), CalculateExtent(), G4ParameterisationConsPhi::ComputeDimensions(), G4ParameterisationConsZ::ComputeDimensions(), G4GDMLWriteParamvol::Cone_dimensionsWrite(), G4GDMLWriteSolids::ConeWrite(), G4ParameterisationConsRho::G4ParameterisationConsRho(), G4VParameterisationCons::G4VParameterisationCons(), G4ParameterisationConsRho::GetMaxParameter(), and G4tgbGeometryDumper::GetSolidParams().

◆ GetOuterRadiusPlusZ()

G4double G4Cons::GetOuterRadiusPlusZ ( ) const
inline

Referenced by G4HepRepFileSceneHandler::AddSolid(), G4HepRepSceneHandler::AddSolid(), BoundingLimits(), G4tgbVolume::BuildSolidForDivision(), CalculateExtent(), G4ParameterisationConsRho::ComputeDimensions(), G4ParameterisationConsPhi::ComputeDimensions(), G4ParameterisationConsZ::ComputeDimensions(), G4GDMLWriteParamvol::Cone_dimensionsWrite(), G4GDMLWriteSolids::ConeWrite(), G4VParameterisationCons::G4VParameterisationCons(), and G4tgbGeometryDumper::GetSolidParams().

◆ GetPointOnSurface()

G4ThreeVector G4Cons::GetPointOnSurface ( ) const
virtual

Reimplemented from G4VSolid.

Definition at line 2139 of file G4Cons.cc.

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

◆ GetRmax1()

G4double G4Cons::GetRmax1 ( ) const
inline

◆ GetRmax2()

G4double G4Cons::GetRmax2 ( ) const
inline

◆ GetRmin1()

G4double G4Cons::GetRmin1 ( ) const
inline

◆ GetRmin2()

G4double G4Cons::GetRmin2 ( ) const
inline

◆ GetSinEndPhi()

G4double G4Cons::GetSinEndPhi ( ) const
inline

Referenced by BoundingLimits(), and CalculateExtent().

◆ GetSinStartPhi()

G4double G4Cons::GetSinStartPhi ( ) const
inline

Referenced by BoundingLimits(), and CalculateExtent().

◆ GetSPhi()

G4double G4Cons::GetSPhi ( ) const
inline

◆ GetStartPhiAngle()

G4double G4Cons::GetStartPhiAngle ( ) const
inline

Referenced by G4tgbVolume::BuildSolidForDivision(), G4ParameterisationConsRho::ComputeDimensions(), G4ParameterisationConsPhi::ComputeDimensions(), G4ParameterisationConsZ::ComputeDimensions(), G4GDMLWriteParamvol::Cone_dimensionsWrite(), G4GDMLWriteSolids::ConeWrite(), G4VParameterisationCons::G4VParameterisationCons(), and G4tgbGeometryDumper::GetSolidParams().

◆ GetSurfaceArea()

G4double G4Cons::GetSurfaceArea ( )
inlinevirtual

Reimplemented from G4VSolid.

◆ GetZHalfLength()

G4double G4Cons::GetZHalfLength ( ) const
inline

Referenced by G4HepRepFileSceneHandler::AddSolid(), G4HepRepSceneHandler::AddSolid(), BoundingLimits(), G4tgbVolume::BuildSolidForDivision(), CalculateExtent(), G4ParameterisationConsRho::ComputeDimensions(), G4ParameterisationConsPhi::ComputeDimensions(), G4ParameterisationConsZ::ComputeDimensions(), G4ParameterisationConsZ::ComputeTransformation(), G4GDMLWriteParamvol::Cone_dimensionsWrite(), G4GDMLWriteSolids::ConeWrite(), G4ParameterisationConsZ::G4ParameterisationConsZ(), G4VParameterisationCons::G4VParameterisationCons(), G4ParameterisationConsZ::GetMaxParameter(), and G4tgbGeometryDumper::GetSolidParams().

◆ Inside()

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

Implements G4VSolid.

Definition at line 195 of file G4Cons.cc.

196{
197 G4double r2, rl, rh, pPhi, tolRMin, tolRMax; // rh2, rl2 ;
198 EInside in;
199
200 if (std::fabs(p.z()) > fDz + halfCarTolerance ) { return in = kOutside; }
201 else if(std::fabs(p.z()) >= fDz - halfCarTolerance ) { in = kSurface; }
202 else { in = kInside; }
203
204 r2 = p.x()*p.x() + p.y()*p.y() ;
205 rl = 0.5*(fRmin2*(p.z() + fDz) + fRmin1*(fDz - p.z()))/fDz ;
206 rh = 0.5*(fRmax2*(p.z()+fDz)+fRmax1*(fDz-p.z()))/fDz;
207
208 // rh2 = rh*rh;
209
210 tolRMin = rl - halfRadTolerance;
211 if ( tolRMin < 0 ) { tolRMin = 0; }
212 tolRMax = rh + halfRadTolerance;
213
214 if ( (r2<tolRMin*tolRMin) || (r2>tolRMax*tolRMax) ) { return in = kOutside; }
215
216 if (rl) { tolRMin = rl + halfRadTolerance; }
217 else { tolRMin = 0.0; }
218 tolRMax = rh - halfRadTolerance;
219
220 if (in == kInside) // else it's kSurface already
221 {
222 if ( (r2 < tolRMin*tolRMin) || (r2 >= tolRMax*tolRMax) ) { in = kSurface; }
223 }
224 if ( !fPhiFullCone && ((p.x() != 0.0) || (p.y() != 0.0)) )
225 {
226 pPhi = std::atan2(p.y(),p.x()) ;
227
228 if ( pPhi < fSPhi - halfAngTolerance ) { pPhi += twopi; }
229 else if ( pPhi > fSPhi + fDPhi + halfAngTolerance ) { pPhi -= twopi; }
230
231 if ( (pPhi < fSPhi - halfAngTolerance) ||
232 (pPhi > fSPhi + fDPhi + halfAngTolerance) ) { return in = kOutside; }
233
234 else if (in == kInside) // else it's kSurface anyway already
235 {
236 if ( (pPhi < fSPhi + halfAngTolerance) ||
237 (pPhi > fSPhi + fDPhi - halfAngTolerance) ) { in = kSurface; }
238 }
239 }
240 else if ( !fPhiFullCone ) { in = kSurface; }
241
242 return in ;
243}
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 162 of file G4Cons.cc.

163{
164 // Check assignment to self
165 //
166 if (this == &rhs) { return *this; }
167
168 // Copy base class data
169 //
170 G4CSGSolid::operator=(rhs);
171
172 // Copy data
173 //
174 kRadTolerance = rhs.kRadTolerance;
175 kAngTolerance = rhs.kAngTolerance;
176 fRmin1 = rhs.fRmin1; fRmin2 = rhs.fRmin2;
177 fRmax1 = rhs.fRmax1; fRmax2 = rhs.fRmax2;
178 fDz = rhs.fDz; fSPhi = rhs.fSPhi; fDPhi = rhs.fDPhi;
179 sinCPhi = rhs.sinCPhi; cosCPhi = rhs.cosCPhi; cosHDPhi = rhs.cosHDPhi;
180 cosHDPhiOT = rhs.cosHDPhiOT; cosHDPhiIT = rhs.cosHDPhiIT;
181 sinSPhi = rhs.sinSPhi; cosSPhi = rhs.cosSPhi;
182 sinEPhi = rhs.sinEPhi; cosEPhi = rhs.cosEPhi;
183 fPhiFullCone = rhs.fPhiFullCone;
184 halfCarTolerance = rhs.halfCarTolerance;
185 halfRadTolerance = rhs.halfRadTolerance;
186 halfAngTolerance = rhs.halfAngTolerance;
187
188 return *this;
189}
G4CSGSolid::operator=
G4CSGSolid & operator=(const G4CSGSolid &rhs)
Definition: G4CSGSolid.cc:89

◆ SetDeltaPhiAngle()

void G4Cons::SetDeltaPhiAngle ( G4double  newDPhi)
inline

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

◆ SetInnerRadiusMinusZ()

void G4Cons::SetInnerRadiusMinusZ ( G4double  Rmin1)
inline

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

◆ SetInnerRadiusPlusZ()

void G4Cons::SetInnerRadiusPlusZ ( G4double  Rmin2)
inline

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

◆ SetOuterRadiusMinusZ()

void G4Cons::SetOuterRadiusMinusZ ( G4double  Rmax1)
inline

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

◆ SetOuterRadiusPlusZ()

void G4Cons::SetOuterRadiusPlusZ ( G4double  Rmax2)
inline

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

◆ SetStartPhiAngle()

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

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

◆ SetZHalfLength()

void G4Cons::SetZHalfLength ( G4double  newDz)
inline

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

◆ StreamInfo()

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

Reimplemented from G4CSGSolid.

Definition at line 2112 of file G4Cons.cc.

2113{
2114 G4int oldprc = os.precision(16);
2115 os << "-----------------------------------------------------------\n"
2116 << " *** Dump for solid - " << GetName() << " ***\n"
2117 << " ===================================================\n"
2118 << " Solid type: G4Cons\n"
2119 << " Parameters: \n"
2120 << " inside -fDz radius: " << fRmin1/mm << " mm \n"
2121 << " outside -fDz radius: " << fRmax1/mm << " mm \n"
2122 << " inside +fDz radius: " << fRmin2/mm << " mm \n"
2123 << " outside +fDz radius: " << fRmax2/mm << " mm \n"
2124 << " half length in Z : " << fDz/mm << " mm \n"
2125 << " starting angle of segment: " << fSPhi/degree << " degrees \n"
2126 << " delta angle of segment : " << fDPhi/degree << " degrees \n"
2127 << "-----------------------------------------------------------\n";
2128 os.precision(oldprc);
2129
2130 return os;
2131}

◆ SurfaceNormal()

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

Implements G4VSolid.

Definition at line 419 of file G4Cons.cc.

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