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

#include <G4Cons.hh>

Inheritance diagram for G4Cons:
G4CSGSolid G4VSolid

Public Member Functions

 G4Cons (const G4String &pName, G4double pRmin1, G4double pRmax1, G4double pRmin2, G4double pRmax2, G4double pDz, G4double pSPhi, G4double pDPhi)
 
 ~G4Cons ()
 
G4double GetInnerRadiusMinusZ () const
 
G4double GetOuterRadiusMinusZ () const
 
G4double GetInnerRadiusPlusZ () const
 
G4double GetOuterRadiusPlusZ () const
 
G4double GetZHalfLength () const
 
G4double GetStartPhiAngle () const
 
G4double GetDeltaPhiAngle () const
 
void SetInnerRadiusMinusZ (G4double Rmin1)
 
void SetOuterRadiusMinusZ (G4double Rmax1)
 
void SetInnerRadiusPlusZ (G4double Rmin2)
 
void SetOuterRadiusPlusZ (G4double Rmax2)
 
void SetZHalfLength (G4double newDz)
 
void SetStartPhiAngle (G4double newSPhi, G4bool trig=true)
 
void SetDeltaPhiAngle (G4double newDPhi)
 
G4double GetCubicVolume ()
 
G4double GetSurfaceArea ()
 
void ComputeDimensions (G4VPVParameterisation *p, const G4int n, const G4VPhysicalVolume *pRep)
 
G4bool CalculateExtent (const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pmin, G4double &pmax) const
 
EInside Inside (const G4ThreeVector &p) const
 
G4ThreeVector SurfaceNormal (const G4ThreeVector &p) const
 
G4double DistanceToIn (const G4ThreeVector &p, const G4ThreeVector &v) const
 
G4double DistanceToIn (const G4ThreeVector &p) const
 
G4double DistanceToOut (const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=G4bool(false), G4bool *validNorm=0, G4ThreeVector *n=0) const
 
G4double DistanceToOut (const G4ThreeVector &p) const
 
G4GeometryType GetEntityType () const
 
G4ThreeVector GetPointOnSurface () const
 
G4VSolid * Clone () const
 
std::ostream & StreamInfo (std::ostream &os) const
 
void DescribeYourselfTo (G4VGraphicsScene &scene) const
 
G4Polyhedron * CreatePolyhedron () const
 
 G4Cons (__void__ &)
 
 G4Cons (const G4Cons &rhs)
 
G4Cons & operator= (const G4Cons &rhs)
 
G4double GetRmin1 () const
 
G4double GetRmax1 () const
 
G4double GetRmin2 () const
 
G4double GetRmax2 () 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
 
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)
 

Additional Inherited Members

- 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
 
- Protected Attributes inherited from G4CSGSolid
G4double fCubicVolume
 
G4double fSurfaceArea
 
G4Polyhedron * fpPolyhedron
 
- Protected Attributes inherited from G4VSolid
G4double kCarTolerance
 

Detailed Description

Definition at line 82 of file G4Cons.hh.

Constructor & Destructor Documentation

G4Cons::G4Cons ( const G4String &  pName,
G4double  pRmin1,
G4double  pRmax1,
G4double  pRmin2,
G4double  pRmax2,
G4double  pDz,
G4double  pSPhi,
G4double  pDPhi 
)

Definition at line 79 of file G4Cons.cc.

References FatalException, G4endl, G4Exception(), G4GeometryTolerance::GetAngularTolerance(), G4GeometryTolerance::GetInstance(), G4VSolid::GetName(), G4GeometryTolerance::GetRadialTolerance(), and G4VSolid::kCarTolerance.

Referenced by Clone().

84  : G4CSGSolid(pName), fRmin1(pRmin1), fRmin2(pRmin2),
85  fRmax1(pRmax1), fRmax2(pRmax2), fDz(pDz), fSPhi(0.), fDPhi(0.)
86 {
87  kRadTolerance = G4GeometryTolerance::GetInstance()->GetRadialTolerance();
88  kAngTolerance = G4GeometryTolerance::GetInstance()->GetAngularTolerance();
89 
90  halfCarTolerance=kCarTolerance*0.5;
91  halfRadTolerance=kRadTolerance*0.5;
92  halfAngTolerance=kAngTolerance*0.5;
93 
94  // Check z-len
95  //
96  if ( pDz < 0 )
97  {
98  std::ostringstream message;
99  message << "Invalid Z half-length for Solid: " << GetName() << G4endl
100  << " hZ = " << pDz;
101  G4Exception("G4Cons::G4Cons()", "GeomSolids0002",
102  FatalException, message);
103  }
104 
105  // Check radii
106  //
107  if (((pRmin1>=pRmax1) || (pRmin2>=pRmax2) || (pRmin1<0)) && (pRmin2<0))
108  {
109  std::ostringstream message;
110  message << "Invalid values of radii for Solid: " << GetName() << G4endl
111  << " pRmin1 = " << pRmin1 << ", pRmin2 = " << pRmin2
112  << ", pRmax1 = " << pRmax1 << ", pRmax2 = " << pRmax2;
113  G4Exception("G4Cons::G4Cons()", "GeomSolids0002",
114  FatalException, message) ;
115  }
116  if( (pRmin1 == 0.0) && (pRmin2 > 0.0) ) { fRmin1 = 1e3*kRadTolerance ; }
117  if( (pRmin2 == 0.0) && (pRmin1 > 0.0) ) { fRmin2 = 1e3*kRadTolerance ; }
118 
119  // Check angles
120  //
121  CheckPhiAngles(pSPhi, pDPhi);
122 }
G4VSolid::GetName
G4String GetName() const
G4GeometryTolerance::GetRadialTolerance
G4double GetRadialTolerance() const
Definition: G4GeometryTolerance.cc:97
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4CSGSolid::G4CSGSolid
G4CSGSolid(const G4String &pName)
Definition: G4CSGSolid.cc:42
FatalException
Definition: G4ExceptionSeverity.hh:60
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
G4Cons::~G4Cons ( )

Definition at line 144 of file G4Cons.cc.

145 {
146 }
G4Cons::G4Cons ( __void__ &  a)

Definition at line 129 of file G4Cons.cc.

130  : G4CSGSolid(a), kRadTolerance(0.), kAngTolerance(0.),
131  fRmin1(0.), fRmin2(0.), fRmax1(0.), fRmax2(0.), fDz(0.),
132  fSPhi(0.), fDPhi(0.), sinCPhi(0.), cosCPhi(0.), cosHDPhiOT(0.),
133  cosHDPhiIT(0.), sinSPhi(0.), cosSPhi(0.), sinEPhi(0.), cosEPhi(0.),
134  fPhiFullCone(false), halfCarTolerance(0.), halfRadTolerance(0.),
135  halfAngTolerance(0.)
136 
137 {
138 }
G4CSGSolid::G4CSGSolid
G4CSGSolid(const G4String &pName)
Definition: G4CSGSolid.cc:42
G4Cons::G4Cons ( const G4Cons &  rhs)

Definition at line 152 of file G4Cons.cc.

153  : G4CSGSolid(rhs), kRadTolerance(rhs.kRadTolerance),
154  kAngTolerance(rhs.kAngTolerance), fRmin1(rhs.fRmin1), fRmin2(rhs.fRmin2),
155  fRmax1(rhs.fRmax1), fRmax2(rhs.fRmax2), fDz(rhs.fDz), fSPhi(rhs.fSPhi),
156  fDPhi(rhs.fDPhi), sinCPhi(rhs.sinCPhi), cosCPhi(rhs.cosCPhi),
157  cosHDPhiOT(rhs.cosHDPhiOT), cosHDPhiIT(rhs.cosHDPhiIT),
158  sinSPhi(rhs.sinSPhi), cosSPhi(rhs.cosSPhi), sinEPhi(rhs.sinEPhi),
159  cosEPhi(rhs.cosEPhi), fPhiFullCone(rhs.fPhiFullCone),
160  halfCarTolerance(rhs.halfCarTolerance),
161  halfRadTolerance(rhs.halfRadTolerance),
162  halfAngTolerance(rhs.halfAngTolerance)
163 {
164 }
G4CSGSolid::G4CSGSolid
G4CSGSolid(const G4String &pName)
Definition: G4CSGSolid.cc:42

Member Function Documentation

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

Implements G4VSolid.

Definition at line 270 of file G4Cons.cc.

References G4VSolid::ClipBetweenSections(), G4VSolid::ClipCrossSection(), G4VoxelLimits::GetMaxExtent(), 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().

275 {
276  if ( !pTransform.IsRotated() && (fDPhi == twopi)
277  && (fRmin1 == 0) && (fRmin2 == 0) )
278  {
279  // Special case handling for unrotated solid cones
280  // Compute z/x/y mins and maxs for bounding box respecting limits,
281  // with early returns if outside limits. Then switch() on pAxis,
282  // and compute exact x and y limit for x/y case
283 
284  G4double xoffset, xMin, xMax ;
285  G4double yoffset, yMin, yMax ;
286  G4double zoffset, zMin, zMax ;
287 
288  G4double diff1, diff2, delta, maxDiff, newMin, newMax, RMax ;
289  G4double xoff1, xoff2, yoff1, yoff2 ;
290 
291  zoffset = pTransform.NetTranslation().z();
292  zMin = zoffset - fDz ;
293  zMax = zoffset + fDz ;
294 
295  if (pVoxelLimit.IsZLimited())
296  {
297  if( (zMin > pVoxelLimit.GetMaxZExtent() + kCarTolerance) ||
298  (zMax < pVoxelLimit.GetMinZExtent() - kCarTolerance) )
299  {
300  return false ;
301  }
302  else
303  {
304  if ( zMin < pVoxelLimit.GetMinZExtent() )
305  {
306  zMin = pVoxelLimit.GetMinZExtent() ;
307  }
308  if ( zMax > pVoxelLimit.GetMaxZExtent() )
309  {
310  zMax = pVoxelLimit.GetMaxZExtent() ;
311  }
312  }
313  }
314  xoffset = pTransform.NetTranslation().x() ;
315  RMax = (fRmax2 >= fRmax1) ? zMax : zMin ;
316  xMax = xoffset + (fRmax1 + fRmax2)*0.5 +
317  (RMax - zoffset)*(fRmax2 - fRmax1)/(2*fDz) ;
318  xMin = 2*xoffset-xMax ;
319 
320  if (pVoxelLimit.IsXLimited())
321  {
322  if ( (xMin > pVoxelLimit.GetMaxXExtent() + kCarTolerance) ||
323  (xMax < pVoxelLimit.GetMinXExtent() - kCarTolerance) )
324  {
325  return false ;
326  }
327  else
328  {
329  if ( xMin < pVoxelLimit.GetMinXExtent() )
330  {
331  xMin = pVoxelLimit.GetMinXExtent() ;
332  }
333  if ( xMax > pVoxelLimit.GetMaxXExtent() )
334  {
335  xMax=pVoxelLimit.GetMaxXExtent() ;
336  }
337  }
338  }
339  yoffset = pTransform.NetTranslation().y() ;
340  yMax = yoffset + (fRmax1 + fRmax2)*0.5 +
341  (RMax - zoffset)*(fRmax2 - fRmax1)/(2*fDz) ;
342  yMin = 2*yoffset-yMax ;
343  RMax = yMax - yoffset ; // = max radius due to Zmax/Zmin cuttings
344 
345  if (pVoxelLimit.IsYLimited())
346  {
347  if ( (yMin > pVoxelLimit.GetMaxYExtent() + kCarTolerance) ||
348  (yMax < pVoxelLimit.GetMinYExtent() - kCarTolerance) )
349  {
350  return false ;
351  }
352  else
353  {
354  if ( yMin < pVoxelLimit.GetMinYExtent() )
355  {
356  yMin = pVoxelLimit.GetMinYExtent() ;
357  }
358  if ( yMax > pVoxelLimit.GetMaxYExtent() )
359  {
360  yMax = pVoxelLimit.GetMaxYExtent() ;
361  }
362  }
363  }
364  switch (pAxis) // Known to cut cones
365  {
366  case kXAxis:
367  yoff1 = yoffset - yMin ;
368  yoff2 = yMax - yoffset ;
369 
370  if ((yoff1 >= 0) && (yoff2 >= 0)) // Y limits cross max/min x
371  { // => no change
372  pMin = xMin ;
373  pMax = xMax ;
374  }
375  else
376  {
377  // Y limits don't cross max/min x => compute max delta x,
378  // hence new mins/maxs
379  delta=RMax*RMax-yoff1*yoff1;
380  diff1=(delta>0.) ? std::sqrt(delta) : 0.;
381  delta=RMax*RMax-yoff2*yoff2;
382  diff2=(delta>0.) ? std::sqrt(delta) : 0.;
383  maxDiff = (diff1>diff2) ? diff1:diff2 ;
384  newMin = xoffset - maxDiff ;
385  newMax = xoffset + maxDiff ;
386  pMin = ( newMin < xMin ) ? xMin : newMin ;
387  pMax = ( newMax > xMax) ? xMax : newMax ;
388  }
389  break ;
390 
391  case kYAxis:
392  xoff1 = xoffset - xMin ;
393  xoff2 = xMax - xoffset ;
394 
395  if ((xoff1 >= 0) && (xoff2 >= 0) ) // X limits cross max/min y
396  { // => no change
397  pMin = yMin ;
398  pMax = yMax ;
399  }
400  else
401  {
402  // X limits don't cross max/min y => compute max delta y,
403  // hence new mins/maxs
404  delta=RMax*RMax-xoff1*xoff1;
405  diff1=(delta>0.) ? std::sqrt(delta) : 0.;
406  delta=RMax*RMax-xoff2*xoff2;
407  diff2=(delta>0.) ? std::sqrt(delta) : 0.;
408  maxDiff = (diff1 > diff2) ? diff1:diff2 ;
409  newMin = yoffset - maxDiff ;
410  newMax = yoffset + maxDiff ;
411  pMin = (newMin < yMin) ? yMin : newMin ;
412  pMax = (newMax > yMax) ? yMax : newMax ;
413  }
414  break ;
415 
416  case kZAxis:
417  pMin = zMin ;
418  pMax = zMax ;
419  break ;
420 
421  default:
422  break ;
423  }
424  pMin -= kCarTolerance ;
425  pMax += kCarTolerance ;
426 
427  return true ;
428  }
429  else // Calculate rotated vertex coordinates
430  {
431  G4int i, noEntries, noBetweenSections4 ;
432  G4bool existsAfterClip = false ;
433  G4ThreeVectorList* vertices = CreateRotatedVertices(pTransform) ;
434 
435  pMin = +kInfinity ;
436  pMax = -kInfinity ;
437 
438  noEntries = vertices->size() ;
439  noBetweenSections4 = noEntries-4 ;
440 
441  for ( i = 0 ; i < noEntries ; i += 4 )
442  {
443  ClipCrossSection(vertices, i, pVoxelLimit, pAxis, pMin, pMax) ;
444  }
445  for ( i = 0 ; i < noBetweenSections4 ; i += 4 )
446  {
447  ClipBetweenSections(vertices, i, pVoxelLimit, pAxis, pMin, pMax) ;
448  }
449  if ( (pMin != kInfinity) || (pMax != -kInfinity) )
450  {
451  existsAfterClip = true ;
452 
453  // Add 2*tolerance to avoid precision troubles
454 
455  pMin -= kCarTolerance ;
456  pMax += kCarTolerance ;
457  }
458  else
459  {
460  // Check for case where completely enveloping clipping volume
461  // If point inside then we are confident that the solid completely
462  // envelopes the clipping volume. Hence set min/max extents according
463  // to clipping volume extents along the specified axis.
464 
465  G4ThreeVector clipCentre(
466  (pVoxelLimit.GetMinXExtent() + pVoxelLimit.GetMaxXExtent())*0.5,
467  (pVoxelLimit.GetMinYExtent() + pVoxelLimit.GetMaxYExtent())*0.5,
468  (pVoxelLimit.GetMinZExtent() + pVoxelLimit.GetMaxZExtent())*0.5 ) ;
469 
470  if (Inside(pTransform.Inverse().TransformPoint(clipCentre)) != kOutside)
471  {
472  existsAfterClip = true ;
473  pMin = pVoxelLimit.GetMinExtent(pAxis) ;
474  pMax = pVoxelLimit.GetMaxExtent(pAxis) ;
475  }
476  }
477  delete vertices ;
478  return existsAfterClip ;
479  }
480 }
G4VSolid::ClipCrossSection
void ClipCrossSection(G4ThreeVectorList *pVertices, const G4int pSectionIndex, const G4VoxelLimits &pVoxelLimit, const EAxis pAxis, G4double &pMin, G4double &pMax) const
Definition: G4VSolid.cc:347
G4VoxelLimits::GetMinYExtent
G4double GetMinYExtent() const
CLHEP::Hep3Vector::x
double x() const
G4AffineTransform::Inverse
G4AffineTransform Inverse() const
G4VoxelLimits::IsYLimited
G4bool IsYLimited() const
G4Cons::Inside
EInside Inside(const G4ThreeVector &p) const
Definition: G4Cons.cc:203
G4AffineTransform::IsRotated
G4bool IsRotated() const
G4AffineTransform::NetTranslation
G4ThreeVector NetTranslation() const
G4VoxelLimits::IsXLimited
G4bool IsXLimited() const
G4int
int G4int
Definition: G4Types.hh:78
CLHEP::Hep3Vector::z
double z() const
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
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 * G4Cons::Clone ( ) const
virtual

Reimplemented from G4VSolid.

Definition at line 2254 of file G4Cons.cc.

References G4Cons().

2255 {
2256  return new G4Cons(*this);
2257 }
G4Cons::G4Cons
G4Cons(const G4String &pName, G4double pRmin1, G4double pRmax1, G4double pRmin2, G4double pRmax2, G4double pDz, G4double pSPhi, G4double pDPhi)
Definition: G4Cons.cc:79
void G4Cons::ComputeDimensions ( G4VPVParameterisation *  p,
const G4int  n,
const G4VPhysicalVolume *  pRep 
)
virtual

Reimplemented from G4VSolid.

Definition at line 258 of file G4Cons.cc.

References G4VPVParameterisation::ComputeDimensions().

261 {
262  p->ComputeDimensions(*this,n,pRep) ;
263 }
n
const G4int n
Definition: G4UrQMD1_3Interface.hh:144
G4VPVParameterisation::ComputeDimensions
virtual void ComputeDimensions(G4Box &, const G4int, const G4VPhysicalVolume *) const
Definition: G4VPVParameterisation.hh:95
G4Polyhedron * G4Cons::CreatePolyhedron ( ) const
virtual

Reimplemented from G4VSolid.

Definition at line 2382 of file G4Cons.cc.

2383 {
2384  return new G4PolyhedronCons(fRmin1,fRmax1,fRmin2,fRmax2,fDz,fSPhi,fDPhi);
2385 }
void G4Cons::DescribeYourselfTo ( G4VGraphicsScene &  scene) const
virtual

Implements G4VSolid.

Definition at line 2377 of file G4Cons.cc.

References G4VGraphicsScene::AddSolid().

2378 {
2379  scene.AddSolid (*this);
2380 }
G4VGraphicsScene::AddSolid
virtual void AddSolid(const G4Box &)=0
G4double G4Cons::DistanceToIn ( const G4ThreeVector &  p,
const G4ThreeVector &  v 
) const
virtual

Implements G4VSolid.

Definition at line 718 of file G4Cons.cc.

References test::b, test::c, CLHEP::Hep3Vector::dot(), plottest35::t1, CLHEP::Hep3Vector::x(), CLHEP::Hep3Vector::y(), and CLHEP::Hep3Vector::z().

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

Implements G4VSolid.

Definition at line 1388 of file G4Cons.cc.

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

1389 {
1390  G4double safe=0.0, rho, safeR1, safeR2, safeZ, safePhi, cosPsi ;
1391  G4double tanRMin, secRMin, pRMin ;
1392  G4double tanRMax, secRMax, pRMax ;
1393 
1394  rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
1395  safeZ = std::fabs(p.z()) - fDz ;
1396 
1397  if ( fRmin1 || fRmin2 )
1398  {
1399  tanRMin = (fRmin2 - fRmin1)*0.5/fDz ;
1400  secRMin = std::sqrt(1.0 + tanRMin*tanRMin) ;
1401  pRMin = tanRMin*p.z() + (fRmin1 + fRmin2)*0.5 ;
1402  safeR1 = (pRMin - rho)/secRMin ;
1403 
1404  tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;
1405  secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
1406  pRMax = tanRMax*p.z() + (fRmax1 + fRmax2)*0.5 ;
1407  safeR2 = (rho - pRMax)/secRMax ;
1408 
1409  if ( safeR1 > safeR2) { safe = safeR1; }
1410  else { safe = safeR2; }
1411  }
1412  else
1413  {
1414  tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;
1415  secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
1416  pRMax = tanRMax*p.z() + (fRmax1 + fRmax2)*0.5 ;
1417  safe = (rho - pRMax)/secRMax ;
1418  }
1419  if ( safeZ > safe ) { safe = safeZ; }
1420 
1421  if ( !fPhiFullCone && rho )
1422  {
1423  // Psi=angle from central phi to point
1424 
1425  cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/rho ;
1426 
1427  if ( cosPsi < std::cos(fDPhi*0.5) ) // Point lies outside phi range
1428  {
1429  if ( (p.y()*cosCPhi - p.x()*sinCPhi) <= 0.0 )
1430  {
1431  safePhi = std::fabs(p.x()*std::sin(fSPhi)-p.y()*std::cos(fSPhi));
1432  }
1433  else
1434  {
1435  safePhi = std::fabs(p.x()*sinEPhi-p.y()*cosEPhi);
1436  }
1437  if ( safePhi > safe ) { safe = safePhi; }
1438  }
1439  }
1440  if ( safe < 0.0 ) { safe = 0.0; }
1441 
1442  return safe ;
1443 }
CLHEP::Hep3Vector::x
double x() const
CLHEP::Hep3Vector::z
double z() const
CLHEP::Hep3Vector::y
double y() const
G4double
double G4double
Definition: G4Types.hh:76
G4double G4Cons::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 1450 of file G4Cons.cc.

References test::b, test::c, python.hepunit::degree, CLHEP::Hep3Vector::dot(), G4VSolid::DumpInfo(), G4cout, G4endl, G4Exception(), JustWarning, G4VSolid::kCarTolerance, python.hepunit::mm, python.hepunit::pi, plottest35::t1, python.hepunit::twopi, CLHEP::Hep3Vector::unit(), CLHEP::Hep3Vector::x(), CLHEP::Hep3Vector::y(), and CLHEP::Hep3Vector::z().

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

Implements G4VSolid.

Definition at line 2073 of file G4Cons.cc.

References python.hepunit::degree, G4VSolid::DumpInfo(), G4cout, G4endl, G4Exception(), Inside(), JustWarning, kOutside, python.hepunit::mm, CLHEP::Hep3Vector::x(), CLHEP::Hep3Vector::y(), and CLHEP::Hep3Vector::z().

2074 {
2075  G4double safe=0.0, rho, safeR1, safeR2, safeZ, safePhi;
2076  G4double tanRMin, secRMin, pRMin;
2077  G4double tanRMax, secRMax, pRMax;
2078 
2079 #ifdef G4CSGDEBUG
2080  if( Inside(p) == kOutside )
2081  {
2082  G4int oldprc=G4cout.precision(16) ;
2083  G4cout << G4endl ;
2084  DumpInfo();
2085  G4cout << "Position:" << G4endl << G4endl ;
2086  G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl ;
2087  G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ;
2088  G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ;
2089  G4cout << "pho at z = " << std::sqrt( p.x()*p.x()+p.y()*p.y() )/mm
2090  << " mm" << G4endl << G4endl ;
2091  if( (p.x() != 0.) || (p.x() != 0.) )
2092  {
2093  G4cout << "point phi = " << std::atan2(p.y(),p.x())/degree
2094  << " degree" << G4endl << G4endl ;
2095  }
2096  G4cout.precision(oldprc) ;
2097  G4Exception("G4Cons::DistanceToOut(p)", "GeomSolids1002",
2098  JustWarning, "Point p is outside !?" );
2099  }
2100 #endif
2101 
2102  rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
2103  safeZ = fDz - std::fabs(p.z()) ;
2104 
2105  if (fRmin1 || fRmin2)
2106  {
2107  tanRMin = (fRmin2 - fRmin1)*0.5/fDz ;
2108  secRMin = std::sqrt(1.0 + tanRMin*tanRMin) ;
2109  pRMin = tanRMin*p.z() + (fRmin1 + fRmin2)*0.5 ;
2110  safeR1 = (rho - pRMin)/secRMin ;
2111  }
2112  else
2113  {
2114  safeR1 = kInfinity ;
2115  }
2116 
2117  tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;
2118  secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
2119  pRMax = tanRMax*p.z() + (fRmax1+fRmax2)*0.5 ;
2120  safeR2 = (pRMax - rho)/secRMax ;
2121 
2122  if (safeR1 < safeR2) { safe = safeR1; }
2123  else { safe = safeR2; }
2124  if (safeZ < safe) { safe = safeZ ; }
2125 
2126  // Check if phi divided, Calc distances closest phi plane
2127 
2128  if (!fPhiFullCone)
2129  {
2130  // Above/below central phi of G4Cons?
2131 
2132  if ( (p.y()*cosCPhi - p.x()*sinCPhi) <= 0 )
2133  {
2134  safePhi = -(p.x()*sinSPhi - p.y()*cosSPhi) ;
2135  }
2136  else
2137  {
2138  safePhi = (p.x()*sinEPhi - p.y()*cosEPhi) ;
2139  }
2140  if (safePhi < safe) { safe = safePhi; }
2141  }
2142  if ( safe < 0 ) { safe = 0; }
2143 
2144  return safe ;
2145 }
CLHEP::Hep3Vector::x
double x() const
G4Cons::Inside
EInside Inside(const G4ThreeVector &p) const
Definition: G4Cons.cc:203
JustWarning
Definition: G4ExceptionSeverity.hh:64
G4int
int G4int
Definition: G4Types.hh:78
CLHEP::Hep3Vector::z
double z() const
G4VSolid::DumpInfo
void DumpInfo() const
G4cout
G4GLOB_DLL std::ostream G4cout
python.hepunit.degree
tuple degree
Definition: hepunit.py:69
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
CLHEP::Hep3Vector::y
double y() const
G4endl
#define G4endl
Definition: G4ios.hh:61
G4double
double G4double
Definition: G4Types.hh:76
G4double G4Cons::GetCubicVolume ( )
inlinevirtual

Reimplemented from G4VSolid.

G4double G4Cons::GetDeltaPhiAngle ( ) const
inline

Referenced by G4HepRepFileSceneHandler::AddSolid(), G4HepRepSceneHandler::AddSolid(), G4tgbVolume::BuildSolidForDivision(), G4ParameterisationConsRho::ComputeDimensions(), G4ParameterisationConsZ::ComputeDimensions(), G4GDMLWriteParamvol::Cone_dimensionsWrite(), G4GDMLWriteSolids::ConeWrite(), export_G4Cons(), G4ParameterisationConsPhi::G4ParameterisationConsPhi(), G4VParameterisationCons::G4VParameterisationCons(), G4ParameterisationConsPhi::GetMaxParameter(), and G4tgbGeometryDumper::GetSolidParams().

G4double G4Cons::GetDPhi ( ) const
inline
G4double G4Cons::GetDz ( ) const
inline
G4GeometryType G4Cons::GetEntityType ( ) const
virtual

Implements G4VSolid.

Definition at line 2245 of file G4Cons.cc.

2246 {
2247  return G4String("G4Cons");
2248 }
G4double G4Cons::GetInnerRadiusMinusZ ( ) const
inline

Referenced by G4HepRepFileSceneHandler::AddSolid(), G4HepRepSceneHandler::AddSolid(), G4tgbVolume::BuildSolidForDivision(), G4ParameterisationConsRho::ComputeDimensions(), G4ParameterisationConsPhi::ComputeDimensions(), G4ParameterisationConsZ::ComputeDimensions(), G4GDMLWriteParamvol::Cone_dimensionsWrite(), G4GDMLWriteSolids::ConeWrite(), export_G4Cons(), G4ParameterisationConsRho::G4ParameterisationConsRho(), G4VParameterisationCons::G4VParameterisationCons(), G4ParameterisationConsRho::GetMaxParameter(), and G4tgbGeometryDumper::GetSolidParams().

G4double G4Cons::GetInnerRadiusPlusZ ( ) const
inline

Referenced by G4HepRepFileSceneHandler::AddSolid(), G4HepRepSceneHandler::AddSolid(), G4tgbVolume::BuildSolidForDivision(), G4ParameterisationConsRho::ComputeDimensions(), G4ParameterisationConsPhi::ComputeDimensions(), G4ParameterisationConsZ::ComputeDimensions(), G4GDMLWriteParamvol::Cone_dimensionsWrite(), G4GDMLWriteSolids::ConeWrite(), export_G4Cons(), G4ParameterisationConsRho::G4ParameterisationConsRho(), G4VParameterisationCons::G4VParameterisationCons(), and G4tgbGeometryDumper::GetSolidParams().

G4double G4Cons::GetOuterRadiusMinusZ ( ) const
inline

Referenced by G4HepRepFileSceneHandler::AddSolid(), G4HepRepSceneHandler::AddSolid(), G4tgbVolume::BuildSolidForDivision(), G4ParameterisationConsPhi::ComputeDimensions(), G4ParameterisationConsZ::ComputeDimensions(), G4GDMLWriteParamvol::Cone_dimensionsWrite(), G4GDMLWriteSolids::ConeWrite(), export_G4Cons(), G4ParameterisationConsRho::G4ParameterisationConsRho(), G4VParameterisationCons::G4VParameterisationCons(), G4ParameterisationConsRho::GetMaxParameter(), and G4tgbGeometryDumper::GetSolidParams().

G4double G4Cons::GetOuterRadiusPlusZ ( ) const
inline

Referenced by G4HepRepFileSceneHandler::AddSolid(), G4HepRepSceneHandler::AddSolid(), G4tgbVolume::BuildSolidForDivision(), G4ParameterisationConsRho::ComputeDimensions(), G4ParameterisationConsPhi::ComputeDimensions(), G4ParameterisationConsZ::ComputeDimensions(), G4GDMLWriteParamvol::Cone_dimensionsWrite(), G4GDMLWriteSolids::ConeWrite(), export_G4Cons(), G4VParameterisationCons::G4VParameterisationCons(), and G4tgbGeometryDumper::GetSolidParams().

G4ThreeVector G4Cons::GetPointOnSurface ( ) const
virtual

Reimplemented from G4VSolid.

Definition at line 2290 of file G4Cons.cc.

References G4CSGSolid::GetRadiusInRing(), G4INCL::DeJongSpin::shoot(), and sqr().

2291 {
2292  // declare working variables
2293  //
2294  G4double Aone, Atwo, Athree, Afour, Afive, slin, slout, phi;
2295  G4double zRand, cosu, sinu, rRand1, rRand2, chose, rone, rtwo, qone, qtwo;
2296  rone = (fRmax1-fRmax2)/(2.*fDz);
2297  rtwo = (fRmin1-fRmin2)/(2.*fDz);
2298  qone=0.; qtwo=0.;
2299  if(fRmax1!=fRmax2) { qone = fDz*(fRmax1+fRmax2)/(fRmax1-fRmax2); }
2300  if(fRmin1!=fRmin2) { qtwo = fDz*(fRmin1+fRmin2)/(fRmin1-fRmin2); }
2301  slin = std::sqrt(sqr(fRmin1-fRmin2)+sqr(2.*fDz));
2302  slout = std::sqrt(sqr(fRmax1-fRmax2)+sqr(2.*fDz));
2303  Aone = 0.5*fDPhi*(fRmax2 + fRmax1)*slout;
2304  Atwo = 0.5*fDPhi*(fRmin2 + fRmin1)*slin;
2305  Athree = 0.5*fDPhi*(fRmax1*fRmax1-fRmin1*fRmin1);
2306  Afour = 0.5*fDPhi*(fRmax2*fRmax2-fRmin2*fRmin2);
2307  Afive = fDz*(fRmax1-fRmin1+fRmax2-fRmin2);
2308 
2309  phi = RandFlat::shoot(fSPhi,fSPhi+fDPhi);
2310  cosu = std::cos(phi); sinu = std::sin(phi);
2311  rRand1 = GetRadiusInRing(fRmin1, fRmin2);
2312  rRand2 = GetRadiusInRing(fRmax1, fRmax2);
2313 
2314  if ( (fSPhi == 0.) && fPhiFullCone ) { Afive = 0.; }
2315  chose = RandFlat::shoot(0.,Aone+Atwo+Athree+Afour+2.*Afive);
2316 
2317  if( (chose >= 0.) && (chose < Aone) )
2318  {
2319  if(fRmin1 != fRmin2)
2320  {
2321  zRand = RandFlat::shoot(-1.*fDz,fDz);
2322  return G4ThreeVector (rtwo*cosu*(qtwo-zRand),
2323  rtwo*sinu*(qtwo-zRand), zRand);
2324  }
2325  else
2326  {
2327  return G4ThreeVector(fRmin1*cosu, fRmin2*sinu,
2328  RandFlat::shoot(-1.*fDz,fDz));
2329  }
2330  }
2331  else if( (chose >= Aone) && (chose <= Aone + Atwo) )
2332  {
2333  if(fRmax1 != fRmax2)
2334  {
2335  zRand = RandFlat::shoot(-1.*fDz,fDz);
2336  return G4ThreeVector (rone*cosu*(qone-zRand),
2337  rone*sinu*(qone-zRand), zRand);
2338  }
2339  else
2340  {
2341  return G4ThreeVector(fRmax1*cosu, fRmax2*sinu,
2342  RandFlat::shoot(-1.*fDz,fDz));
2343  }
2344  }
2345  else if( (chose >= Aone + Atwo) && (chose < Aone + Atwo + Athree) )
2346  {
2347  return G4ThreeVector (rRand1*cosu, rRand1*sinu, -1*fDz);
2348  }
2349  else if( (chose >= Aone + Atwo + Athree)
2350  && (chose < Aone + Atwo + Athree + Afour) )
2351  {
2352  return G4ThreeVector (rRand2*cosu,rRand2*sinu,fDz);
2353  }
2354  else if( (chose >= Aone + Atwo + Athree + Afour)
2355  && (chose < Aone + Atwo + Athree + Afour + Afive) )
2356  {
2357  zRand = RandFlat::shoot(-1.*fDz,fDz);
2358  rRand1 = RandFlat::shoot(fRmin2-((zRand-fDz)/(2.*fDz))*(fRmin1-fRmin2),
2359  fRmax2-((zRand-fDz)/(2.*fDz))*(fRmax1-fRmax2));
2360  return G4ThreeVector (rRand1*std::cos(fSPhi),
2361  rRand1*std::sin(fSPhi), zRand);
2362  }
2363  else
2364  {
2365  zRand = RandFlat::shoot(-1.*fDz,fDz);
2366  rRand1 = RandFlat::shoot(fRmin2-((zRand-fDz)/(2.*fDz))*(fRmin1-fRmin2),
2367  fRmax2-((zRand-fDz)/(2.*fDz))*(fRmax1-fRmax2));
2368  return G4ThreeVector (rRand1*std::cos(fSPhi+fDPhi),
2369  rRand1*std::sin(fSPhi+fDPhi), zRand);
2370  }
2371 }
G4INCL::DeJongSpin::shoot
ThreeVector shoot(const G4int Ap, const G4int Af)
Definition: G4INCLDeJongSpin.cc:67
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:42
G4CSGSolid::GetRadiusInRing
G4double GetRadiusInRing(G4double rmin, G4double rmax) const
Definition: G4CSGSolid.cc:101
sqr
T sqr(const T &x)
Definition: templates.hh:145
G4double
double G4double
Definition: G4Types.hh:76
G4double G4Cons::GetRmax1 ( ) const
inline
G4double G4Cons::GetRmax2 ( ) const
inline
G4double G4Cons::GetRmin1 ( ) const
inline
G4double G4Cons::GetRmin2 ( ) const
inline
G4double G4Cons::GetSPhi ( ) const
inline
G4double G4Cons::GetStartPhiAngle ( ) const
inline

Referenced by G4tgbVolume::BuildSolidForDivision(), G4ParameterisationConsRho::ComputeDimensions(), G4ParameterisationConsPhi::ComputeDimensions(), G4ParameterisationConsZ::ComputeDimensions(), G4GDMLWriteParamvol::Cone_dimensionsWrite(), G4GDMLWriteSolids::ConeWrite(), export_G4Cons(), G4VParameterisationCons::G4VParameterisationCons(), and G4tgbGeometryDumper::GetSolidParams().

G4double G4Cons::GetSurfaceArea ( )
inlinevirtual

Reimplemented from G4VSolid.

G4double G4Cons::GetZHalfLength ( ) const
inline

Referenced by G4HepRepFileSceneHandler::AddSolid(), G4HepRepSceneHandler::AddSolid(), G4tgbVolume::BuildSolidForDivision(), G4ParameterisationConsRho::ComputeDimensions(), G4ParameterisationConsPhi::ComputeDimensions(), G4ParameterisationConsZ::ComputeDimensions(), G4ParameterisationConsZ::ComputeTransformation(), G4GDMLWriteParamvol::Cone_dimensionsWrite(), G4GDMLWriteSolids::ConeWrite(), export_G4Cons(), G4ParameterisationConsZ::G4ParameterisationConsZ(), G4VParameterisationCons::G4VParameterisationCons(), G4ParameterisationConsZ::GetMaxParameter(), and G4tgbGeometryDumper::GetSolidParams().

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

Implements G4VSolid.

Definition at line 203 of file G4Cons.cc.

References kInside, kOutside, kSurface, python.hepunit::twopi, CLHEP::Hep3Vector::x(), CLHEP::Hep3Vector::y(), and CLHEP::Hep3Vector::z().

Referenced by CalculateExtent(), and DistanceToOut().

204 {
205  G4double r2, rl, rh, pPhi, tolRMin, tolRMax; // rh2, rl2 ;
206  EInside in;
207 
208  if (std::fabs(p.z()) > fDz + halfCarTolerance ) { return in = kOutside; }
209  else if(std::fabs(p.z()) >= fDz - halfCarTolerance ) { in = kSurface; }
210  else { in = kInside; }
211 
212  r2 = p.x()*p.x() + p.y()*p.y() ;
213  rl = 0.5*(fRmin2*(p.z() + fDz) + fRmin1*(fDz - p.z()))/fDz ;
214  rh = 0.5*(fRmax2*(p.z()+fDz)+fRmax1*(fDz-p.z()))/fDz;
215 
216  // rh2 = rh*rh;
217 
218  tolRMin = rl - halfRadTolerance;
219  if ( tolRMin < 0 ) { tolRMin = 0; }
220  tolRMax = rh + halfRadTolerance;
221 
222  if ( (r2<tolRMin*tolRMin) || (r2>tolRMax*tolRMax) ) { return in = kOutside; }
223 
224  if (rl) { tolRMin = rl + halfRadTolerance; }
225  else { tolRMin = 0.0; }
226  tolRMax = rh - halfRadTolerance;
227 
228  if (in == kInside) // else it's kSurface already
229  {
230  if ( (r2 < tolRMin*tolRMin) || (r2 >= tolRMax*tolRMax) ) { in = kSurface; }
231  }
232  if ( !fPhiFullCone && ((p.x() != 0.0) || (p.y() != 0.0)) )
233  {
234  pPhi = std::atan2(p.y(),p.x()) ;
235 
236  if ( pPhi < fSPhi - halfAngTolerance ) { pPhi += twopi; }
237  else if ( pPhi > fSPhi + fDPhi + halfAngTolerance ) { pPhi -= twopi; }
238 
239  if ( (pPhi < fSPhi - halfAngTolerance) ||
240  (pPhi > fSPhi + fDPhi + halfAngTolerance) ) { return in = kOutside; }
241 
242  else if (in == kInside) // else it's kSurface anyway already
243  {
244  if ( (pPhi < fSPhi + halfAngTolerance) ||
245  (pPhi > fSPhi + fDPhi - halfAngTolerance) ) { in = kSurface; }
246  }
247  }
248  else if ( !fPhiFullCone ) { in = kSurface; }
249 
250  return in ;
251 }
CLHEP::Hep3Vector::x
double x() const
CLHEP::Hep3Vector::z
double z() const
EInside
EInside
Definition: geomdefs.hh:58
CLHEP::Hep3Vector::y
double y() const
G4double
double G4double
Definition: G4Types.hh:76
G4Cons & G4Cons::operator= ( const G4Cons &  rhs)

Definition at line 170 of file G4Cons.cc.

References G4CSGSolid::operator=().

171 {
172  // Check assignment to self
173  //
174  if (this == &rhs) { return *this; }
175 
176  // Copy base class data
177  //
178  G4CSGSolid::operator=(rhs);
179 
180  // Copy data
181  //
182  kRadTolerance = rhs.kRadTolerance;
183  kAngTolerance = rhs.kAngTolerance;
184  fRmin1 = rhs.fRmin1; fRmin2 = rhs.fRmin2;
185  fRmax1 = rhs.fRmax1; fRmax2 = rhs.fRmax2;
186  fDz = rhs.fDz; fSPhi = rhs.fSPhi; fDPhi = rhs.fDPhi;
187  sinCPhi = rhs.sinCPhi; cosCPhi = rhs.cosCPhi;
188  cosHDPhiOT = rhs.cosHDPhiOT; cosHDPhiIT = rhs.cosHDPhiIT;
189  sinSPhi = rhs.sinSPhi; cosSPhi = rhs.cosSPhi;
190  sinEPhi = rhs.sinEPhi; cosEPhi = rhs.cosEPhi;
191  fPhiFullCone = rhs.fPhiFullCone;
192  halfCarTolerance = rhs.halfCarTolerance;
193  halfRadTolerance = rhs.halfRadTolerance;
194  halfAngTolerance = rhs.halfAngTolerance;
195 
196  return *this;
197 }
G4CSGSolid::operator=
G4CSGSolid & operator=(const G4CSGSolid &rhs)
Definition: G4CSGSolid.cc:82
void G4Cons::SetDeltaPhiAngle ( G4double  newDPhi)
inline

Referenced by G4ParameterisationConsRho::ComputeDimensions(), G4ParameterisationConsPhi::ComputeDimensions(), G4ParameterisationConsZ::ComputeDimensions(), and export_G4Cons().

void G4Cons::SetInnerRadiusMinusZ ( G4double  Rmin1)
inline

Referenced by UltraFresnelLensParameterisation::ComputeDimensions(), G4ParameterisationConsRho::ComputeDimensions(), G4ParameterisationConsPhi::ComputeDimensions(), G4ParameterisationConsZ::ComputeDimensions(), and export_G4Cons().

void G4Cons::SetInnerRadiusPlusZ ( G4double  Rmin2)
inline

Referenced by UltraFresnelLensParameterisation::ComputeDimensions(), G4ParameterisationConsRho::ComputeDimensions(), G4ParameterisationConsPhi::ComputeDimensions(), G4ParameterisationConsZ::ComputeDimensions(), and export_G4Cons().

void G4Cons::SetOuterRadiusMinusZ ( G4double  Rmax1)
inline

Referenced by UltraFresnelLensParameterisation::ComputeDimensions(), G4ParameterisationConsRho::ComputeDimensions(), G4ParameterisationConsPhi::ComputeDimensions(), G4ParameterisationConsZ::ComputeDimensions(), and export_G4Cons().

void G4Cons::SetOuterRadiusPlusZ ( G4double  Rmax2)
inline

Referenced by UltraFresnelLensParameterisation::ComputeDimensions(), G4ParameterisationConsRho::ComputeDimensions(), G4ParameterisationConsPhi::ComputeDimensions(), G4ParameterisationConsZ::ComputeDimensions(), and export_G4Cons().

void G4Cons::SetStartPhiAngle ( G4double  newSPhi,
G4bool  trig = true 
)
inline

Referenced by G4ParameterisationConsRho::ComputeDimensions(), G4ParameterisationConsPhi::ComputeDimensions(), G4ParameterisationConsZ::ComputeDimensions(), and export_G4Cons().

void G4Cons::SetZHalfLength ( G4double  newDz)
inline

Referenced by UltraFresnelLensParameterisation::ComputeDimensions(), G4ParameterisationConsRho::ComputeDimensions(), G4ParameterisationConsPhi::ComputeDimensions(), G4ParameterisationConsZ::ComputeDimensions(), and export_G4Cons().

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

Reimplemented from G4CSGSolid.

Definition at line 2263 of file G4Cons.cc.

References python.hepunit::degree, G4VSolid::GetName(), and python.hepunit::mm.

2264 {
2265  G4int oldprc = os.precision(16);
2266  os << "-----------------------------------------------------------\n"
2267  << " *** Dump for solid - " << GetName() << " ***\n"
2268  << " ===================================================\n"
2269  << " Solid type: G4Cons\n"
2270  << " Parameters: \n"
2271  << " inside -fDz radius: " << fRmin1/mm << " mm \n"
2272  << " outside -fDz radius: " << fRmax1/mm << " mm \n"
2273  << " inside +fDz radius: " << fRmin2/mm << " mm \n"
2274  << " outside +fDz radius: " << fRmax2/mm << " mm \n"
2275  << " half length in Z : " << fDz/mm << " mm \n"
2276  << " starting angle of segment: " << fSPhi/degree << " degrees \n"
2277  << " delta angle of segment : " << fDPhi/degree << " degrees \n"
2278  << "-----------------------------------------------------------\n";
2279  os.precision(oldprc);
2280 
2281  return os;
2282 }
G4VSolid::GetName
G4String GetName() const
G4int
int G4int
Definition: G4Types.hh:78
python.hepunit.degree
tuple degree
Definition: hepunit.py:69
G4ThreeVector G4Cons::SurfaceNormal ( const G4ThreeVector &  p) const
virtual

Implements G4VSolid.

Definition at line 488 of file G4Cons.cc.

References G4Exception(), JustWarning, python.hepunit::twopi, CLHEP::Hep3Vector::unit(), CLHEP::Hep3Vector::x(), CLHEP::Hep3Vector::y(), and CLHEP::Hep3Vector::z().

489 {
490  G4int noSurfaces = 0;
491  G4double rho, pPhi;
492  G4double distZ, distRMin, distRMax;
493  G4double distSPhi = kInfinity, distEPhi = kInfinity;
494  G4double tanRMin, secRMin, pRMin, widRMin;
495  G4double tanRMax, secRMax, pRMax, widRMax;
496 
497  G4ThreeVector norm, sumnorm(0.,0.,0.), nZ = G4ThreeVector(0.,0.,1.);
498  G4ThreeVector nR, nr(0.,0.,0.), nPs, nPe;
499 
500  distZ = std::fabs(std::fabs(p.z()) - fDz);
501  rho = std::sqrt(p.x()*p.x() + p.y()*p.y());
502 
503  tanRMin = (fRmin2 - fRmin1)*0.5/fDz;
504  secRMin = std::sqrt(1 + tanRMin*tanRMin);
505  pRMin = rho - p.z()*tanRMin;
506  widRMin = fRmin2 - fDz*tanRMin;
507  distRMin = std::fabs(pRMin - widRMin)/secRMin;
508 
509  tanRMax = (fRmax2 - fRmax1)*0.5/fDz;
510  secRMax = std::sqrt(1+tanRMax*tanRMax);
511  pRMax = rho - p.z()*tanRMax;
512  widRMax = fRmax2 - fDz*tanRMax;
513  distRMax = std::fabs(pRMax - widRMax)/secRMax;
514 
515  if (!fPhiFullCone) // Protected against (0,0,z)
516  {
517  if ( rho )
518  {
519  pPhi = std::atan2(p.y(),p.x());
520 
521  if (pPhi < fSPhi-halfCarTolerance) { pPhi += twopi; }
522  else if (pPhi > fSPhi+fDPhi+halfCarTolerance) { pPhi -= twopi; }
523 
524  distSPhi = std::fabs( pPhi - fSPhi );
525  distEPhi = std::fabs( pPhi - fSPhi - fDPhi );
526  }
527  else if( !(fRmin1) || !(fRmin2) )
528  {
529  distSPhi = 0.;
530  distEPhi = 0.;
531  }
532  nPs = G4ThreeVector(std::sin(fSPhi), -std::cos(fSPhi), 0);
533  nPe = G4ThreeVector(-std::sin(fSPhi+fDPhi), std::cos(fSPhi+fDPhi), 0);
534  }
535  if ( rho > halfCarTolerance )
536  {
537  nR = G4ThreeVector(p.x()/rho/secRMax, p.y()/rho/secRMax, -tanRMax/secRMax);
538  if (fRmin1 || fRmin2)
539  {
540  nr = G4ThreeVector(-p.x()/rho/secRMin,-p.y()/rho/secRMin,tanRMin/secRMin);
541  }
542  }
543 
544  if( distRMax <= halfCarTolerance )
545  {
546  noSurfaces ++;
547  sumnorm += nR;
548  }
549  if( (fRmin1 || fRmin2) && (distRMin <= halfCarTolerance) )
550  {
551  noSurfaces ++;
552  sumnorm += nr;
553  }
554  if( !fPhiFullCone )
555  {
556  if (distSPhi <= halfAngTolerance)
557  {
558  noSurfaces ++;
559  sumnorm += nPs;
560  }
561  if (distEPhi <= halfAngTolerance)
562  {
563  noSurfaces ++;
564  sumnorm += nPe;
565  }
566  }
567  if (distZ <= halfCarTolerance)
568  {
569  noSurfaces ++;
570  if ( p.z() >= 0.) { sumnorm += nZ; }
571  else { sumnorm -= nZ; }
572  }
573  if ( noSurfaces == 0 )
574  {
575 #ifdef G4CSGDEBUG
576  G4Exception("G4Cons::SurfaceNormal(p)", "GeomSolids1002",
577  JustWarning, "Point p is not on surface !?" );
578 #endif
579  norm = ApproxSurfaceNormal(p);
580  }
581  else if ( noSurfaces == 1 ) { norm = sumnorm; }
582  else { norm = sumnorm.unit(); }
583 
584  return norm ;
585 }
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:42
CLHEP::Hep3Vector::x
double x() const
JustWarning
Definition: G4ExceptionSeverity.hh:64
G4int
int G4int
Definition: G4Types.hh:78
CLHEP::Hep3Vector::z
double z() const
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
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:17 for Geant4.10 by   doxygen 1.8.7