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 default: // Should never reach this case...
645 DumpInfo();
646 G4Exception("G4ShiftedCone::ApproxSurfaceNormal()",
647 "GeomSolids1002", JustWarning,
648 "Undefined side for valid surface normal to solid.");
649 break ;
650 }
651 return norm ;
652}
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:17
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

◆ Clone()

G4VSolid * G4ShiftedCone::Clone ( ) const

Definition at line 2130 of file G4ShiftedCone.cxx.

2131{
2132 return new G4ShiftedCone(*this);
2133}
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 2259 of file G4ShiftedCone.cxx.

2260{
2261 G4double rmin[2] = { GetRmin1(), GetRmin2() };
2262 G4double rmax[2] = { GetRmax1(), GetRmax2() };
2263 G4double z[2] = { GetZ1(), GetZ2() };
2264 return new G4PolyhedronPcon(
2265 GetStartPhiAngle(), GetDeltaPhiAngle(), 2, z, rmin, rmax
2266 );
2267}
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 2254 of file G4ShiftedCone.cxx.

2255{
2256 scene.AddSolid (*this);
2257}

◆ DistanceToIn() [1/2]

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

Definition at line 1351 of file G4ShiftedCone.cxx.

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

◆ DistanceToIn() [2/2]

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

Definition at line 678 of file G4ShiftedCone.cxx.

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

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

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

◆ 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 2121 of file G4ShiftedCone.cxx.

2122{
2123 return G4String("G4ShiftedCone");
2124}

◆ 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 2167 of file G4ShiftedCone.cxx.

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

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

◆ 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.14.0