Geant4  10.03.p01
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros 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
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 583 of file G4CutTubs.cc.

584 {
585  ENorm side ;
586  G4ThreeVector norm ;
587  G4double rho, phi ;
588  G4double distZLow,distZHigh,distZ;
589  G4double distRMin, distRMax, distSPhi, distEPhi, distMin ;
590  G4ThreeVector vZ=G4ThreeVector(0,0,fDz);
591 
592  rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
593 
594  distRMin = std::fabs(rho - fRMin) ;
595  distRMax = std::fabs(rho - fRMax) ;
596 
597  //dist to Low Cut
598  //
599  distZLow =std::fabs((p+vZ).dot(fLowNorm));
600 
601  //dist to High Cut
602  //
603  distZHigh = std::fabs((p-vZ).dot(fHighNorm));
604  distZ=std::min(distZLow,distZHigh);
605 
606  if (distRMin < distRMax) // First minimum
607  {
608  if ( distZ < distRMin )
609  {
610  distMin = distZ ;
611  side = kNZ ;
612  }
613  else
614  {
615  distMin = distRMin ;
616  side = kNRMin ;
617  }
618  }
619  else
620  {
621  if ( distZ < distRMax )
622  {
623  distMin = distZ ;
624  side = kNZ ;
625  }
626  else
627  {
628  distMin = distRMax ;
629  side = kNRMax ;
630  }
631  }
632  if (!fPhiFullCutTube && rho ) // Protected against (0,0,z)
633  {
634  phi = std::atan2(p.y(),p.x()) ;
635 
636  if ( phi < 0 ) { phi += twopi; }
637 
638  if ( fSPhi < 0 )
639  {
640  distSPhi = std::fabs(phi - (fSPhi + twopi))*rho ;
641  }
642  else
643  {
644  distSPhi = std::fabs(phi - fSPhi)*rho ;
645  }
646  distEPhi = std::fabs(phi - fSPhi - fDPhi)*rho ;
647 
648  if (distSPhi < distEPhi) // Find new minimum
649  {
650  if ( distSPhi < distMin )
651  {
652  side = kNSPhi ;
653  }
654  }
655  else
656  {
657  if ( distEPhi < distMin )
658  {
659  side = kNEPhi ;
660  }
661  }
662  }
663  switch ( side )
664  {
665  case kNRMin : // Inner radius
666  {
667  norm = G4ThreeVector(-p.x()/rho, -p.y()/rho, 0) ;
668  break ;
669  }
670  case kNRMax : // Outer radius
671  {
672  norm = G4ThreeVector(p.x()/rho, p.y()/rho, 0) ;
673  break ;
674  }
675  case kNZ : // + or - dz
676  {
677  if ( distZHigh > distZLow ) { norm = fHighNorm ; }
678  else { norm = fLowNorm; }
679  break ;
680  }
681  case kNSPhi:
682  {
683  norm = G4ThreeVector(std::sin(fSPhi), -std::cos(fSPhi), 0) ;
684  break ;
685  }
686  case kNEPhi:
687  {
688  norm = G4ThreeVector(-std::sin(fSPhi+fDPhi), std::cos(fSPhi+fDPhi), 0) ;
689  break;
690  }
691  default: // Should never reach this case ...
692  {
693  DumpInfo();
694  G4Exception("G4CutTubs::ApproxSurfaceNormal()",
695  "GeomSolids1002", JustWarning,
696  "Undefined side for valid surface normal to solid.");
697  break ;
698  }
699  }
700  return norm;
701 }
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 = (360/NSTEPS)*deg; // 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

Here is the call graph for this function:

G4VSolid * G4CutTubs::Clone ( ) const
virtual

Reimplemented from G4VSolid.

Definition at line 1725 of file G4CutTubs.cc.

1726 {
1727  return new G4CutTubs(*this);
1728 }
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 1838 of file G4CutTubs.cc.

1839 {
1840  typedef G4double G4double3[3];
1841  typedef G4int G4int4[4];
1842 
1843  G4Polyhedron *ph = new G4Polyhedron;
1844  G4Polyhedron *ph1 = G4OTubs::CreatePolyhedron();
1845  G4int nn=ph1->GetNoVertices();
1846  G4int nf=ph1->GetNoFacets();
1847  G4double3* xyz = new G4double3[nn]; // number of nodes
1848  G4int4* faces = new G4int4[nf] ; // number of faces
1849 
1850  for(G4int i=0;i<nn;i++)
1851  {
1852  xyz[i][0]=ph1->GetVertex(i+1).x();
1853  xyz[i][1]=ph1->GetVertex(i+1).y();
1854  G4double tmpZ=ph1->GetVertex(i+1).z();
1855  if(tmpZ>=fDz-kCarTolerance)
1856  {
1857  xyz[i][2]=GetCutZ(G4ThreeVector(xyz[i][0],xyz[i][1],fDz));
1858  }
1859  else if(tmpZ<=-fDz+kCarTolerance)
1860  {
1861  xyz[i][2]=GetCutZ(G4ThreeVector(xyz[i][0],xyz[i][1],-fDz));
1862  }
1863  else
1864  {
1865  xyz[i][2]=tmpZ;
1866  }
1867  }
1868  G4int iNodes[4];
1869  G4int *iEdge=0;
1870  G4int n;
1871  for(G4int i=0;i<nf;i++)
1872  {
1873  ph1->GetFacet(i+1,n,iNodes,iEdge);
1874  for(G4int k=0;k<n;k++)
1875  {
1876  faces[i][k]=iNodes[k];
1877  }
1878  for(G4int k=n;k<4;k++)
1879  {
1880  faces[i][k]=0;
1881  }
1882  }
1883  ph->createPolyhedron(nn,nf,xyz,faces);
1884 
1885  delete [] xyz;
1886  delete [] faces;
1887  delete ph1;
1888 
1889  return ph;
1890 }
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
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:1921

Here is the call graph for this function:

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

Implements G4VSolid.

Definition at line 1833 of file G4CutTubs.cc.

1834 {
1835  scene.AddSolid (*this) ;
1836 }
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 725 of file G4CutTubs.cc.

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

Here is the call graph for this function:

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

Implements G4VSolid.

Definition at line 1178 of file G4CutTubs.cc.

1179 {
1180  G4double safRMin,safRMax,safZLow,safZHigh,safePhi,safe,rho,cosPsi;
1181  G4ThreeVector vZ=G4ThreeVector(0,0,fDz);
1182 
1183  // Distance to R
1184  //
1185  rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
1186 
1187  safRMin = fRMin- rho ;
1188  safRMax = rho - fRMax ;
1189 
1190  // Distances to ZCut(Low/High)
1191 
1192  // Dist to Low Cut
1193  //
1194  safZLow = (p+vZ).dot(fLowNorm);
1195 
1196  // Dist to High Cut
1197  //
1198  safZHigh = (p-vZ).dot(fHighNorm);
1199 
1200  safe = std::max(safZLow,safZHigh);
1201 
1202  if ( safRMin > safe ) { safe = safRMin; }
1203  if ( safRMax> safe ) { safe = safRMax; }
1204 
1205  // Distance to Phi
1206  //
1207  if ( (!fPhiFullCutTube) && (rho) )
1208  {
1209  // Psi=angle from central phi to point
1210  //
1211  cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/rho ;
1212 
1213  if ( cosPsi < std::cos(fDPhi*0.5) )
1214  {
1215  // Point lies outside phi range
1216 
1217  if ( (p.y()*cosCPhi - p.x()*sinCPhi) <= 0 )
1218  {
1219  safePhi = std::fabs(p.x()*sinSPhi - p.y()*cosSPhi) ;
1220  }
1221  else
1222  {
1223  safePhi = std::fabs(p.x()*sinEPhi - p.y()*cosEPhi) ;
1224  }
1225  if ( safePhi > safe ) { safe = safePhi; }
1226  }
1227  }
1228  if ( safe < 0 ) { safe = 0; }
1229 
1230  return safe ;
1231 }
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 1238 of file G4CutTubs.cc.

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

1670 {
1671  G4double safRMin,safRMax,safZLow,safZHigh,safePhi,safe,rho;
1672  G4ThreeVector vZ=G4ThreeVector(0,0,fDz);
1673 
1674  rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ; // Distance to R
1675 
1676  safRMin = rho - fRMin ;
1677  safRMax = fRMax - rho ;
1678 
1679  // Distances to ZCut(Low/High)
1680 
1681  // Dist to Low Cut
1682  //
1683  safZLow = std::fabs((p+vZ).dot(fLowNorm));
1684 
1685  // Dist to High Cut
1686  //
1687  safZHigh = std::fabs((p-vZ).dot(fHighNorm));
1688  safe = std::min(safZLow,safZHigh);
1689 
1690  if ( safRMin < safe ) { safe = safRMin; }
1691  if ( safRMax< safe ) { safe = safRMax; }
1692 
1693  // Check if phi divided, Calc distances closest phi plane
1694  //
1695  if ( !fPhiFullCutTube )
1696  {
1697  if ( p.y()*cosCPhi-p.x()*sinCPhi <= 0 )
1698  {
1699  safePhi = -(p.x()*sinSPhi - p.y()*cosSPhi) ;
1700  }
1701  else
1702  {
1703  safePhi = (p.x()*sinEPhi - p.y()*cosEPhi) ;
1704  }
1705  if (safePhi < safe) { safe = safePhi ; }
1706  }
1707  if ( safe < 0 ) { safe = 0; }
1708 
1709  return safe ;
1710 }
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 1921 of file G4CutTubs.cc.

1922 {
1923  G4double newz = p.z(); // p.z() should be either +fDz or -fDz
1924  if (p.z()<0)
1925  {
1926  if(fLowNorm.z()!=0.)
1927  {
1928  newz = -fDz-(p.x()*fLowNorm.x()+p.y()*fLowNorm.y())/fLowNorm.z();
1929  }
1930  }
1931  else
1932  {
1933  if(fHighNorm.z()!=0.)
1934  {
1935  newz = fDz-(p.x()*fHighNorm.x()+p.y()*fHighNorm.y())/fHighNorm.z();
1936  }
1937  }
1938  return newz;
1939 }
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 1716 of file G4CutTubs.cc.

1717 {
1718  return G4String("G4CutTubs");
1719 }
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 1945 of file G4CutTubs.cc.

1947 {
1948  G4double phiLow = std::atan2(fLowNorm.y(),fLowNorm.x());
1949  G4double phiHigh= std::atan2(fHighNorm.y(),fHighNorm.x());
1950 
1951  G4double xc=0, yc=0,z1;
1952  G4double z[8];
1953  G4bool in_range_low = false;
1954  G4bool in_range_hi = false;
1955 
1956  G4int i;
1957  for (i=0; i<2; i++)
1958  {
1959  if (phiLow<0) { phiLow+=twopi; }
1960  G4double ddp = phiLow-fSPhi;
1961  if (ddp<0) { ddp += twopi; }
1962  if (ddp <= fDPhi)
1963  {
1964  xc = fRMin*std::cos(phiLow);
1965  yc = fRMin*std::sin(phiLow);
1966  z1 = GetCutZ(G4ThreeVector(xc, yc, -fDz));
1967  xc = fRMax*std::cos(phiLow);
1968  yc = fRMax*std::sin(phiLow);
1969  z1 = std::min(z1, GetCutZ(G4ThreeVector(xc, yc, -fDz)));
1970  if (in_range_low) { zmin = std::min(zmin, z1); }
1971  else { zmin = z1; }
1972  in_range_low = true;
1973  }
1974  phiLow += pi;
1975  if (phiLow>twopi) { phiLow-=twopi; }
1976  }
1977  for (i=0; i<2; i++)
1978  {
1979  if (phiHigh<0) { phiHigh+=twopi; }
1980  G4double ddp = phiHigh-fSPhi;
1981  if (ddp<0) { ddp += twopi; }
1982  if (ddp <= fDPhi)
1983  {
1984  xc = fRMin*std::cos(phiHigh);
1985  yc = fRMin*std::sin(phiHigh);
1986  z1 = GetCutZ(G4ThreeVector(xc, yc, fDz));
1987  xc = fRMax*std::cos(phiHigh);
1988  yc = fRMax*std::sin(phiHigh);
1989  z1 = std::min(z1, GetCutZ(G4ThreeVector(xc, yc, fDz)));
1990  if (in_range_hi) { zmax = std::min(zmax, z1); }
1991  else { zmax = z1; }
1992  in_range_hi = true;
1993  }
1994  phiHigh += pi;
1995  if (phiHigh>twopi) { phiHigh-=twopi; }
1996  }
1997 
1998  xc = fRMin*std::cos(fSPhi);
1999  yc = fRMin*std::sin(fSPhi);
2000  z[0] = GetCutZ(G4ThreeVector(xc, yc, -fDz));
2001  z[4] = GetCutZ(G4ThreeVector(xc, yc, fDz));
2002 
2003  xc = fRMin*std::cos(fSPhi+fDPhi);
2004  yc = fRMin*std::sin(fSPhi+fDPhi);
2005  z[1] = GetCutZ(G4ThreeVector(xc, yc, -fDz));
2006  z[5] = GetCutZ(G4ThreeVector(xc, yc, fDz));
2007 
2008  xc = fRMax*std::cos(fSPhi);
2009  yc = fRMax*std::sin(fSPhi);
2010  z[2] = GetCutZ(G4ThreeVector(xc, yc, -fDz));
2011  z[6] = GetCutZ(G4ThreeVector(xc, yc, fDz));
2012 
2013  xc = fRMax*std::cos(fSPhi+fDPhi);
2014  yc = fRMax*std::sin(fSPhi+fDPhi);
2015  z[3] = GetCutZ(G4ThreeVector(xc, yc, -fDz));
2016  z[7] = GetCutZ(G4ThreeVector(xc, yc, fDz));
2017 
2018  // Find min/max
2019 
2020  z1=z[0];
2021  for (i = 1; i < 4; i++)
2022  {
2023  if(z[i] < z[i-1])z1=z[i];
2024  }
2025 
2026  if (in_range_low)
2027  {
2028  zmin = std::min(zmin, z1);
2029  }
2030  else
2031  {
2032  zmin = z1;
2033  }
2034  z1=z[4];
2035  for (i = 1; i < 4; i++)
2036  {
2037  if(z[4+i] > z[4+i-1]) { z1=z[4+i]; }
2038  }
2039 
2040  if (in_range_hi) { zmax = std::max(zmax, z1); }
2041  else { zmax = z1; }
2042 }
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
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:1921

Here is the call graph for this function:

G4ThreeVector G4CutTubs::GetPointOnSurface ( ) const
virtual

Reimplemented from G4VSolid.

Definition at line 1759 of file G4CutTubs.cc.

1760 {
1761  G4double xRand, yRand, zRand, phi, cosphi, sinphi, chose,
1762  aOne, aTwo, aThr, aFou;
1763  G4double rRand;
1764 
1765  aOne = 2.*fDz*fDPhi*fRMax;
1766  aTwo = 2.*fDz*fDPhi*fRMin;
1767  aThr = 0.5*fDPhi*(fRMax*fRMax-fRMin*fRMin);
1768  aFou = 2.*fDz*(fRMax-fRMin);
1769 
1770  phi = G4RandFlat::shoot(fSPhi, fSPhi+fDPhi);
1771  cosphi = std::cos(phi);
1772  sinphi = std::sin(phi);
1773 
1774  rRand = GetRadiusInRing(fRMin,fRMax);
1775 
1776  if( (fSPhi == 0) && (fDPhi == twopi) ) { aFou = 0; }
1777 
1778  chose = G4RandFlat::shoot(0.,aOne+aTwo+2.*aThr+2.*aFou);
1779 
1780  if( (chose >=0) && (chose < aOne) )
1781  {
1782  xRand = fRMax*cosphi;
1783  yRand = fRMax*sinphi;
1784  zRand = G4RandFlat::shoot(GetCutZ(G4ThreeVector(xRand,yRand,-fDz)),
1785  GetCutZ(G4ThreeVector(xRand,yRand,fDz)));
1786  return G4ThreeVector (xRand, yRand, zRand);
1787  }
1788  else if( (chose >= aOne) && (chose < aOne + aTwo) )
1789  {
1790  xRand = fRMin*cosphi;
1791  yRand = fRMin*sinphi;
1792  zRand = G4RandFlat::shoot(GetCutZ(G4ThreeVector(xRand,yRand,-fDz)),
1793  GetCutZ(G4ThreeVector(xRand,yRand,fDz)));
1794  return G4ThreeVector (xRand, yRand, zRand);
1795  }
1796  else if( (chose >= aOne + aTwo) && (chose < aOne + aTwo + aThr) )
1797  {
1798  xRand = rRand*cosphi;
1799  yRand = rRand*sinphi;
1800  zRand = GetCutZ(G4ThreeVector(xRand,yRand,fDz));
1801  return G4ThreeVector (xRand, yRand, zRand);
1802  }
1803  else if( (chose >= aOne + aTwo + aThr) && (chose < aOne + aTwo + 2.*aThr) )
1804  {
1805  xRand = rRand*cosphi;
1806  yRand = rRand*sinphi;
1807  zRand = GetCutZ(G4ThreeVector(xRand,yRand,-fDz));
1808  return G4ThreeVector (xRand, yRand, zRand);
1809  }
1810  else if( (chose >= aOne + aTwo + 2.*aThr)
1811  && (chose < aOne + aTwo + 2.*aThr + aFou) )
1812  {
1813  xRand = rRand*std::cos(fSPhi);
1814  yRand = rRand*std::sin(fSPhi);
1815  zRand = G4RandFlat::shoot(GetCutZ(G4ThreeVector(xRand,yRand,-fDz)),
1816  GetCutZ(G4ThreeVector(xRand,yRand,fDz)));
1817  return G4ThreeVector (xRand, yRand, zRand);
1818  }
1819  else
1820  {
1821  xRand = rRand*std::cos(fSPhi+fDPhi);
1822  yRand = rRand*std::sin(fSPhi+fDPhi);
1823  zRand = G4RandFlat::shoot(GetCutZ(G4ThreeVector(xRand,yRand,-fDz)),
1824  GetCutZ(G4ThreeVector(xRand,yRand,fDz)));
1825  return G4ThreeVector (xRand, yRand, zRand);
1826  }
1827 }
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:1921

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  G4double zinLow,zinHigh,r2,pPhi=0.;
376  G4double tolRMin,tolRMax;
377  G4ThreeVector vZ=G4ThreeVector(0,0,fDz);
378  EInside in = kInside;
379 
380  // Check if point is contained in the cut plane in -/+ Z
381 
382  // Check the lower cut plane
383  //
384  zinLow =(p+vZ).dot(fLowNorm);
385  if (zinLow>halfCarTolerance) { return kOutside; }
386 
387  // Check the higher cut plane
388  //
389  zinHigh = (p-vZ).dot(fHighNorm);
390  if (zinHigh>halfCarTolerance) { return kOutside; }
391 
392  // Check radius
393  //
394  r2 = p.x()*p.x() + p.y()*p.y() ;
395 
396  // First check 'generous' boundaries R+tolerance
397  //
398  tolRMin = fRMin - halfRadTolerance ;
399  tolRMax = fRMax + halfRadTolerance ;
400  if ( tolRMin < 0 ) { tolRMin = 0; }
401 
402  if ( ((r2 < tolRMin*tolRMin) || (r2 > tolRMax*tolRMax))
403  && (r2 >=halfRadTolerance*halfRadTolerance) ) { return kOutside; }
404 
405  // Check Phi
406  //
407  if(!fPhiFullCutTube)
408  {
409  // Try outer tolerant phi boundaries only
410 
411  if ( (tolRMin==0) && (std::fabs(p.x())<=halfCarTolerance)
412  && (std::fabs(p.y())<=halfCarTolerance) )
413  {
414  return kSurface;
415  }
416 
417  pPhi = std::atan2(p.y(),p.x()) ;
418 
419  if ( pPhi < -halfAngTolerance) { pPhi += twopi; } // 0<=pPhi<2pi
420  if ( fSPhi >= 0 )
421  {
422  if ( (std::fabs(pPhi) < halfAngTolerance)
423  && (std::fabs(fSPhi + fDPhi - twopi) < halfAngTolerance) )
424  {
425  pPhi += twopi ; // 0 <= pPhi < 2pi
426  }
427  if ( (pPhi <= fSPhi - halfAngTolerance)
428  || (pPhi >= fSPhi + fDPhi + halfAngTolerance) )
429  {
430  in = kOutside ;
431  }
432  else if ( (pPhi <= fSPhi + halfAngTolerance)
433  || (pPhi >= fSPhi + fDPhi - halfAngTolerance) )
434  {
435  in=kSurface;
436  }
437  }
438  else // fSPhi < 0
439  {
440  if ( (pPhi <= fSPhi + twopi - halfAngTolerance)
441  && (pPhi >= fSPhi + fDPhi + halfAngTolerance) )
442  {
443  in = kOutside;
444  }
445  else
446  {
447  in = kSurface ;
448  }
449  }
450  }
451 
452  // Check on the Surface for Z
453  //
454  if ((zinLow>=-halfCarTolerance)
455  || (zinHigh>=-halfCarTolerance))
456  {
457  in=kSurface;
458  }
459 
460  // Check on the Surface for R
461  //
462  if (fRMin) { tolRMin = fRMin + halfRadTolerance ; }
463  else { tolRMin = 0 ; }
464  tolRMax = fRMax - halfRadTolerance ;
465  if ( ((r2 <= tolRMin*tolRMin) || (r2 >= tolRMax*tolRMax))&&
466  (r2 >=halfRadTolerance*halfRadTolerance) )
467  {
468  return kSurface;
469  }
470 
471  return in;
472 }
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
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 1898 of file G4CutTubs.cc.

1899 {
1900  G4double zXLow1,zXLow2,zYLow1,zYLow2;
1901  G4double zXHigh1,zXHigh2,zYHigh1,zYHigh2;
1902 
1903  zXLow1 = GetCutZ(G4ThreeVector(-fRMax, 0,-fDz));
1904  zXLow2 = GetCutZ(G4ThreeVector( fRMax, 0,-fDz));
1905  zYLow1 = GetCutZ(G4ThreeVector( 0,-fRMax,-fDz));
1906  zYLow2 = GetCutZ(G4ThreeVector( 0, fRMax,-fDz));
1907  zXHigh1 = GetCutZ(G4ThreeVector(-fRMax, 0, fDz));
1908  zXHigh2 = GetCutZ(G4ThreeVector( fRMax, 0, fDz));
1909  zYHigh1 = GetCutZ(G4ThreeVector( 0,-fRMax, fDz));
1910  zYHigh2 = GetCutZ(G4ThreeVector( 0, fRMax, fDz));
1911  if ( (zXLow1>zXHigh1) ||(zXLow2>zXHigh2)
1912  || (zYLow1>zYHigh1) ||(zYLow2>zYHigh2)) { return true; }
1913 
1914  return false;
1915 }
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:1921

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

1735 {
1736  G4int oldprc = os.precision(16);
1737  os << "-----------------------------------------------------------\n"
1738  << " *** Dump for solid - " << GetName() << " ***\n"
1739  << " ===================================================\n"
1740  << " Solid type: G4CutTubs\n"
1741  << " Parameters: \n"
1742  << " inner radius : " << fRMin/mm << " mm \n"
1743  << " outer radius : " << fRMax/mm << " mm \n"
1744  << " half length Z: " << fDz/mm << " mm \n"
1745  << " starting phi : " << fSPhi/degree << " degrees \n"
1746  << " delta phi : " << fDPhi/degree << " degrees \n"
1747  << " low Norm : " << fLowNorm << " \n"
1748  << " high Norm : " <<fHighNorm << " \n"
1749  << "-----------------------------------------------------------\n";
1750  os.precision(oldprc);
1751 
1752  return os;
1753 }
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 480 of file G4CutTubs.cc.

481 {
482  G4int noSurfaces = 0;
483  G4double rho, pPhi;
484  G4double distZLow,distZHigh, distRMin, distRMax;
485  G4double distSPhi = kInfinity, distEPhi = kInfinity;
486  G4ThreeVector vZ=G4ThreeVector(0,0,fDz);
487 
488  G4ThreeVector norm, sumnorm(0.,0.,0.);
489  G4ThreeVector nZ = G4ThreeVector(0, 0, 1.0);
490  G4ThreeVector nR, nPs, nPe;
491 
492  rho = std::sqrt(p.x()*p.x() + p.y()*p.y());
493 
494  distRMin = std::fabs(rho - fRMin);
495  distRMax = std::fabs(rho - fRMax);
496 
497  // dist to Low Cut
498  //
499  distZLow =std::fabs((p+vZ).dot(fLowNorm));
500 
501  // dist to High Cut
502  //
503  distZHigh = std::fabs((p-vZ).dot(fHighNorm));
504 
505  if (!fPhiFullCutTube) // Protected against (0,0,z)
506  {
507  if ( rho > halfCarTolerance )
508  {
509  pPhi = std::atan2(p.y(),p.x());
510 
511  if(pPhi < fSPhi- halfCarTolerance) { pPhi += twopi; }
512  else if(pPhi > fSPhi+fDPhi+ halfCarTolerance) { pPhi -= twopi; }
513 
514  distSPhi = std::fabs(pPhi - fSPhi);
515  distEPhi = std::fabs(pPhi - fSPhi - fDPhi);
516  }
517  else if( !fRMin )
518  {
519  distSPhi = 0.;
520  distEPhi = 0.;
521  }
522  nPs = G4ThreeVector(std::sin(fSPhi),-std::cos(fSPhi),0);
523  nPe = G4ThreeVector(-std::sin(fSPhi+fDPhi),std::cos(fSPhi+fDPhi),0);
524  }
525  if ( rho > halfCarTolerance ) { nR = G4ThreeVector(p.x()/rho,p.y()/rho,0); }
526 
527  if( distRMax <= halfCarTolerance )
528  {
529  noSurfaces ++;
530  sumnorm += nR;
531  }
532  if( fRMin && (distRMin <= halfCarTolerance) )
533  {
534  noSurfaces ++;
535  sumnorm -= nR;
536  }
537  if( fDPhi < twopi )
538  {
539  if (distSPhi <= halfAngTolerance)
540  {
541  noSurfaces ++;
542  sumnorm += nPs;
543  }
544  if (distEPhi <= halfAngTolerance)
545  {
546  noSurfaces ++;
547  sumnorm += nPe;
548  }
549  }
550  if (distZLow <= halfCarTolerance)
551  {
552  noSurfaces ++;
553  sumnorm += fLowNorm;
554  }
555  if (distZHigh <= halfCarTolerance)
556  {
557  noSurfaces ++;
558  sumnorm += fHighNorm;
559  }
560  if ( noSurfaces == 0 )
561  {
562 #ifdef G4CSGDEBUG
563  G4Exception("G4CutTubs::SurfaceNormal(p)", "GeomSolids1002",
564  JustWarning, "Point p is not on surface !?" );
565  G4int oldprc = G4cout.precision(20);
566  G4cout<< "G4CutTubs::SN ( "<<p.x()<<", "<<p.y()<<", "<<p.z()<<" ); "
567  << G4endl << G4endl;
568  G4cout.precision(oldprc) ;
569 #endif
570  norm = ApproxSurfaceNormal(p);
571  }
572  else if ( noSurfaces == 1 ) { norm = sumnorm; }
573  else { norm = sumnorm.unit(); }
574 
575  return norm;
576 }
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:583
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: