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

#include <G4CutTubs.hh>

+ Inheritance diagram for G4CutTubs:

Public Member Functions

 G4CutTubs (const G4String &pName, G4double pRMin, G4double pRMax, G4double pDz, G4double pSPhi, G4double pDPhi, G4ThreeVector pLowNorm, G4ThreeVector pHighNorm)
 
 ~G4CutTubs ()
 
G4ThreeVector GetLowNorm () const
 
G4ThreeVector GetHighNorm () const
 
G4double GetCubicVolume ()
 
G4double GetSurfaceArea ()
 
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
 
 G4CutTubs (__void__ &)
 
 G4CutTubs (const G4CutTubs &rhs)
 
G4CutTubs & operator= (const G4CutTubs &rhs)
 
- Public Member Functions inherited from G4Tubs
 G4Tubs (const G4String &pName, G4double pRMin, G4double pRMax, G4double pDz, G4double pSPhi, G4double pDPhi)
 
virtual ~G4Tubs ()
 
G4double GetInnerRadius () const
 
G4double GetOuterRadius () const
 
G4double GetZHalfLength () const
 
G4double GetStartPhiAngle () const
 
G4double GetDeltaPhiAngle () const
 
void SetInnerRadius (G4double newRMin)
 
void SetOuterRadius (G4double newRMax)
 
void SetZHalfLength (G4double newDz)
 
void SetStartPhiAngle (G4double newSPhi, G4bool trig=true)
 
void SetDeltaPhiAngle (G4double newDPhi)
 
G4double GetCubicVolume ()
 
G4double GetSurfaceArea ()
 
void ComputeDimensions (G4VPVParameterisation *p, const G4int n, const G4VPhysicalVolume *pRep)
 
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
 
G4NURBS * CreateNURBS () const
 
 G4Tubs (__void__ &)
 
 G4Tubs (const G4Tubs &rhs)
 
G4Tubs & operator= (const G4Tubs &rhs)
 
G4double GetRMin () const
 
G4double GetRMax () const
 
G4double GetDz () const
 
G4double GetSPhi () const
 
G4double GetDPhi () const
 
- Public Member Functions inherited from G4CSGSolid
 G4CSGSolid (const G4String &pName)
 
virtual ~G4CSGSolid ()
 
virtual std::ostream & StreamInfo (std::ostream &os) const
 
virtual G4Polyhedron * GetPolyhedron () const
 
 G4CSGSolid (__void__ &)
 
 G4CSGSolid (const G4CSGSolid &rhs)
 
G4CSGSolid & operator= (const G4CSGSolid &rhs)
 
- Public Member Functions inherited from G4VSolid
 G4VSolid (const G4String &name)
 
virtual ~G4VSolid ()
 
G4bool operator== (const G4VSolid &s) const
 
G4String GetName () const
 
void SetName (const G4String &name)
 
G4double GetTolerance () const
 
virtual G4bool CalculateExtent (const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pMin, G4double &pMax) const =0
 
virtual EInside Inside (const G4ThreeVector &p) const =0
 
virtual G4ThreeVector SurfaceNormal (const G4ThreeVector &p) const =0
 
virtual G4double DistanceToIn (const G4ThreeVector &p, const G4ThreeVector &v) const =0
 
virtual G4double DistanceToIn (const G4ThreeVector &p) const =0
 
virtual G4double DistanceToOut (const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=false, G4bool *validNorm=0, G4ThreeVector *n=0) const =0
 
virtual G4double DistanceToOut (const G4ThreeVector &p) const =0
 
virtual void ComputeDimensions (G4VPVParameterisation *p, const G4int n, const G4VPhysicalVolume *pRep)
 
virtual G4double GetCubicVolume ()
 
virtual G4double GetSurfaceArea ()
 
virtual G4GeometryType GetEntityType () const =0
 
virtual G4ThreeVector GetPointOnSurface () const
 
virtual G4VSolid * Clone () const
 
virtual std::ostream & StreamInfo (std::ostream &os) const =0
 
void DumpInfo () const
 
virtual void DescribeYourselfTo (G4VGraphicsScene &scene) const =0
 
virtual G4VisExtent GetExtent () const
 
virtual G4Polyhedron * CreatePolyhedron () const
 
virtual G4NURBS * CreateNURBS () const
 
virtual G4Polyhedron * GetPolyhedron () const
 
virtual const G4VSolid * GetConstituentSolid (G4int no) const
 
virtual G4VSolid * GetConstituentSolid (G4int no)
 
virtual const G4DisplacedSolid * GetDisplacedSolidPtr () const
 
virtual G4DisplacedSolid * GetDisplacedSolidPtr ()
 
 G4VSolid (__void__ &)
 
 G4VSolid (const G4VSolid &rhs)
 
G4VSolid & operator= (const G4VSolid &rhs)
 

Protected Member Functions

G4ThreeVectorList * CreateRotatedVertices (const G4AffineTransform &pTransform) const
 
G4ThreeVector ApproxSurfaceNormal (const G4ThreeVector &p) const
 
G4bool IsCrossingCutPlanes () const
 
G4double GetCutZ (const G4ThreeVector &p) const
 
void GetMaxMinZ (G4double &zmin, G4double &zmax) const
 
- Protected Member Functions inherited from G4Tubs
G4ThreeVectorList * CreateRotatedVertices (const G4AffineTransform &pTransform) const
 
void Initialize ()
 
void CheckSPhiAngle (G4double sPhi)
 
void CheckDPhiAngle (G4double dPhi)
 
void CheckPhiAngles (G4double sPhi, G4double dPhi)
 
void InitializeTrigonometry ()
 
virtual G4ThreeVector ApproxSurfaceNormal (const G4ThreeVector &p) const
 
- Protected Member Functions inherited from G4CSGSolid
G4double GetRadiusInRing (G4double rmin, G4double rmax) const
 
- Protected Member Functions inherited from G4VSolid
void CalculateClippedPolygonExtent (G4ThreeVectorList &pPolygon, const G4VoxelLimits &pVoxelLimit, const EAxis pAxis, G4double &pMin, G4double &pMax) const
 
void ClipCrossSection (G4ThreeVectorList *pVertices, const G4int pSectionIndex, const G4VoxelLimits &pVoxelLimit, const EAxis pAxis, G4double &pMin, G4double &pMax) const
 
void ClipBetweenSections (G4ThreeVectorList *pVertices, const G4int pSectionIndex, const G4VoxelLimits &pVoxelLimit, const EAxis pAxis, G4double &pMin, G4double &pMax) const
 
void ClipPolygon (G4ThreeVectorList &pPolygon, const G4VoxelLimits &pVoxelLimit, const EAxis pAxis) const
 
G4double EstimateCubicVolume (G4int nStat, G4double epsilon) const
 
G4double EstimateSurfaceArea (G4int nStat, G4double ell) const
 

Additional Inherited Members

- Protected Types inherited from G4Tubs
enum  ESide {
  kNull , kRMin , kRMax , kSPhi ,
  kEPhi , kPZ , kMZ
}
 
enum  ENorm {
  kNRMin , kNRMax , kNSPhi , kNEPhi ,
  kNZ
}
 
- Protected Attributes inherited from G4Tubs
G4double kRadTolerance
 
G4double kAngTolerance
 
G4double fRMin
 
G4double fRMax
 
G4double fDz
 
G4double fSPhi
 
G4double fDPhi
 
G4double sinCPhi
 
G4double cosCPhi
 
G4double cosHDPhiOT
 
G4double cosHDPhiIT
 
G4double sinSPhi
 
G4double cosSPhi
 
G4double sinEPhi
 
G4double cosEPhi
 
G4bool fPhiFullTube
 
- Protected Attributes inherited from G4CSGSolid
G4double fCubicVolume
 
G4double fSurfaceArea
 
G4Polyhedron * fpPolyhedron
 
- Protected Attributes inherited from G4VSolid
G4double kCarTolerance
 

Detailed Description

Definition at line 52 of file G4CutTubs.hh.

Constructor & Destructor Documentation

◆ G4CutTubs() [1/3]

G4CutTubs::G4CutTubs ( const G4String &  pName,
G4double  pRMin,
G4double  pRMax,
G4double  pDz,
G4double  pSPhi,
G4double  pDPhi,
G4ThreeVector  pLowNorm,
G4ThreeVector  pHighNorm 
)

Definition at line 66 of file G4CutTubs.cc.

71 : G4Tubs(pName, pRMin, pRMax, pDz, pSPhi, pDPhi),
72 fPhiFullCutTube(true)
73{
74 kRadTolerance = G4GeometryTolerance::GetInstance()->GetRadialTolerance();
75 kAngTolerance = G4GeometryTolerance::GetInstance()->GetAngularTolerance();
76
77 // Check on Cutted Planes Normals
78 // If there is NO CUT, propose to use G4Tubs instead
79 //
80 if(pDPhi<twopi) { fPhiFullCutTube=false; }
81
82 if ( ( !pLowNorm.x()) && ( !pLowNorm.y())
83 && ( !pHighNorm.x()) && (!pHighNorm.y()) )
84 {
85 std::ostringstream message;
86 message << "Inexisting Low/High Normal to Z plane or Parallel to Z."
87 << G4endl
88 << "Normals to Z plane are (" << pLowNorm <<" and "
89 << pHighNorm << ") in solid: " << GetName();
90 G4Exception("G4CutTubs::G4CutTubs()", "GeomSolids1001",
91 JustWarning, message, "Should use G4Tubs!");
92 }
93
94 // If Normal is (0,0,0),means parallel to R, give it value of (0,0,+/-1)
95 //
96 if (pLowNorm.mag2() == 0.) { pLowNorm.setZ(-1.); }
97 if (pHighNorm.mag2()== 0.) { pHighNorm.setZ(1.); }
98
99 // Given Normals to Cut Planes have to be an unit vectors.
100 // Normalize if it is needed.
101 //
102 if (pLowNorm.mag2() != 1.) { pLowNorm = pLowNorm.unit(); }
103 if (pHighNorm.mag2()!= 1.) { pHighNorm = pHighNorm.unit(); }
104
105 // Normals to cutted planes have to point outside Solid
106 //
107 if( (pLowNorm.mag2() != 0.) && (pHighNorm.mag2()!= 0. ) )
108 {
109 if( ( pLowNorm.z()>= 0. ) || ( pHighNorm.z() <= 0.))
110 {
111 std::ostringstream message;
112 message << "Invalid Low or High Normal to Z plane; "
113 "has to point outside Solid." << G4endl
114 << "Invalid Norm to Z plane (" << pLowNorm << " or "
115 << pHighNorm << ") in solid: " << GetName();
116 G4Exception("G4CutTubs::G4CutTubs()", "GeomSolids0002",
117 FatalException, message);
118 }
119 }
120 fLowNorm = pLowNorm;
121 fHighNorm = pHighNorm;
122
123 // Check Intersection of Cutted planes. They MUST NOT Intersect
124 //
125 if(IsCrossingCutPlanes())
126 {
127 std::ostringstream message;
128 message << "Invalid Low or High Normal to Z plane; "
129 << "Crossing Cutted Planes." << G4endl
130 << "Invalid Norm to Z plane (" << pLowNorm << " and "
131 << pHighNorm << ") in solid: " << GetName();
132 G4Exception("G4CutTubs::G4CutTubs()", "GeomSolids0002",
133 FatalException, message);
134 }
135}
JustWarning
@ JustWarning
Definition: G4ExceptionSeverity.hh:64
FatalException
@ FatalException
Definition: G4ExceptionSeverity.hh:60
G4endl
#define G4endl
Definition: G4ios.hh:52
CLHEP::Hep3Vector::z
double z() const
CLHEP::Hep3Vector::unit
Hep3Vector unit() const
CLHEP::Hep3Vector::x
double x() const
CLHEP::Hep3Vector::mag2
double mag2() const
CLHEP::Hep3Vector::y
double y() const
CLHEP::Hep3Vector::setZ
void setZ(double)
G4CutTubs::IsCrossingCutPlanes
G4bool IsCrossingCutPlanes() const
Definition: G4CutTubs.cc:2035
G4GeometryTolerance::GetRadialTolerance
G4double GetRadialTolerance() const
Definition: G4GeometryTolerance.cc:97
G4GeometryTolerance::GetInstance
static G4GeometryTolerance * GetInstance()
Definition: G4GeometryTolerance.cc:73
G4GeometryTolerance::GetAngularTolerance
G4double GetAngularTolerance() const
Definition: G4GeometryTolerance.cc:92
G4Tubs
Definition: G4Tubs.hh:77
G4Tubs::kAngTolerance
G4double kAngTolerance
Definition: G4Tubs.hh:205
G4Tubs::kRadTolerance
G4double kRadTolerance
Definition: G4Tubs.hh:205
G4VSolid::GetName
G4String GetName() const
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41

◆ ~G4CutTubs()

G4CutTubs::~G4CutTubs ( )

Definition at line 153 of file G4CutTubs.cc.

154{
155}

◆ G4CutTubs() [2/3]

G4CutTubs::G4CutTubs ( __void__ &  a)

Definition at line 142 of file G4CutTubs.cc.

143 : G4Tubs(a),
144 fLowNorm(G4ThreeVector()), fHighNorm(G4ThreeVector()),
145 fPhiFullCutTube(false)
146{
147}
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:42

◆ G4CutTubs() [3/3]

G4CutTubs::G4CutTubs ( const G4CutTubs &  rhs)

Definition at line 161 of file G4CutTubs.cc.

162 : G4Tubs(rhs),
163 fLowNorm(rhs.fLowNorm), fHighNorm(rhs.fHighNorm),
164 fPhiFullCutTube(rhs.fPhiFullCutTube)
165{
166}

Member Function Documentation

◆ ApproxSurfaceNormal()

G4ThreeVector G4CutTubs::ApproxSurfaceNormal ( const G4ThreeVector &  p) const
protectedvirtual

Reimplemented from G4Tubs.

Definition at line 623 of file G4CutTubs.cc.

624{
625 ENorm side ;
626 G4ThreeVector norm ;
627 G4double rho, phi ;
628 G4double distZLow,distZHigh,distZ;
629 G4double distRMin, distRMax, distSPhi, distEPhi, distMin ;
630 G4ThreeVector vZ=G4ThreeVector(0,0,fDz);
631
632 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
633
634 distRMin = std::fabs(rho - fRMin) ;
635 distRMax = std::fabs(rho - fRMax) ;
636
637 //dist to Low Cut
638 //
639 distZLow =std::fabs((p+vZ).dot(fLowNorm));
640
641 //dist to High Cut
642 //
643 distZHigh = std::fabs((p-vZ).dot(fHighNorm));
644 distZ=std::min(distZLow,distZHigh);
645
646 if (distRMin < distRMax) // First minimum
647 {
648 if ( distZ < distRMin )
649 {
650 distMin = distZ ;
651 side = kNZ ;
652 }
653 else
654 {
655 distMin = distRMin ;
656 side = kNRMin ;
657 }
658 }
659 else
660 {
661 if ( distZ < distRMax )
662 {
663 distMin = distZ ;
664 side = kNZ ;
665 }
666 else
667 {
668 distMin = distRMax ;
669 side = kNRMax ;
670 }
671 }
672 if (!fPhiFullCutTube && rho ) // Protected against (0,0,z)
673 {
674 phi = std::atan2(p.y(),p.x()) ;
675
676 if ( phi < 0 ) { phi += twopi; }
677
678 if ( fSPhi < 0 )
679 {
680 distSPhi = std::fabs(phi - (fSPhi + twopi))*rho ;
681 }
682 else
683 {
684 distSPhi = std::fabs(phi - fSPhi)*rho ;
685 }
686 distEPhi = std::fabs(phi - fSPhi - fDPhi)*rho ;
687
688 if (distSPhi < distEPhi) // Find new minimum
689 {
690 if ( distSPhi < distMin )
691 {
692 side = kNSPhi ;
693 }
694 }
695 else
696 {
697 if ( distEPhi < distMin )
698 {
699 side = kNEPhi ;
700 }
701 }
702 }
703 switch ( side )
704 {
705 case kNRMin : // Inner radius
706 {
707 norm = G4ThreeVector(-p.x()/rho, -p.y()/rho, 0) ;
708 break ;
709 }
710 case kNRMax : // Outer radius
711 {
712 norm = G4ThreeVector(p.x()/rho, p.y()/rho, 0) ;
713 break ;
714 }
715 case kNZ : // + or - dz
716 {
717 if ( distZHigh > distZLow ) { norm = fHighNorm ; }
718 else { norm = fLowNorm; }
719 break ;
720 }
721 case kNSPhi:
722 {
723 norm = G4ThreeVector(std::sin(fSPhi), -std::cos(fSPhi), 0) ;
724 break ;
725 }
726 case kNEPhi:
727 {
728 norm = G4ThreeVector(-std::sin(fSPhi+fDPhi), std::cos(fSPhi+fDPhi), 0) ;
729 break;
730 }
731 default: // Should never reach this case ...
732 {
733 DumpInfo();
734 G4Exception("G4CutTubs::ApproxSurfaceNormal()",
735 "GeomSolids1002", JustWarning,
736 "Undefined side for valid surface normal to solid.");
737 break ;
738 }
739 }
740 return norm;
741}
ENorm
ENorm
Definition: G4Cons.cc:72
G4double
double G4double
Definition: G4Types.hh:64
G4Tubs::fRMin
G4double fRMin
Definition: G4Tubs.hh:209
G4Tubs::fDPhi
G4double fDPhi
Definition: G4Tubs.hh:209
G4Tubs::fRMax
G4double fRMax
Definition: G4Tubs.hh:209
G4Tubs::kNRMax
@ kNRMax
Definition: G4Tubs.hh:203
G4Tubs::kNZ
@ kNZ
Definition: G4Tubs.hh:203
G4Tubs::kNRMin
@ kNRMin
Definition: G4Tubs.hh:203
G4Tubs::kNEPhi
@ kNEPhi
Definition: G4Tubs.hh:203
G4Tubs::kNSPhi
@ kNSPhi
Definition: G4Tubs.hh:203
G4Tubs::fDz
G4double fDz
Definition: G4Tubs.hh:209
G4Tubs::fSPhi
G4double fSPhi
Definition: G4Tubs.hh:209
G4VSolid::DumpInfo
void DumpInfo() const

Referenced by SurfaceNormal().

◆ CalculateExtent()

G4bool G4CutTubs::CalculateExtent ( const EAxis  pAxis,
const G4VoxelLimits &  pVoxelLimit,
const G4AffineTransform &  pTransform,
G4double &  pmin,
G4double &  pmax 
) const
virtual

Implements G4VSolid.

Definition at line 194 of file G4CutTubs.cc.

199{
200 if ( (!pTransform.IsRotated()) && (fDPhi == twopi) && (fRMin == 0) )
201 {
202 // Special case handling for unrotated solid tubes
203 // Compute x/y/z mins and maxs fro bounding box respecting limits,
204 // with early returns if outside limits. Then switch() on pAxis,
205 // and compute exact x and y limit for x/y case
206
207 G4double xoffset, xMin, xMax;
208 G4double yoffset, yMin, yMax;
209 G4double zoffset, zMin, zMax;
210
211 G4double diff1, diff2, maxDiff, newMin, newMax;
212 G4double xoff1, xoff2, yoff1, yoff2, delta;
213
214 xoffset = pTransform.NetTranslation().x();
215 xMin = xoffset - fRMax;
216 xMax = xoffset + fRMax;
217
218 if (pVoxelLimit.IsXLimited())
219 {
220 if ( (xMin > pVoxelLimit.GetMaxXExtent())
221 || (xMax < pVoxelLimit.GetMinXExtent()) )
222 {
223 return false;
224 }
225 else
226 {
227 if (xMin < pVoxelLimit.GetMinXExtent())
228 {
229 xMin = pVoxelLimit.GetMinXExtent();
230 }
231 if (xMax > pVoxelLimit.GetMaxXExtent())
232 {
233 xMax = pVoxelLimit.GetMaxXExtent();
234 }
235 }
236 }
237 yoffset = pTransform.NetTranslation().y();
238 yMin = yoffset - fRMax;
239 yMax = yoffset + fRMax;
240
241 if ( pVoxelLimit.IsYLimited() )
242 {
243 if ( (yMin > pVoxelLimit.GetMaxYExtent())
244 || (yMax < pVoxelLimit.GetMinYExtent()) )
245 {
246 return false;
247 }
248 else
249 {
250 if (yMin < pVoxelLimit.GetMinYExtent())
251 {
252 yMin = pVoxelLimit.GetMinYExtent();
253 }
254 if (yMax > pVoxelLimit.GetMaxYExtent())
255 {
256 yMax=pVoxelLimit.GetMaxYExtent();
257 }
258 }
259 }
260 zoffset = pTransform.NetTranslation().z();
261 GetMaxMinZ(zMin,zMax);
262 zMin += zoffset;
263 zMax += zoffset;
264
265 if ( pVoxelLimit.IsZLimited() )
266 {
267 if ( (zMin > pVoxelLimit.GetMaxZExtent())
268 || (zMax < pVoxelLimit.GetMinZExtent()) )
269 {
270 return false;
271 }
272 else
273 {
274 if (zMin < pVoxelLimit.GetMinZExtent())
275 {
276 zMin = pVoxelLimit.GetMinZExtent();
277 }
278 if (zMax > pVoxelLimit.GetMaxZExtent())
279 {
280 zMax = pVoxelLimit.GetMaxZExtent();
281 }
282 }
283 }
284 switch ( pAxis ) // Known to cut cylinder
285 {
286 case kXAxis :
287 {
288 yoff1 = yoffset - yMin;
289 yoff2 = yMax - yoffset;
290
291 if ( (yoff1 >= 0) && (yoff2 >= 0) ) // Y limits cross max/min x
292 { // => no change
293 pMin = xMin;
294 pMax = xMax;
295 }
296 else
297 {
298 // Y limits don't cross max/min x => compute max delta x,
299 // hence new mins/maxs
300
301 delta = fRMax*fRMax - yoff1*yoff1;
302 diff1 = (delta>0.) ? std::sqrt(delta) : 0.;
303 delta = fRMax*fRMax - yoff2*yoff2;
304 diff2 = (delta>0.) ? std::sqrt(delta) : 0.;
305 maxDiff = (diff1 > diff2) ? diff1:diff2;
306 newMin = xoffset - maxDiff;
307 newMax = xoffset + maxDiff;
308 pMin = (newMin < xMin) ? xMin : newMin;
309 pMax = (newMax > xMax) ? xMax : newMax;
310 }
311 break;
312 }
313 case kYAxis :
314 {
315 xoff1 = xoffset - xMin;
316 xoff2 = xMax - xoffset;
317
318 if ( (xoff1 >= 0) && (xoff2 >= 0) ) // X limits cross max/min y
319 { // => no change
320 pMin = yMin;
321 pMax = yMax;
322 }
323 else
324 {
325 // X limits don't cross max/min y => compute max delta y,
326 // hence new mins/maxs
327
328 delta = fRMax*fRMax - xoff1*xoff1;
329 diff1 = (delta>0.) ? std::sqrt(delta) : 0.;
330 delta = fRMax*fRMax - xoff2*xoff2;
331 diff2 = (delta>0.) ? std::sqrt(delta) : 0.;
332 maxDiff = (diff1 > diff2) ? diff1 : diff2;
333 newMin = yoffset - maxDiff;
334 newMax = yoffset + maxDiff;
335 pMin = (newMin < yMin) ? yMin : newMin;
336 pMax = (newMax > yMax) ? yMax : newMax;
337 }
338 break;
339 }
340 case kZAxis:
341 {
342 pMin = zMin;
343 pMax = zMax;
344 break;
345 }
346 default:
347 break;
348 }
349 pMin -= kCarTolerance;
350 pMax += kCarTolerance;
351 return true;
352 }
353 else // Calculate rotated vertex coordinates
354 {
355 G4int i, noEntries, noBetweenSections4;
356 G4bool existsAfterClip = false;
357 G4ThreeVectorList* vertices = CreateRotatedVertices(pTransform);
358
359 pMin = kInfinity;
360 pMax = -kInfinity;
361
362 noEntries = vertices->size();
363 noBetweenSections4 = noEntries - 4;
364
365 for ( i = 0 ; i < noEntries ; i += 4 )
366 {
367 ClipCrossSection(vertices, i, pVoxelLimit, pAxis, pMin, pMax);
368 }
369 for ( i = 0 ; i < noBetweenSections4 ; i += 4 )
370 {
371 ClipBetweenSections(vertices, i, pVoxelLimit, pAxis, pMin, pMax);
372 }
373 if ( (pMin != kInfinity) || (pMax != -kInfinity) )
374 {
375 existsAfterClip = true;
376 pMin -= kCarTolerance; // Add 2*tolerance to avoid precision troubles
377 pMax += kCarTolerance;
378 }
379 else
380 {
381 // Check for case where completely enveloping clipping volume
382 // If point inside then we are confident that the solid completely
383 // envelopes the clipping volume. Hence set min/max extents according
384 // to clipping volume extents along the specified axis.
385
386 G4ThreeVector clipCentre(
387 (pVoxelLimit.GetMinXExtent()+pVoxelLimit.GetMaxXExtent())*0.5,
388 (pVoxelLimit.GetMinYExtent()+pVoxelLimit.GetMaxYExtent())*0.5,
389 (pVoxelLimit.GetMinZExtent()+pVoxelLimit.GetMaxZExtent())*0.5 );
390
391 if ( Inside(pTransform.Inverse().TransformPoint(clipCentre)) != kOutside )
392 {
393 existsAfterClip = true;
394 pMin = pVoxelLimit.GetMinExtent(pAxis);
395 pMax = pVoxelLimit.GetMaxExtent(pAxis);
396 }
397 }
398 delete vertices;
399 return existsAfterClip;
400 }
401}
G4int
int G4int
Definition: G4Types.hh:66
G4bool
bool G4bool
Definition: G4Types.hh:67
G4ThreeVectorList
std::vector< G4ThreeVector > G4ThreeVectorList
Definition: G4VSolid.hh:85
G4AffineTransform::IsRotated
G4bool IsRotated() const
G4AffineTransform::Inverse
G4AffineTransform Inverse() const
G4AffineTransform::NetTranslation
G4ThreeVector NetTranslation() const
G4AffineTransform::TransformPoint
G4ThreeVector TransformPoint(const G4ThreeVector &vec) const
G4CutTubs::Inside
EInside Inside(const G4ThreeVector &p) const
Definition: G4CutTubs.cc:407
G4CutTubs::GetMaxMinZ
void GetMaxMinZ(G4double &zmin, G4double &zmax) const
Definition: G4CutTubs.cc:2082
G4CutTubs::CreateRotatedVertices
G4ThreeVectorList * CreateRotatedVertices(const G4AffineTransform &pTransform) const
Definition: G4CutTubs.cc:1767
G4VSolid::ClipBetweenSections
void ClipBetweenSections(G4ThreeVectorList *pVertices, const G4int pSectionIndex, const G4VoxelLimits &pVoxelLimit, const EAxis pAxis, G4double &pMin, G4double &pMax) const
Definition: G4VSolid.cc:376
G4VSolid::kCarTolerance
G4double kCarTolerance
Definition: G4VSolid.hh:307
G4VSolid::ClipCrossSection
void ClipCrossSection(G4ThreeVectorList *pVertices, const G4int pSectionIndex, const G4VoxelLimits &pVoxelLimit, const EAxis pAxis, G4double &pMin, G4double &pMax) const
Definition: G4VSolid.cc:345
G4VoxelLimits::GetMinExtent
G4double GetMinExtent(const EAxis pAxis) const
G4VoxelLimits::IsYLimited
G4bool IsYLimited() const
G4VoxelLimits::GetMinZExtent
G4double GetMinZExtent() const
G4VoxelLimits::IsXLimited
G4bool IsXLimited() const
G4VoxelLimits::GetMaxExtent
G4double GetMaxExtent(const EAxis pAxis) const
G4VoxelLimits::GetMaxYExtent
G4double GetMaxYExtent() const
G4VoxelLimits::GetMaxZExtent
G4double GetMaxZExtent() const
G4VoxelLimits::GetMinYExtent
G4double GetMinYExtent() const
G4VoxelLimits::GetMinXExtent
G4double GetMinXExtent() const
G4VoxelLimits::IsZLimited
G4bool IsZLimited() const
G4VoxelLimits::GetMaxXExtent
G4double GetMaxXExtent() const
kYAxis
@ kYAxis
Definition: geomdefs.hh:54
kXAxis
@ kXAxis
Definition: geomdefs.hh:54
kZAxis
@ kZAxis
Definition: geomdefs.hh:54
kOutside
@ kOutside
Definition: geomdefs.hh:58

◆ Clone()

G4VSolid * G4CutTubs::Clone ( ) const
virtual

Reimplemented from G4VSolid.

Definition at line 1862 of file G4CutTubs.cc.

1863{
1864 return new G4CutTubs(*this);
1865}

◆ CreatePolyhedron()

G4Polyhedron * G4CutTubs::CreatePolyhedron ( ) const
virtual

Reimplemented from G4VSolid.

Definition at line 1975 of file G4CutTubs.cc.

1976{
1977 typedef G4double G4double3[3];
1978 typedef G4int G4int4[4];
1979
1980 G4Polyhedron *ph = new G4Polyhedron;
1981 G4Polyhedron *ph1 = G4Tubs::CreatePolyhedron();
1982 G4int nn=ph1->GetNoVertices();
1983 G4int nf=ph1->GetNoFacets();
1984 G4double3* xyz = new G4double3[nn]; // number of nodes
1985 G4int4* faces = new G4int4[nf] ; // number of faces
1986
1987 for(G4int i=0;i<nn;i++)
1988 {
1989 xyz[i][0]=ph1->GetVertex(i+1).x();
1990 xyz[i][1]=ph1->GetVertex(i+1).y();
1991 G4double tmpZ=ph1->GetVertex(i+1).z();
1992 if(tmpZ>=fDz-kCarTolerance)
1993 {
1994 xyz[i][2]=GetCutZ(G4ThreeVector(xyz[i][0],xyz[i][1],fDz));
1995 }
1996 else if(tmpZ<=-fDz+kCarTolerance)
1997 {
1998 xyz[i][2]=GetCutZ(G4ThreeVector(xyz[i][0],xyz[i][1],-fDz));
1999 }
2000 else
2001 {
2002 xyz[i][2]=tmpZ;
2003 }
2004 }
2005 G4int iNodes[4];
2006 G4int *iEdge=0;
2007 G4int n;
2008 for(G4int i=0;i<nf;i++)
2009 {
2010 ph1->GetFacet(i+1,n,iNodes,iEdge);
2011 for(G4int k=0;k<n;k++)
2012 {
2013 faces[i][k]=iNodes[k];
2014 }
2015 for(G4int k=n;k<4;k++)
2016 {
2017 faces[i][k]=0;
2018 }
2019 }
2020 ph->createPolyhedron(nn,nf,xyz,faces);
2021
2022 delete [] xyz;
2023 delete [] faces;
2024 delete ph1;
2025
2026 return ph;
2027}
G4CutTubs::GetCutZ
G4double GetCutZ(const G4ThreeVector &p) const
Definition: G4CutTubs.cc:2058
G4Tubs::CreatePolyhedron
G4Polyhedron * CreatePolyhedron() const
Definition: G4Tubs.cc:1923
HepPolyhedron::GetFacet
void GetFacet(G4int iFace, G4int &n, G4int *iNodes, G4int *edgeFlags=0, G4int *iFaces=0) const
Definition: HepPolyhedron.cc:1063
HepPolyhedron::GetVertex
G4Point3D GetVertex(G4int index) const
Definition: HepPolyhedron.cc:883
HepPolyhedron::createPolyhedron
G4int createPolyhedron(G4int Nnodes, G4int Nfaces, const G4double xyz[][3], const G4int faces[][4])
Definition: HepPolyhedron.cc:1356
HepPolyhedron::GetNoFacets
G4int GetNoFacets() const
Definition: HepPolyhedron.h:255
HepPolyhedron::GetNoVertices
G4int GetNoVertices() const
Definition: HepPolyhedron.h:252

◆ CreateRotatedVertices()

G4ThreeVectorList * G4CutTubs::CreateRotatedVertices ( const G4AffineTransform &  pTransform) const
protected

Definition at line 1767 of file G4CutTubs.cc.

1768{
1769 G4ThreeVectorList* vertices ;
1770 G4ThreeVector vertex0, vertex1, vertex2, vertex3 ;
1771 G4double meshAngle, meshRMax, crossAngle,
1772 cosCrossAngle, sinCrossAngle, sAngle;
1773 G4double rMaxX, rMaxY, rMinX, rMinY, meshRMin ;
1774 G4int crossSection, noCrossSections;
1775
1776 // Compute no of cross-sections necessary to mesh tube
1777 //
1778 noCrossSections = G4int(fDPhi/kMeshAngleDefault) + 1 ;
1779
1780 if ( noCrossSections < kMinMeshSections )
1781 {
1782 noCrossSections = kMinMeshSections ;
1783 }
1784 else if (noCrossSections>kMaxMeshSections)
1785 {
1786 noCrossSections = kMaxMeshSections ;
1787 }
1788 // noCrossSections = 4 ;
1789
1790 meshAngle = fDPhi/(noCrossSections - 1) ;
1791 // meshAngle = fDPhi/(noCrossSections) ;
1792
1793 meshRMax = (fRMax+100*kCarTolerance)/std::cos(meshAngle*0.5) ;
1794 meshRMin = fRMin - 100*kCarTolerance ;
1795
1796 // If complete in phi, set start angle such that mesh will be at fRMax
1797 // on the x axis. Will give better extent calculations when not rotated.
1798
1799 if (fPhiFullCutTube && (fSPhi == 0) ) { sAngle = -meshAngle*0.5 ; }
1800 else { sAngle = fSPhi ; }
1801
1802 vertices = new G4ThreeVectorList();
1803
1804 if ( vertices )
1805 {
1806 vertices->reserve(noCrossSections*4);
1807 for (crossSection = 0 ; crossSection < noCrossSections ; crossSection++ )
1808 {
1809 // Compute coordinates of cross section at section crossSection
1810
1811 crossAngle = sAngle + crossSection*meshAngle ;
1812 cosCrossAngle = std::cos(crossAngle) ;
1813 sinCrossAngle = std::sin(crossAngle) ;
1814
1815 rMaxX = meshRMax*cosCrossAngle ;
1816 rMaxY = meshRMax*sinCrossAngle ;
1817
1818 if(meshRMin <= 0.0)
1819 {
1820 rMinX = 0.0 ;
1821 rMinY = 0.0 ;
1822 }
1823 else
1824 {
1825 rMinX = meshRMin*cosCrossAngle ;
1826 rMinY = meshRMin*sinCrossAngle ;
1827 }
1828 vertex0 = G4ThreeVector(rMinX,rMinY,GetCutZ(G4ThreeVector(rMinX,rMinY,-fDz))) ;
1829 vertex1 = G4ThreeVector(rMaxX,rMaxY,GetCutZ(G4ThreeVector(rMaxX,rMaxY,-fDz))) ;
1830 vertex2 = G4ThreeVector(rMaxX,rMaxY,GetCutZ(G4ThreeVector(rMaxX,rMaxY,+fDz))) ;
1831 vertex3 = G4ThreeVector(rMinX,rMinY,GetCutZ(G4ThreeVector(rMinX,rMinY,+fDz))) ;
1832
1833 vertices->push_back(pTransform.TransformPoint(vertex0)) ;
1834 vertices->push_back(pTransform.TransformPoint(vertex1)) ;
1835 vertices->push_back(pTransform.TransformPoint(vertex2)) ;
1836 vertices->push_back(pTransform.TransformPoint(vertex3)) ;
1837 }
1838 }
1839 else
1840 {
1841 DumpInfo();
1842 G4Exception("G4CutTubs::CreateRotatedVertices()",
1843 "GeomSolids0003", FatalException,
1844 "Error in allocation of vertices. Out of memory !");
1845 }
1846 return vertices ;
1847}
kMaxMeshSections
const G4int kMaxMeshSections
Definition: meshdefs.hh:46
kMinMeshSections
const G4int kMinMeshSections
Definition: meshdefs.hh:45
kMeshAngleDefault
const G4double kMeshAngleDefault
Definition: meshdefs.hh:42

Referenced by CalculateExtent().

◆ DescribeYourselfTo()

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

Implements G4VSolid.

Definition at line 1970 of file G4CutTubs.cc.

1971{
1972 scene.AddSolid (*this) ;
1973}
G4VGraphicsScene::AddSolid
virtual void AddSolid(const G4Box &)=0

◆ DistanceToIn() [1/2]

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

Implements G4VSolid.

Definition at line 1221 of file G4CutTubs.cc.

1222{
1223 G4double safRMin,safRMax,safZLow,safZHigh,safePhi,safe,rho,cosPsi;
1224 G4ThreeVector vZ=G4ThreeVector(0,0,fDz);
1225
1226 // Distance to R
1227 //
1228 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
1229
1230 safRMin = fRMin- rho ;
1231 safRMax = rho - fRMax ;
1232
1233 // Distances to ZCut(Low/High)
1234
1235 // Dist to Low Cut
1236 //
1237 safZLow = (p+vZ).dot(fLowNorm);
1238
1239 // Dist to High Cut
1240 //
1241 safZHigh = (p-vZ).dot(fHighNorm);
1242
1243 safe = std::max(safZLow,safZHigh);
1244
1245 if ( safRMin > safe ) { safe = safRMin; }
1246 if ( safRMax> safe ) { safe = safRMax; }
1247
1248 // Distance to Phi
1249 //
1250 if ( (!fPhiFullCutTube) && (rho) )
1251 {
1252 // Psi=angle from central phi to point
1253 //
1254 cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/rho ;
1255
1256 if ( cosPsi < std::cos(fDPhi*0.5) )
1257 {
1258 // Point lies outside phi range
1259
1260 if ( (p.y()*cosCPhi - p.x()*sinCPhi) <= 0 )
1261 {
1262 safePhi = std::fabs(p.x()*sinSPhi - p.y()*cosSPhi) ;
1263 }
1264 else
1265 {
1266 safePhi = std::fabs(p.x()*sinEPhi - p.y()*cosEPhi) ;
1267 }
1268 if ( safePhi > safe ) { safe = safePhi; }
1269 }
1270 }
1271 if ( safe < 0 ) { safe = 0; }
1272
1273 return safe ;
1274}
G4Tubs::sinCPhi
G4double sinCPhi
Definition: G4Tubs.hh:213
G4Tubs::cosEPhi
G4double cosEPhi
Definition: G4Tubs.hh:214
G4Tubs::sinSPhi
G4double sinSPhi
Definition: G4Tubs.hh:214
G4Tubs::cosCPhi
G4double cosCPhi
Definition: G4Tubs.hh:213
G4Tubs::cosSPhi
G4double cosSPhi
Definition: G4Tubs.hh:214
G4Tubs::sinEPhi
G4double sinEPhi
Definition: G4Tubs.hh:214

◆ DistanceToIn() [2/2]

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

Implements G4VSolid.

Definition at line 765 of file G4CutTubs.cc.

767{
768 G4double snxt = kInfinity ; // snxt = default return value
769 G4double tolORMin2, tolIRMax2 ; // 'generous' radii squared
770 G4double tolORMax2, tolIRMin2;
771 const G4double dRmax = 100.*fRMax;
772 G4ThreeVector vZ=G4ThreeVector(0,0,fDz);
773
774 static const G4double halfCarTolerance = 0.5*kCarTolerance;
775 static const G4double halfRadTolerance = 0.5*kRadTolerance;
776
777 // Intersection point variables
778 //
779 G4double Dist, sd=0, xi, yi, zi, rho2, inum, iden, cosPsi, Comp,calf ;
780 G4double t1, t2, t3, b, c, d ; // Quadratic solver variables
781 G4double distZLow,distZHigh;
782 // Calculate tolerant rmin and rmax
783
784 if (fRMin > kRadTolerance)
785 {
786 tolORMin2 = (fRMin - halfRadTolerance)*(fRMin - halfRadTolerance) ;
787 tolIRMin2 = (fRMin + halfRadTolerance)*(fRMin + halfRadTolerance) ;
788 }
789 else
790 {
791 tolORMin2 = 0.0 ;
792 tolIRMin2 = 0.0 ;
793 }
794 tolORMax2 = (fRMax + halfRadTolerance)*(fRMax + halfRadTolerance) ;
795 tolIRMax2 = (fRMax - halfRadTolerance)*(fRMax - halfRadTolerance) ;
796
797 // Intersection with ZCut surfaces
798
799 // dist to Low Cut
800 //
801 distZLow =(p+vZ).dot(fLowNorm);
802
803 // dist to High Cut
804 //
805 distZHigh = (p-vZ).dot(fHighNorm);
806
807 if ( distZLow >= -halfCarTolerance )
808 {
809 calf = v.dot(fLowNorm);
810 if (calf<0)
811 {
812 sd = -distZLow/calf;
813 if(sd < 0.0) { sd = 0.0; }
814
815 xi = p.x() + sd*v.x() ; // Intersection coords
816 yi = p.y() + sd*v.y() ;
817 rho2 = xi*xi + yi*yi ;
818
819 // Check validity of intersection
820
821 if ((tolIRMin2 <= rho2) && (rho2 <= tolIRMax2))
822 {
823 if (!fPhiFullCutTube && rho2)
824 {
825 // Psi = angle made with central (average) phi of shape
826 //
827 inum = xi*cosCPhi + yi*sinCPhi ;
828 iden = std::sqrt(rho2) ;
829 cosPsi = inum/iden ;
830 if (cosPsi >= cosHDPhiIT) { return sd ; }
831 }
832 else
833 {
834 return sd ;
835 }
836 }
837 }
838 else
839 {
840 if ( sd<halfCarTolerance )
841 {
842 if(calf>=0) { sd=kInfinity; }
843 return sd ; // On/outside extent, and heading away
844 } // -> cannot intersect
845 }
846 }
847
848 if(distZHigh >= -halfCarTolerance )
849 {
850 calf = v.dot(fHighNorm);
851 if (calf<0)
852 {
853 sd = -distZHigh/calf;
854
855 if(sd < 0.0) { sd = 0.0; }
856
857 xi = p.x() + sd*v.x() ; // Intersection coords
858 yi = p.y() + sd*v.y() ;
859 rho2 = xi*xi + yi*yi ;
860
861 // Check validity of intersection
862
863 if ((tolIRMin2 <= rho2) && (rho2 <= tolIRMax2))
864 {
865 if (!fPhiFullCutTube && rho2)
866 {
867 // Psi = angle made with central (average) phi of shape
868 //
869 inum = xi*cosCPhi + yi*sinCPhi ;
870 iden = std::sqrt(rho2) ;
871 cosPsi = inum/iden ;
872 if (cosPsi >= cosHDPhiIT) { return sd ; }
873 }
874 else
875 {
876 return sd ;
877 }
878 }
879 }
880 else
881 {
882 if ( sd<halfCarTolerance )
883 {
884 if(calf>=0) { sd=kInfinity; }
885 return sd ; // On/outside extent, and heading away
886 } // -> cannot intersect
887 }
888 }
889
890 // -> Can not intersect z surfaces
891 //
892 // Intersection with rmax (possible return) and rmin (must also check phi)
893 //
894 // Intersection point (xi,yi,zi) on line x=p.x+t*v.x etc.
895 //
896 // Intersects with x^2+y^2=R^2
897 //
898 // Hence (v.x^2+v.y^2)t^2+ 2t(p.x*v.x+p.y*v.y)+p.x^2+p.y^2-R^2=0
899 // t1 t2 t3
900
901 t1 = 1.0 - v.z()*v.z() ;
902 t2 = p.x()*v.x() + p.y()*v.y() ;
903 t3 = p.x()*p.x() + p.y()*p.y() ;
904 if ( t1 > 0 ) // Check not || to z axis
905 {
906 b = t2/t1 ;
907 c = t3 - fRMax*fRMax ;
908
909 if ((t3 >= tolORMax2) && (t2<0)) // This also handles the tangent case
910 {
911 // Try outer cylinder intersection, c=(t3-fRMax*fRMax)/t1;
912
913 c /= t1 ;
914 d = b*b - c ;
915
916 if (d >= 0) // If real root
917 {
918 sd = c/(-b+std::sqrt(d));
919 if (sd >= 0) // If 'forwards'
920 {
921 if ( sd>dRmax ) // Avoid rounding errors due to precision issues on
922 { // 64 bits systems. Split long distances and recompute
923 G4double fTerm = sd-std::fmod(sd,dRmax);
924 sd = fTerm + DistanceToIn(p+fTerm*v,v);
925 }
926 // Check z intersection
927 //
928 zi = p.z() + sd*v.z() ;
929 xi = p.x() + sd*v.x() ;
930 yi = p.y() + sd*v.y() ;
931 if ((-xi*fLowNorm.x()-yi*fLowNorm.y()
932 -(zi+fDz)*fLowNorm.z())>-halfCarTolerance)
933 {
934 if ((-xi*fHighNorm.x()-yi*fHighNorm.y()
935 +(fDz-zi)*fHighNorm.z())>-halfCarTolerance)
936 {
937 // Z ok. Check phi intersection if reqd
938 //
939 if (fPhiFullCutTube)
940 {
941 return sd ;
942 }
943 else
944 {
945 xi = p.x() + sd*v.x() ;
946 yi = p.y() + sd*v.y() ;
947 cosPsi = (xi*cosCPhi + yi*sinCPhi)/fRMax ;
948 if (cosPsi >= cosHDPhiIT) { return sd ; }
949 }
950 } // end if std::fabs(zi)
951 }
952 } // end if (sd>=0)
953 } // end if (d>=0)
954 } // end if (r>=fRMax)
955 else
956 {
957 // Inside outer radius :
958 // check not inside, and heading through tubs (-> 0 to in)
959 if ((t3 > tolIRMin2) && (t2 < 0)
960 && (std::fabs(p.z()) <= std::fabs(GetCutZ(p))-halfCarTolerance ))
961 {
962 // Inside both radii, delta r -ve, inside z extent
963
964 if (!fPhiFullCutTube)
965 {
966 inum = p.x()*cosCPhi + p.y()*sinCPhi ;
967 iden = std::sqrt(t3) ;
968 cosPsi = inum/iden ;
969 if (cosPsi >= cosHDPhiIT)
970 {
971 // In the old version, the small negative tangent for the point
972 // on surface was not taken in account, and returning 0.0 ...
973 // New version: check the tangent for the point on surface and
974 // if no intersection, return kInfinity, if intersection instead
975 // return sd.
976 //
977 c = t3-fRMax*fRMax;
978 if ( c<=0.0 )
979 {
980 return 0.0;
981 }
982 else
983 {
984 c = c/t1 ;
985 d = b*b-c;
986 if ( d>=0.0 )
987 {
988 snxt = c/(-b+std::sqrt(d)); // using safe solution
989 // for quadratic equation
990 if ( snxt < halfCarTolerance ) { snxt=0; }
991 return snxt ;
992 }
993 else
994 {
995 return kInfinity;
996 }
997 }
998 }
999 }
1000 else
1001 {
1002 // In the old version, the small negative tangent for the point
1003 // on surface was not taken in account, and returning 0.0 ...
1004 // New version: check the tangent for the point on surface and
1005 // if no intersection, return kInfinity, if intersection instead
1006 // return sd.
1007 //
1008 c = t3 - fRMax*fRMax;
1009 if ( c<=0.0 )
1010 {
1011 return 0.0;
1012 }
1013 else
1014 {
1015 c = c/t1 ;
1016 d = b*b-c;
1017 if ( d>=0.0 )
1018 {
1019 snxt= c/(-b+std::sqrt(d)); // using safe solution
1020 // for quadratic equation
1021 if ( snxt < halfCarTolerance ) { snxt=0; }
1022 return snxt ;
1023 }
1024 else
1025 {
1026 return kInfinity;
1027 }
1028 }
1029 } // end if (!fPhiFullCutTube)
1030 } // end if (t3>tolIRMin2)
1031 } // end if (Inside Outer Radius)
1032
1033 if ( fRMin ) // Try inner cylinder intersection
1034 {
1035 c = (t3 - fRMin*fRMin)/t1 ;
1036 d = b*b - c ;
1037 if ( d >= 0.0 ) // If real root
1038 {
1039 // Always want 2nd root - we are outside and know rmax Hit was bad
1040 // - If on surface of rmin also need farthest root
1041
1042 sd =( b > 0. )? c/(-b - std::sqrt(d)) : (-b + std::sqrt(d));
1043 if (sd >= -10*halfCarTolerance) // check forwards
1044 {
1045 // Check z intersection
1046 //
1047 if (sd < 0.0) { sd = 0.0; }
1048 if (sd>dRmax) // Avoid rounding errors due to precision issues seen
1049 { // 64 bits systems. Split long distances and recompute
1050 G4double fTerm = sd-std::fmod(sd,dRmax);
1051 sd = fTerm + DistanceToIn(p+fTerm*v,v);
1052 }
1053 zi = p.z() + sd*v.z() ;
1054 xi = p.x() + sd*v.x() ;
1055 yi = p.y() + sd*v.y() ;
1056 if ((-xi*fLowNorm.x()-yi*fLowNorm.y()
1057 -(zi+fDz)*fLowNorm.z())>-halfCarTolerance)
1058 {
1059 if ((-xi*fHighNorm.x()-yi*fHighNorm.y()
1060 +(fDz-zi)*fHighNorm.z())>-halfCarTolerance)
1061 {
1062 // Z ok. Check phi
1063 //
1064 if ( fPhiFullCutTube )
1065 {
1066 return sd ;
1067 }
1068 else
1069 {
1070 cosPsi = (xi*cosCPhi + yi*sinCPhi)/fRMin ;
1071 if (cosPsi >= cosHDPhiIT)
1072 {
1073 // Good inner radius isect
1074 // - but earlier phi isect still possible
1075 //
1076 snxt = sd ;
1077 }
1078 }
1079 } // end if std::fabs(zi)
1080 }
1081 } // end if (sd>=0)
1082 } // end if (d>=0)
1083 } // end if (fRMin)
1084 }
1085
1086 // Phi segment intersection
1087 //
1088 // o Tolerant of points inside phi planes by up to kCarTolerance*0.5
1089 //
1090 // o NOTE: Large duplication of code between sphi & ephi checks
1091 // -> only diffs: sphi -> ephi, Comp -> -Comp and half-plane
1092 // intersection check <=0 -> >=0
1093 // -> use some form of loop Construct ?
1094 //
1095 if ( !fPhiFullCutTube )
1096 {
1097 // First phi surface (Starting phi)
1098 //
1099 Comp = v.x()*sinSPhi - v.y()*cosSPhi ;
1100
1101 if ( Comp < 0 ) // Component in outwards normal dirn
1102 {
1103 Dist = (p.y()*cosSPhi - p.x()*sinSPhi) ;
1104
1105 if ( Dist < halfCarTolerance )
1106 {
1107 sd = Dist/Comp ;
1108
1109 if (sd < snxt)
1110 {
1111 if ( sd < 0 ) { sd = 0.0; }
1112 zi = p.z() + sd*v.z() ;
1113 xi = p.x() + sd*v.x() ;
1114 yi = p.y() + sd*v.y() ;
1115 if ((-xi*fLowNorm.x()-yi*fLowNorm.y()
1116 -(zi+fDz)*fLowNorm.z())>-halfCarTolerance)
1117 {
1118 if ((-xi*fHighNorm.x()-yi*fHighNorm.y()
1119 +(fDz-zi)*fHighNorm.z())>-halfCarTolerance)
1120 {
1121 rho2 = xi*xi + yi*yi ;
1122 if ( ( (rho2 >= tolIRMin2) && (rho2 <= tolIRMax2) )
1123 || ( (rho2 > tolORMin2) && (rho2 < tolIRMin2)
1124 && ( v.y()*cosSPhi - v.x()*sinSPhi > 0 )
1125 && ( v.x()*cosSPhi + v.y()*sinSPhi >= 0 ) )
1126 || ( (rho2 > tolIRMax2) && (rho2 < tolORMax2)
1127 && (v.y()*cosSPhi - v.x()*sinSPhi > 0)
1128 && (v.x()*cosSPhi + v.y()*sinSPhi < 0) ) )
1129 {
1130 // z and r intersections good
1131 // - check intersecting with correct half-plane
1132 //
1133 if ((yi*cosCPhi-xi*sinCPhi) <= halfCarTolerance) { snxt = sd; }
1134 }
1135 } //two Z conditions
1136 }
1137 }
1138 }
1139 }
1140
1141 // Second phi surface (Ending phi)
1142 //
1143 Comp = -(v.x()*sinEPhi - v.y()*cosEPhi) ;
1144
1145 if (Comp < 0 ) // Component in outwards normal dirn
1146 {
1147 Dist = -(p.y()*cosEPhi - p.x()*sinEPhi) ;
1148
1149 if ( Dist < halfCarTolerance )
1150 {
1151 sd = Dist/Comp ;
1152
1153 if (sd < snxt)
1154 {
1155 if ( sd < 0 ) { sd = 0; }
1156 zi = p.z() + sd*v.z() ;
1157 xi = p.x() + sd*v.x() ;
1158 yi = p.y() + sd*v.y() ;
1159 if ((-xi*fLowNorm.x()-yi*fLowNorm.y()
1160 -(zi+fDz)*fLowNorm.z())>-halfCarTolerance)
1161 {
1162 if ((-xi*fHighNorm.x()-yi*fHighNorm.y()
1163 +(fDz-zi)*fHighNorm.z())>-halfCarTolerance)
1164 {
1165 xi = p.x() + sd*v.x() ;
1166 yi = p.y() + sd*v.y() ;
1167 rho2 = xi*xi + yi*yi ;
1168 if ( ( (rho2 >= tolIRMin2) && (rho2 <= tolIRMax2) )
1169 || ( (rho2 > tolORMin2) && (rho2 < tolIRMin2)
1170 && (v.x()*sinEPhi - v.y()*cosEPhi > 0)
1171 && (v.x()*cosEPhi + v.y()*sinEPhi >= 0) )
1172 || ( (rho2 > tolIRMax2) && (rho2 < tolORMax2)
1173 && (v.x()*sinEPhi - v.y()*cosEPhi > 0)
1174 && (v.x()*cosEPhi + v.y()*sinEPhi < 0) ) )
1175 {
1176 // z and r intersections good
1177 // - check intersecting with correct half-plane
1178 //
1179 if ( (yi*cosCPhi-xi*sinCPhi) >= -halfCarTolerance )
1180 {
1181 snxt = sd;
1182 }
1183 } //?? >=-halfCarTolerance
1184 }
1185 } // two Z conditions
1186 }
1187 }
1188 } // Comp < 0
1189 } // !fPhiFullTube
1190 if ( snxt<halfCarTolerance ) { snxt=0; }
1191
1192 return snxt ;
1193}
CLHEP::Hep3Vector::dot
double dot(const Hep3Vector &) const
G4CutTubs::DistanceToIn
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
Definition: G4CutTubs.cc:765
G4Tubs::cosHDPhiIT
G4double cosHDPhiIT
Definition: G4Tubs.hh:213

Referenced by DistanceToIn().

◆ DistanceToOut() [1/2]

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

Implements G4VSolid.

Definition at line 1714 of file G4CutTubs.cc.

1715{
1716 G4double safRMin,safRMax,safZLow,safZHigh,safePhi,safe,rho;
1717 G4ThreeVector vZ=G4ThreeVector(0,0,fDz);
1718
1719 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ; // Distance to R
1720
1721 safRMin = rho - fRMin ;
1722 safRMax = fRMax - rho ;
1723
1724 // Distances to ZCut(Low/High)
1725
1726 // Dist to Low Cut
1727 //
1728 safZLow = std::fabs((p+vZ).dot(fLowNorm));
1729
1730 // Dist to High Cut
1731 //
1732 safZHigh = std::fabs((p-vZ).dot(fHighNorm));
1733 safe = std::min(safZLow,safZHigh);
1734
1735 if ( safRMin < safe ) { safe = safRMin; }
1736 if ( safRMax< safe ) { safe = safRMax; }
1737
1738 // Check if phi divided, Calc distances closest phi plane
1739 //
1740 if ( !fPhiFullCutTube )
1741 {
1742 if ( p.y()*cosCPhi-p.x()*sinCPhi <= 0 )
1743 {
1744 safePhi = -(p.x()*sinSPhi - p.y()*cosSPhi) ;
1745 }
1746 else
1747 {
1748 safePhi = (p.x()*sinEPhi - p.y()*cosEPhi) ;
1749 }
1750 if (safePhi < safe) { safe = safePhi ; }
1751 }
1752 if ( safe < 0 ) { safe = 0; }
1753
1754 return safe ;
1755}

◆ DistanceToOut() [2/2]

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

Implements G4VSolid.

Definition at line 1281 of file G4CutTubs.cc.

1286{
1287 ESide side=kNull , sider=kNull, sidephi=kNull ;
1288 G4double snxt=kInfinity, srd=kInfinity,sz=kInfinity, sphi=kInfinity ;
1289 G4double deltaR, t1, t2, t3, b, c, d2, roMin2 ;
1290 G4double distZLow,distZHigh,calfH,calfL;
1291 G4ThreeVector vZ=G4ThreeVector(0,0,fDz);
1292 static const G4double halfCarTolerance = kCarTolerance*0.5;
1293 static const G4double halfAngTolerance = kAngTolerance*0.5;
1294
1295 // Vars for phi intersection:
1296 //
1297 G4double pDistS, compS, pDistE, compE, sphi2, xi, yi, vphi, roi2 ;
1298
1299 // Z plane intersection
1300 // Distances to ZCut(Low/High)
1301
1302 // dist to Low Cut
1303 //
1304 distZLow =(p+vZ).dot(fLowNorm);
1305
1306 // dist to High Cut
1307 //
1308 distZHigh = (p-vZ).dot(fHighNorm);
1309
1310 calfH = v.dot(fHighNorm);
1311 calfL = v.dot(fLowNorm);
1312
1313 if (calfH > 0 )
1314 {
1315 if ( distZHigh < halfCarTolerance )
1316 {
1317 snxt = -distZHigh/calfH ;
1318 side = kPZ ;
1319 }
1320 else
1321 {
1322 if (calcNorm)
1323 {
1324 *n = G4ThreeVector(0,0,1) ;
1325 *validNorm = true ;
1326 }
1327 return snxt = 0 ;
1328 }
1329 }
1330 if ( calfL>0)
1331 {
1332
1333 if ( distZLow < halfCarTolerance )
1334 {
1335 sz = -distZLow/calfL ;
1336 if(sz<snxt){
1337 snxt=sz;
1338 side = kMZ ;
1339 }
1340
1341 }
1342 else
1343 {
1344 if (calcNorm)
1345 {
1346 *n = G4ThreeVector(0,0,-1) ;
1347 *validNorm = true ;
1348 }
1349 return snxt = 0.0 ;
1350 }
1351 }
1352 if((calfH<=0)&&(calfL<=0))
1353 {
1354 snxt = kInfinity ; // Travel perpendicular to z axis
1355 side = kNull;
1356 }
1357 // Radial Intersections
1358 //
1359 // Find intersection with cylinders at rmax/rmin
1360 // Intersection point (xi,yi,zi) on line x=p.x+t*v.x etc.
1361 //
1362 // Intersects with x^2+y^2=R^2
1363 //
1364 // Hence (v.x^2+v.y^2)t^2+ 2t(p.x*v.x+p.y*v.y)+p.x^2+p.y^2-R^2=0
1365 //
1366 // t1 t2 t3
1367
1368 t1 = 1.0 - v.z()*v.z() ; // since v normalised
1369 t2 = p.x()*v.x() + p.y()*v.y() ;
1370 t3 = p.x()*p.x() + p.y()*p.y() ;
1371
1372 if ( snxt > 10*(fDz+fRMax) ) { roi2 = 2*fRMax*fRMax; }
1373 else { roi2 = snxt*snxt*t1 + 2*snxt*t2 + t3; } // radius^2 on +-fDz
1374
1375 if ( t1 > 0 ) // Check not parallel
1376 {
1377 // Calculate srd, r exit distance
1378
1379 if ( (t2 >= 0.0) && (roi2 > fRMax*(fRMax + kRadTolerance)) )
1380 {
1381 // Delta r not negative => leaving via rmax
1382
1383 deltaR = t3 - fRMax*fRMax ;
1384
1385 // NOTE: Should use rho-fRMax<-kRadTolerance*0.5
1386 // - avoid sqrt for efficiency
1387
1388 if ( deltaR < -kRadTolerance*fRMax )
1389 {
1390 b = t2/t1 ;
1391 c = deltaR/t1 ;
1392 d2 = b*b-c;
1393 if( d2 >= 0 ) { srd = c/( -b - std::sqrt(d2)); }
1394 else { srd = 0.; }
1395 sider = kRMax ;
1396 }
1397 else
1398 {
1399 // On tolerant boundary & heading outwards (or perpendicular to)
1400 // outer radial surface -> leaving immediately
1401
1402 if ( calcNorm )
1403 {
1404 *n = G4ThreeVector(p.x()/fRMax,p.y()/fRMax,0) ;
1405 *validNorm = true ;
1406 }
1407 return snxt = 0 ; // Leaving by rmax immediately
1408 }
1409 }
1410 else if ( t2 < 0. ) // i.e. t2 < 0; Possible rmin intersection
1411 {
1412 roMin2 = t3 - t2*t2/t1 ; // min ro2 of the plane of movement
1413
1414 if ( fRMin && (roMin2 < fRMin*(fRMin - kRadTolerance)) )
1415 {
1416 deltaR = t3 - fRMin*fRMin ;
1417 b = t2/t1 ;
1418 c = deltaR/t1 ;
1419 d2 = b*b - c ;
1420
1421 if ( d2 >= 0 ) // Leaving via rmin
1422 {
1423 // NOTE: SHould use rho-rmin>kRadTolerance*0.5
1424 // - avoid sqrt for efficiency
1425
1426 if (deltaR > kRadTolerance*fRMin)
1427 {
1428 srd = c/(-b+std::sqrt(d2));
1429 sider = kRMin ;
1430 }
1431 else
1432 {
1433 if ( calcNorm ) { *validNorm = false; } // Concave side
1434 return snxt = 0.0;
1435 }
1436 }
1437 else // No rmin intersect -> must be rmax intersect
1438 {
1439 deltaR = t3 - fRMax*fRMax ;
1440 c = deltaR/t1 ;
1441 d2 = b*b-c;
1442 if( d2 >=0. )
1443 {
1444 srd = -b + std::sqrt(d2) ;
1445 sider = kRMax ;
1446 }
1447 else // Case: On the border+t2<kRadTolerance
1448 // (v is perpendicular to the surface)
1449 {
1450 if (calcNorm)
1451 {
1452 *n = G4ThreeVector(p.x()/fRMax,p.y()/fRMax,0) ;
1453 *validNorm = true ;
1454 }
1455 return snxt = 0.0;
1456 }
1457 }
1458 }
1459 else if ( roi2 > fRMax*(fRMax + kRadTolerance) )
1460 // No rmin intersect -> must be rmax intersect
1461 {
1462 deltaR = t3 - fRMax*fRMax ;
1463 b = t2/t1 ;
1464 c = deltaR/t1;
1465 d2 = b*b-c;
1466 if( d2 >= 0 )
1467 {
1468 srd = -b + std::sqrt(d2) ;
1469 sider = kRMax ;
1470 }
1471 else // Case: On the border+t2<kRadTolerance
1472 // (v is perpendicular to the surface)
1473 {
1474 if (calcNorm)
1475 {
1476 *n = G4ThreeVector(p.x()/fRMax,p.y()/fRMax,0) ;
1477 *validNorm = true ;
1478 }
1479 return snxt = 0.0;
1480 }
1481 }
1482 }
1483 // Phi Intersection
1484
1485 if ( !fPhiFullCutTube )
1486 {
1487 // add angle calculation with correction
1488 // of the difference in domain of atan2 and Sphi
1489 //
1490 vphi = std::atan2(v.y(),v.x()) ;
1491
1492 if ( vphi < fSPhi - halfAngTolerance ) { vphi += twopi; }
1493 else if ( vphi > fSPhi + fDPhi + halfAngTolerance ) { vphi -= twopi; }
1494
1495
1496 if ( p.x() || p.y() ) // Check if on z axis (rho not needed later)
1497 {
1498 // pDist -ve when inside
1499
1500 pDistS = p.x()*sinSPhi - p.y()*cosSPhi ;
1501 pDistE = -p.x()*sinEPhi + p.y()*cosEPhi ;
1502
1503 // Comp -ve when in direction of outwards normal
1504
1505 compS = -sinSPhi*v.x() + cosSPhi*v.y() ;
1506 compE = sinEPhi*v.x() - cosEPhi*v.y() ;
1507
1508 sidephi = kNull;
1509
1510 if( ( (fDPhi <= pi) && ( (pDistS <= halfCarTolerance)
1511 && (pDistE <= halfCarTolerance) ) )
1512 || ( (fDPhi > pi) && !((pDistS > halfCarTolerance)
1513 && (pDistE > halfCarTolerance) ) ) )
1514 {
1515 // Inside both phi *full* planes
1516
1517 if ( compS < 0 )
1518 {
1519 sphi = pDistS/compS ;
1520
1521 if (sphi >= -halfCarTolerance)
1522 {
1523 xi = p.x() + sphi*v.x() ;
1524 yi = p.y() + sphi*v.y() ;
1525
1526 // Check intersecting with correct half-plane
1527 // (if not -> no intersect)
1528 //
1529 if( (std::fabs(xi)<=kCarTolerance)
1530 && (std::fabs(yi)<=kCarTolerance) )
1531 {
1532 sidephi = kSPhi;
1533 if (((fSPhi-halfAngTolerance)<=vphi)
1534 &&((fSPhi+fDPhi+halfAngTolerance)>=vphi))
1535 {
1536 sphi = kInfinity;
1537 }
1538 }
1539 else if ( yi*cosCPhi-xi*sinCPhi >=0 )
1540 {
1541 sphi = kInfinity ;
1542 }
1543 else
1544 {
1545 sidephi = kSPhi ;
1546 if ( pDistS > -halfCarTolerance )
1547 {
1548 sphi = 0.0 ; // Leave by sphi immediately
1549 }
1550 }
1551 }
1552 else
1553 {
1554 sphi = kInfinity ;
1555 }
1556 }
1557 else
1558 {
1559 sphi = kInfinity ;
1560 }
1561
1562 if ( compE < 0 )
1563 {
1564 sphi2 = pDistE/compE ;
1565
1566 // Only check further if < starting phi intersection
1567 //
1568 if ( (sphi2 > -halfCarTolerance) && (sphi2 < sphi) )
1569 {
1570 xi = p.x() + sphi2*v.x() ;
1571 yi = p.y() + sphi2*v.y() ;
1572
1573 if ((std::fabs(xi)<=kCarTolerance)&&(std::fabs(yi)<=kCarTolerance))
1574 {
1575 // Leaving via ending phi
1576 //
1577 if( !((fSPhi-halfAngTolerance <= vphi)
1578 &&(fSPhi+fDPhi+halfAngTolerance >= vphi)) )
1579 {
1580 sidephi = kEPhi ;
1581 if ( pDistE <= -halfCarTolerance ) { sphi = sphi2 ; }
1582 else { sphi = 0.0 ; }
1583 }
1584 }
1585 else // Check intersecting with correct half-plane
1586
1587 if ( (yi*cosCPhi-xi*sinCPhi) >= 0)
1588 {
1589 // Leaving via ending phi
1590 //
1591 sidephi = kEPhi ;
1592 if ( pDistE <= -halfCarTolerance ) { sphi = sphi2 ; }
1593 else { sphi = 0.0 ; }
1594 }
1595 }
1596 }
1597 }
1598 else
1599 {
1600 sphi = kInfinity ;
1601 }
1602 }
1603 else
1604 {
1605 // On z axis + travel not || to z axis -> if phi of vector direction
1606 // within phi of shape, Step limited by rmax, else Step =0
1607
1608 if ( (fSPhi - halfAngTolerance <= vphi)
1609 && (vphi <= fSPhi + fDPhi + halfAngTolerance ) )
1610 {
1611 sphi = kInfinity ;
1612 }
1613 else
1614 {
1615 sidephi = kSPhi ; // arbitrary
1616 sphi = 0.0 ;
1617 }
1618 }
1619 if (sphi < snxt) // Order intersecttions
1620 {
1621 snxt = sphi ;
1622 side = sidephi ;
1623 }
1624 }
1625 if (srd < snxt) // Order intersections
1626 {
1627 snxt = srd ;
1628 side = sider ;
1629 }
1630 }
1631 if (calcNorm)
1632 {
1633 switch(side)
1634 {
1635 case kRMax:
1636 // Note: returned vector not normalised
1637 // (divide by fRMax for unit vector)
1638 //
1639 xi = p.x() + snxt*v.x() ;
1640 yi = p.y() + snxt*v.y() ;
1641 *n = G4ThreeVector(xi/fRMax,yi/fRMax,0) ;
1642 *validNorm = true ;
1643 break ;
1644
1645 case kRMin:
1646 *validNorm = false ; // Rmin is inconvex
1647 break ;
1648
1649 case kSPhi:
1650 if ( fDPhi <= pi )
1651 {
1652 *n = G4ThreeVector(sinSPhi,-cosSPhi,0) ;
1653 *validNorm = true ;
1654 }
1655 else
1656 {
1657 *validNorm = false ;
1658 }
1659 break ;
1660
1661 case kEPhi:
1662 if (fDPhi <= pi)
1663 {
1664 *n = G4ThreeVector(-sinEPhi,cosEPhi,0) ;
1665 *validNorm = true ;
1666 }
1667 else
1668 {
1669 *validNorm = false ;
1670 }
1671 break ;
1672
1673 case kPZ:
1674 *n = fHighNorm ;
1675 *validNorm = true ;
1676 break ;
1677
1678 case kMZ:
1679 *n = fLowNorm ;
1680 *validNorm = true ;
1681 break ;
1682
1683 default:
1684 G4cout << G4endl ;
1685 DumpInfo();
1686 std::ostringstream message;
1687 G4int oldprc = message.precision(16);
1688 message << "Undefined side for valid surface normal to solid."
1689 << G4endl
1690 << "Position:" << G4endl << G4endl
1691 << "p.x() = " << p.x()/mm << " mm" << G4endl
1692 << "p.y() = " << p.y()/mm << " mm" << G4endl
1693 << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl
1694 << "Direction:" << G4endl << G4endl
1695 << "v.x() = " << v.x() << G4endl
1696 << "v.y() = " << v.y() << G4endl
1697 << "v.z() = " << v.z() << G4endl << G4endl
1698 << "Proposed distance :" << G4endl << G4endl
1699 << "snxt = " << snxt/mm << " mm" << G4endl ;
1700 message.precision(oldprc) ;
1701 G4Exception("G4CutTubs::DistanceToOut(p,v,..)", "GeomSolids1002",
1702 JustWarning, message);
1703 break ;
1704 }
1705 }
1706 if ( snxt<halfCarTolerance ) { snxt=0 ; }
1707 return snxt ;
1708}
ESide
ESide
Definition: G4Cons.cc:68
G4cout
G4DLLIMPORT std::ostream G4cout
G4Tubs::kEPhi
@ kEPhi
Definition: G4Tubs.hh:199
G4Tubs::kRMax
@ kRMax
Definition: G4Tubs.hh:199
G4Tubs::kPZ
@ kPZ
Definition: G4Tubs.hh:199
G4Tubs::kMZ
@ kMZ
Definition: G4Tubs.hh:199
G4Tubs::kRMin
@ kRMin
Definition: G4Tubs.hh:199
G4Tubs::kSPhi
@ kSPhi
Definition: G4Tubs.hh:199
G4Tubs::kNull
@ kNull
Definition: G4Tubs.hh:199

◆ GetCubicVolume()

G4double G4CutTubs::GetCubicVolume ( )
inlinevirtual

Reimplemented from G4VSolid.

◆ GetCutZ()

G4double G4CutTubs::GetCutZ ( const G4ThreeVector &  p) const
protected

Definition at line 2058 of file G4CutTubs.cc.

2059{
2060 G4double newz = p.z(); // p.z() should be either +fDz or -fDz
2061 if (p.z()<0)
2062 {
2063 if(fLowNorm.z()!=0.)
2064 {
2065 newz = -fDz-(p.x()*fLowNorm.x()+p.y()*fLowNorm.y())/fLowNorm.z();
2066 }
2067 }
2068 else
2069 {
2070 if(fHighNorm.z()!=0.)
2071 {
2072 newz = fDz-(p.x()*fHighNorm.x()+p.y()*fHighNorm.y())/fHighNorm.z();
2073 }
2074 }
2075 return newz;
2076}

Referenced by CreatePolyhedron(), CreateRotatedVertices(), DistanceToIn(), GetMaxMinZ(), GetPointOnSurface(), and IsCrossingCutPlanes().

◆ GetEntityType()

G4GeometryType G4CutTubs::GetEntityType ( ) const
virtual

Implements G4VSolid.

Definition at line 1853 of file G4CutTubs.cc.

1854{
1855 return G4String("G4CutTubs");
1856}

◆ GetHighNorm()

G4ThreeVector G4CutTubs::GetHighNorm ( ) const
inline

Referenced by G4GDMLWriteSolids::CutTubeWrite().

◆ GetLowNorm()

G4ThreeVector G4CutTubs::GetLowNorm ( ) const
inline

Referenced by G4GDMLWriteSolids::CutTubeWrite().

◆ GetMaxMinZ()

void G4CutTubs::GetMaxMinZ ( G4double &  zmin,
G4double &  zmax 
) const
protected

Definition at line 2082 of file G4CutTubs.cc.

2084{
2085 G4double phiLow = std::atan2(fLowNorm.y(),fLowNorm.x());
2086 G4double phiHigh= std::atan2(fHighNorm.y(),fHighNorm.x());
2087
2088 G4double xc=0, yc=0,z1;
2089 G4double z[8];
2090 G4bool in_range_low = false;
2091 G4bool in_range_hi = false;
2092
2093 G4int i;
2094 for (i=0; i<2; i++)
2095 {
2096 if (phiLow<0) { phiLow+=twopi; }
2097 G4double ddp = phiLow-fSPhi;
2098 if (ddp<0) { ddp += twopi; }
2099 if (ddp <= fDPhi)
2100 {
2101 xc = fRMin*std::cos(phiLow);
2102 yc = fRMin*std::sin(phiLow);
2103 z1 = GetCutZ(G4ThreeVector(xc, yc, -fDz));
2104 xc = fRMax*std::cos(phiLow);
2105 yc = fRMax*std::sin(phiLow);
2106 z1 = std::min(z1, GetCutZ(G4ThreeVector(xc, yc, -fDz)));
2107 if (in_range_low) { zmin = std::min(zmin, z1); }
2108 else { zmin = z1; }
2109 in_range_low = true;
2110 }
2111 phiLow += pi;
2112 if (phiLow>twopi) { phiLow-=twopi; }
2113 }
2114 for (i=0; i<2; i++)
2115 {
2116 if (phiHigh<0) { phiHigh+=twopi; }
2117 G4double ddp = phiHigh-fSPhi;
2118 if (ddp<0) { ddp += twopi; }
2119 if (ddp <= fDPhi)
2120 {
2121 xc = fRMin*std::cos(phiHigh);
2122 yc = fRMin*std::sin(phiHigh);
2123 z1 = GetCutZ(G4ThreeVector(xc, yc, fDz));
2124 xc = fRMax*std::cos(phiHigh);
2125 yc = fRMax*std::sin(phiHigh);
2126 z1 = std::min(z1, GetCutZ(G4ThreeVector(xc, yc, fDz)));
2127 if (in_range_hi) { zmax = std::min(zmax, z1); }
2128 else { zmax = z1; }
2129 in_range_hi = true;
2130 }
2131 phiHigh += pi;
2132 if (phiLow>twopi) { phiHigh-=twopi; }
2133 }
2134
2135 xc = fRMin*std::cos(fSPhi);
2136 yc = fRMin*std::sin(fSPhi);
2137 z[0] = GetCutZ(G4ThreeVector(xc, yc, -fDz));
2138 z[4] = GetCutZ(G4ThreeVector(xc, yc, fDz));
2139
2140 xc = fRMin*std::cos(fSPhi+fDPhi);
2141 yc = fRMin*std::sin(fSPhi+fDPhi);
2142 z[1] = GetCutZ(G4ThreeVector(xc, yc, -fDz));
2143 z[5] = GetCutZ(G4ThreeVector(xc, yc, fDz));
2144
2145 xc = fRMax*std::cos(fSPhi);
2146 yc = fRMax*std::sin(fSPhi);
2147 z[2] = GetCutZ(G4ThreeVector(xc, yc, -fDz));
2148 z[6] = GetCutZ(G4ThreeVector(xc, yc, fDz));
2149
2150 xc = fRMax*std::cos(fSPhi+fDPhi);
2151 yc = fRMax*std::sin(fSPhi+fDPhi);
2152 z[3] = GetCutZ(G4ThreeVector(xc, yc, -fDz));
2153 z[7] = GetCutZ(G4ThreeVector(xc, yc, fDz));
2154
2155 // Find min/max
2156
2157 z1=z[0];
2158 for (i = 1; i < 4; i++)
2159 {
2160 if(z[i] < z[i-1])z1=z[i];
2161 }
2162
2163 if (in_range_low)
2164 {
2165 zmin = std::min(zmin, z1);
2166 }
2167 else
2168 {
2169 zmin = z1;
2170 }
2171 z1=z[4];
2172 for (i = 1; i < 4; i++)
2173 {
2174 if(z[4+i] > z[4+i-1]) { z1=z[4+i]; }
2175 }
2176
2177 if (in_range_hi) { zmax = std::max(zmax, z1); }
2178 else { zmax = z1; }
2179}
G4INCL::Math::pi
const G4double pi
Definition: G4INCLGlobals.hh:70

Referenced by CalculateExtent().

◆ GetPointOnSurface()

G4ThreeVector G4CutTubs::GetPointOnSurface ( ) const
virtual

Reimplemented from G4VSolid.

Definition at line 1896 of file G4CutTubs.cc.

1897{
1898 G4double xRand, yRand, zRand, phi, cosphi, sinphi, chose,
1899 aOne, aTwo, aThr, aFou;
1900 G4double rRand;
1901
1902 aOne = 2.*fDz*fDPhi*fRMax;
1903 aTwo = 2.*fDz*fDPhi*fRMin;
1904 aThr = 0.5*fDPhi*(fRMax*fRMax-fRMin*fRMin);
1905 aFou = 2.*fDz*(fRMax-fRMin);
1906
1907 phi = RandFlat::shoot(fSPhi, fSPhi+fDPhi);
1908 cosphi = std::cos(phi);
1909 sinphi = std::sin(phi);
1910
1911 rRand = GetRadiusInRing(fRMin,fRMax);
1912
1913 if( (fSPhi == 0) && (fDPhi == twopi) ) { aFou = 0; }
1914
1915 chose = RandFlat::shoot(0.,aOne+aTwo+2.*aThr+2.*aFou);
1916
1917 if( (chose >=0) && (chose < aOne) )
1918 {
1919 xRand = fRMax*cosphi;
1920 yRand = fRMax*sinphi;
1921 zRand = RandFlat::shoot(GetCutZ(G4ThreeVector(xRand,yRand,-fDz)),
1922 GetCutZ(G4ThreeVector(xRand,yRand,fDz)));
1923 return G4ThreeVector (xRand, yRand, zRand);
1924 }
1925 else if( (chose >= aOne) && (chose < aOne + aTwo) )
1926 {
1927 xRand = fRMin*cosphi;
1928 yRand = fRMin*sinphi;
1929 zRand = RandFlat::shoot(GetCutZ(G4ThreeVector(xRand,yRand,-fDz)),
1930 GetCutZ(G4ThreeVector(xRand,yRand,fDz)));
1931 return G4ThreeVector (xRand, yRand, zRand);
1932 }
1933 else if( (chose >= aOne + aTwo) && (chose < aOne + aTwo + aThr) )
1934 {
1935 xRand = rRand*cosphi;
1936 yRand = rRand*sinphi;
1937 zRand = GetCutZ(G4ThreeVector(xRand,yRand,fDz));
1938 return G4ThreeVector (xRand, yRand, zRand);
1939 }
1940 else if( (chose >= aOne + aTwo + aThr) && (chose < aOne + aTwo + 2.*aThr) )
1941 {
1942 xRand = rRand*cosphi;
1943 yRand = rRand*sinphi;
1944 zRand = GetCutZ(G4ThreeVector(xRand,yRand,-fDz));
1945 return G4ThreeVector (xRand, yRand, zRand);
1946 }
1947 else if( (chose >= aOne + aTwo + 2.*aThr)
1948 && (chose < aOne + aTwo + 2.*aThr + aFou) )
1949 {
1950 xRand = rRand*std::cos(fSPhi);
1951 yRand = rRand*std::sin(fSPhi);
1952 zRand = RandFlat::shoot(GetCutZ(G4ThreeVector(xRand,yRand,-fDz)),
1953 GetCutZ(G4ThreeVector(xRand,yRand,fDz)));
1954 return G4ThreeVector (xRand, yRand, zRand);
1955 }
1956 else
1957 {
1958 xRand = rRand*std::cos(fSPhi+fDPhi);
1959 yRand = rRand*std::sin(fSPhi+fDPhi);
1960 zRand = RandFlat::shoot(GetCutZ(G4ThreeVector(xRand,yRand,-fDz)),
1961 GetCutZ(G4ThreeVector(xRand,yRand,fDz)));
1962 return G4ThreeVector (xRand, yRand, zRand);
1963 }
1964}
CLHEP::RandFlat::shoot
static double shoot()
Definition: RandFlat.cc:59
G4CSGSolid::GetRadiusInRing
G4double GetRadiusInRing(G4double rmin, G4double rmax) const
Definition: G4CSGSolid.cc:101

◆ GetSurfaceArea()

G4double G4CutTubs::GetSurfaceArea ( )
inlinevirtual

Reimplemented from G4VSolid.

◆ Inside()

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

Implements G4VSolid.

Definition at line 407 of file G4CutTubs.cc.

408{
409 G4double zinLow,zinHigh,r2,pPhi=0.;
410 G4double tolRMin,tolRMax;
411 G4ThreeVector vZ=G4ThreeVector(0,0,fDz);
412 EInside in = kInside;
413 static const G4double halfCarTolerance=kCarTolerance*0.5;
414 static const G4double halfRadTolerance=kRadTolerance*0.5;
415 static const G4double halfAngTolerance=kAngTolerance*0.5;
416
417 // Check if point is contained in the cut plane in -/+ Z
418
419 // Check the lower cut plane
420 //
421 zinLow =(p+vZ).dot(fLowNorm);
422 if (zinLow>halfCarTolerance) { return kOutside; }
423
424 // Check the higher cut plane
425 //
426 zinHigh = (p-vZ).dot(fHighNorm);
427 if (zinHigh>halfCarTolerance) { return kOutside; }
428
429 // Check radius
430 //
431 r2 = p.x()*p.x() + p.y()*p.y() ;
432
433 // First check 'generous' boundaries R+tolerance
434 //
435 tolRMin = fRMin - halfRadTolerance ;
436 tolRMax = fRMax + halfRadTolerance ;
437 if ( tolRMin < 0 ) { tolRMin = 0; }
438
439 if ( ((r2 < tolRMin*tolRMin) || (r2 > tolRMax*tolRMax))
440 && (r2 >=halfRadTolerance*halfRadTolerance) ) { return kOutside; }
441
442 // Check Phi
443 //
444 if(!fPhiFullCutTube)
445 {
446 // Try outer tolerant phi boundaries only
447
448 if ( (tolRMin==0) && (std::fabs(p.x())<=halfCarTolerance)
449 && (std::fabs(p.y())<=halfCarTolerance) )
450 {
451 return kSurface;
452 }
453
454 pPhi = std::atan2(p.y(),p.x()) ;
455
456 if ( pPhi < -halfAngTolerance) { pPhi += twopi; } // 0<=pPhi<2pi
457 if ( fSPhi >= 0 )
458 {
459 if ( (std::fabs(pPhi) < halfAngTolerance)
460 && (std::fabs(fSPhi + fDPhi - twopi) < halfAngTolerance) )
461 {
462 pPhi += twopi ; // 0 <= pPhi < 2pi
463 }
464 if ( (pPhi <= fSPhi - halfAngTolerance)
465 || (pPhi >= fSPhi + fDPhi + halfAngTolerance) )
466 {
467 in = kOutside ;
468 }
469 else if ( (pPhi <= fSPhi + halfAngTolerance)
470 || (pPhi >= fSPhi + fDPhi - halfAngTolerance) )
471 {
472 in=kSurface;
473 }
474 }
475 else // fSPhi < 0
476 {
477 if ( (pPhi <= fSPhi + twopi - halfAngTolerance)
478 && (pPhi >= fSPhi + fDPhi + halfAngTolerance) )
479 {
480 in = kOutside;
481 }
482 else
483 {
484 in = kSurface ;
485 }
486 }
487 }
488
489 // Check on the Surface for Z
490 //
491 if ((zinLow>=-halfCarTolerance)
492 || (zinHigh>=-halfCarTolerance))
493 {
494 in=kSurface;
495 }
496
497 // Check on the Surface for R
498 //
499 if (fRMin) { tolRMin = fRMin + halfRadTolerance ; }
500 else { tolRMin = 0 ; }
501 tolRMax = fRMax - halfRadTolerance ;
502 if ( ((r2 <= tolRMin*tolRMin) || (r2 >= tolRMax*tolRMax))&&
503 (r2 >=halfRadTolerance*halfRadTolerance) )
504 {
505 return kSurface;
506 }
507
508 return in;
509}
EInside
EInside
Definition: geomdefs.hh:58
kInside
@ kInside
Definition: geomdefs.hh:58
kSurface
@ kSurface
Definition: geomdefs.hh:58

Referenced by CalculateExtent().

◆ IsCrossingCutPlanes()

G4bool G4CutTubs::IsCrossingCutPlanes ( ) const
protected

Definition at line 2035 of file G4CutTubs.cc.

2036{
2037 G4double zXLow1,zXLow2,zYLow1,zYLow2;
2038 G4double zXHigh1,zXHigh2,zYHigh1,zYHigh2;
2039
2040 zXLow1 = GetCutZ(G4ThreeVector(-fRMax, 0,-fDz));
2041 zXLow2 = GetCutZ(G4ThreeVector( fRMax, 0,-fDz));
2042 zYLow1 = GetCutZ(G4ThreeVector( 0,-fRMax,-fDz));
2043 zYLow2 = GetCutZ(G4ThreeVector( 0, fRMax,-fDz));
2044 zXHigh1 = GetCutZ(G4ThreeVector(-fRMax, 0, fDz));
2045 zXHigh2 = GetCutZ(G4ThreeVector( fRMax, 0, fDz));
2046 zYHigh1 = GetCutZ(G4ThreeVector( 0,-fRMax, fDz));
2047 zYHigh2 = GetCutZ(G4ThreeVector( 0, fRMax, fDz));
2048 if ( (zXLow1>zXHigh1) ||(zXLow2>zXHigh2)
2049 || (zYLow1>zYHigh1) ||(zYLow2>zYHigh2)) { return true; }
2050
2051 return false;
2052}

Referenced by G4CutTubs().

◆ operator=()

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

Definition at line 172 of file G4CutTubs.cc.

173{
174 // Check assignment to self
175 //
176 if (this == &rhs) { return *this; }
177
178 // Copy base class data
179 //
180 G4Tubs::operator=(rhs);
181
182 // Copy data
183 //
184 fLowNorm = rhs.fLowNorm; fHighNorm = rhs.fHighNorm;
185 fPhiFullCutTube = rhs.fPhiFullCutTube;
186
187 return *this;
188}
G4Tubs::operator=
G4Tubs & operator=(const G4Tubs &rhs)
Definition: G4Tubs.cc:160

◆ StreamInfo()

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

Reimplemented from G4CSGSolid.

Definition at line 1871 of file G4CutTubs.cc.

1872{
1873 G4int oldprc = os.precision(16);
1874 os << "-----------------------------------------------------------\n"
1875 << " *** Dump for solid - " << GetName() << " ***\n"
1876 << " ===================================================\n"
1877 << " Solid type: G4CutTubs\n"
1878 << " Parameters: \n"
1879 << " inner radius : " << fRMin/mm << " mm \n"
1880 << " outer radius : " << fRMax/mm << " mm \n"
1881 << " half length Z: " << fDz/mm << " mm \n"
1882 << " starting phi : " << fSPhi/degree << " degrees \n"
1883 << " delta phi : " << fDPhi/degree << " degrees \n"
1884 << " low Norm : " << fLowNorm << " \n"
1885 << " high Norm : " <<fHighNorm << " \n"
1886 << "-----------------------------------------------------------\n";
1887 os.precision(oldprc);
1888
1889 return os;
1890}

◆ SurfaceNormal()

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

Implements G4VSolid.

Definition at line 517 of file G4CutTubs.cc.

518{
519 G4int noSurfaces = 0;
520 G4double rho, pPhi;
521 G4double distZLow,distZHigh, distRMin, distRMax;
522 G4double distSPhi = kInfinity, distEPhi = kInfinity;
523 G4ThreeVector vZ=G4ThreeVector(0,0,fDz);
524
525 static const G4double halfCarTolerance = 0.5*kCarTolerance;
526 static const G4double halfAngTolerance = 0.5*kAngTolerance;
527
528 G4ThreeVector norm, sumnorm(0.,0.,0.);
529 G4ThreeVector nZ = G4ThreeVector(0, 0, 1.0);
530 G4ThreeVector nR, nPs, nPe;
531
532 rho = std::sqrt(p.x()*p.x() + p.y()*p.y());
533
534 distRMin = std::fabs(rho - fRMin);
535 distRMax = std::fabs(rho - fRMax);
536
537 // dist to Low Cut
538 //
539 distZLow =std::fabs((p+vZ).dot(fLowNorm));
540
541 // dist to High Cut
542 //
543 distZHigh = std::fabs((p-vZ).dot(fHighNorm));
544
545 if (!fPhiFullCutTube) // Protected against (0,0,z)
546 {
547 if ( rho > halfCarTolerance )
548 {
549 pPhi = std::atan2(p.y(),p.x());
550
551 if(pPhi < fSPhi- halfCarTolerance) { pPhi += twopi; }
552 else if(pPhi > fSPhi+fDPhi+ halfCarTolerance) { pPhi -= twopi; }
553
554 distSPhi = std::fabs(pPhi - fSPhi);
555 distEPhi = std::fabs(pPhi - fSPhi - fDPhi);
556 }
557 else if( !fRMin )
558 {
559 distSPhi = 0.;
560 distEPhi = 0.;
561 }
562 nPs = G4ThreeVector(std::sin(fSPhi),-std::cos(fSPhi),0);
563 nPe = G4ThreeVector(-std::sin(fSPhi+fDPhi),std::cos(fSPhi+fDPhi),0);
564 }
565 if ( rho > halfCarTolerance ) { nR = G4ThreeVector(p.x()/rho,p.y()/rho,0); }
566
567 if( distRMax <= halfCarTolerance )
568 {
569 noSurfaces ++;
570 sumnorm += nR;
571 }
572 if( fRMin && (distRMin <= halfCarTolerance) )
573 {
574 noSurfaces ++;
575 sumnorm -= nR;
576 }
577 if( fDPhi < twopi )
578 {
579 if (distSPhi <= halfAngTolerance)
580 {
581 noSurfaces ++;
582 sumnorm += nPs;
583 }
584 if (distEPhi <= halfAngTolerance)
585 {
586 noSurfaces ++;
587 sumnorm += nPe;
588 }
589 }
590 if (distZLow <= halfCarTolerance)
591 {
592 noSurfaces ++;
593 sumnorm += fLowNorm;
594 }
595 if (distZHigh <= halfCarTolerance)
596 {
597 noSurfaces ++;
598 sumnorm += fHighNorm;
599 }
600 if ( noSurfaces == 0 )
601 {
602#ifdef G4CSGDEBUG
603 G4Exception("G4CutTubs::SurfaceNormal(p)", "GeomSolids1002",
604 JustWarning, "Point p is not on surface !?" );
605 G4int oldprc = G4cout.precision(20);
606 G4cout<< "G4CutTubs::SN ( "<<p.x()<<", "<<p.y()<<", "<<p.z()<<" ); "
607 << G4endl << G4endl;
608 G4cout.precision(oldprc) ;
609#endif
610 norm = ApproxSurfaceNormal(p);
611 }
612 else if ( noSurfaces == 1 ) { norm = sumnorm; }
613 else { norm = sumnorm.unit(); }
614
615 return norm;
616}
G4CutTubs::ApproxSurfaceNormal
G4ThreeVector ApproxSurfaceNormal(const G4ThreeVector &p) const
Definition: G4CutTubs.cc:623
The documentation for this class was generated from the following files: