Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Public Member Functions | Protected Member Functions | List of all members
G4CutTubs Class Reference

#include <G4CutTubs.hh>

Inheritance diagram for G4CutTubs:
Inheritance graph
[legend]
Collaboration diagram for G4CutTubs:
Collaboration graph
[legend]

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 ()
 
void Extent (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
 
 G4CutTubs (__void__ &)
 
 G4CutTubs (const G4CutTubs &rhs)
 
G4CutTubs & operator= (const G4CutTubs &rhs)
 
- Public Member Functions inherited from G4OTubs
 G4OTubs (const G4String &pName, G4double pRMin, G4double pRMax, G4double pDz, G4double pSPhi, G4double pDPhi)
 
virtual ~G4OTubs ()
 
G4double GetInnerRadius () const
 
G4double GetOuterRadius () const
 
G4double GetZHalfLength () const
 
G4double GetStartPhiAngle () const
 
G4double GetDeltaPhiAngle () const
 
G4double GetSinStartPhi () const
 
G4double GetCosStartPhi () const
 
G4double GetSinEndPhi () const
 
G4double GetCosEndPhi () 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 Extent (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
 
 G4OTubs (__void__ &)
 
 G4OTubs (const G4OTubs &rhs)
 
G4OTubs & operator= (const G4OTubs &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 G4Polyhedron * GetPolyhedron () const
 
 G4CSGSolid (__void__ &)
 
 G4CSGSolid (const G4CSGSolid &rhs)
 
G4CSGSolid & operator= (const G4CSGSolid &rhs)
 
- Public Member Functions inherited from G4VSolid
 G4VSolid (const G4String &name)
 
virtual ~G4VSolid ()
 
G4bool operator== (const G4VSolid &s) const
 
G4String GetName () const
 
void SetName (const G4String &name)
 
G4double GetTolerance () const
 
virtual void ComputeDimensions (G4VPVParameterisation *p, const G4int n, const G4VPhysicalVolume *pRep)
 
void DumpInfo () const
 
virtual G4VisExtent GetExtent () const
 
virtual const G4VSolid * GetConstituentSolid (G4int no) const
 
virtual G4VSolid * GetConstituentSolid (G4int no)
 
virtual const G4DisplacedSolid * GetDisplacedSolidPtr () const
 
virtual G4DisplacedSolid * GetDisplacedSolidPtr ()
 
 G4VSolid (__void__ &)
 
 G4VSolid (const G4VSolid &rhs)
 
G4VSolid & operator= (const G4VSolid &rhs)
 
G4double EstimateCubicVolume (G4int nStat, G4double epsilon) const
 
G4double EstimateSurfaceArea (G4int nStat, G4double ell) const
 

Protected Member Functions

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 G4OTubs
void Initialize ()
 
void CheckSPhiAngle (G4double sPhi)
 
void CheckDPhiAngle (G4double dPhi)
 
void CheckPhiAngles (G4double sPhi, G4double dPhi)
 
void InitializeTrigonometry ()
 
- 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
 

Additional Inherited Members

- Protected Types inherited from G4OTubs
enum  ESide {
  kNull, kRMin, kRMax, kSPhi,
  kEPhi, kPZ, kMZ
}
 
enum  ENorm {
  kNRMin, kNRMax, kNSPhi, kNEPhi,
  kNZ
}
 
- Protected Attributes inherited from G4OTubs
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
 
G4double halfCarTolerance
 
G4double halfRadTolerance
 
G4double halfAngTolerance
 
- Protected Attributes inherited from G4CSGSolid
G4double fCubicVolume
 
G4double fSurfaceArea
 
G4bool fRebuildPolyhedron
 
G4Polyhedron * fpPolyhedron
 
- Protected Attributes inherited from G4VSolid
G4double kCarTolerance
 

Detailed Description

Definition at line 51 of file G4CutTubs.hh.

Constructor & Destructor Documentation

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

Definition at line 65 of file G4CutTubs.cc.

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

Here is the call graph for this function:

Here is the caller graph for this function:

G4CutTubs::~G4CutTubs ( )

Definition at line 159 of file G4CutTubs.cc.

160 {
161 }
G4CutTubs::G4CutTubs ( __void__ &  a)

Definition at line 148 of file G4CutTubs.cc.

149  : G4OTubs(a), fLowNorm(G4ThreeVector()),
150  fHighNorm(G4ThreeVector()), fPhiFullCutTube(false),
151  halfCarTolerance(0.), halfRadTolerance(0.), halfAngTolerance(0.)
152 {
153 }
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:42
a
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
G4OTubs::G4OTubs
G4OTubs(const G4String &pName, G4double pRMin, G4double pRMax, G4double pDz, G4double pSPhi, G4double pDPhi)
Definition: G4OTubs.cc:59
G4CutTubs::G4CutTubs ( const G4CutTubs &  rhs)

Definition at line 167 of file G4CutTubs.cc.

168  : G4OTubs(rhs), fLowNorm(rhs.fLowNorm), fHighNorm(rhs.fHighNorm),
169  fPhiFullCutTube(rhs.fPhiFullCutTube),
170  halfCarTolerance(rhs.halfCarTolerance),
171  halfRadTolerance(rhs.halfRadTolerance),
172  halfAngTolerance(rhs.halfAngTolerance)
173 {
174 }
G4OTubs::G4OTubs
G4OTubs(const G4String &pName, G4double pRMin, G4double pRMax, G4double pDz, G4double pSPhi, G4double pDPhi)
Definition: G4OTubs.cc:59

Member Function Documentation

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

Reimplemented from G4OTubs.

Definition at line 558 of file G4CutTubs.cc.

559 {
560  ENorm side ;
561  G4ThreeVector norm ;
562  G4double rho, phi ;
563  G4double distZLow,distZHigh,distZ;
564  G4double distRMin, distRMax, distSPhi, distEPhi, distMin ;
565  G4ThreeVector vZ=G4ThreeVector(0,0,fDz);
566 
567  rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
568 
569  distRMin = std::fabs(rho - fRMin) ;
570  distRMax = std::fabs(rho - fRMax) ;
571 
572  //dist to Low Cut
573  //
574  distZLow =std::fabs((p+vZ).dot(fLowNorm));
575 
576  //dist to High Cut
577  //
578  distZHigh = std::fabs((p-vZ).dot(fHighNorm));
579  distZ=std::min(distZLow,distZHigh);
580 
581  if (distRMin < distRMax) // First minimum
582  {
583  if ( distZ < distRMin )
584  {
585  distMin = distZ ;
586  side = kNZ ;
587  }
588  else
589  {
590  distMin = distRMin ;
591  side = kNRMin ;
592  }
593  }
594  else
595  {
596  if ( distZ < distRMax )
597  {
598  distMin = distZ ;
599  side = kNZ ;
600  }
601  else
602  {
603  distMin = distRMax ;
604  side = kNRMax ;
605  }
606  }
607  if (!fPhiFullCutTube && rho ) // Protected against (0,0,z)
608  {
609  phi = std::atan2(p.y(),p.x()) ;
610 
611  if ( phi < 0 ) { phi += twopi; }
612 
613  if ( fSPhi < 0 )
614  {
615  distSPhi = std::fabs(phi - (fSPhi + twopi))*rho ;
616  }
617  else
618  {
619  distSPhi = std::fabs(phi - fSPhi)*rho ;
620  }
621  distEPhi = std::fabs(phi - fSPhi - fDPhi)*rho ;
622 
623  if (distSPhi < distEPhi) // Find new minimum
624  {
625  if ( distSPhi < distMin )
626  {
627  side = kNSPhi ;
628  }
629  }
630  else
631  {
632  if ( distEPhi < distMin )
633  {
634  side = kNEPhi ;
635  }
636  }
637  }
638  switch ( side )
639  {
640  case kNRMin : // Inner radius
641  {
642  norm = G4ThreeVector(-p.x()/rho, -p.y()/rho, 0) ;
643  break ;
644  }
645  case kNRMax : // Outer radius
646  {
647  norm = G4ThreeVector(p.x()/rho, p.y()/rho, 0) ;
648  break ;
649  }
650  case kNZ : // + or - dz
651  {
652  if ( distZHigh > distZLow ) { norm = fHighNorm ; }
653  else { norm = fLowNorm; }
654  break ;
655  }
656  case kNSPhi:
657  {
658  norm = G4ThreeVector(std::sin(fSPhi), -std::cos(fSPhi), 0) ;
659  break ;
660  }
661  case kNEPhi:
662  {
663  norm = G4ThreeVector(-std::sin(fSPhi+fDPhi), std::cos(fSPhi+fDPhi), 0) ;
664  break;
665  }
666  default: // Should never reach this case ...
667  {
668  DumpInfo();
669  G4Exception("G4CutTubs::ApproxSurfaceNormal()",
670  "GeomSolids1002", JustWarning,
671  "Undefined side for valid surface normal to solid.");
672  break ;
673  }
674  }
675  return norm;
676 }
G4OTubs::fDPhi
G4double fDPhi
Definition: G4OTubs.hh:177
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:42
CLHEP::Hep3Vector::x
double x() const
G4OTubs::fRMax
G4double fRMax
Definition: G4OTubs.hh:177
JustWarning
Definition: G4ExceptionSeverity.hh:64
G4VSolid::DumpInfo
void DumpInfo() const
ENorm
ENorm
Definition: G4Cons.cc:80
twopi
static constexpr double twopi
Definition: G4SIunits.hh:76
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4OTubs::fDz
G4double fDz
Definition: G4OTubs.hh:177
CLHEP::Hep3Vector::y
double y() const
G4OTubs::fSPhi
G4double fSPhi
Definition: G4OTubs.hh:177
G4INCL::Math::min
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
Definition: G4INCLGlobals.hh:117
G4OTubs::fRMin
G4double fRMin
Definition: G4OTubs.hh:177
G4double
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

Here is the caller graph for this function:

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

Implements G4VSolid.

Definition at line 262 of file G4CutTubs.cc.

267 {
268  G4ThreeVector bmin, bmax;
269  G4bool exist;
270 
271  // Get bounding box
272  Extent(bmin,bmax);
273 
274  // Check bounding box
275  G4BoundingEnvelope bbox(bmin,bmax);
276 #ifdef G4BBOX_EXTENT
277  if (true) return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
278 #endif
279  if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVoxelLimit,pTransform,pMin,pMax))
280  {
281  return exist = (pMin < pMax) ? true : false;
282  }
283 
284  // Get parameters of the solid
285  G4double rmin = GetInnerRadius();
286  G4double rmax = GetOuterRadius();
287  G4double dphi = GetDeltaPhiAngle();
288  G4double zmin = bmin.z();
289  G4double zmax = bmax.z();
290 
291  // Find bounding envelope and calculate extent
292  //
293  const G4int NSTEPS = 24; // number of steps for whole circle
294  G4double astep = twopi/NSTEPS; // max angle for one step
295  G4int ksteps = (dphi <= astep) ? 1 : (G4int)((dphi-deg)/astep) + 1;
296  G4double ang = dphi/ksteps;
297 
298  G4double sinHalf = std::sin(0.5*ang);
299  G4double cosHalf = std::cos(0.5*ang);
300  G4double sinStep = 2.*sinHalf*cosHalf;
301  G4double cosStep = 1. - 2.*sinHalf*sinHalf;
302  G4double rext = rmax/cosHalf;
303 
304  // bounding envelope for full cylinder consists of two polygons,
305  // in other cases it is a sequence of quadrilaterals
306  if (rmin == 0 && dphi == twopi)
307  {
308  G4double sinCur = sinHalf;
309  G4double cosCur = cosHalf;
310 
311  G4ThreeVectorList baseA(NSTEPS),baseB(NSTEPS);
312  for (G4int k=0; k<NSTEPS; ++k)
313  {
314  baseA[k].set(rext*cosCur,rext*sinCur,zmin);
315  baseB[k].set(rext*cosCur,rext*sinCur,zmax);
316 
317  G4double sinTmp = sinCur;
318  sinCur = sinCur*cosStep + cosCur*sinStep;
319  cosCur = cosCur*cosStep - sinTmp*sinStep;
320  }
321  std::vector<const G4ThreeVectorList *> polygons(2);
322  polygons[0] = &baseA;
323  polygons[1] = &baseB;
324  G4BoundingEnvelope benv(bmin,bmax,polygons);
325  exist = benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
326  }
327  else
328  {
329  G4double sinStart = GetSinStartPhi();
330  G4double cosStart = GetCosStartPhi();
331  G4double sinEnd = GetSinEndPhi();
332  G4double cosEnd = GetCosEndPhi();
333  G4double sinCur = sinStart*cosHalf + cosStart*sinHalf;
334  G4double cosCur = cosStart*cosHalf - sinStart*sinHalf;
335 
336  // set quadrilaterals
337  G4ThreeVectorList pols[NSTEPS+2];
338  for (G4int k=0; k<ksteps+2; ++k) pols[k].resize(4);
339  pols[0][0].set(rmin*cosStart,rmin*sinStart,zmax);
340  pols[0][1].set(rmin*cosStart,rmin*sinStart,zmin);
341  pols[0][2].set(rmax*cosStart,rmax*sinStart,zmin);
342  pols[0][3].set(rmax*cosStart,rmax*sinStart,zmax);
343  for (G4int k=1; k<ksteps+1; ++k)
344  {
345  pols[k][0].set(rmin*cosCur,rmin*sinCur,zmax);
346  pols[k][1].set(rmin*cosCur,rmin*sinCur,zmin);
347  pols[k][2].set(rext*cosCur,rext*sinCur,zmin);
348  pols[k][3].set(rext*cosCur,rext*sinCur,zmax);
349 
350  G4double sinTmp = sinCur;
351  sinCur = sinCur*cosStep + cosCur*sinStep;
352  cosCur = cosCur*cosStep - sinTmp*sinStep;
353  }
354  pols[ksteps+1][0].set(rmin*cosEnd,rmin*sinEnd,zmax);
355  pols[ksteps+1][1].set(rmin*cosEnd,rmin*sinEnd,zmin);
356  pols[ksteps+1][2].set(rmax*cosEnd,rmax*sinEnd,zmin);
357  pols[ksteps+1][3].set(rmax*cosEnd,rmax*sinEnd,zmax);
358 
359  // set envelope and calculate extent
360  std::vector<const G4ThreeVectorList *> polygons;
361  polygons.resize(ksteps+2);
362  for (G4int k=0; k<ksteps+2; ++k) polygons[k] = &pols[k];
363  G4BoundingEnvelope benv(bmin,bmax,polygons);
364  exist = benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
365  }
366  return exist;
367 }
G4OTubs::GetCosStartPhi
G4double GetCosStartPhi() const
G4OTubs::GetSinStartPhi
G4double GetSinStartPhi() const
G4int
int G4int
Definition: G4Types.hh:78
CLHEP::Hep3Vector::z
double z() const
twopi
static constexpr double twopi
Definition: G4SIunits.hh:76
G4bool
bool G4bool
Definition: G4Types.hh:79
G4OTubs::GetSinEndPhi
G4double GetSinEndPhi() const
G4ThreeVectorList
std::vector< G4ThreeVector > G4ThreeVectorList
Definition: G4BoundingEnvelope.hh:61
G4OTubs::GetOuterRadius
G4double GetOuterRadius() const
G4OTubs::GetCosEndPhi
G4double GetCosEndPhi() const
G4OTubs::GetInnerRadius
G4double GetInnerRadius() const
G4CutTubs::Extent
void Extent(G4ThreeVector &pMin, G4ThreeVector &pMax) const
Definition: G4CutTubs.cc:205
G4OTubs::GetDeltaPhiAngle
G4double GetDeltaPhiAngle() const
G4double
double G4double
Definition: G4Types.hh:76
deg
static constexpr double deg
Definition: G4SIunits.hh:152
test1.astep
tuple astep
Definition: test1.py:13

Here is the call graph for this function:

G4VSolid * G4CutTubs::Clone ( ) const
virtual

Reimplemented from G4VSolid.

Definition at line 1700 of file G4CutTubs.cc.

1701 {
1702  return new G4CutTubs(*this);
1703 }
G4CutTubs::G4CutTubs
G4CutTubs(const G4String &pName, G4double pRMin, G4double pRMax, G4double pDz, G4double pSPhi, G4double pDPhi, G4ThreeVector pLowNorm, G4ThreeVector pHighNorm)
Definition: G4CutTubs.cc:65

Here is the call graph for this function:

G4Polyhedron * G4CutTubs::CreatePolyhedron ( ) const
virtual

Reimplemented from G4VSolid.

Definition at line 1813 of file G4CutTubs.cc.

1814 {
1815  typedef G4double G4double3[3];
1816  typedef G4int G4int4[4];
1817 
1818  G4Polyhedron *ph = new G4Polyhedron;
1819  G4Polyhedron *ph1 = G4OTubs::CreatePolyhedron();
1820  G4int nn=ph1->GetNoVertices();
1821  G4int nf=ph1->GetNoFacets();
1822  G4double3* xyz = new G4double3[nn]; // number of nodes
1823  G4int4* faces = new G4int4[nf] ; // number of faces
1824 
1825  for(G4int i=0;i<nn;i++)
1826  {
1827  xyz[i][0]=ph1->GetVertex(i+1).x();
1828  xyz[i][1]=ph1->GetVertex(i+1).y();
1829  G4double tmpZ=ph1->GetVertex(i+1).z();
1830  if(tmpZ>=fDz-kCarTolerance)
1831  {
1832  xyz[i][2]=GetCutZ(G4ThreeVector(xyz[i][0],xyz[i][1],fDz));
1833  }
1834  else if(tmpZ<=-fDz+kCarTolerance)
1835  {
1836  xyz[i][2]=GetCutZ(G4ThreeVector(xyz[i][0],xyz[i][1],-fDz));
1837  }
1838  else
1839  {
1840  xyz[i][2]=tmpZ;
1841  }
1842  }
1843  G4int iNodes[4];
1844  G4int *iEdge=0;
1845  G4int n;
1846  for(G4int i=0;i<nf;i++)
1847  {
1848  ph1->GetFacet(i+1,n,iNodes,iEdge);
1849  for(G4int k=0;k<n;k++)
1850  {
1851  faces[i][k]=iNodes[k];
1852  }
1853  for(G4int k=n;k<4;k++)
1854  {
1855  faces[i][k]=0;
1856  }
1857  }
1858  ph->createPolyhedron(nn,nf,xyz,faces);
1859 
1860  delete [] xyz;
1861  delete [] faces;
1862  delete ph1;
1863 
1864  return ph;
1865 }
HepPolyhedron::GetFacet
void GetFacet(G4int iFace, G4int &n, G4int *iNodes, G4int *edgeFlags=0, G4int *iFaces=0) const
Definition: HepPolyhedron.cc:1078
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:42
HepPolyhedron::createPolyhedron
G4int createPolyhedron(G4int Nnodes, G4int Nfaces, const G4double xyz[][3], const G4int faces[][4])
Definition: HepPolyhedron.cc:1371
G4int
int G4int
Definition: G4Types.hh:78
n
const G4int n
Definition: G4UrQMD1_3Interface.hh:144
G4OTubs::CreatePolyhedron
G4Polyhedron * CreatePolyhedron() const
Definition: G4OTubs.cc:1728
HepPolyhedron::GetNoVertices
G4int GetNoVertices() const
Definition: HepPolyhedron.h:252
G4OTubs::fDz
G4double fDz
Definition: G4OTubs.hh:177
G4VSolid::kCarTolerance
G4double kCarTolerance
Definition: G4VSolid.hh:307
G4double
double G4double
Definition: G4Types.hh:76
HepPolyhedron::GetNoFacets
G4int GetNoFacets() const
Definition: HepPolyhedron.h:256
HepPolyhedron::GetVertex
G4Point3D GetVertex(G4int index) const
Definition: HepPolyhedron.cc:898
G4CutTubs::GetCutZ
G4double GetCutZ(const G4ThreeVector &p) const
Definition: G4CutTubs.cc:1896

Here is the call graph for this function:

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

Implements G4VSolid.

Definition at line 1808 of file G4CutTubs.cc.

1809 {
1810  scene.AddSolid (*this) ;
1811 }
G4VGraphicsScene::AddSolid
virtual void AddSolid(const G4Box &)=0

Here is the call graph for this function:

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

Implements G4VSolid.

Definition at line 700 of file G4CutTubs.cc.

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

Here is the call graph for this function:

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

Implements G4VSolid.

Definition at line 1153 of file G4CutTubs.cc.

1154 {
1155  G4double safRMin,safRMax,safZLow,safZHigh,safePhi,safe,rho,cosPsi;
1156  G4ThreeVector vZ=G4ThreeVector(0,0,fDz);
1157 
1158  // Distance to R
1159  //
1160  rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
1161 
1162  safRMin = fRMin- rho ;
1163  safRMax = rho - fRMax ;
1164 
1165  // Distances to ZCut(Low/High)
1166 
1167  // Dist to Low Cut
1168  //
1169  safZLow = (p+vZ).dot(fLowNorm);
1170 
1171  // Dist to High Cut
1172  //
1173  safZHigh = (p-vZ).dot(fHighNorm);
1174 
1175  safe = std::max(safZLow,safZHigh);
1176 
1177  if ( safRMin > safe ) { safe = safRMin; }
1178  if ( safRMax> safe ) { safe = safRMax; }
1179 
1180  // Distance to Phi
1181  //
1182  if ( (!fPhiFullCutTube) && (rho) )
1183  {
1184  // Psi=angle from central phi to point
1185  //
1186  cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/rho ;
1187 
1188  if ( cosPsi < std::cos(fDPhi*0.5) )
1189  {
1190  // Point lies outside phi range
1191 
1192  if ( (p.y()*cosCPhi - p.x()*sinCPhi) <= 0 )
1193  {
1194  safePhi = std::fabs(p.x()*sinSPhi - p.y()*cosSPhi) ;
1195  }
1196  else
1197  {
1198  safePhi = std::fabs(p.x()*sinEPhi - p.y()*cosEPhi) ;
1199  }
1200  if ( safePhi > safe ) { safe = safePhi; }
1201  }
1202  }
1203  if ( safe < 0 ) { safe = 0; }
1204 
1205  return safe ;
1206 }
G4OTubs::fDPhi
G4double fDPhi
Definition: G4OTubs.hh:177
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:42
CLHEP::Hep3Vector::x
double x() const
G4OTubs::cosSPhi
G4double cosSPhi
Definition: G4OTubs.hh:181
G4OTubs::fRMax
G4double fRMax
Definition: G4OTubs.hh:177
G4OTubs::cosEPhi
G4double cosEPhi
Definition: G4OTubs.hh:181
G4OTubs::sinSPhi
G4double sinSPhi
Definition: G4OTubs.hh:181
G4OTubs::cosCPhi
G4double cosCPhi
Definition: G4OTubs.hh:181
G4OTubs::sinEPhi
G4double sinEPhi
Definition: G4OTubs.hh:181
G4INCL::Math::max
T max(const T t1, const T t2)
brief Return the largest of the two arguments
Definition: G4INCLGlobals.hh:112
G4OTubs::fDz
G4double fDz
Definition: G4OTubs.hh:177
CLHEP::Hep3Vector::y
double y() const
G4OTubs::fRMin
G4double fRMin
Definition: G4OTubs.hh:177
G4double
double G4double
Definition: G4Types.hh:76
G4OTubs::sinCPhi
G4double sinCPhi
Definition: G4OTubs.hh:181

Here is the call graph for this function:

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 1213 of file G4CutTubs.cc.

1218 {
1219  ESide side=kNull , sider=kNull, sidephi=kNull ;
1220  G4double snxt=kInfinity, srd=kInfinity,sz=kInfinity, sphi=kInfinity ;
1221  G4double deltaR, t1, t2, t3, b, c, d2, roMin2 ;
1222  G4double distZLow,distZHigh,calfH,calfL;
1223  G4ThreeVector vZ=G4ThreeVector(0,0,fDz);
1224 
1225  // Vars for phi intersection:
1226  //
1227  G4double pDistS, compS, pDistE, compE, sphi2, xi, yi, vphi, roi2 ;
1228 
1229  // Z plane intersection
1230  // Distances to ZCut(Low/High)
1231 
1232  // dist to Low Cut
1233  //
1234  distZLow =(p+vZ).dot(fLowNorm);
1235 
1236  // dist to High Cut
1237  //
1238  distZHigh = (p-vZ).dot(fHighNorm);
1239 
1240  calfH = v.dot(fHighNorm);
1241  calfL = v.dot(fLowNorm);
1242 
1243  if (calfH > 0 )
1244  {
1245  if ( distZHigh < halfCarTolerance )
1246  {
1247  snxt = -distZHigh/calfH ;
1248  side = kPZ ;
1249  }
1250  else
1251  {
1252  if (calcNorm)
1253  {
1254  *n = G4ThreeVector(0,0,1) ;
1255  *validNorm = true ;
1256  }
1257  return snxt = 0 ;
1258  }
1259  }
1260  if ( calfL>0)
1261  {
1262 
1263  if ( distZLow < halfCarTolerance )
1264  {
1265  sz = -distZLow/calfL ;
1266  if(sz<snxt){
1267  snxt=sz;
1268  side = kMZ ;
1269  }
1270 
1271  }
1272  else
1273  {
1274  if (calcNorm)
1275  {
1276  *n = G4ThreeVector(0,0,-1) ;
1277  *validNorm = true ;
1278  }
1279  return snxt = 0.0 ;
1280  }
1281  }
1282  if((calfH<=0)&&(calfL<=0))
1283  {
1284  snxt = kInfinity ; // Travel perpendicular to z axis
1285  side = kNull;
1286  }
1287  // Radial Intersections
1288  //
1289  // Find intersection with cylinders at rmax/rmin
1290  // Intersection point (xi,yi,zi) on line x=p.x+t*v.x etc.
1291  //
1292  // Intersects with x^2+y^2=R^2
1293  //
1294  // 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
1295  //
1296  // t1 t2 t3
1297 
1298  t1 = 1.0 - v.z()*v.z() ; // since v normalised
1299  t2 = p.x()*v.x() + p.y()*v.y() ;
1300  t3 = p.x()*p.x() + p.y()*p.y() ;
1301 
1302  if ( snxt > 10*(fDz+fRMax) ) { roi2 = 2*fRMax*fRMax; }
1303  else { roi2 = snxt*snxt*t1 + 2*snxt*t2 + t3; } // radius^2 on +-fDz
1304 
1305  if ( t1 > 0 ) // Check not parallel
1306  {
1307  // Calculate srd, r exit distance
1308 
1309  if ( (t2 >= 0.0) && (roi2 > fRMax*(fRMax + kRadTolerance)) )
1310  {
1311  // Delta r not negative => leaving via rmax
1312 
1313  deltaR = t3 - fRMax*fRMax ;
1314 
1315  // NOTE: Should use rho-fRMax<-kRadTolerance*0.5
1316  // - avoid sqrt for efficiency
1317 
1318  if ( deltaR < -kRadTolerance*fRMax )
1319  {
1320  b = t2/t1 ;
1321  c = deltaR/t1 ;
1322  d2 = b*b-c;
1323  if( d2 >= 0 ) { srd = c/( -b - std::sqrt(d2)); }
1324  else { srd = 0.; }
1325  sider = kRMax ;
1326  }
1327  else
1328  {
1329  // On tolerant boundary & heading outwards (or perpendicular to)
1330  // outer radial surface -> leaving immediately
1331 
1332  if ( calcNorm )
1333  {
1334  *n = G4ThreeVector(p.x()/fRMax,p.y()/fRMax,0) ;
1335  *validNorm = true ;
1336  }
1337  return snxt = 0 ; // Leaving by rmax immediately
1338  }
1339  }
1340  else if ( t2 < 0. ) // i.e. t2 < 0; Possible rmin intersection
1341  {
1342  roMin2 = t3 - t2*t2/t1 ; // min ro2 of the plane of movement
1343 
1344  if ( fRMin && (roMin2 < fRMin*(fRMin - kRadTolerance)) )
1345  {
1346  deltaR = t3 - fRMin*fRMin ;
1347  b = t2/t1 ;
1348  c = deltaR/t1 ;
1349  d2 = b*b - c ;
1350 
1351  if ( d2 >= 0 ) // Leaving via rmin
1352  {
1353  // NOTE: SHould use rho-rmin>kRadTolerance*0.5
1354  // - avoid sqrt for efficiency
1355 
1356  if (deltaR > kRadTolerance*fRMin)
1357  {
1358  srd = c/(-b+std::sqrt(d2));
1359  sider = kRMin ;
1360  }
1361  else
1362  {
1363  if ( calcNorm ) { *validNorm = false; } // Concave side
1364  return snxt = 0.0;
1365  }
1366  }
1367  else // No rmin intersect -> must be rmax intersect
1368  {
1369  deltaR = t3 - fRMax*fRMax ;
1370  c = deltaR/t1 ;
1371  d2 = b*b-c;
1372  if( d2 >=0. )
1373  {
1374  srd = -b + std::sqrt(d2) ;
1375  sider = kRMax ;
1376  }
1377  else // Case: On the border+t2<kRadTolerance
1378  // (v is perpendicular to the surface)
1379  {
1380  if (calcNorm)
1381  {
1382  *n = G4ThreeVector(p.x()/fRMax,p.y()/fRMax,0) ;
1383  *validNorm = true ;
1384  }
1385  return snxt = 0.0;
1386  }
1387  }
1388  }
1389  else if ( roi2 > fRMax*(fRMax + kRadTolerance) )
1390  // No rmin intersect -> must be rmax intersect
1391  {
1392  deltaR = t3 - fRMax*fRMax ;
1393  b = t2/t1 ;
1394  c = deltaR/t1;
1395  d2 = b*b-c;
1396  if( d2 >= 0 )
1397  {
1398  srd = -b + std::sqrt(d2) ;
1399  sider = kRMax ;
1400  }
1401  else // Case: On the border+t2<kRadTolerance
1402  // (v is perpendicular to the surface)
1403  {
1404  if (calcNorm)
1405  {
1406  *n = G4ThreeVector(p.x()/fRMax,p.y()/fRMax,0) ;
1407  *validNorm = true ;
1408  }
1409  return snxt = 0.0;
1410  }
1411  }
1412  }
1413  // Phi Intersection
1414 
1415  if ( !fPhiFullCutTube )
1416  {
1417  // add angle calculation with correction
1418  // of the difference in domain of atan2 and Sphi
1419  //
1420  vphi = std::atan2(v.y(),v.x()) ;
1421 
1422  if ( vphi < fSPhi - halfAngTolerance ) { vphi += twopi; }
1423  else if ( vphi > fSPhi + fDPhi + halfAngTolerance ) { vphi -= twopi; }
1424 
1425 
1426  if ( p.x() || p.y() ) // Check if on z axis (rho not needed later)
1427  {
1428  // pDist -ve when inside
1429 
1430  pDistS = p.x()*sinSPhi - p.y()*cosSPhi ;
1431  pDistE = -p.x()*sinEPhi + p.y()*cosEPhi ;
1432 
1433  // Comp -ve when in direction of outwards normal
1434 
1435  compS = -sinSPhi*v.x() + cosSPhi*v.y() ;
1436  compE = sinEPhi*v.x() - cosEPhi*v.y() ;
1437 
1438  sidephi = kNull;
1439 
1440  if( ( (fDPhi <= pi) && ( (pDistS <= halfCarTolerance)
1441  && (pDistE <= halfCarTolerance) ) )
1442  || ( (fDPhi > pi) && !((pDistS > halfCarTolerance)
1443  && (pDistE > halfCarTolerance) ) ) )
1444  {
1445  // Inside both phi *full* planes
1446 
1447  if ( compS < 0 )
1448  {
1449  sphi = pDistS/compS ;
1450 
1451  if (sphi >= -halfCarTolerance)
1452  {
1453  xi = p.x() + sphi*v.x() ;
1454  yi = p.y() + sphi*v.y() ;
1455 
1456  // Check intersecting with correct half-plane
1457  // (if not -> no intersect)
1458  //
1459  if( (std::fabs(xi)<=kCarTolerance)
1460  && (std::fabs(yi)<=kCarTolerance) )
1461  {
1462  sidephi = kSPhi;
1463  if (((fSPhi-halfAngTolerance)<=vphi)
1464  &&((fSPhi+fDPhi+halfAngTolerance)>=vphi))
1465  {
1466  sphi = kInfinity;
1467  }
1468  }
1469  else if ( yi*cosCPhi-xi*sinCPhi >=0 )
1470  {
1471  sphi = kInfinity ;
1472  }
1473  else
1474  {
1475  sidephi = kSPhi ;
1476  if ( pDistS > -halfCarTolerance )
1477  {
1478  sphi = 0.0 ; // Leave by sphi immediately
1479  }
1480  }
1481  }
1482  else
1483  {
1484  sphi = kInfinity ;
1485  }
1486  }
1487  else
1488  {
1489  sphi = kInfinity ;
1490  }
1491 
1492  if ( compE < 0 )
1493  {
1494  sphi2 = pDistE/compE ;
1495 
1496  // Only check further if < starting phi intersection
1497  //
1498  if ( (sphi2 > -halfCarTolerance) && (sphi2 < sphi) )
1499  {
1500  xi = p.x() + sphi2*v.x() ;
1501  yi = p.y() + sphi2*v.y() ;
1502 
1503  if ((std::fabs(xi)<=kCarTolerance)&&(std::fabs(yi)<=kCarTolerance))
1504  {
1505  // Leaving via ending phi
1506  //
1507  if( !((fSPhi-halfAngTolerance <= vphi)
1508  &&(fSPhi+fDPhi+halfAngTolerance >= vphi)) )
1509  {
1510  sidephi = kEPhi ;
1511  if ( pDistE <= -halfCarTolerance ) { sphi = sphi2 ; }
1512  else { sphi = 0.0 ; }
1513  }
1514  }
1515  else // Check intersecting with correct half-plane
1516 
1517  if ( (yi*cosCPhi-xi*sinCPhi) >= 0)
1518  {
1519  // Leaving via ending phi
1520  //
1521  sidephi = kEPhi ;
1522  if ( pDistE <= -halfCarTolerance ) { sphi = sphi2 ; }
1523  else { sphi = 0.0 ; }
1524  }
1525  }
1526  }
1527  }
1528  else
1529  {
1530  sphi = kInfinity ;
1531  }
1532  }
1533  else
1534  {
1535  // On z axis + travel not || to z axis -> if phi of vector direction
1536  // within phi of shape, Step limited by rmax, else Step =0
1537 
1538  if ( (fSPhi - halfAngTolerance <= vphi)
1539  && (vphi <= fSPhi + fDPhi + halfAngTolerance ) )
1540  {
1541  sphi = kInfinity ;
1542  }
1543  else
1544  {
1545  sidephi = kSPhi ; // arbitrary
1546  sphi = 0.0 ;
1547  }
1548  }
1549  if (sphi < snxt) // Order intersecttions
1550  {
1551  snxt = sphi ;
1552  side = sidephi ;
1553  }
1554  }
1555  if (srd < snxt) // Order intersections
1556  {
1557  snxt = srd ;
1558  side = sider ;
1559  }
1560  }
1561  if (calcNorm)
1562  {
1563  switch(side)
1564  {
1565  case kRMax:
1566  // Note: returned vector not normalised
1567  // (divide by fRMax for unit vector)
1568  //
1569  xi = p.x() + snxt*v.x() ;
1570  yi = p.y() + snxt*v.y() ;
1571  *n = G4ThreeVector(xi/fRMax,yi/fRMax,0) ;
1572  *validNorm = true ;
1573  break ;
1574 
1575  case kRMin:
1576  *validNorm = false ; // Rmin is inconvex
1577  break ;
1578 
1579  case kSPhi:
1580  if ( fDPhi <= pi )
1581  {
1582  *n = G4ThreeVector(sinSPhi,-cosSPhi,0) ;
1583  *validNorm = true ;
1584  }
1585  else
1586  {
1587  *validNorm = false ;
1588  }
1589  break ;
1590 
1591  case kEPhi:
1592  if (fDPhi <= pi)
1593  {
1594  *n = G4ThreeVector(-sinEPhi,cosEPhi,0) ;
1595  *validNorm = true ;
1596  }
1597  else
1598  {
1599  *validNorm = false ;
1600  }
1601  break ;
1602 
1603  case kPZ:
1604  *n = fHighNorm ;
1605  *validNorm = true ;
1606  break ;
1607 
1608  case kMZ:
1609  *n = fLowNorm ;
1610  *validNorm = true ;
1611  break ;
1612 
1613  default:
1614  G4cout << G4endl ;
1615  DumpInfo();
1616  std::ostringstream message;
1617  G4int oldprc = message.precision(16);
1618  message << "Undefined side for valid surface normal to solid."
1619  << G4endl
1620  << "Position:" << G4endl << G4endl
1621  << "p.x() = " << p.x()/mm << " mm" << G4endl
1622  << "p.y() = " << p.y()/mm << " mm" << G4endl
1623  << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl
1624  << "Direction:" << G4endl << G4endl
1625  << "v.x() = " << v.x() << G4endl
1626  << "v.y() = " << v.y() << G4endl
1627  << "v.z() = " << v.z() << G4endl << G4endl
1628  << "Proposed distance :" << G4endl << G4endl
1629  << "snxt = " << snxt/mm << " mm" << G4endl ;
1630  message.precision(oldprc) ;
1631  G4Exception("G4CutTubs::DistanceToOut(p,v,..)", "GeomSolids1002",
1632  JustWarning, message);
1633  break ;
1634  }
1635  }
1636  if ( snxt<halfCarTolerance ) { snxt=0 ; }
1637  return snxt ;
1638 }
mm
static constexpr double mm
Definition: G4SIunits.hh:115
G4OTubs::fDPhi
G4double fDPhi
Definition: G4OTubs.hh:177
kInfinity
static const G4double kInfinity
Definition: geomdefs.hh:42
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:42
CLHEP::Hep3Vector::x
double x() const
CLHEP::Hep3Vector::dot
double dot(const Hep3Vector &) const
G4OTubs::cosSPhi
G4double cosSPhi
Definition: G4OTubs.hh:181
G4OTubs::fRMax
G4double fRMax
Definition: G4OTubs.hh:177
d2
static const G4double d2
Definition: G4ElectroNuclearCrossSection.cc:89
G4OTubs::cosEPhi
G4double cosEPhi
Definition: G4OTubs.hh:181
JustWarning
Definition: G4ExceptionSeverity.hh:64
G4int
int G4int
Definition: G4Types.hh:78
G4OTubs::sinSPhi
G4double sinSPhi
Definition: G4OTubs.hh:181
CLHEP::Hep3Vector::z
double z() const
G4VSolid::DumpInfo
void DumpInfo() const
G4OTubs::cosCPhi
G4double cosCPhi
Definition: G4OTubs.hh:181
twopi
static constexpr double twopi
Definition: G4SIunits.hh:76
test.b
tuple b
Definition: test.py:12
G4cout
G4GLOB_DLL std::ostream G4cout
G4OTubs::sinEPhi
G4double sinEPhi
Definition: G4OTubs.hh:181
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
plottest35.t1
tuple t1
Definition: plottest35.py:33
G4OTubs::fDz
G4double fDz
Definition: G4OTubs.hh:177
CLHEP::Hep3Vector::y
double y() const
G4OTubs::fSPhi
G4double fSPhi
Definition: G4OTubs.hh:177
G4OTubs::kRadTolerance
G4double kRadTolerance
Definition: G4OTubs.hh:173
G4endl
#define G4endl
Definition: G4ios.hh:61
G4VSolid::kCarTolerance
G4double kCarTolerance
Definition: G4VSolid.hh:307
G4OTubs::fRMin
G4double fRMin
Definition: G4OTubs.hh:177
pi
static constexpr double pi
Definition: G4SIunits.hh:75
ESide
ESide
Definition: G4Cons.cc:76
G4double
double G4double
Definition: G4Types.hh:76
test.c
tuple c
Definition: test.py:13
G4OTubs::sinCPhi
G4double sinCPhi
Definition: G4OTubs.hh:181

Here is the call graph for this function:

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

Implements G4VSolid.

Definition at line 1644 of file G4CutTubs.cc.

1645 {
1646  G4double safRMin,safRMax,safZLow,safZHigh,safePhi,safe,rho;
1647  G4ThreeVector vZ=G4ThreeVector(0,0,fDz);
1648 
1649  rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ; // Distance to R
1650 
1651  safRMin = rho - fRMin ;
1652  safRMax = fRMax - rho ;
1653 
1654  // Distances to ZCut(Low/High)
1655 
1656  // Dist to Low Cut
1657  //
1658  safZLow = std::fabs((p+vZ).dot(fLowNorm));
1659 
1660  // Dist to High Cut
1661  //
1662  safZHigh = std::fabs((p-vZ).dot(fHighNorm));
1663  safe = std::min(safZLow,safZHigh);
1664 
1665  if ( safRMin < safe ) { safe = safRMin; }
1666  if ( safRMax< safe ) { safe = safRMax; }
1667 
1668  // Check if phi divided, Calc distances closest phi plane
1669  //
1670  if ( !fPhiFullCutTube )
1671  {
1672  if ( p.y()*cosCPhi-p.x()*sinCPhi <= 0 )
1673  {
1674  safePhi = -(p.x()*sinSPhi - p.y()*cosSPhi) ;
1675  }
1676  else
1677  {
1678  safePhi = (p.x()*sinEPhi - p.y()*cosEPhi) ;
1679  }
1680  if (safePhi < safe) { safe = safePhi ; }
1681  }
1682  if ( safe < 0 ) { safe = 0; }
1683 
1684  return safe ;
1685 }
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:42
CLHEP::Hep3Vector::x
double x() const
G4OTubs::cosSPhi
G4double cosSPhi
Definition: G4OTubs.hh:181
G4OTubs::fRMax
G4double fRMax
Definition: G4OTubs.hh:177
G4OTubs::cosEPhi
G4double cosEPhi
Definition: G4OTubs.hh:181
G4OTubs::sinSPhi
G4double sinSPhi
Definition: G4OTubs.hh:181
G4OTubs::cosCPhi
G4double cosCPhi
Definition: G4OTubs.hh:181
G4OTubs::sinEPhi
G4double sinEPhi
Definition: G4OTubs.hh:181
G4OTubs::fDz
G4double fDz
Definition: G4OTubs.hh:177
CLHEP::Hep3Vector::y
double y() const
G4INCL::Math::min
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
Definition: G4INCLGlobals.hh:117
G4OTubs::fRMin
G4double fRMin
Definition: G4OTubs.hh:177
G4double
double G4double
Definition: G4Types.hh:76
G4OTubs::sinCPhi
G4double sinCPhi
Definition: G4OTubs.hh:181

Here is the call graph for this function:

void G4CutTubs::Extent ( G4ThreeVector &  pMin,
G4ThreeVector &  pMax 
) const
virtual

Reimplemented from G4VSolid.

Definition at line 205 of file G4CutTubs.cc.

206 {
207  G4double rmin = GetInnerRadius();
208  G4double rmax = GetOuterRadius();
209  G4double dz = GetZHalfLength();
210 
211  G4ThreeVector norm;
212  G4double xynorm, znorm;
213 
214  // get zmin
215  norm = GetLowNorm();
216  xynorm = std::sqrt(norm.x()*norm.x()+norm.y()*norm.y());
217  znorm = std::abs(norm.z());
218  G4double zmin = -(dz + rmax*xynorm/znorm);
219 
220  // get zmax
221  norm = GetHighNorm();
222  xynorm = std::sqrt(norm.x()*norm.x()+norm.y()*norm.y());
223  znorm = std::abs(norm.z());
224  G4double zmax = dz + rmax*xynorm/znorm;
225 
226  // Find bounding box
227  //
228  if (GetDeltaPhiAngle() < twopi)
229  {
230  G4TwoVector vmin,vmax;
231  G4GeomTools::DiskExtent(rmin,rmax,
232  GetSinStartPhi(),GetCosStartPhi(),
233  GetSinEndPhi(),GetCosEndPhi(),
234  vmin,vmax);
235  pMin.set(vmin.x(),vmin.y(), zmin);
236  pMax.set(vmax.x(),vmax.y(), zmax);
237  }
238  else
239  {
240  pMin.set(-rmax,-rmax, zmin);
241  pMax.set( rmax, rmax, zmax);
242  }
243 
244  // Check correctness of the bounding box
245  //
246  if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
247  {
248  std::ostringstream message;
249  message << "Bad bounding box (min >= max) for solid: "
250  << GetName() << " !"
251  << "\npMin = " << pMin
252  << "\npMax = " << pMax;
253  G4Exception("G4CutTubs::Extent()", "GeomMgt0001", JustWarning, message);
254  DumpInfo();
255  }
256 }
G4OTubs::GetCosStartPhi
G4double GetCosStartPhi() const
CLHEP::Hep3Vector::set
void set(double x, double y, double z)
G4VSolid::GetName
G4String GetName() const
CLHEP::Hep2Vector::y
double y() const
CLHEP::Hep2Vector::x
double x() const
CLHEP::Hep3Vector::x
double x() const
G4OTubs::GetSinStartPhi
G4double GetSinStartPhi() const
G4CutTubs::GetLowNorm
G4ThreeVector GetLowNorm() const
JustWarning
Definition: G4ExceptionSeverity.hh:64
CLHEP::Hep3Vector::z
double z() const
G4VSolid::DumpInfo
void DumpInfo() const
twopi
static constexpr double twopi
Definition: G4SIunits.hh:76
G4OTubs::GetSinEndPhi
G4double GetSinEndPhi() const
G4OTubs::GetZHalfLength
G4double GetZHalfLength() const
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4OTubs::GetOuterRadius
G4double GetOuterRadius() const
G4OTubs::GetCosEndPhi
G4double GetCosEndPhi() const
G4CutTubs::GetHighNorm
G4ThreeVector GetHighNorm() const
CLHEP::Hep3Vector::y
double y() const
G4OTubs::GetInnerRadius
G4double GetInnerRadius() const
G4OTubs::GetDeltaPhiAngle
G4double GetDeltaPhiAngle() const
G4double
double G4double
Definition: G4Types.hh:76
G4GeomTools::DiskExtent
static G4bool DiskExtent(G4double rmin, G4double rmax, G4double startPhi, G4double delPhi, G4TwoVector &pmin, G4TwoVector &pmax)
Definition: G4GeomTools.cc:378

Here is the call graph for this function:

Here is the caller graph for this function:

G4double G4CutTubs::GetCubicVolume ( )
inlinevirtual

Reimplemented from G4VSolid.

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

Definition at line 1896 of file G4CutTubs.cc.

1897 {
1898  G4double newz = p.z(); // p.z() should be either +fDz or -fDz
1899  if (p.z()<0)
1900  {
1901  if(fLowNorm.z()!=0.)
1902  {
1903  newz = -fDz-(p.x()*fLowNorm.x()+p.y()*fLowNorm.y())/fLowNorm.z();
1904  }
1905  }
1906  else
1907  {
1908  if(fHighNorm.z()!=0.)
1909  {
1910  newz = fDz-(p.x()*fHighNorm.x()+p.y()*fHighNorm.y())/fHighNorm.z();
1911  }
1912  }
1913  return newz;
1914 }
CLHEP::Hep3Vector::x
double x() const
CLHEP::Hep3Vector::z
double z() const
G4OTubs::fDz
G4double fDz
Definition: G4OTubs.hh:177
CLHEP::Hep3Vector::y
double y() const
G4double
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

Here is the caller graph for this function:

G4GeometryType G4CutTubs::GetEntityType ( ) const
virtual

Implements G4VSolid.

Definition at line 1691 of file G4CutTubs.cc.

1692 {
1693  return G4String("G4CutTubs");
1694 }
G4ThreeVector G4CutTubs::GetHighNorm ( ) const
inline

Here is the caller graph for this function:

G4ThreeVector G4CutTubs::GetLowNorm ( ) const
inline

Here is the caller graph for this function:

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

Definition at line 1920 of file G4CutTubs.cc.

1922 {
1923  G4double phiLow = std::atan2(fLowNorm.y(),fLowNorm.x());
1924  G4double phiHigh= std::atan2(fHighNorm.y(),fHighNorm.x());
1925 
1926  G4double xc=0, yc=0,z1;
1927  G4double z[8];
1928  G4bool in_range_low = false;
1929  G4bool in_range_hi = false;
1930 
1931  G4int i;
1932  for (i=0; i<2; i++)
1933  {
1934  if (phiLow<0) { phiLow+=twopi; }
1935  G4double ddp = phiLow-fSPhi;
1936  if (ddp<0) { ddp += twopi; }
1937  if (ddp <= fDPhi)
1938  {
1939  xc = fRMin*std::cos(phiLow);
1940  yc = fRMin*std::sin(phiLow);
1941  z1 = GetCutZ(G4ThreeVector(xc, yc, -fDz));
1942  xc = fRMax*std::cos(phiLow);
1943  yc = fRMax*std::sin(phiLow);
1944  z1 = std::min(z1, GetCutZ(G4ThreeVector(xc, yc, -fDz)));
1945  if (in_range_low) { zmin = std::min(zmin, z1); }
1946  else { zmin = z1; }
1947  in_range_low = true;
1948  }
1949  phiLow += pi;
1950  if (phiLow>twopi) { phiLow-=twopi; }
1951  }
1952  for (i=0; i<2; i++)
1953  {
1954  if (phiHigh<0) { phiHigh+=twopi; }
1955  G4double ddp = phiHigh-fSPhi;
1956  if (ddp<0) { ddp += twopi; }
1957  if (ddp <= fDPhi)
1958  {
1959  xc = fRMin*std::cos(phiHigh);
1960  yc = fRMin*std::sin(phiHigh);
1961  z1 = GetCutZ(G4ThreeVector(xc, yc, fDz));
1962  xc = fRMax*std::cos(phiHigh);
1963  yc = fRMax*std::sin(phiHigh);
1964  z1 = std::min(z1, GetCutZ(G4ThreeVector(xc, yc, fDz)));
1965  if (in_range_hi) { zmax = std::min(zmax, z1); }
1966  else { zmax = z1; }
1967  in_range_hi = true;
1968  }
1969  phiHigh += pi;
1970  if (phiHigh>twopi) { phiHigh-=twopi; }
1971  }
1972 
1973  xc = fRMin*std::cos(fSPhi);
1974  yc = fRMin*std::sin(fSPhi);
1975  z[0] = GetCutZ(G4ThreeVector(xc, yc, -fDz));
1976  z[4] = GetCutZ(G4ThreeVector(xc, yc, fDz));
1977 
1978  xc = fRMin*std::cos(fSPhi+fDPhi);
1979  yc = fRMin*std::sin(fSPhi+fDPhi);
1980  z[1] = GetCutZ(G4ThreeVector(xc, yc, -fDz));
1981  z[5] = GetCutZ(G4ThreeVector(xc, yc, fDz));
1982 
1983  xc = fRMax*std::cos(fSPhi);
1984  yc = fRMax*std::sin(fSPhi);
1985  z[2] = GetCutZ(G4ThreeVector(xc, yc, -fDz));
1986  z[6] = GetCutZ(G4ThreeVector(xc, yc, fDz));
1987 
1988  xc = fRMax*std::cos(fSPhi+fDPhi);
1989  yc = fRMax*std::sin(fSPhi+fDPhi);
1990  z[3] = GetCutZ(G4ThreeVector(xc, yc, -fDz));
1991  z[7] = GetCutZ(G4ThreeVector(xc, yc, fDz));
1992 
1993  // Find min/max
1994 
1995  z1=z[0];
1996  for (i = 1; i < 4; i++)
1997  {
1998  if(z[i] < z[i-1])z1=z[i];
1999  }
2000 
2001  if (in_range_low)
2002  {
2003  zmin = std::min(zmin, z1);
2004  }
2005  else
2006  {
2007  zmin = z1;
2008  }
2009  z1=z[4];
2010  for (i = 1; i < 4; i++)
2011  {
2012  if(z[4+i] > z[4+i-1]) { z1=z[4+i]; }
2013  }
2014 
2015  if (in_range_hi) { zmax = std::max(zmax, z1); }
2016  else { zmax = z1; }
2017 }
G4OTubs::fDPhi
G4double fDPhi
Definition: G4OTubs.hh:177
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:42
CLHEP::Hep3Vector::x
double x() const
G4OTubs::fRMax
G4double fRMax
Definition: G4OTubs.hh:177
G4int
int G4int
Definition: G4Types.hh:78
twopi
static constexpr double twopi
Definition: G4SIunits.hh:76
G4bool
bool G4bool
Definition: G4Types.hh:79
G4INCL::Math::max
T max(const T t1, const T t2)
brief Return the largest of the two arguments
Definition: G4INCLGlobals.hh:112
G4OTubs::fDz
G4double fDz
Definition: G4OTubs.hh:177
CLHEP::Hep3Vector::y
double y() const
G4OTubs::fSPhi
G4double fSPhi
Definition: G4OTubs.hh:177
G4INCL::Math::min
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
Definition: G4INCLGlobals.hh:117
test.z
tuple z
Definition: test.py:28
G4OTubs::fRMin
G4double fRMin
Definition: G4OTubs.hh:177
pi
static constexpr double pi
Definition: G4SIunits.hh:75
G4double
double G4double
Definition: G4Types.hh:76
G4CutTubs::GetCutZ
G4double GetCutZ(const G4ThreeVector &p) const
Definition: G4CutTubs.cc:1896

Here is the call graph for this function:

G4ThreeVector G4CutTubs::GetPointOnSurface ( ) const
virtual

Reimplemented from G4VSolid.

Definition at line 1734 of file G4CutTubs.cc.

1735 {
1736  G4double xRand, yRand, zRand, phi, cosphi, sinphi, chose,
1737  aOne, aTwo, aThr, aFou;
1738  G4double rRand;
1739 
1740  aOne = 2.*fDz*fDPhi*fRMax;
1741  aTwo = 2.*fDz*fDPhi*fRMin;
1742  aThr = 0.5*fDPhi*(fRMax*fRMax-fRMin*fRMin);
1743  aFou = 2.*fDz*(fRMax-fRMin);
1744 
1745  phi = G4RandFlat::shoot(fSPhi, fSPhi+fDPhi);
1746  cosphi = std::cos(phi);
1747  sinphi = std::sin(phi);
1748 
1749  rRand = GetRadiusInRing(fRMin,fRMax);
1750 
1751  if( (fSPhi == 0) && (fDPhi == twopi) ) { aFou = 0; }
1752 
1753  chose = G4RandFlat::shoot(0.,aOne+aTwo+2.*aThr+2.*aFou);
1754 
1755  if( (chose >=0) && (chose < aOne) )
1756  {
1757  xRand = fRMax*cosphi;
1758  yRand = fRMax*sinphi;
1759  zRand = G4RandFlat::shoot(GetCutZ(G4ThreeVector(xRand,yRand,-fDz)),
1760  GetCutZ(G4ThreeVector(xRand,yRand,fDz)));
1761  return G4ThreeVector (xRand, yRand, zRand);
1762  }
1763  else if( (chose >= aOne) && (chose < aOne + aTwo) )
1764  {
1765  xRand = fRMin*cosphi;
1766  yRand = fRMin*sinphi;
1767  zRand = G4RandFlat::shoot(GetCutZ(G4ThreeVector(xRand,yRand,-fDz)),
1768  GetCutZ(G4ThreeVector(xRand,yRand,fDz)));
1769  return G4ThreeVector (xRand, yRand, zRand);
1770  }
1771  else if( (chose >= aOne + aTwo) && (chose < aOne + aTwo + aThr) )
1772  {
1773  xRand = rRand*cosphi;
1774  yRand = rRand*sinphi;
1775  zRand = GetCutZ(G4ThreeVector(xRand,yRand,fDz));
1776  return G4ThreeVector (xRand, yRand, zRand);
1777  }
1778  else if( (chose >= aOne + aTwo + aThr) && (chose < aOne + aTwo + 2.*aThr) )
1779  {
1780  xRand = rRand*cosphi;
1781  yRand = rRand*sinphi;
1782  zRand = GetCutZ(G4ThreeVector(xRand,yRand,-fDz));
1783  return G4ThreeVector (xRand, yRand, zRand);
1784  }
1785  else if( (chose >= aOne + aTwo + 2.*aThr)
1786  && (chose < aOne + aTwo + 2.*aThr + aFou) )
1787  {
1788  xRand = rRand*std::cos(fSPhi);
1789  yRand = rRand*std::sin(fSPhi);
1790  zRand = G4RandFlat::shoot(GetCutZ(G4ThreeVector(xRand,yRand,-fDz)),
1791  GetCutZ(G4ThreeVector(xRand,yRand,fDz)));
1792  return G4ThreeVector (xRand, yRand, zRand);
1793  }
1794  else
1795  {
1796  xRand = rRand*std::cos(fSPhi+fDPhi);
1797  yRand = rRand*std::sin(fSPhi+fDPhi);
1798  zRand = G4RandFlat::shoot(GetCutZ(G4ThreeVector(xRand,yRand,-fDz)),
1799  GetCutZ(G4ThreeVector(xRand,yRand,fDz)));
1800  return G4ThreeVector (xRand, yRand, zRand);
1801  }
1802 }
G4INCL::DeJongSpin::shoot
ThreeVector shoot(const G4int Ap, const G4int Af)
Definition: G4INCLDeJongSpin.cc:68
G4OTubs::fDPhi
G4double fDPhi
Definition: G4OTubs.hh:177
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:42
G4OTubs::fRMax
G4double fRMax
Definition: G4OTubs.hh:177
G4CSGSolid::GetRadiusInRing
G4double GetRadiusInRing(G4double rmin, G4double rmax) const
Definition: G4CSGSolid.cc:111
twopi
static constexpr double twopi
Definition: G4SIunits.hh:76
G4OTubs::fDz
G4double fDz
Definition: G4OTubs.hh:177
G4OTubs::fSPhi
G4double fSPhi
Definition: G4OTubs.hh:177
G4OTubs::fRMin
G4double fRMin
Definition: G4OTubs.hh:177
G4double
double G4double
Definition: G4Types.hh:76
G4CutTubs::GetCutZ
G4double GetCutZ(const G4ThreeVector &p) const
Definition: G4CutTubs.cc:1896

Here is the call graph for this function:

G4double G4CutTubs::GetSurfaceArea ( )
inlinevirtual

Reimplemented from G4VSolid.

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

Implements G4VSolid.

Definition at line 373 of file G4CutTubs.cc.

374 {
375  G4ThreeVector vZ = G4ThreeVector(0,0,fDz);
376  EInside in = kInside;
377 
378  // Check the lower cut plane
379  //
380  G4double zinLow =(p+vZ).dot(fLowNorm);
381  if (zinLow > halfCarTolerance) { return kOutside; }
382 
383  // Check the higher cut plane
384  //
385  G4double zinHigh = (p-vZ).dot(fHighNorm);
386  if (zinHigh > halfCarTolerance) { return kOutside; }
387 
388  // Check radius
389  //
390  G4double r2 = p.x()*p.x() + p.y()*p.y() ;
391 
392  G4double tolRMin = fRMin - halfRadTolerance;
393  G4double tolRMax = fRMax + halfRadTolerance;
394  if ( tolRMin < 0 ) { tolRMin = 0; }
395 
396  if (r2 < tolRMin*tolRMin || r2 > tolRMax*tolRMax) { return kOutside; }
397 
398  // Check Phi cut
399  //
400  if(!fPhiFullCutTube)
401  {
402  if ((tolRMin == 0) && (std::fabs(p.x()) <= halfCarTolerance)
403  && (std::fabs(p.y()) <= halfCarTolerance))
404  {
405  return kSurface;
406  }
407 
408  G4double phi0 = std::atan2(p.y(),p.x());
409  G4double phi1 = phi0 - twopi;
410  G4double phi2 = phi0 + twopi;
411 
412  in = kOutside;
413  G4double sphi = fSPhi - halfAngTolerance;
414  G4double ephi = sphi + fDPhi + kAngTolerance;
415  if ((phi0 >= sphi && phi0 <= ephi) ||
416  (phi1 >= sphi && phi1 <= ephi) ||
417  (phi2 >= sphi && phi2 <= ephi)) in = kSurface;
418  if (in == kOutside) { return kOutside; }
419 
420  sphi += kAngTolerance;
421  ephi -= kAngTolerance;
422  if ((phi0 >= sphi && phi0 <= ephi) ||
423  (phi1 >= sphi && phi1 <= ephi) ||
424  (phi2 >= sphi && phi2 <= ephi)) in = kInside;
425  if (in == kSurface) { return kSurface; }
426  }
427 
428  // Check on the Surface for Z
429  //
430  if ((zinLow >= -halfCarTolerance) || (zinHigh >= -halfCarTolerance))
431  {
432  return kSurface;
433  }
434 
435  // Check on the Surface for R
436  //
437  if (fRMin) { tolRMin = fRMin + halfRadTolerance; }
438  else { tolRMin = 0; }
439  tolRMax = fRMax - halfRadTolerance;
440  if (((r2 <= tolRMin*tolRMin) || (r2 >= tolRMax*tolRMax)) &&
441  (r2 >= halfRadTolerance*halfRadTolerance))
442  {
443  return kSurface;
444  }
445 
446  return in;
447 }
G4OTubs::fDPhi
G4double fDPhi
Definition: G4OTubs.hh:177
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:42
CLHEP::Hep3Vector::x
double x() const
G4OTubs::fRMax
G4double fRMax
Definition: G4OTubs.hh:177
G4OTubs::kAngTolerance
G4double kAngTolerance
Definition: G4OTubs.hh:173
twopi
static constexpr double twopi
Definition: G4SIunits.hh:76
EInside
EInside
Definition: geomdefs.hh:58
G4OTubs::fDz
G4double fDz
Definition: G4OTubs.hh:177
CLHEP::Hep3Vector::y
double y() const
G4OTubs::fSPhi
G4double fSPhi
Definition: G4OTubs.hh:177
G4OTubs::fRMin
G4double fRMin
Definition: G4OTubs.hh:177
G4double
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

G4bool G4CutTubs::IsCrossingCutPlanes ( ) const
protected

Definition at line 1873 of file G4CutTubs.cc.

1874 {
1875  G4double zXLow1,zXLow2,zYLow1,zYLow2;
1876  G4double zXHigh1,zXHigh2,zYHigh1,zYHigh2;
1877 
1878  zXLow1 = GetCutZ(G4ThreeVector(-fRMax, 0,-fDz));
1879  zXLow2 = GetCutZ(G4ThreeVector( fRMax, 0,-fDz));
1880  zYLow1 = GetCutZ(G4ThreeVector( 0,-fRMax,-fDz));
1881  zYLow2 = GetCutZ(G4ThreeVector( 0, fRMax,-fDz));
1882  zXHigh1 = GetCutZ(G4ThreeVector(-fRMax, 0, fDz));
1883  zXHigh2 = GetCutZ(G4ThreeVector( fRMax, 0, fDz));
1884  zYHigh1 = GetCutZ(G4ThreeVector( 0,-fRMax, fDz));
1885  zYHigh2 = GetCutZ(G4ThreeVector( 0, fRMax, fDz));
1886  if ( (zXLow1>zXHigh1) ||(zXLow2>zXHigh2)
1887  || (zYLow1>zYHigh1) ||(zYLow2>zYHigh2)) { return true; }
1888 
1889  return false;
1890 }
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:42
G4OTubs::fRMax
G4double fRMax
Definition: G4OTubs.hh:177
G4OTubs::fDz
G4double fDz
Definition: G4OTubs.hh:177
G4double
double G4double
Definition: G4Types.hh:76
G4CutTubs::GetCutZ
G4double GetCutZ(const G4ThreeVector &p) const
Definition: G4CutTubs.cc:1896

Here is the call graph for this function:

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

Definition at line 180 of file G4CutTubs.cc.

181 {
182  // Check assignment to self
183  //
184  if (this == &rhs) { return *this; }
185 
186  // Copy base class data
187  //
188  G4OTubs::operator=(rhs);
189 
190  // Copy data
191  //
192  fLowNorm = rhs.fLowNorm; fHighNorm = rhs.fHighNorm;
193  fPhiFullCutTube = rhs.fPhiFullCutTube;
194  halfCarTolerance = rhs.halfCarTolerance;
195  halfRadTolerance = rhs.halfRadTolerance;
196  halfAngTolerance = rhs.halfAngTolerance;
197 
198  return *this;
199 }
G4OTubs::operator=
G4OTubs & operator=(const G4OTubs &rhs)
Definition: G4OTubs.cc:140

Here is the call graph for this function:

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

Reimplemented from G4CSGSolid.

Definition at line 1709 of file G4CutTubs.cc.

1710 {
1711  G4int oldprc = os.precision(16);
1712  os << "-----------------------------------------------------------\n"
1713  << " *** Dump for solid - " << GetName() << " ***\n"
1714  << " ===================================================\n"
1715  << " Solid type: G4CutTubs\n"
1716  << " Parameters: \n"
1717  << " inner radius : " << fRMin/mm << " mm \n"
1718  << " outer radius : " << fRMax/mm << " mm \n"
1719  << " half length Z: " << fDz/mm << " mm \n"
1720  << " starting phi : " << fSPhi/degree << " degrees \n"
1721  << " delta phi : " << fDPhi/degree << " degrees \n"
1722  << " low Norm : " << fLowNorm << " \n"
1723  << " high Norm : " <<fHighNorm << " \n"
1724  << "-----------------------------------------------------------\n";
1725  os.precision(oldprc);
1726 
1727  return os;
1728 }
G4VSolid::GetName
G4String GetName() const
mm
static constexpr double mm
Definition: G4SIunits.hh:115
G4OTubs::fDPhi
G4double fDPhi
Definition: G4OTubs.hh:177
G4OTubs::fRMax
G4double fRMax
Definition: G4OTubs.hh:177
G4int
int G4int
Definition: G4Types.hh:78
degree
static constexpr double degree
Definition: G4SIunits.hh:144
G4OTubs::fDz
G4double fDz
Definition: G4OTubs.hh:177
G4OTubs::fSPhi
G4double fSPhi
Definition: G4OTubs.hh:177
G4OTubs::fRMin
G4double fRMin
Definition: G4OTubs.hh:177

Here is the call graph for this function:

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

Implements G4VSolid.

Definition at line 455 of file G4CutTubs.cc.

456 {
457  G4int noSurfaces = 0;
458  G4double rho, pPhi;
459  G4double distZLow,distZHigh, distRMin, distRMax;
460  G4double distSPhi = kInfinity, distEPhi = kInfinity;
461  G4ThreeVector vZ=G4ThreeVector(0,0,fDz);
462 
463  G4ThreeVector norm, sumnorm(0.,0.,0.);
464  G4ThreeVector nZ = G4ThreeVector(0, 0, 1.0);
465  G4ThreeVector nR, nPs, nPe;
466 
467  rho = std::sqrt(p.x()*p.x() + p.y()*p.y());
468 
469  distRMin = std::fabs(rho - fRMin);
470  distRMax = std::fabs(rho - fRMax);
471 
472  // dist to Low Cut
473  //
474  distZLow =std::fabs((p+vZ).dot(fLowNorm));
475 
476  // dist to High Cut
477  //
478  distZHigh = std::fabs((p-vZ).dot(fHighNorm));
479 
480  if (!fPhiFullCutTube) // Protected against (0,0,z)
481  {
482  if ( rho > halfCarTolerance )
483  {
484  pPhi = std::atan2(p.y(),p.x());
485 
486  if(pPhi < fSPhi- halfCarTolerance) { pPhi += twopi; }
487  else if(pPhi > fSPhi+fDPhi+ halfCarTolerance) { pPhi -= twopi; }
488 
489  distSPhi = std::fabs(pPhi - fSPhi);
490  distEPhi = std::fabs(pPhi - fSPhi - fDPhi);
491  }
492  else if( !fRMin )
493  {
494  distSPhi = 0.;
495  distEPhi = 0.;
496  }
497  nPs = G4ThreeVector(std::sin(fSPhi),-std::cos(fSPhi),0);
498  nPe = G4ThreeVector(-std::sin(fSPhi+fDPhi),std::cos(fSPhi+fDPhi),0);
499  }
500  if ( rho > halfCarTolerance ) { nR = G4ThreeVector(p.x()/rho,p.y()/rho,0); }
501 
502  if( distRMax <= halfCarTolerance )
503  {
504  noSurfaces ++;
505  sumnorm += nR;
506  }
507  if( fRMin && (distRMin <= halfCarTolerance) )
508  {
509  noSurfaces ++;
510  sumnorm -= nR;
511  }
512  if( fDPhi < twopi )
513  {
514  if (distSPhi <= halfAngTolerance)
515  {
516  noSurfaces ++;
517  sumnorm += nPs;
518  }
519  if (distEPhi <= halfAngTolerance)
520  {
521  noSurfaces ++;
522  sumnorm += nPe;
523  }
524  }
525  if (distZLow <= halfCarTolerance)
526  {
527  noSurfaces ++;
528  sumnorm += fLowNorm;
529  }
530  if (distZHigh <= halfCarTolerance)
531  {
532  noSurfaces ++;
533  sumnorm += fHighNorm;
534  }
535  if ( noSurfaces == 0 )
536  {
537 #ifdef G4CSGDEBUG
538  G4Exception("G4CutTubs::SurfaceNormal(p)", "GeomSolids1002",
539  JustWarning, "Point p is not on surface !?" );
540  G4int oldprc = G4cout.precision(20);
541  G4cout<< "G4CutTubs::SN ( "<<p.x()<<", "<<p.y()<<", "<<p.z()<<" ); "
542  << G4endl << G4endl;
543  G4cout.precision(oldprc) ;
544 #endif
545  norm = ApproxSurfaceNormal(p);
546  }
547  else if ( noSurfaces == 1 ) { norm = sumnorm; }
548  else { norm = sumnorm.unit(); }
549 
550  return norm;
551 }
G4OTubs::fDPhi
G4double fDPhi
Definition: G4OTubs.hh:177
kInfinity
static const G4double kInfinity
Definition: geomdefs.hh:42
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:42
CLHEP::Hep3Vector::x
double x() const
G4OTubs::fRMax
G4double fRMax
Definition: G4OTubs.hh:177
JustWarning
Definition: G4ExceptionSeverity.hh:64
G4int
int G4int
Definition: G4Types.hh:78
CLHEP::Hep3Vector::z
double z() const
twopi
static constexpr double twopi
Definition: G4SIunits.hh:76
G4cout
G4GLOB_DLL std::ostream G4cout
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4CutTubs::ApproxSurfaceNormal
G4ThreeVector ApproxSurfaceNormal(const G4ThreeVector &p) const
Definition: G4CutTubs.cc:558
CLHEP::Hep3Vector::unit
Hep3Vector unit() const
G4OTubs::fDz
G4double fDz
Definition: G4OTubs.hh:177
CLHEP::Hep3Vector::y
double y() const
G4OTubs::fSPhi
G4double fSPhi
Definition: G4OTubs.hh:177
G4endl
#define G4endl
Definition: G4ios.hh:61
G4OTubs::fRMin
G4double fRMin
Definition: G4OTubs.hh:177
G4double
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

The documentation for this class was generated from the following files: