Geant4.10
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Public Member Functions | Protected Member Functions
G4CutTubs Class Reference

#include <G4CutTubs.hh>

Inheritance diagram for G4CutTubs:
G4OTubs G4CSGSolid G4VSolid

Public Member Functions

 G4CutTubs (const G4String &pName, G4double pRMin, G4double pRMax, G4double pDz, G4double pSPhi, G4double pDPhi, G4ThreeVector pLowNorm, G4ThreeVector pHighNorm)
 
 ~G4CutTubs ()
 
G4ThreeVector GetLowNorm () const
 
G4ThreeVector GetHighNorm () const
 
G4double GetCubicVolume ()
 
G4double GetSurfaceArea ()
 
G4bool CalculateExtent (const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pmin, G4double &pmax) const
 
EInside Inside (const G4ThreeVector &p) const
 
G4ThreeVector SurfaceNormal (const G4ThreeVector &p) const
 
G4double DistanceToIn (const G4ThreeVector &p, const G4ThreeVector &v) const
 
G4double DistanceToIn (const G4ThreeVector &p) const
 
G4double DistanceToOut (const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=G4bool(false), G4bool *validNorm=0, G4ThreeVector *n=0) const
 
G4double DistanceToOut (const G4ThreeVector &p) const
 
G4GeometryType GetEntityType () const
 
G4ThreeVector GetPointOnSurface () const
 
G4VSolid * Clone () const
 
std::ostream & StreamInfo (std::ostream &os) const
 
void DescribeYourselfTo (G4VGraphicsScene &scene) const
 
G4Polyhedron * CreatePolyhedron () const
 
 G4CutTubs (__void__ &)
 
 G4CutTubs (const G4CutTubs &rhs)
 
G4CutTubs & operator= (const G4CutTubs &rhs)
 
- Public Member Functions inherited from 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
 
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 ()
 
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)
 

Protected Member Functions

G4ThreeVectorList * CreateRotatedVertices (const G4AffineTransform &pTransform) const
 
G4ThreeVector ApproxSurfaceNormal (const G4ThreeVector &p) const
 
G4bool IsCrossingCutPlanes () const
 
G4double GetCutZ (const G4ThreeVector &p) const
 
void GetMaxMinZ (G4double &zmin, G4double &zmax) const
 
- Protected Member Functions inherited from G4OTubs
G4ThreeVectorList * CreateRotatedVertices (const G4AffineTransform &pTransform) const
 
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
 
G4double EstimateCubicVolume (G4int nStat, G4double epsilon) const
 
G4double EstimateSurfaceArea (G4int nStat, G4double ell) 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
 
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 61 of file G4CutTubs.cc.

References FatalException, G4endl, G4Exception(), G4GeometryTolerance::GetAngularTolerance(), G4GeometryTolerance::GetInstance(), G4VSolid::GetName(), G4GeometryTolerance::GetRadialTolerance(), IsCrossingCutPlanes(), JustWarning, G4OTubs::kAngTolerance, G4VSolid::kCarTolerance, G4OTubs::kRadTolerance, CLHEP::Hep3Vector::mag2(), CLHEP::Hep3Vector::setZ(), python.hepunit::twopi, CLHEP::Hep3Vector::unit(), CLHEP::Hep3Vector::x(), CLHEP::Hep3Vector::y(), and CLHEP::Hep3Vector::z().

Referenced by Clone().

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

Definition at line 152 of file G4CutTubs.cc.

153 {
154 }
G4CutTubs::G4CutTubs ( __void__ &  a)

Definition at line 141 of file G4CutTubs.cc.

142  : G4OTubs(a), fLowNorm(G4ThreeVector()),
143  fHighNorm(G4ThreeVector()), fPhiFullCutTube(false),
144  halfCarTolerance(0.), halfRadTolerance(0.), halfAngTolerance(0.)
145 {
146 }
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:60
G4CutTubs::G4CutTubs ( const G4CutTubs &  rhs)

Definition at line 160 of file G4CutTubs.cc.

161  : G4OTubs(rhs), fLowNorm(rhs.fLowNorm), fHighNorm(rhs.fHighNorm),
162  fPhiFullCutTube(rhs.fPhiFullCutTube),
163  halfCarTolerance(rhs.halfCarTolerance),
164  halfRadTolerance(rhs.halfRadTolerance),
165  halfAngTolerance(rhs.halfAngTolerance)
166 {
167 }
G4OTubs::G4OTubs
G4OTubs(const G4String &pName, G4double pRMin, G4double pRMax, G4double pDz, G4double pSPhi, G4double pDPhi)
Definition: G4OTubs.cc:60

Member Function Documentation

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

Reimplemented from G4OTubs.

Definition at line 621 of file G4CutTubs.cc.

References G4VSolid::DumpInfo(), G4OTubs::fDPhi, G4OTubs::fDz, G4OTubs::fRMax, G4OTubs::fRMin, G4OTubs::fSPhi, G4Exception(), JustWarning, G4OTubs::kNEPhi, G4OTubs::kNRMax, G4OTubs::kNRMin, G4OTubs::kNSPhi, G4OTubs::kNZ, G4INCL::Math::min(), python.hepunit::twopi, CLHEP::Hep3Vector::x(), and CLHEP::Hep3Vector::y().

Referenced by SurfaceNormal().

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

Implements G4VSolid.

Definition at line 198 of file G4CutTubs.cc.

References G4VSolid::ClipBetweenSections(), G4VSolid::ClipCrossSection(), CreateRotatedVertices(), G4OTubs::fDPhi, G4OTubs::fRMax, G4OTubs::fRMin, G4VoxelLimits::GetMaxExtent(), GetMaxMinZ(), G4VoxelLimits::GetMaxXExtent(), G4VoxelLimits::GetMaxYExtent(), G4VoxelLimits::GetMaxZExtent(), G4VoxelLimits::GetMinExtent(), G4VoxelLimits::GetMinXExtent(), G4VoxelLimits::GetMinYExtent(), G4VoxelLimits::GetMinZExtent(), Inside(), G4AffineTransform::Inverse(), G4AffineTransform::IsRotated(), G4VoxelLimits::IsXLimited(), G4VoxelLimits::IsYLimited(), G4VoxelLimits::IsZLimited(), G4VSolid::kCarTolerance, kOutside, kXAxis, kYAxis, kZAxis, G4AffineTransform::NetTranslation(), G4AffineTransform::TransformPoint(), python.hepunit::twopi, CLHEP::Hep3Vector::x(), CLHEP::Hep3Vector::y(), and CLHEP::Hep3Vector::z().

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

Reimplemented from G4VSolid.

Definition at line 1855 of file G4CutTubs.cc.

References G4CutTubs().

1856 {
1857  return new G4CutTubs(*this);
1858 }
G4CutTubs::G4CutTubs
G4CutTubs(const G4String &pName, G4double pRMin, G4double pRMax, G4double pDz, G4double pSPhi, G4double pDPhi, G4ThreeVector pLowNorm, G4ThreeVector pHighNorm)
Definition: G4CutTubs.cc:61
G4Polyhedron * G4CutTubs::CreatePolyhedron ( ) const
virtual

Reimplemented from G4VSolid.

Definition at line 1968 of file G4CutTubs.cc.

References G4OTubs::CreatePolyhedron(), HepPolyhedron::createPolyhedron(), G4OTubs::fDz, GetCutZ(), HepPolyhedron::GetFacet(), HepPolyhedron::GetNoFacets(), HepPolyhedron::GetNoVertices(), HepPolyhedron::GetVertex(), G4VSolid::kCarTolerance, n, G4InuclParticleNames::nn, HepGeom::BasicVector3D< T >::x(), HepGeom::BasicVector3D< T >::y(), and HepGeom::BasicVector3D< T >::z().

1969 {
1970  typedef G4double G4double3[3];
1971  typedef G4int G4int4[4];
1972 
1973  G4Polyhedron *ph = new G4Polyhedron;
1974  G4Polyhedron *ph1 = G4OTubs::CreatePolyhedron();
1975  G4int nn=ph1->GetNoVertices();
1976  G4int nf=ph1->GetNoFacets();
1977  G4double3* xyz = new G4double3[nn]; // number of nodes
1978  G4int4* faces = new G4int4[nf] ; // number of faces
1979 
1980  for(G4int i=0;i<nn;i++)
1981  {
1982  xyz[i][0]=ph1->GetVertex(i+1).x();
1983  xyz[i][1]=ph1->GetVertex(i+1).y();
1984  G4double tmpZ=ph1->GetVertex(i+1).z();
1985  if(tmpZ>=fDz-kCarTolerance)
1986  {
1987  xyz[i][2]=GetCutZ(G4ThreeVector(xyz[i][0],xyz[i][1],fDz));
1988  }
1989  else if(tmpZ<=-fDz+kCarTolerance)
1990  {
1991  xyz[i][2]=GetCutZ(G4ThreeVector(xyz[i][0],xyz[i][1],-fDz));
1992  }
1993  else
1994  {
1995  xyz[i][2]=tmpZ;
1996  }
1997  }
1998  G4int iNodes[4];
1999  G4int *iEdge=0;
2000  G4int n;
2001  for(G4int i=0;i<nf;i++)
2002  {
2003  ph1->GetFacet(i+1,n,iNodes,iEdge);
2004  for(G4int k=0;k<n;k++)
2005  {
2006  faces[i][k]=iNodes[k];
2007  }
2008  for(G4int k=n;k<4;k++)
2009  {
2010  faces[i][k]=0;
2011  }
2012  }
2013  ph->createPolyhedron(nn,nf,xyz,faces);
2014 
2015  delete [] xyz;
2016  delete [] faces;
2017  delete ph1;
2018 
2019  return ph;
2020 }
HepPolyhedron::GetFacet
void GetFacet(G4int iFace, G4int &n, G4int *iNodes, G4int *edgeFlags=0, G4int *iFaces=0) const
Definition: HepPolyhedron.cc:1065
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:1358
G4int
int G4int
Definition: G4Types.hh:78
n
const G4int n
Definition: G4UrQMD1_3Interface.hh:144
G4OTubs::CreatePolyhedron
G4Polyhedron * CreatePolyhedron() const
Definition: G4OTubs.cc:1883
HepPolyhedron::GetNoVertices
G4int GetNoVertices() const
Definition: HepPolyhedron.h:252
G4OTubs::fDz
G4double fDz
Definition: G4OTubs.hh:176
G4VSolid::kCarTolerance
G4double kCarTolerance
Definition: G4VSolid.hh:305
G4double
double G4double
Definition: G4Types.hh:76
HepPolyhedron::GetNoFacets
G4int GetNoFacets() const
Definition: HepPolyhedron.h:255
HepPolyhedron::GetVertex
G4Point3D GetVertex(G4int index) const
Definition: HepPolyhedron.cc:885
G4CutTubs::GetCutZ
G4double GetCutZ(const G4ThreeVector &p) const
Definition: G4CutTubs.cc:2051
G4ThreeVectorList * G4CutTubs::CreateRotatedVertices ( const G4AffineTransform &  pTransform) const
protected

Definition at line 1760 of file G4CutTubs.cc.

References G4VSolid::DumpInfo(), FatalException, G4OTubs::fDPhi, G4OTubs::fDz, G4OTubs::fRMax, G4OTubs::fRMin, G4OTubs::fSPhi, G4Exception(), GetCutZ(), G4VSolid::kCarTolerance, kMaxMeshSections, kMeshAngleDefault, kMinMeshSections, and G4AffineTransform::TransformPoint().

Referenced by CalculateExtent().

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

Implements G4VSolid.

Definition at line 1963 of file G4CutTubs.cc.

References G4VGraphicsScene::AddSolid().

1964 {
1965  scene.AddSolid (*this) ;
1966 }
G4VGraphicsScene::AddSolid
virtual void AddSolid(const G4Box &)=0
G4double G4CutTubs::DistanceToIn ( const G4ThreeVector &  p,
const G4ThreeVector &  v 
) const
virtual

Implements G4VSolid.

Definition at line 763 of file G4CutTubs.cc.

References test::b, test::c, G4OTubs::cosCPhi, G4OTubs::cosEPhi, G4OTubs::cosHDPhiIT, G4OTubs::cosSPhi, CLHEP::Hep3Vector::dot(), G4OTubs::fDz, G4OTubs::fRMax, G4OTubs::fRMin, GetCutZ(), G4OTubs::kRadTolerance, G4OTubs::sinCPhi, G4OTubs::sinEPhi, G4OTubs::sinSPhi, plottest35::t1, CLHEP::Hep3Vector::x(), CLHEP::Hep3Vector::y(), and CLHEP::Hep3Vector::z().

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

Implements G4VSolid.

Definition at line 1216 of file G4CutTubs.cc.

References G4OTubs::cosCPhi, G4OTubs::cosEPhi, G4OTubs::cosSPhi, G4OTubs::fDPhi, G4OTubs::fDz, G4OTubs::fRMax, G4OTubs::fRMin, G4INCL::Math::max(), G4OTubs::sinCPhi, G4OTubs::sinEPhi, G4OTubs::sinSPhi, CLHEP::Hep3Vector::x(), and CLHEP::Hep3Vector::y().

1217 {
1218  G4double safRMin,safRMax,safZLow,safZHigh,safePhi,safe,rho,cosPsi;
1219  G4ThreeVector vZ=G4ThreeVector(0,0,fDz);
1220 
1221  // Distance to R
1222  //
1223  rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
1224 
1225  safRMin = fRMin- rho ;
1226  safRMax = rho - fRMax ;
1227 
1228  // Distances to ZCut(Low/High)
1229 
1230  // Dist to Low Cut
1231  //
1232  safZLow = (p+vZ).dot(fLowNorm);
1233 
1234  // Dist to High Cut
1235  //
1236  safZHigh = (p-vZ).dot(fHighNorm);
1237 
1238  safe = std::max(safZLow,safZHigh);
1239 
1240  if ( safRMin > safe ) { safe = safRMin; }
1241  if ( safRMax> safe ) { safe = safRMax; }
1242 
1243  // Distance to Phi
1244  //
1245  if ( (!fPhiFullCutTube) && (rho) )
1246  {
1247  // Psi=angle from central phi to point
1248  //
1249  cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/rho ;
1250 
1251  if ( cosPsi < std::cos(fDPhi*0.5) )
1252  {
1253  // Point lies outside phi range
1254 
1255  if ( (p.y()*cosCPhi - p.x()*sinCPhi) <= 0 )
1256  {
1257  safePhi = std::fabs(p.x()*sinSPhi - p.y()*cosSPhi) ;
1258  }
1259  else
1260  {
1261  safePhi = std::fabs(p.x()*sinEPhi - p.y()*cosEPhi) ;
1262  }
1263  if ( safePhi > safe ) { safe = safePhi; }
1264  }
1265  }
1266  if ( safe < 0 ) { safe = 0; }
1267 
1268  return safe ;
1269 }
G4OTubs::fDPhi
G4double fDPhi
Definition: G4OTubs.hh:176
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:42
CLHEP::Hep3Vector::x
double x() const
G4OTubs::cosSPhi
G4double cosSPhi
Definition: G4OTubs.hh:180
G4OTubs::fRMax
G4double fRMax
Definition: G4OTubs.hh:176
G4OTubs::cosEPhi
G4double cosEPhi
Definition: G4OTubs.hh:180
G4OTubs::sinSPhi
G4double sinSPhi
Definition: G4OTubs.hh:180
G4OTubs::cosCPhi
G4double cosCPhi
Definition: G4OTubs.hh:180
G4OTubs::sinEPhi
G4double sinEPhi
Definition: G4OTubs.hh:180
G4INCL::Math::max
T max(const T t1, const T t2)
brief Return the largest of the two arguments
Definition: G4INCLGlobals.hh:116
G4OTubs::fDz
G4double fDz
Definition: G4OTubs.hh:176
CLHEP::Hep3Vector::y
double y() const
G4OTubs::fRMin
G4double fRMin
Definition: G4OTubs.hh:176
G4double
double G4double
Definition: G4Types.hh:76
G4OTubs::sinCPhi
G4double sinCPhi
Definition: G4OTubs.hh:180
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 1276 of file G4CutTubs.cc.

References test::b, test::c, G4OTubs::cosCPhi, G4OTubs::cosEPhi, G4OTubs::cosSPhi, CLHEP::Hep3Vector::dot(), G4VSolid::DumpInfo(), G4OTubs::fDPhi, G4OTubs::fDz, G4OTubs::fRMax, G4OTubs::fRMin, G4OTubs::fSPhi, G4cout, G4endl, G4Exception(), JustWarning, G4VSolid::kCarTolerance, G4OTubs::kEPhi, G4OTubs::kMZ, G4OTubs::kNull, G4OTubs::kPZ, G4OTubs::kRadTolerance, G4OTubs::kRMax, G4OTubs::kRMin, G4OTubs::kSPhi, python.hepunit::mm, python.hepunit::pi, G4OTubs::sinCPhi, G4OTubs::sinEPhi, G4OTubs::sinSPhi, plottest35::t1, python.hepunit::twopi, CLHEP::Hep3Vector::x(), CLHEP::Hep3Vector::y(), and CLHEP::Hep3Vector::z().

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

Implements G4VSolid.

Definition at line 1707 of file G4CutTubs.cc.

References G4OTubs::cosCPhi, G4OTubs::cosEPhi, G4OTubs::cosSPhi, G4OTubs::fDz, G4OTubs::fRMax, G4OTubs::fRMin, G4INCL::Math::min(), G4OTubs::sinCPhi, G4OTubs::sinEPhi, G4OTubs::sinSPhi, CLHEP::Hep3Vector::x(), and CLHEP::Hep3Vector::y().

1708 {
1709  G4double safRMin,safRMax,safZLow,safZHigh,safePhi,safe,rho;
1710  G4ThreeVector vZ=G4ThreeVector(0,0,fDz);
1711 
1712  rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ; // Distance to R
1713 
1714  safRMin = rho - fRMin ;
1715  safRMax = fRMax - rho ;
1716 
1717  // Distances to ZCut(Low/High)
1718 
1719  // Dist to Low Cut
1720  //
1721  safZLow = std::fabs((p+vZ).dot(fLowNorm));
1722 
1723  // Dist to High Cut
1724  //
1725  safZHigh = std::fabs((p-vZ).dot(fHighNorm));
1726  safe = std::min(safZLow,safZHigh);
1727 
1728  if ( safRMin < safe ) { safe = safRMin; }
1729  if ( safRMax< safe ) { safe = safRMax; }
1730 
1731  // Check if phi divided, Calc distances closest phi plane
1732  //
1733  if ( !fPhiFullCutTube )
1734  {
1735  if ( p.y()*cosCPhi-p.x()*sinCPhi <= 0 )
1736  {
1737  safePhi = -(p.x()*sinSPhi - p.y()*cosSPhi) ;
1738  }
1739  else
1740  {
1741  safePhi = (p.x()*sinEPhi - p.y()*cosEPhi) ;
1742  }
1743  if (safePhi < safe) { safe = safePhi ; }
1744  }
1745  if ( safe < 0 ) { safe = 0; }
1746 
1747  return safe ;
1748 }
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:42
CLHEP::Hep3Vector::x
double x() const
G4OTubs::cosSPhi
G4double cosSPhi
Definition: G4OTubs.hh:180
G4OTubs::fRMax
G4double fRMax
Definition: G4OTubs.hh:176
G4OTubs::cosEPhi
G4double cosEPhi
Definition: G4OTubs.hh:180
G4OTubs::sinSPhi
G4double sinSPhi
Definition: G4OTubs.hh:180
G4OTubs::cosCPhi
G4double cosCPhi
Definition: G4OTubs.hh:180
G4OTubs::sinEPhi
G4double sinEPhi
Definition: G4OTubs.hh:180
G4OTubs::fDz
G4double fDz
Definition: G4OTubs.hh:176
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:121
G4OTubs::fRMin
G4double fRMin
Definition: G4OTubs.hh:176
G4double
double G4double
Definition: G4Types.hh:76
G4OTubs::sinCPhi
G4double sinCPhi
Definition: G4OTubs.hh:180
G4double G4CutTubs::GetCubicVolume ( )
inlinevirtual

Reimplemented from G4VSolid.

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

Definition at line 2051 of file G4CutTubs.cc.

References G4OTubs::fDz, CLHEP::Hep3Vector::x(), CLHEP::Hep3Vector::y(), and CLHEP::Hep3Vector::z().

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

2052 {
2053  G4double newz = p.z(); // p.z() should be either +fDz or -fDz
2054  if (p.z()<0)
2055  {
2056  if(fLowNorm.z()!=0.)
2057  {
2058  newz = -fDz-(p.x()*fLowNorm.x()+p.y()*fLowNorm.y())/fLowNorm.z();
2059  }
2060  }
2061  else
2062  {
2063  if(fHighNorm.z()!=0.)
2064  {
2065  newz = fDz-(p.x()*fHighNorm.x()+p.y()*fHighNorm.y())/fHighNorm.z();
2066  }
2067  }
2068  return newz;
2069 }
CLHEP::Hep3Vector::x
double x() const
CLHEP::Hep3Vector::z
double z() const
G4OTubs::fDz
G4double fDz
Definition: G4OTubs.hh:176
CLHEP::Hep3Vector::y
double y() const
G4double
double G4double
Definition: G4Types.hh:76
G4GeometryType G4CutTubs::GetEntityType ( ) const
virtual

Implements G4VSolid.

Definition at line 1846 of file G4CutTubs.cc.

1847 {
1848  return G4String("G4CutTubs");
1849 }
G4ThreeVector G4CutTubs::GetHighNorm ( ) const
inline

Referenced by G4GDMLWriteSolids::CutTubeWrite().

G4ThreeVector G4CutTubs::GetLowNorm ( ) const
inline

Referenced by G4GDMLWriteSolids::CutTubeWrite().

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

Definition at line 2075 of file G4CutTubs.cc.

References G4OTubs::fDPhi, G4OTubs::fDz, G4OTubs::fRMax, G4OTubs::fRMin, G4OTubs::fSPhi, GetCutZ(), G4INCL::Math::max(), G4INCL::Math::min(), python.hepunit::pi, python.hepunit::twopi, CLHEP::Hep3Vector::x(), CLHEP::Hep3Vector::y(), and z.

Referenced by CalculateExtent().

2077 {
2078  G4double phiLow = std::atan2(fLowNorm.y(),fLowNorm.x());
2079  G4double phiHigh= std::atan2(fHighNorm.y(),fHighNorm.x());
2080 
2081  G4double xc=0, yc=0,z1;
2082  G4double z[8];
2083  G4bool in_range_low = false;
2084  G4bool in_range_hi = false;
2085 
2086  G4int i;
2087  for (i=0; i<2; i++)
2088  {
2089  if (phiLow<0) { phiLow+=twopi; }
2090  G4double ddp = phiLow-fSPhi;
2091  if (ddp<0) { ddp += twopi; }
2092  if (ddp <= fDPhi)
2093  {
2094  xc = fRMin*std::cos(phiLow);
2095  yc = fRMin*std::sin(phiLow);
2096  z1 = GetCutZ(G4ThreeVector(xc, yc, -fDz));
2097  xc = fRMax*std::cos(phiLow);
2098  yc = fRMax*std::sin(phiLow);
2099  z1 = std::min(z1, GetCutZ(G4ThreeVector(xc, yc, -fDz)));
2100  if (in_range_low) { zmin = std::min(zmin, z1); }
2101  else { zmin = z1; }
2102  in_range_low = true;
2103  }
2104  phiLow += pi;
2105  if (phiLow>twopi) { phiLow-=twopi; }
2106  }
2107  for (i=0; i<2; i++)
2108  {
2109  if (phiHigh<0) { phiHigh+=twopi; }
2110  G4double ddp = phiHigh-fSPhi;
2111  if (ddp<0) { ddp += twopi; }
2112  if (ddp <= fDPhi)
2113  {
2114  xc = fRMin*std::cos(phiHigh);
2115  yc = fRMin*std::sin(phiHigh);
2116  z1 = GetCutZ(G4ThreeVector(xc, yc, fDz));
2117  xc = fRMax*std::cos(phiHigh);
2118  yc = fRMax*std::sin(phiHigh);
2119  z1 = std::min(z1, GetCutZ(G4ThreeVector(xc, yc, fDz)));
2120  if (in_range_hi) { zmax = std::min(zmax, z1); }
2121  else { zmax = z1; }
2122  in_range_hi = true;
2123  }
2124  phiHigh += pi;
2125  if (phiLow>twopi) { phiHigh-=twopi; }
2126  }
2127 
2128  xc = fRMin*std::cos(fSPhi);
2129  yc = fRMin*std::sin(fSPhi);
2130  z[0] = GetCutZ(G4ThreeVector(xc, yc, -fDz));
2131  z[4] = GetCutZ(G4ThreeVector(xc, yc, fDz));
2132 
2133  xc = fRMin*std::cos(fSPhi+fDPhi);
2134  yc = fRMin*std::sin(fSPhi+fDPhi);
2135  z[1] = GetCutZ(G4ThreeVector(xc, yc, -fDz));
2136  z[5] = GetCutZ(G4ThreeVector(xc, yc, fDz));
2137 
2138  xc = fRMax*std::cos(fSPhi);
2139  yc = fRMax*std::sin(fSPhi);
2140  z[2] = GetCutZ(G4ThreeVector(xc, yc, -fDz));
2141  z[6] = GetCutZ(G4ThreeVector(xc, yc, fDz));
2142 
2143  xc = fRMax*std::cos(fSPhi+fDPhi);
2144  yc = fRMax*std::sin(fSPhi+fDPhi);
2145  z[3] = GetCutZ(G4ThreeVector(xc, yc, -fDz));
2146  z[7] = GetCutZ(G4ThreeVector(xc, yc, fDz));
2147 
2148  // Find min/max
2149 
2150  z1=z[0];
2151  for (i = 1; i < 4; i++)
2152  {
2153  if(z[i] < z[i-1])z1=z[i];
2154  }
2155 
2156  if (in_range_low)
2157  {
2158  zmin = std::min(zmin, z1);
2159  }
2160  else
2161  {
2162  zmin = z1;
2163  }
2164  z1=z[4];
2165  for (i = 1; i < 4; i++)
2166  {
2167  if(z[4+i] > z[4+i-1]) { z1=z[4+i]; }
2168  }
2169 
2170  if (in_range_hi) { zmax = std::max(zmax, z1); }
2171  else { zmax = z1; }
2172 }
G4OTubs::fDPhi
G4double fDPhi
Definition: G4OTubs.hh:176
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:42
CLHEP::Hep3Vector::x
double x() const
G4OTubs::fRMax
G4double fRMax
Definition: G4OTubs.hh:176
z
G4double z
Definition: TRTMaterials.hh:39
G4int
int G4int
Definition: G4Types.hh:78
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:116
G4OTubs::fDz
G4double fDz
Definition: G4OTubs.hh:176
CLHEP::Hep3Vector::y
double y() const
G4OTubs::fSPhi
G4double fSPhi
Definition: G4OTubs.hh:176
G4INCL::Math::min
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
Definition: G4INCLGlobals.hh:121
G4OTubs::fRMin
G4double fRMin
Definition: G4OTubs.hh:176
G4double
double G4double
Definition: G4Types.hh:76
G4CutTubs::GetCutZ
G4double GetCutZ(const G4ThreeVector &p) const
Definition: G4CutTubs.cc:2051
G4ThreeVector G4CutTubs::GetPointOnSurface ( ) const
virtual

Reimplemented from G4VSolid.

Definition at line 1889 of file G4CutTubs.cc.

References G4OTubs::fDPhi, G4OTubs::fDz, G4OTubs::fRMax, G4OTubs::fRMin, G4OTubs::fSPhi, GetCutZ(), G4CSGSolid::GetRadiusInRing(), G4INCL::DeJongSpin::shoot(), and python.hepunit::twopi.

1890 {
1891  G4double xRand, yRand, zRand, phi, cosphi, sinphi, chose,
1892  aOne, aTwo, aThr, aFou;
1893  G4double rRand;
1894 
1895  aOne = 2.*fDz*fDPhi*fRMax;
1896  aTwo = 2.*fDz*fDPhi*fRMin;
1897  aThr = 0.5*fDPhi*(fRMax*fRMax-fRMin*fRMin);
1898  aFou = 2.*fDz*(fRMax-fRMin);
1899 
1900  phi = RandFlat::shoot(fSPhi, fSPhi+fDPhi);
1901  cosphi = std::cos(phi);
1902  sinphi = std::sin(phi);
1903 
1904  rRand = GetRadiusInRing(fRMin,fRMax);
1905 
1906  if( (fSPhi == 0) && (fDPhi == twopi) ) { aFou = 0; }
1907 
1908  chose = RandFlat::shoot(0.,aOne+aTwo+2.*aThr+2.*aFou);
1909 
1910  if( (chose >=0) && (chose < aOne) )
1911  {
1912  xRand = fRMax*cosphi;
1913  yRand = fRMax*sinphi;
1914  zRand = RandFlat::shoot(GetCutZ(G4ThreeVector(xRand,yRand,-fDz)),
1915  GetCutZ(G4ThreeVector(xRand,yRand,fDz)));
1916  return G4ThreeVector (xRand, yRand, zRand);
1917  }
1918  else if( (chose >= aOne) && (chose < aOne + aTwo) )
1919  {
1920  xRand = fRMin*cosphi;
1921  yRand = fRMin*sinphi;
1922  zRand = RandFlat::shoot(GetCutZ(G4ThreeVector(xRand,yRand,-fDz)),
1923  GetCutZ(G4ThreeVector(xRand,yRand,fDz)));
1924  return G4ThreeVector (xRand, yRand, zRand);
1925  }
1926  else if( (chose >= aOne + aTwo) && (chose < aOne + aTwo + aThr) )
1927  {
1928  xRand = rRand*cosphi;
1929  yRand = rRand*sinphi;
1930  zRand = GetCutZ(G4ThreeVector(xRand,yRand,fDz));
1931  return G4ThreeVector (xRand, yRand, zRand);
1932  }
1933  else if( (chose >= aOne + aTwo + aThr) && (chose < aOne + aTwo + 2.*aThr) )
1934  {
1935  xRand = rRand*cosphi;
1936  yRand = rRand*sinphi;
1937  zRand = GetCutZ(G4ThreeVector(xRand,yRand,-fDz));
1938  return G4ThreeVector (xRand, yRand, zRand);
1939  }
1940  else if( (chose >= aOne + aTwo + 2.*aThr)
1941  && (chose < aOne + aTwo + 2.*aThr + aFou) )
1942  {
1943  xRand = rRand*std::cos(fSPhi);
1944  yRand = rRand*std::sin(fSPhi);
1945  zRand = RandFlat::shoot(GetCutZ(G4ThreeVector(xRand,yRand,-fDz)),
1946  GetCutZ(G4ThreeVector(xRand,yRand,fDz)));
1947  return G4ThreeVector (xRand, yRand, zRand);
1948  }
1949  else
1950  {
1951  xRand = rRand*std::cos(fSPhi+fDPhi);
1952  yRand = rRand*std::sin(fSPhi+fDPhi);
1953  zRand = RandFlat::shoot(GetCutZ(G4ThreeVector(xRand,yRand,-fDz)),
1954  GetCutZ(G4ThreeVector(xRand,yRand,fDz)));
1955  return G4ThreeVector (xRand, yRand, zRand);
1956  }
1957 }
G4INCL::DeJongSpin::shoot
ThreeVector shoot(const G4int Ap, const G4int Af)
Definition: G4INCLDeJongSpin.cc:67
G4OTubs::fDPhi
G4double fDPhi
Definition: G4OTubs.hh:176
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:42
G4OTubs::fRMax
G4double fRMax
Definition: G4OTubs.hh:176
G4CSGSolid::GetRadiusInRing
G4double GetRadiusInRing(G4double rmin, G4double rmax) const
Definition: G4CSGSolid.cc:101
G4OTubs::fDz
G4double fDz
Definition: G4OTubs.hh:176
G4OTubs::fSPhi
G4double fSPhi
Definition: G4OTubs.hh:176
G4OTubs::fRMin
G4double fRMin
Definition: G4OTubs.hh:176
G4double
double G4double
Definition: G4Types.hh:76
G4CutTubs::GetCutZ
G4double GetCutZ(const G4ThreeVector &p) const
Definition: G4CutTubs.cc:2051
G4double G4CutTubs::GetSurfaceArea ( )
inlinevirtual

Reimplemented from G4VSolid.

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

Implements G4VSolid.

Definition at line 411 of file G4CutTubs.cc.

References G4OTubs::fDPhi, G4OTubs::fDz, G4OTubs::fRMax, G4OTubs::fRMin, G4OTubs::fSPhi, kInside, kOutside, kSurface, python.hepunit::twopi, CLHEP::Hep3Vector::x(), and CLHEP::Hep3Vector::y().

Referenced by CalculateExtent().

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

Definition at line 2028 of file G4CutTubs.cc.

References G4OTubs::fDz, G4OTubs::fRMax, and GetCutZ().

Referenced by G4CutTubs().

2029 {
2030  G4double zXLow1,zXLow2,zYLow1,zYLow2;
2031  G4double zXHigh1,zXHigh2,zYHigh1,zYHigh2;
2032 
2033  zXLow1 = GetCutZ(G4ThreeVector(-fRMax, 0,-fDz));
2034  zXLow2 = GetCutZ(G4ThreeVector( fRMax, 0,-fDz));
2035  zYLow1 = GetCutZ(G4ThreeVector( 0,-fRMax,-fDz));
2036  zYLow2 = GetCutZ(G4ThreeVector( 0, fRMax,-fDz));
2037  zXHigh1 = GetCutZ(G4ThreeVector(-fRMax, 0, fDz));
2038  zXHigh2 = GetCutZ(G4ThreeVector( fRMax, 0, fDz));
2039  zYHigh1 = GetCutZ(G4ThreeVector( 0,-fRMax, fDz));
2040  zYHigh2 = GetCutZ(G4ThreeVector( 0, fRMax, fDz));
2041  if ( (zXLow1>zXHigh1) ||(zXLow2>zXHigh2)
2042  || (zYLow1>zYHigh1) ||(zYLow2>zYHigh2)) { return true; }
2043 
2044  return false;
2045 }
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:42
G4OTubs::fRMax
G4double fRMax
Definition: G4OTubs.hh:176
G4OTubs::fDz
G4double fDz
Definition: G4OTubs.hh:176
G4double
double G4double
Definition: G4Types.hh:76
G4CutTubs::GetCutZ
G4double GetCutZ(const G4ThreeVector &p) const
Definition: G4CutTubs.cc:2051
G4CutTubs & G4CutTubs::operator= ( const G4CutTubs &  rhs)

Definition at line 173 of file G4CutTubs.cc.

References G4OTubs::operator=().

174 {
175  // Check assignment to self
176  //
177  if (this == &rhs) { return *this; }
178 
179  // Copy base class data
180  //
181  G4OTubs::operator=(rhs);
182 
183  // Copy data
184  //
185  fLowNorm = rhs.fLowNorm; fHighNorm = rhs.fHighNorm;
186  fPhiFullCutTube = rhs.fPhiFullCutTube;
187  halfCarTolerance = rhs.halfCarTolerance;
188  halfRadTolerance = rhs.halfRadTolerance;
189  halfAngTolerance = rhs.halfAngTolerance;
190 
191  return *this;
192 }
G4OTubs::operator=
G4OTubs & operator=(const G4OTubs &rhs)
Definition: G4OTubs.cc:141
std::ostream & G4CutTubs::StreamInfo ( std::ostream &  os) const
virtual

Reimplemented from G4CSGSolid.

Definition at line 1864 of file G4CutTubs.cc.

References python.hepunit::degree, G4OTubs::fDPhi, G4OTubs::fDz, G4OTubs::fRMax, G4OTubs::fRMin, G4OTubs::fSPhi, G4VSolid::GetName(), and python.hepunit::mm.

1865 {
1866  G4int oldprc = os.precision(16);
1867  os << "-----------------------------------------------------------\n"
1868  << " *** Dump for solid - " << GetName() << " ***\n"
1869  << " ===================================================\n"
1870  << " Solid type: G4CutTubs\n"
1871  << " Parameters: \n"
1872  << " inner radius : " << fRMin/mm << " mm \n"
1873  << " outer radius : " << fRMax/mm << " mm \n"
1874  << " half length Z: " << fDz/mm << " mm \n"
1875  << " starting phi : " << fSPhi/degree << " degrees \n"
1876  << " delta phi : " << fDPhi/degree << " degrees \n"
1877  << " low Norm : " << fLowNorm << " \n"
1878  << " high Norm : " <<fHighNorm << " \n"
1879  << "-----------------------------------------------------------\n";
1880  os.precision(oldprc);
1881 
1882  return os;
1883 }
G4VSolid::GetName
G4String GetName() const
G4OTubs::fDPhi
G4double fDPhi
Definition: G4OTubs.hh:176
G4OTubs::fRMax
G4double fRMax
Definition: G4OTubs.hh:176
G4int
int G4int
Definition: G4Types.hh:78
python.hepunit.degree
tuple degree
Definition: hepunit.py:69
G4OTubs::fDz
G4double fDz
Definition: G4OTubs.hh:176
G4OTubs::fSPhi
G4double fSPhi
Definition: G4OTubs.hh:176
G4OTubs::fRMin
G4double fRMin
Definition: G4OTubs.hh:176
G4ThreeVector G4CutTubs::SurfaceNormal ( const G4ThreeVector &  p) const
virtual

Implements G4VSolid.

Definition at line 518 of file G4CutTubs.cc.

References ApproxSurfaceNormal(), G4OTubs::fDPhi, G4OTubs::fDz, G4OTubs::fRMax, G4OTubs::fRMin, G4OTubs::fSPhi, G4cout, G4endl, G4Exception(), JustWarning, python.hepunit::twopi, CLHEP::Hep3Vector::unit(), CLHEP::Hep3Vector::x(), CLHEP::Hep3Vector::y(), and CLHEP::Hep3Vector::z().

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