ATLAS Offline Software
Loading...
Searching...
No Matches
Public Member Functions | Private Types | Private Member Functions | Private Attributes | List of all members
G4ShiftedCone Class Reference

#include <G4ShiftedCone.h>

Inheritance diagram for G4ShiftedCone:
Collaboration diagram for G4ShiftedCone:

Public Member Functions

 G4ShiftedCone (const G4String &pName, G4double pZ1, G4double pZ2, G4double pRmin1, G4double pRmax1, G4double pRmin2, G4double pRmax2)
 ~G4ShiftedCone ()
G4double GetInnerRadiusMinusZ () const
G4double GetOuterRadiusMinusZ () const
G4double GetInnerRadiusPlusZ () const
G4double GetOuterRadiusPlusZ () const
G4double GetZHalfLength () const
G4double GetZ1 () const
G4double GetZ2 () 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)
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=G4bool(false), G4bool *validNorm=0, G4ThreeVector *n=0) const
G4double DistanceToOut (const G4ThreeVector &p) const
G4GeometryType GetEntityType () const
G4ThreeVector GetPointOnSurface () const
G4VSolid * Clone () const
std::ostream & StreamInfo (std::ostream &os) const
void DescribeYourselfTo (G4VGraphicsScene &scene) const
G4Polyhedron * CreatePolyhedron () const
 G4ShiftedCone (__void__ &)
 G4ShiftedCone (const G4ShiftedCone &rhs)
G4ShiftedCone & operator= (const G4ShiftedCone &rhs)
G4double GetRmin1 () const
G4double GetRmax1 () const
G4double GetRmin2 () const
G4double GetRmax2 () const

Private Types

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

Private Member Functions

void Initialize ()
void InitializeTrigonometry ()
G4ThreeVector ApproxSurfaceNormal (const G4ThreeVector &p) const

Private Attributes

G4double kRadTolerance
G4double kAngTolerance
G4double fRmin1
G4double fRmin2
G4double fRmax1
G4double fRmax2
G4double fDz
G4double fZshift
G4double halfCarTolerance
G4double halfRadTolerance
G4double halfAngTolerance

Detailed Description

Definition at line 86 of file G4ShiftedCone.h.

Member Enumeration Documentation

◆ ENorm

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

Definition at line 220 of file G4ShiftedCone.h.

220{kNRMin,kNRMax,kNSPhi,kNEPhi,kNZ};

◆ ESide

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

Definition at line 216 of file G4ShiftedCone.h.

216{kNull,kRMin,kRMax,kSPhi,kEPhi,kPZ,kMZ};

Constructor & Destructor Documentation

◆ G4ShiftedCone() [1/3]

G4ShiftedCone::G4ShiftedCone ( const G4String & pName,
G4double pZ1,
G4double pZ2,
G4double pRmin1,
G4double pRmax1,
G4double pRmin2,
G4double pRmax2 )

Definition at line 82 of file G4ShiftedCone.cxx.

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

◆ ~G4ShiftedCone()

G4ShiftedCone::~G4ShiftedCone ( )

Definition at line 149 of file G4ShiftedCone.cxx.

150{
151}

◆ G4ShiftedCone() [2/3]

G4ShiftedCone::G4ShiftedCone ( __void__ & a)

Definition at line 134 of file G4ShiftedCone.cxx.

135 : G4CSGSolid(a), kRadTolerance(0.), kAngTolerance(0.),
136 fRmin1(0.), fRmin2(0.), fRmax1(0.), fRmax2(0.),
137 fDz(0.), fZshift(0.),
138// fSPhi(0.), fDPhi(0.), sinCPhi(0.), cosCPhi(0.), cosHDPhiOT(0.),
139// cosHDPhiIT(0.), sinSPhi(0.), cosSPhi(0.), sinEPhi(0.), cosEPhi(0.),
140// fPhiFullCone(false),
141 halfCarTolerance(0.), halfRadTolerance(0.), halfAngTolerance(0.)
142{
143}
a
static Double_t a
Definition LArPhysWaveHECTool.cxx:38

◆ G4ShiftedCone() [3/3]

G4ShiftedCone::G4ShiftedCone ( const G4ShiftedCone & rhs)

Definition at line 157 of file G4ShiftedCone.cxx.

158 : G4CSGSolid(rhs), kRadTolerance(rhs.kRadTolerance),
159 kAngTolerance(rhs.kAngTolerance), fRmin1(rhs.fRmin1), fRmin2(rhs.fRmin2),
160 fRmax1(rhs.fRmax1), fRmax2(rhs.fRmax2), fDz(rhs.fDz), fZshift(rhs.fZshift),
161// fSPhi(rhs.fSPhi),
162// fDPhi(rhs.fDPhi), sinCPhi(rhs.sinCPhi), cosCPhi(rhs.cosCPhi),
163// cosHDPhiOT(rhs.cosHDPhiOT), cosHDPhiIT(rhs.cosHDPhiIT),
164// sinSPhi(rhs.sinSPhi), cosSPhi(rhs.cosSPhi), sinEPhi(rhs.sinEPhi),
165// cosEPhi(rhs.cosEPhi), fPhiFullCone(rhs.fPhiFullCone),
166 halfCarTolerance(rhs.halfCarTolerance),
167 halfRadTolerance(rhs.halfRadTolerance),
168 halfAngTolerance(rhs.halfAngTolerance)
169{
170}

Member Function Documentation

◆ ApproxSurfaceNormal()

G4ThreeVector G4ShiftedCone::ApproxSurfaceNormal ( const G4ThreeVector & p) const
private

Definition at line 549 of file G4ShiftedCone.cxx.

550{
551 ENorm side ;
552 G4ThreeVector norm ;
553 G4double rho;//, phi ;
554 G4double distZ, distRMin, distRMax;//, distSPhi, distEPhi, distMin ;
555 G4double tanRMin, secRMin, pRMin, widRMin ;
556 G4double tanRMax, secRMax, pRMax, widRMax ;
557
558 G4double z = p.z() - fZshift;
559
560 distZ = std::fabs(std::fabs(z) - fDz) ;
561 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
562
563 tanRMin = (fRmin2 - fRmin1)*0.5/fDz ;
564 secRMin = std::sqrt(1 + tanRMin*tanRMin) ;
565 pRMin = rho - z*tanRMin ;
566 widRMin = fRmin2 - fDz*tanRMin ;
567 distRMin = std::fabs(pRMin - widRMin)/secRMin ;
568
569 tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;
570 secRMax = std::sqrt(1+tanRMax*tanRMax) ;
571 pRMax = rho - z*tanRMax ;
572 widRMax = fRmax2 - fDz*tanRMax ;
573 distRMax = std::fabs(pRMax - widRMax)/secRMax ;
574
575 if (distRMin < distRMax) // First minimum
576 {
577 if (distZ < distRMin)
578 {
579// distMin = distZ ;
580 side = kNZ ;
581 }
582 else
583 {
584// distMin = distRMin ;
585 side = kNRMin ;
586 }
587 }
588 else
589 {
590 if (distZ < distRMax)
591 {
592// distMin = distZ ;
593 side = kNZ ;
594 }
595 else
596 {
597// distMin = distRMax ;
598 side = kNRMax ;
599 }
600 }
601/*
602 if ( !fPhiFullCone && rho ) // Protected against (0,0,z)
603 {
604 phi = std::atan2(p.y(),p.x()) ;
605
606 if (phi < 0) { phi += twopi; }
607
608 if (fSPhi < 0) { distSPhi = std::fabs(phi - (fSPhi + twopi))*rho; }
609 else { distSPhi = std::fabs(phi - fSPhi)*rho; }
610
611 distEPhi = std::fabs(phi - fSPhi - fDPhi)*rho ;
612
613 // Find new minimum
614
615 if (distSPhi < distEPhi)
616 {
617 if (distSPhi < distMin) { side = kNSPhi; }
618 }
619 else
620 {
621 if (distEPhi < distMin) { side = kNEPhi; }
622 }
623 }*/
624 switch (side)
625 {
626 case kNRMin: // Inner radius
627 rho *= secRMin ;
628 norm = G4ThreeVector(-p.x()/rho, -p.y()/rho, tanRMin/secRMin) ;
629 break ;
630 case kNRMax: // Outer radius
631 rho *= secRMax ;
632 norm = G4ThreeVector(p.x()/rho, p.y()/rho, -tanRMax/secRMax) ;
633 break ;
634 case kNZ: // +/- dz
635 if (z > 0) { norm = G4ThreeVector(0,0,1); }
636 else { norm = G4ThreeVector(0,0,-1); }
637 break ;
638// case kNSPhi:
639// norm = G4ThreeVector(std::sin(fSPhi), -std::cos(fSPhi), 0) ;
640// break ;
641// case kNEPhi:
642// norm=G4ThreeVector(-std::sin(fSPhi+fDPhi), std::cos(fSPhi+fDPhi), 0) ;
643// break ;
644 //coverity[DEADCODE]
645 default: // Should never reach this case...
646 DumpInfo();
647 G4Exception("G4ShiftedCone::ApproxSurfaceNormal()",
648 "GeomSolids1002", JustWarning,
649 "Undefined side for valid surface normal to solid.");
650 break ;
651 }
652 return norm ;
653}
z
#define z

◆ BoundingLimits()

void G4ShiftedCone::BoundingLimits ( G4ThreeVector & pMin,
G4ThreeVector & pMax ) const

Definition at line 283 of file G4ShiftedCone.cxx.

284{
285// G4double rmin = std::min(GetInnerRadiusMinusZ(),GetInnerRadiusPlusZ());
286 G4double rmax = std::max(GetOuterRadiusMinusZ(),GetOuterRadiusPlusZ());
287
288 // Find bounding box
289 //
290/* if (GetDeltaPhiAngle() < twopi)
291 {
292 G4double dz = GetZHalfLength();
293 G4TwoVector vmin,vmax;
294 G4GeomTools::DiskExtent(rmin,rmax,
295 GetSinStartPhi(),GetCosStartPhi(),
296 GetSinEndPhi(),GetCosEndPhi(),
297 vmin,vmax);
298 pMin.set(vmin.x(),vmin.y(),-dz);
299 pMax.set(vmax.x(),vmax.y(), dz);
300 }
301 else*/
302 {
303 pMin.set(-rmax,-rmax, fZshift - fDz);
304 pMax.set( rmax, rmax, fZshift + fDz);
305 }
306
307 // Check correctness of the bounding box
308 //
309 if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
310 {
311 std::ostringstream message;
312 message << "Bad bounding box (min >= max) for solid: "
313 << GetName() << " !"
314 << "\npMin = " << pMin
315 << "\npMax = " << pMax;
316 G4Exception("G4ShiftedCone::BoundingLimits()", "GeomMgt0001",
317 JustWarning, message);
318 DumpInfo();
319 }
320}
G4ShiftedCone::GetOuterRadiusMinusZ
G4double GetOuterRadiusMinusZ() const
G4ShiftedCone::GetOuterRadiusPlusZ
G4double GetOuterRadiusPlusZ() const

◆ CalculateExtent()

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

Definition at line 326 of file G4ShiftedCone.cxx.

331{
332 G4ThreeVector bmin, bmax;
333 G4bool exist;
334
335 // Get bounding box
336 BoundingLimits(bmin,bmax);
337
338 // Check bounding box
339 G4BoundingEnvelope bbox(bmin,bmax);
340#ifdef G4BBOX_EXTENT
341 if (true) return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
342#endif
343 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVoxelLimit,pTransform,pMin,pMax))
344 {
345 return exist = (pMin < pMax) ? true : false;
346 }
347
348 // Get parameters of the solid
349 G4double rmin1 = GetInnerRadiusMinusZ();
350 G4double rmax1 = GetOuterRadiusMinusZ();
351 G4double rmin2 = GetInnerRadiusPlusZ();
352 G4double rmax2 = GetOuterRadiusPlusZ();
353 G4double z1 = GetZ1();
354 G4double z2 = GetZ2();
355 G4double dphi = GetDeltaPhiAngle();
356
357 // Find bounding envelope and calculate extent
358 //
359 const G4int NSTEPS = 24; // number of steps for whole circle
360 G4double astep = twopi/NSTEPS; // max angle for one step
361 G4int ksteps = (dphi <= astep) ? 1 : (G4int)((dphi-deg)/astep) + 1;
362 G4double ang = dphi/ksteps;
363
364 G4double sinHalf = std::sin(0.5*ang);
365 G4double cosHalf = std::cos(0.5*ang);
366 G4double sinStep = 2.*sinHalf*cosHalf;
367 G4double cosStep = 1. - 2.*sinHalf*sinHalf;
368 G4double rext1 = rmax1/cosHalf;
369 G4double rext2 = rmax2/cosHalf;
370
371 // bounding envelope for full cone without hole consists of two polygons,
372 // in other cases it is a sequence of quadrilaterals
373 if (rmin1 == 0 && rmin2 == 0 && dphi == twopi)
374 {
375 G4double sinCur = sinHalf;
376 G4double cosCur = cosHalf;
377
378 G4ThreeVectorList baseA(NSTEPS),baseB(NSTEPS);
379 for (G4int k=0; k<NSTEPS; ++k)
380 {
381 baseA[k].set(rext1*cosCur,rext1*sinCur, z1);
382 baseB[k].set(rext2*cosCur,rext2*sinCur, z2);
383
384 G4double sinTmp = sinCur;
385 sinCur = sinCur*cosStep + cosCur*sinStep;
386 cosCur = cosCur*cosStep - sinTmp*sinStep;
387 }
388 std::vector<const G4ThreeVectorList *> polygons(2);
389 polygons[0] = &baseA;
390 polygons[1] = &baseB;
391 G4BoundingEnvelope benv(bmin,bmax,polygons);
392 exist = benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
393 }
394 else
395 {
396
397 G4double sinStart = GetSinStartPhi();
398 G4double cosStart = GetCosStartPhi();
399 G4double sinEnd = GetSinEndPhi();
400 G4double cosEnd = GetCosEndPhi();
401 G4double sinCur = sinStart*cosHalf + cosStart*sinHalf;
402 G4double cosCur = cosStart*cosHalf - sinStart*sinHalf;
403
404 // set quadrilaterals
405 G4ThreeVectorList pols[NSTEPS+2];
406 for (G4int k=0; k<ksteps+2; ++k) pols[k].resize(4);
407 pols[0][0].set(rmin2*cosStart,rmin2*sinStart, z2);
408 pols[0][1].set(rmin1*cosStart,rmin1*sinStart, z1);
409 pols[0][2].set(rmax1*cosStart,rmax1*sinStart, z1);
410 pols[0][3].set(rmax2*cosStart,rmax2*sinStart, z2);
411 for (G4int k=1; k<ksteps+1; ++k)
412 {
413 pols[k][0].set(rmin2*cosCur,rmin2*sinCur, z2);
414 pols[k][1].set(rmin1*cosCur,rmin1*sinCur, z1);
415 pols[k][2].set(rext1*cosCur,rext1*sinCur, z1);
416 pols[k][3].set(rext2*cosCur,rext2*sinCur, z2);
417
418 G4double sinTmp = sinCur;
419 sinCur = sinCur*cosStep + cosCur*sinStep;
420 cosCur = cosCur*cosStep - sinTmp*sinStep;
421 }
422 pols[ksteps+1][0].set(rmin2*cosEnd,rmin2*sinEnd, z2);
423 pols[ksteps+1][1].set(rmin1*cosEnd,rmin1*sinEnd, z1);
424 pols[ksteps+1][2].set(rmax1*cosEnd,rmax1*sinEnd, z1);
425 pols[ksteps+1][3].set(rmax2*cosEnd,rmax2*sinEnd, z2);
426
427 // set envelope and calculate extent
428 std::vector<const G4ThreeVectorList *> polygons;
429 polygons.resize(ksteps+2);
430 for (G4int k=0; k<ksteps+2; ++k) polygons[k] = &pols[k];
431 G4BoundingEnvelope benv(bmin,bmax,polygons);
432 exist = benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
433 }
434 return exist;
435}
deg
#define deg
Definition SbPolyhedron.cxx:18
twopi
constexpr double twopi
Definition VertexPointEstimator.cxx:16
G4ShiftedCone::GetDeltaPhiAngle
G4double GetDeltaPhiAngle() const
G4ShiftedCone::GetInnerRadiusMinusZ
G4double GetInnerRadiusMinusZ() const
G4ShiftedCone::GetSinEndPhi
G4double GetSinEndPhi() const
G4ShiftedCone::GetInnerRadiusPlusZ
G4double GetInnerRadiusPlusZ() const
G4ShiftedCone::GetZ1
G4double GetZ1() const
G4ShiftedCone::GetZ2
G4double GetZ2() const
G4ShiftedCone::BoundingLimits
void BoundingLimits(G4ThreeVector &pMin, G4ThreeVector &pMax) const
Definition G4ShiftedCone.cxx:283
G4ShiftedCone::GetCosEndPhi
G4double GetCosEndPhi() const
G4ShiftedCone::GetCosStartPhi
G4double GetCosStartPhi() const
G4ShiftedCone::GetSinStartPhi
G4double GetSinStartPhi() const
SG::resize
virtual bool resize(size_t sz) override
Change the size of all aux data vectors.

◆ Clone()

G4VSolid * G4ShiftedCone::Clone ( ) const

Definition at line 2127 of file G4ShiftedCone.cxx.

2128{
2129 return new G4ShiftedCone(*this);
2130}
G4ShiftedCone::G4ShiftedCone
G4ShiftedCone(const G4String &pName, G4double pZ1, G4double pZ2, G4double pRmin1, G4double pRmax1, G4double pRmin2, G4double pRmax2)
Definition G4ShiftedCone.cxx:82

◆ ComputeDimensions()

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

Definition at line 268 of file G4ShiftedCone.cxx.

271{
272 std::ostringstream message;
273 message << "ComputeDimensions is not implemented for Solid: " << GetName();
274 G4Exception("G4ShiftedCone::ComputeDimensions()", "GeomSolids0002",
275 FatalException, message) ;
276// p->ComputeDimensions(*this,n,pRep) ;
277}

◆ CreatePolyhedron()

G4Polyhedron * G4ShiftedCone::CreatePolyhedron ( ) const

Definition at line 2256 of file G4ShiftedCone.cxx.

2257{
2258 G4double rmin[2] = { GetRmin1(), GetRmin2() };
2259 G4double rmax[2] = { GetRmax1(), GetRmax2() };
2260 G4double z[2] = { GetZ1(), GetZ2() };
2261 return new G4PolyhedronPcon(
2262 GetStartPhiAngle(), GetDeltaPhiAngle(), 2, z, rmin, rmax
2263 );
2264}
G4ShiftedCone::GetRmin2
G4double GetRmin2() const
G4ShiftedCone::GetRmax2
G4double GetRmax2() const
G4ShiftedCone::GetRmax1
G4double GetRmax1() const
G4ShiftedCone::GetStartPhiAngle
G4double GetStartPhiAngle() const
G4ShiftedCone::GetRmin1
G4double GetRmin1() const

◆ DescribeYourselfTo()

void G4ShiftedCone::DescribeYourselfTo ( G4VGraphicsScene & scene) const

Definition at line 2251 of file G4ShiftedCone.cxx.

2252{
2253 scene.AddSolid (*this);
2254}

◆ DistanceToIn() [1/2]

G4double G4ShiftedCone::DistanceToIn ( const G4ThreeVector & p) const

Definition at line 1348 of file G4ShiftedCone.cxx.

1349{
1350 G4double safe=0.0, rho, safeR1, safeR2, safeZ;//, safePhi, cosPsi ;
1351 G4double tanRMin, secRMin, pRMin ;
1352 G4double tanRMax, secRMax, pRMax ;
1353
1354 G4double z = p.z() - fZshift;
1355
1356 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
1357 safeZ = std::fabs(z) - fDz ;
1358
1359 if ( fRmin1 || fRmin2 )
1360 {
1361 tanRMin = (fRmin2 - fRmin1)*0.5/fDz ;
1362 secRMin = std::sqrt(1.0 + tanRMin*tanRMin) ;
1363 pRMin = tanRMin*z + (fRmin1 + fRmin2)*0.5 ;
1364 safeR1 = (pRMin - rho)/secRMin ;
1365
1366 tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;
1367 secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
1368 pRMax = tanRMax*z + (fRmax1 + fRmax2)*0.5 ;
1369 safeR2 = (rho - pRMax)/secRMax ;
1370
1371 if ( safeR1 > safeR2) { safe = safeR1; }
1372 else { safe = safeR2; }
1373 }
1374 else
1375 {
1376 tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;
1377 secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
1378 pRMax = tanRMax*z + (fRmax1 + fRmax2)*0.5 ;
1379 safe = (rho - pRMax)/secRMax ;
1380 }
1381 if ( safeZ > safe ) { safe = safeZ; }
1382
1383/* if ( !fPhiFullCone && rho )
1384 {
1385 // Psi=angle from central phi to point
1386
1387 cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/rho ;
1388
1389 if ( cosPsi < std::cos(fDPhi*0.5) ) // Point lies outside phi range
1390 {
1391 if ( (p.y()*cosCPhi - p.x()*sinCPhi) <= 0.0 )
1392 {
1393 safePhi = std::fabs(p.x()*std::sin(fSPhi)-p.y()*std::cos(fSPhi));
1394 }
1395 else
1396 {
1397 safePhi = std::fabs(p.x()*sinEPhi-p.y()*cosEPhi);
1398 }
1399 if ( safePhi > safe ) { safe = safePhi; }
1400 }
1401 }*/
1402 if ( safe < 0.0 ) { safe = 0.0; }
1403
1404 return safe ;
1405}

◆ DistanceToIn() [2/2]

G4double G4ShiftedCone::DistanceToIn ( const G4ThreeVector & p,
const G4ThreeVector & v ) const

Definition at line 679 of file G4ShiftedCone.cxx.

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

◆ DistanceToOut() [1/2]

G4double G4ShiftedCone::DistanceToOut ( const G4ThreeVector & p) const

Definition at line 2038 of file G4ShiftedCone.cxx.

2039{
2040 G4double safe=0.0, rho, safeR1, safeR2, safeZ;//, safePhi;
2041 G4double tanRMin, secRMin, pRMin;
2042 G4double tanRMax, secRMax, pRMax;
2043
2044#ifdef G4CSGDEBUG
2045 if( Inside(p) == kOutside )
2046 {
2047 G4int oldprc=G4cout.precision(16) ;
2048 G4cout << G4endl ;
2049 DumpInfo();
2050 G4cout << "Position:" << G4endl << G4endl ;
2051 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl ;
2052 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ;
2053 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ;
2054 G4cout << "pho at z = " << std::sqrt( p.x()*p.x()+p.y()*p.y() )/mm
2055 << " mm" << G4endl << G4endl ;
2056 if( (p.x() != 0.) || (p.x() != 0.) )
2057 {
2058 G4cout << "point phi = " << std::atan2(p.y(),p.x())/degree
2059 << " degree" << G4endl << G4endl ;
2060 }
2061 G4cout.precision(oldprc) ;
2062 G4Exception("G4ShiftedCone::DistanceToOut(p)", "GeomSolids1002",
2063 JustWarning, "Point p is outside !?" );
2064 }
2065#endif
2066
2067 G4double z = p.z() - fZshift;
2068
2069 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
2070 safeZ = fDz - std::fabs(z) ;
2071
2072 if (fRmin1 || fRmin2)
2073 {
2074 tanRMin = (fRmin2 - fRmin1)*0.5/fDz ;
2075 secRMin = std::sqrt(1.0 + tanRMin*tanRMin) ;
2076 pRMin = tanRMin*z + (fRmin1 + fRmin2)*0.5 ;
2077 safeR1 = (rho - pRMin)/secRMin ;
2078 }
2079 else
2080 {
2081 safeR1 = kInfinity ;
2082 }
2083
2084 tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;
2085 secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
2086 pRMax = tanRMax*z + (fRmax1+fRmax2)*0.5 ;
2087 safeR2 = (pRMax - rho)/secRMax ;
2088
2089 if (safeR1 < safeR2) { safe = safeR1; }
2090 else { safe = safeR2; }
2091 if (safeZ < safe) { safe = safeZ ; }
2092
2093 // Check if phi divided, Calc distances closest phi plane
2094/*
2095 if (!fPhiFullCone)
2096 {
2097 // Above/below central phi of G4ShiftedCone?
2098
2099 if ( (p.y()*cosCPhi - p.x()*sinCPhi) <= 0 )
2100 {
2101 safePhi = -(p.x()*sinSPhi - p.y()*cosSPhi) ;
2102 }
2103 else
2104 {
2105 safePhi = (p.x()*sinEPhi - p.y()*cosEPhi) ;
2106 }
2107 if (safePhi < safe) { safe = safePhi; }
2108 }*/
2109 if ( safe < 0 ) { safe = 0; }
2110
2111 return safe ;
2112}
G4ShiftedCone::Inside
EInside Inside(const G4ThreeVector &p) const
Definition G4ShiftedCone.cxx:210

◆ DistanceToOut() [2/2]

G4double G4ShiftedCone::DistanceToOut ( const G4ThreeVector & p,
const G4ThreeVector & v,
const G4bool calcNorm = G4bool(false),
G4bool * validNorm = 0,
G4ThreeVector * n = 0 ) const

Definition at line 1412 of file G4ShiftedCone.cxx.

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

◆ GetCosEndPhi()

G4double G4ShiftedCone::GetCosEndPhi ( ) const
inline

◆ GetCosStartPhi()

G4double G4ShiftedCone::GetCosStartPhi ( ) const
inline

◆ GetCubicVolume()

G4double G4ShiftedCone::GetCubicVolume ( )
inline

◆ GetDeltaPhiAngle()

G4double G4ShiftedCone::GetDeltaPhiAngle ( ) const
inline

◆ GetEntityType()

G4GeometryType G4ShiftedCone::GetEntityType ( ) const

Definition at line 2118 of file G4ShiftedCone.cxx.

2119{
2120 return G4String("G4ShiftedCone");
2121}

◆ GetInnerRadiusMinusZ()

G4double G4ShiftedCone::GetInnerRadiusMinusZ ( ) const
inline

◆ GetInnerRadiusPlusZ()

G4double G4ShiftedCone::GetInnerRadiusPlusZ ( ) const
inline

◆ GetOuterRadiusMinusZ()

G4double G4ShiftedCone::GetOuterRadiusMinusZ ( ) const
inline

◆ GetOuterRadiusPlusZ()

G4double G4ShiftedCone::GetOuterRadiusPlusZ ( ) const
inline

◆ GetPointOnSurface()

G4ThreeVector G4ShiftedCone::GetPointOnSurface ( ) const

Definition at line 2164 of file G4ShiftedCone.cxx.

2165{
2166 // declare working variables
2167 //
2168 G4double rone = (fRmax1-fRmax2)/(2.*fDz);
2169 G4double rtwo = (fRmin1-fRmin2)/(2.*fDz);
2170 G4double qone = (fRmax1 == fRmax2) ? 0. : fDz*(fRmax1+fRmax2)/(fRmax1-fRmax2);
2171 G4double qtwo = (fRmin1 == fRmin2) ? 0. : fDz*(fRmin1+fRmin2)/(fRmin1-fRmin2);
2172
2173 G4double slin = std::hypot(fRmin1-fRmin2, 2.*fDz);
2174 G4double slout = std::hypot(fRmax1-fRmax2, 2.*fDz);
2175 G4double Aone = 0.5*GetDeltaPhiAngle()*(fRmax2 + fRmax1)*slout; // outer surface
2176 G4double Atwo = 0.5*GetDeltaPhiAngle()*(fRmin2 + fRmin1)*slin; // inner surface
2177 G4double Athree = 0.5*GetDeltaPhiAngle()*(fRmax1*fRmax1-fRmin1*fRmin1); // base at -Dz
2178 G4double Afour = 0.5*GetDeltaPhiAngle()*(fRmax2*fRmax2-fRmin2*fRmin2); // base at +Dz
2179 G4double Afive = fDz*(fRmax1-fRmin1+fRmax2-fRmin2); // phi section
2180
2181 G4double phi = G4RandFlat::shoot(GetStartPhiAngle(),GetStartPhiAngle() + GetDeltaPhiAngle());
2182 G4double cosu = std::cos(phi);
2183 G4double sinu = std::sin(phi);
2184 G4double rRand1 = GetRadiusInRing(fRmin1, fRmax1);
2185 G4double rRand2 = GetRadiusInRing(fRmin2, fRmax2);
2186
2187 G4bool fPhiFullCone = true;
2188 if ( (GetStartPhiAngle() == 0.) && fPhiFullCone ) { Afive = 0.; }
2189 G4double chose = G4RandFlat::shoot(0.,Aone+Atwo+Athree+Afour+2.*Afive);
2190
2191 if( (chose >= 0.) && (chose < Aone) ) // outer surface
2192 {
2193 if(fRmax1 != fRmax2)
2194 {
2195 G4double zRand = G4RandFlat::shoot(-1.*fDz,fDz);
2196 return G4ThreeVector (rone*cosu*(qone-zRand),
2197 rone*sinu*(qone-zRand), zRand + fZshift);
2198 }
2199 else
2200 {
2201 return G4ThreeVector(fRmax1*cosu, fRmax2*sinu,
2202 G4RandFlat::shoot(-1.*fDz,fDz) + fZshift);
2203 }
2204 }
2205 else if( (chose >= Aone) && (chose < Aone + Atwo) ) // inner surface
2206 {
2207 if(fRmin1 != fRmin2)
2208 {
2209 G4double zRand = G4RandFlat::shoot(-1.*fDz,fDz);
2210 return G4ThreeVector (rtwo*cosu*(qtwo-zRand),
2211 rtwo*sinu*(qtwo-zRand), zRand + fZshift);
2212 }
2213 else
2214 {
2215 return G4ThreeVector(fRmin1*cosu, fRmin2*sinu,
2216 G4RandFlat::shoot(-1.*fDz,fDz) + fZshift);
2217 }
2218 }
2219 else if( (chose >= Aone + Atwo) && (chose < Aone + Atwo + Athree) ) // base at -Dz
2220 {
2221 return G4ThreeVector (rRand1*cosu, rRand1*sinu, -1*fDz + fZshift);
2222 }
2223 else if( (chose >= Aone + Atwo + Athree)
2224 && (chose < Aone + Atwo + Athree + Afour) ) // base at +Dz
2225 {
2226 return G4ThreeVector (rRand2*cosu,rRand2*sinu,fDz + fZshift);
2227 }
2228 else if( (chose >= Aone + Atwo + Athree + Afour) // SPhi section
2229 && (chose < Aone + Atwo + Athree + Afour + Afive) )
2230 {
2231 G4double zRand = G4RandFlat::shoot(-1.*fDz,fDz);
2232 rRand1 = G4RandFlat::shoot(fRmin2-((zRand-fDz)/(2.*fDz))*(fRmin1-fRmin2),
2233 fRmax2-((zRand-fDz)/(2.*fDz))*(fRmax1-fRmax2));
2234 return G4ThreeVector (rRand1*GetCosStartPhi(),
2235 rRand1*GetSinStartPhi(), zRand + fZshift);
2236 }
2237 else // SPhi+DPhi section
2238 {
2239 G4double zRand = G4RandFlat::shoot(-1.*fDz,fDz);
2240 rRand1 = G4RandFlat::shoot(fRmin2-((zRand-fDz)/(2.*fDz))*(fRmin1-fRmin2),
2241 fRmax2-((zRand-fDz)/(2.*fDz))*(fRmax1-fRmax2));
2242 return G4ThreeVector (rRand1*GetCosEndPhi(),
2243 rRand1*GetSinEndPhi(), zRand + fZshift);
2244 }
2245}
phi
Scalar phi() const
phi method
Definition AmgMatrixBasePlugin.h:67

◆ GetRmax1()

G4double G4ShiftedCone::GetRmax1 ( ) const
inline

◆ GetRmax2()

G4double G4ShiftedCone::GetRmax2 ( ) const
inline

◆ GetRmin1()

G4double G4ShiftedCone::GetRmin1 ( ) const
inline

◆ GetRmin2()

G4double G4ShiftedCone::GetRmin2 ( ) const
inline

◆ GetSinEndPhi()

G4double G4ShiftedCone::GetSinEndPhi ( ) const
inline

◆ GetSinStartPhi()

G4double G4ShiftedCone::GetSinStartPhi ( ) const
inline

◆ GetStartPhiAngle()

G4double G4ShiftedCone::GetStartPhiAngle ( ) const
inline

◆ GetSurfaceArea()

G4double G4ShiftedCone::GetSurfaceArea ( )
inline

◆ GetZ1()

G4double G4ShiftedCone::GetZ1 ( ) const
inline

◆ GetZ2()

G4double G4ShiftedCone::GetZ2 ( ) const
inline

◆ GetZHalfLength()

G4double G4ShiftedCone::GetZHalfLength ( ) const
inline

◆ Initialize()

void G4ShiftedCone::Initialize ( )
inlineprivate

◆ InitializeTrigonometry()

void G4ShiftedCone::InitializeTrigonometry ( )
inlineprivate

◆ Inside()

EInside G4ShiftedCone::Inside ( const G4ThreeVector & p) const

Definition at line 210 of file G4ShiftedCone.cxx.

211{
212 G4double r2, rl, rh, /*pPhi,*/ tolRMin, tolRMax; // rh2, rl2 ;
213 EInside in;
214
215 G4double z = p.z() - fZshift;
216
217 if (std::fabs(z) > fDz + halfCarTolerance ) { return in = kOutside; }
218 else if(std::fabs(z) >= fDz - halfCarTolerance ) { in = kSurface; }
219 else { in = kInside; }
220
221 r2 = p.x()*p.x() + p.y()*p.y() ;
222 rl = 0.5*(fRmin2*(z + fDz) + fRmin1*(fDz - z))/fDz ;
223 rh = 0.5*(fRmax2*(z+fDz)+fRmax1*(fDz-z))/fDz;
224
225 // rh2 = rh*rh;
226
227 tolRMin = rl - halfRadTolerance;
228 if ( tolRMin < 0 ) { tolRMin = 0; }
229 tolRMax = rh + halfRadTolerance;
230
231 if ( (r2<tolRMin*tolRMin) || (r2>tolRMax*tolRMax) ) { return in = kOutside; }
232
233 if (rl) { tolRMin = rl + halfRadTolerance; }
234 else { tolRMin = 0.0; }
235 tolRMax = rh - halfRadTolerance;
236
237 if (in == kInside) // else it's kSurface already
238 {
239 if ( (r2 < tolRMin*tolRMin) || (r2 >= tolRMax*tolRMax) ) { in = kSurface; }
240 }
241/*
242 if ( !fPhiFullCone && ((p.x() != 0.0) || (p.y() != 0.0)) )
243 {
244 pPhi = std::atan2(p.y(),p.x()) ;
245
246 if ( pPhi < fSPhi - halfAngTolerance ) { pPhi += twopi; }
247 else if ( pPhi > fSPhi + fDPhi + halfAngTolerance ) { pPhi -= twopi; }
248
249 if ( (pPhi < fSPhi - halfAngTolerance) ||
250 (pPhi > fSPhi + fDPhi + halfAngTolerance) ) { return in = kOutside; }
251
252 else if (in == kInside) // else it's kSurface anyway already
253 {
254 if ( (pPhi < fSPhi + halfAngTolerance) ||
255 (pPhi > fSPhi + fDPhi - halfAngTolerance) ) { in = kSurface; }
256 }
257 }
258 else if ( !fPhiFullCone ) { in = kSurface; }
259*/
260 return in ;
261}

◆ operator=()

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

Definition at line 176 of file G4ShiftedCone.cxx.

177{
178 // Check assignment to self
179 //
180 if (this == &rhs) { return *this; }
181
182 // Copy base class data
183 //
184 G4CSGSolid::operator=(rhs);
185
186 // Copy data
187 //
188 kRadTolerance = rhs.kRadTolerance;
189 kAngTolerance = rhs.kAngTolerance;
190 fRmin1 = rhs.fRmin1; fRmin2 = rhs.fRmin2;
191 fRmax1 = rhs.fRmax1; fRmax2 = rhs.fRmax2;
192 fDz = rhs.fDz; fZshift = rhs.fZshift;
193// fSPhi = rhs.fSPhi; fDPhi = rhs.fDPhi;
194// sinCPhi = rhs.sinCPhi; cosCPhi = rhs.cosCPhi;
195// cosHDPhiOT = rhs.cosHDPhiOT; cosHDPhiIT = rhs.cosHDPhiIT;
196// sinSPhi = rhs.sinSPhi; cosSPhi = rhs.cosSPhi;
197// sinEPhi = rhs.sinEPhi; cosEPhi = rhs.cosEPhi;
198// fPhiFullCone = rhs.fPhiFullCone;
199 halfCarTolerance = rhs.halfCarTolerance;
200 halfRadTolerance = rhs.halfRadTolerance;
201 halfAngTolerance = rhs.halfAngTolerance;
202
203 return *this;
204}

◆ SetInnerRadiusMinusZ()

void G4ShiftedCone::SetInnerRadiusMinusZ ( G4double Rmin1)
inline

◆ SetInnerRadiusPlusZ()

void G4ShiftedCone::SetInnerRadiusPlusZ ( G4double Rmin2)
inline

◆ SetOuterRadiusMinusZ()

void G4ShiftedCone::SetOuterRadiusMinusZ ( G4double Rmax1)
inline

◆ SetOuterRadiusPlusZ()

void G4ShiftedCone::SetOuterRadiusPlusZ ( G4double Rmax2)
inline

◆ StreamInfo()

std::ostream & G4ShiftedCone::StreamInfo ( std::ostream & os) const

Definition at line 2136 of file G4ShiftedCone.cxx.

2137{
2138 G4int oldprc = os.precision(16);
2139 os << "-----------------------------------------------------------\n"
2140 << " *** Dump for solid - " << GetName() << " ***\n"
2141 << " ===================================================\n"
2142 << " Solid type: G4ShiftedCone\n"
2143 << " Parameters: \n"
2144 << " inside -fDz radius: " << fRmin1/mm << " mm \n"
2145 << " outside -fDz radius: " << fRmax1/mm << " mm \n"
2146 << " inside +fDz radius: " << fRmin2/mm << " mm \n"
2147 << " outside +fDz radius: " << fRmax2/mm << " mm \n"
2148 << " Z1 : " << GetZ1()/mm << " mm \n"
2149 << " Z2 : " << GetZ2()/mm << " mm \n"
2150// << " starting angle of segment: " << fSPhi/degree << " degrees \n"
2151// << " delta angle of segment : " << fDPhi/degree << " degrees \n"
2152 << "-----------------------------------------------------------\n";
2153 os.precision(oldprc);
2154
2155 return os;
2156}

◆ SurfaceNormal()

G4ThreeVector G4ShiftedCone::SurfaceNormal ( const G4ThreeVector & p) const

Definition at line 443 of file G4ShiftedCone.cxx.

444{
445 G4int noSurfaces = 0;
446 G4double rho;//, pPhi;
447 G4double distZ, distRMin, distRMax;
448// G4double distSPhi = kInfinity, distEPhi = kInfinity;
449 G4double tanRMin, secRMin, pRMin, widRMin;
450 G4double tanRMax, secRMax, pRMax, widRMax;
451
452 G4ThreeVector norm, sumnorm(0.,0.,0.), nZ = G4ThreeVector(0.,0.,1.);
453 G4ThreeVector nR, nr(0.,0.,0.), nPs, nPe;
454
455 G4double z = p.z() - fZshift;
456
457 distZ = std::fabs(std::fabs(z) - fDz);
458 rho = std::sqrt(p.x()*p.x() + p.y()*p.y());
459
460 tanRMin = (fRmin2 - fRmin1)*0.5/fDz;
461 secRMin = std::sqrt(1 + tanRMin*tanRMin);
462 pRMin = rho - z*tanRMin;
463 widRMin = fRmin2 - fDz*tanRMin;
464 distRMin = std::fabs(pRMin - widRMin)/secRMin;
465
466 tanRMax = (fRmax2 - fRmax1)*0.5/fDz;
467 secRMax = std::sqrt(1+tanRMax*tanRMax);
468 pRMax = rho - z*tanRMax;
469 widRMax = fRmax2 - fDz*tanRMax;
470 distRMax = std::fabs(pRMax - widRMax)/secRMax;
471/*
472 if (!fPhiFullCone) // Protected against (0,0,z)
473 {
474 if ( rho )
475 {
476 pPhi = std::atan2(p.y(),p.x());
477
478 if (pPhi < fSPhi-halfCarTolerance) { pPhi += twopi; }
479 else if (pPhi > fSPhi+fDPhi+halfCarTolerance) { pPhi -= twopi; }
480
481 distSPhi = std::fabs( pPhi - fSPhi );
482 distEPhi = std::fabs( pPhi - fSPhi - fDPhi );
483 }
484 else if( !(fRmin1) || !(fRmin2) )
485 {
486 distSPhi = 0.;
487 distEPhi = 0.;
488 }
489 nPs = G4ThreeVector(std::sin(fSPhi), -std::cos(fSPhi), 0);
490 nPe = G4ThreeVector(-std::sin(fSPhi+fDPhi), std::cos(fSPhi+fDPhi), 0);
491 }*/
492 if ( rho > halfCarTolerance )
493 {
494 nR = G4ThreeVector(p.x()/rho/secRMax, p.y()/rho/secRMax, -tanRMax/secRMax);
495 if (fRmin1 || fRmin2)
496 {
497 nr = G4ThreeVector(-p.x()/rho/secRMin,-p.y()/rho/secRMin,tanRMin/secRMin);
498 }
499 }
500
501 if( distRMax <= halfCarTolerance )
502 {
503 noSurfaces ++;
504 sumnorm += nR;
505 }
506 if( (fRmin1 || fRmin2) && (distRMin <= halfCarTolerance) )
507 {
508 noSurfaces ++;
509 sumnorm += nr;
510 }
511/* if( !fPhiFullCone )
512 {
513 if (distSPhi <= halfAngTolerance)
514 {
515 noSurfaces ++;
516 sumnorm += nPs;
517 }
518 if (distEPhi <= halfAngTolerance)
519 {
520 noSurfaces ++;
521 sumnorm += nPe;
522 }
523 }*/
524 if (distZ <= halfCarTolerance)
525 {
526 noSurfaces ++;
527 if ( z >= 0.) { sumnorm += nZ; }
528 else { sumnorm -= nZ; }
529 }
530 if ( noSurfaces == 0 )
531 {
532#ifdef G4CSGDEBUG
533 G4Exception("G4ShiftedCone::SurfaceNormal(p)", "GeomSolids1002",
534 JustWarning, "Point p is not on surface !?" );
535#endif
536 norm = ApproxSurfaceNormal(p);
537 }
538 else if ( noSurfaces == 1 ) { norm = sumnorm; }
539 else { norm = sumnorm.unit(); }
540
541 return norm ;
542}
G4ShiftedCone::ApproxSurfaceNormal
G4ThreeVector ApproxSurfaceNormal(const G4ThreeVector &p) const
Definition G4ShiftedCone.cxx:549

Member Data Documentation

◆ fDz

G4double G4ShiftedCone::fDz
private

Definition at line 227 of file G4ShiftedCone.h.

◆ fRmax1

G4double G4ShiftedCone::fRmax1
private

Definition at line 226 of file G4ShiftedCone.h.

◆ fRmax2

G4double G4ShiftedCone::fRmax2
private

Definition at line 226 of file G4ShiftedCone.h.

◆ fRmin1

G4double G4ShiftedCone::fRmin1
private

Definition at line 226 of file G4ShiftedCone.h.

◆ fRmin2

G4double G4ShiftedCone::fRmin2
private

Definition at line 226 of file G4ShiftedCone.h.

◆ fZshift

G4double G4ShiftedCone::fZshift
private

Definition at line 227 of file G4ShiftedCone.h.

◆ halfAngTolerance

G4double G4ShiftedCone::halfAngTolerance
private

Definition at line 240 of file G4ShiftedCone.h.

◆ halfCarTolerance

G4double G4ShiftedCone::halfCarTolerance
private

Definition at line 240 of file G4ShiftedCone.h.

◆ halfRadTolerance

G4double G4ShiftedCone::halfRadTolerance
private

Definition at line 240 of file G4ShiftedCone.h.

◆ kAngTolerance

G4double G4ShiftedCone::kAngTolerance
private

Definition at line 222 of file G4ShiftedCone.h.

◆ kRadTolerance

G4double G4ShiftedCone::kRadTolerance
private

Definition at line 222 of file G4ShiftedCone.h.

The documentation for this class was generated from the following files:
Generated on for ATLAS Offline Software by doxygen 1.16.1