Geant4  10.03.p01
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
G4EllipticalCone Class Reference

#include <G4EllipticalCone.hh>

Inheritance diagram for G4EllipticalCone:
Collaboration diagram for G4EllipticalCone:

Public Member Functions

 G4EllipticalCone (const G4String &pName, G4double pxSemiAxis, G4double pySemiAxis, G4double zMax, G4double pzTopCut)
 
virtual ~G4EllipticalCone ()
 
G4double GetSemiAxisMax () const
 
G4double GetSemiAxisX () const
 
G4double GetSemiAxisY () const
 
G4double GetZMax () const
 
G4double GetZTopCut () const
 
void SetSemiAxis (G4double x, G4double y, G4double z)
 
void SetZCut (G4double newzTopCut)
 
G4double GetCubicVolume ()
 
G4double GetSurfaceArea ()
 
void Extent (G4ThreeVector &pMin, G4ThreeVector &pMax) const
 
G4bool CalculateExtent (const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pMin, G4double &pMax) const
 
EInside Inside (const G4ThreeVector &p) const
 
G4ThreeVector SurfaceNormal (const G4ThreeVector &p) const
 
G4double DistanceToIn (const G4ThreeVector &p, const G4ThreeVector &v) const
 
G4double DistanceToIn (const G4ThreeVector &p) const
 
G4double DistanceToOut (const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=G4bool(false), G4bool *validNorm=0, G4ThreeVector *n=0) const
 
G4double DistanceToOut (const G4ThreeVector &p) const
 
G4GeometryType GetEntityType () const
 
G4VSolidClone () const
 
G4ThreeVector GetPointOnSurface () const
 
std::ostream & StreamInfo (std::ostream &os) const
 
G4PolyhedronGetPolyhedron () const
 
void DescribeYourselfTo (G4VGraphicsScene &scene) const
 
G4VisExtent GetExtent () const
 
G4PolyhedronCreatePolyhedron () const
 
 G4EllipticalCone (__void__ &)
 
 G4EllipticalCone (const G4EllipticalCone &rhs)
 
G4EllipticalConeoperator= (const G4EllipticalCone &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 const G4VSolidGetConstituentSolid (G4int no) const
 
virtual G4VSolidGetConstituentSolid (G4int no)
 
virtual const G4DisplacedSolidGetDisplacedSolidPtr () const
 
virtual G4DisplacedSolidGetDisplacedSolidPtr ()
 
 G4VSolid (__void__ &)
 
 G4VSolid (const G4VSolid &rhs)
 
G4VSolidoperator= (const G4VSolid &rhs)
 
G4double EstimateCubicVolume (G4int nStat, G4double epsilon) const
 
G4double EstimateSurfaceArea (G4int nStat, G4double ell) const
 

Protected Attributes

G4bool fRebuildPolyhedron
 
G4PolyhedronfpPolyhedron
 
- Protected Attributes inherited from G4VSolid
G4double kCarTolerance
 

Additional Inherited Members

- 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
 

Detailed Description

Definition at line 85 of file G4EllipticalCone.hh.

Constructor & Destructor Documentation

G4EllipticalCone::G4EllipticalCone ( const G4String pName,
G4double  pxSemiAxis,
G4double  pySemiAxis,
G4double  zMax,
G4double  pzTopCut 
)

Definition at line 70 of file G4EllipticalCone.cc.

75  : G4VSolid(pName), fRebuildPolyhedron(false), fpPolyhedron(0),
76  fCubicVolume(0.), fSurfaceArea(0.), zTopCut(0.)
77 {
78 
80 
81  halfRadTol = 0.5*kRadTolerance;
82  halfCarTol = 0.5*kCarTolerance;
83 
84  // Check Semi-Axis & Z-cut
85  //
86  if ( (pxSemiAxis <= 0.) || (pySemiAxis <= 0.) || (pzMax <= 0.) )
87  {
88  std::ostringstream message;
89  message << "Invalid semi-axis or height - " << GetName();
90  G4Exception("G4EllipticalCone::G4EllipticalCone()", "GeomSolids0002",
91  FatalErrorInArgument, message);
92  }
93  if ( pzTopCut <= 0 )
94  {
95  std::ostringstream message;
96  message << "Invalid z-coordinate for cutting plane - " << GetName();
97  G4Exception("G4EllipticalCone::G4EllipticalCone()", "InvalidSetup",
98  FatalErrorInArgument, message);
99  }
100 
101  SetSemiAxis( pxSemiAxis, pySemiAxis, pzMax );
102  SetZCut(pzTopCut);
103 }
G4String GetName() const
G4Polyhedron * fpPolyhedron
G4double GetRadialTolerance() const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
void SetSemiAxis(G4double x, G4double y, G4double z)
void SetZCut(G4double newzTopCut)
G4VSolid(const G4String &name)
Definition: G4VSolid.cc:61
G4double kCarTolerance
Definition: G4VSolid.hh:307
static G4GeometryTolerance * GetInstance()

Here is the call graph for this function:

Here is the caller graph for this function:

G4EllipticalCone::~G4EllipticalCone ( )
virtual

Definition at line 122 of file G4EllipticalCone.cc.

123 {
124  delete fpPolyhedron; fpPolyhedron = 0;
125 }
G4Polyhedron * fpPolyhedron
G4EllipticalCone::G4EllipticalCone ( __void__ &  a)

Definition at line 110 of file G4EllipticalCone.cc.

111  : G4VSolid(a), fRebuildPolyhedron(false), fpPolyhedron(0),
112  kRadTolerance(0.), halfRadTol(0.), halfCarTol(0.), fCubicVolume(0.),
113  fSurfaceArea(0.), xSemiAxis(0.), ySemiAxis(0.), zheight(0.),
114  semiAxisMax(0.), zTopCut(0.)
115 {
116 }
G4Polyhedron * fpPolyhedron
G4VSolid(const G4String &name)
Definition: G4VSolid.cc:61
G4EllipticalCone::G4EllipticalCone ( const G4EllipticalCone rhs)

Definition at line 131 of file G4EllipticalCone.cc.

132  : G4VSolid(rhs),
133  fRebuildPolyhedron(false), fpPolyhedron(0),
134  kRadTolerance(rhs.kRadTolerance),
135  halfRadTol(rhs.halfRadTol), halfCarTol(rhs.halfCarTol),
136  fCubicVolume(rhs.fCubicVolume), fSurfaceArea(rhs.fSurfaceArea),
137  xSemiAxis(rhs.xSemiAxis), ySemiAxis(rhs.ySemiAxis), zheight(rhs.zheight),
138  semiAxisMax(rhs.semiAxisMax), zTopCut(rhs.zTopCut)
139 {
140 }
G4Polyhedron * fpPolyhedron
G4VSolid(const G4String &name)
Definition: G4VSolid.cc:61

Member Function Documentation

G4bool G4EllipticalCone::CalculateExtent ( const EAxis  pAxis,
const G4VoxelLimits pVoxelLimit,
const G4AffineTransform pTransform,
G4double pMin,
G4double pMax 
) const
virtual

Implements G4VSolid.

Definition at line 202 of file G4EllipticalCone.cc.

206 {
207  G4ThreeVector bmin,bmax;
208  G4bool exist;
209 
210  // Check bounding box (bbox)
211  //
212  Extent(bmin,bmax);
213  G4BoundingEnvelope bbox(bmin,bmax);
214 #ifdef G4BBOX_EXTENT
215  if (true) return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
216 #endif
217  if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVoxelLimit,pTransform,pMin,pMax))
218  {
219  return exist = (pMin < pMax) ? true : false;
220  }
221 
222  // Set bounding envelope (benv) and calculate extent
223  //
224  static const G4int NSTEPS = 48; // number of steps for whole circle
225  static const G4double ang = twopi/NSTEPS;
226  static const G4double sinHalf = std::sin(0.5*ang);
227  static const G4double cosHalf = std::cos(0.5*ang);
228  static const G4double sinStep = 2.*sinHalf*cosHalf;
229  static const G4double cosStep = 1. - 2.*sinHalf*sinHalf;
230  G4double zcut = bmax.z();
231  G4double height = GetZMax();
232  G4double sxmin = GetSemiAxisX()*(height-zcut)/cosHalf;
233  G4double symin = GetSemiAxisY()*(height-zcut)/cosHalf;
234  G4double sxmax = bmax.x()/cosHalf;
235  G4double symax = bmax.y()/cosHalf;
236 
237  G4double sinCur = sinHalf;
238  G4double cosCur = cosHalf;
239  G4ThreeVectorList baseA(NSTEPS),baseB(NSTEPS);
240  for (G4int k=0; k<NSTEPS; ++k)
241  {
242  baseA[k].set(sxmax*cosCur,symax*sinCur,-zcut);
243  baseB[k].set(sxmin*cosCur,symin*sinCur, zcut);
244 
245  G4double sinTmp = sinCur;
246  sinCur = sinCur*cosStep + cosCur*sinStep;
247  cosCur = cosCur*cosStep - sinTmp*sinStep;
248  }
249 
250  std::vector<const G4ThreeVectorList *> polygons(2);
251  polygons[0] = &baseA;
252  polygons[1] = &baseB;
253  G4BoundingEnvelope benv(bmin,bmax,polygons);
254  exist = benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
255  return exist;
256 }
G4double GetSemiAxisX() const
double x() const
int G4int
Definition: G4Types.hh:78
double z() const
static constexpr double twopi
Definition: G4SIunits.hh:76
bool G4bool
Definition: G4Types.hh:79
void Extent(G4ThreeVector &pMin, G4ThreeVector &pMax) const
G4double GetZMax() const
std::vector< G4ThreeVector > G4ThreeVectorList
G4double GetSemiAxisY() const
double y() const
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

G4VSolid * G4EllipticalCone::Clone ( ) const
virtual

Reimplemented from G4VSolid.

Definition at line 937 of file G4EllipticalCone.cc.

938 {
939  return new G4EllipticalCone(*this);
940 }
G4EllipticalCone(const G4String &pName, G4double pxSemiAxis, G4double pySemiAxis, G4double zMax, G4double pzTopCut)

Here is the call graph for this function:

G4Polyhedron * G4EllipticalCone::CreatePolyhedron ( ) const
virtual

Reimplemented from G4VSolid.

Definition at line 1046 of file G4EllipticalCone.cc.

1047 {
1048  return new G4PolyhedronEllipticalCone(xSemiAxis, ySemiAxis, zheight, zTopCut);
1049 }

Here is the caller graph for this function:

void G4EllipticalCone::DescribeYourselfTo ( G4VGraphicsScene scene) const
virtual

Implements G4VSolid.

Definition at line 1030 of file G4EllipticalCone.cc.

1031 {
1032  scene.AddSolid(*this);
1033 }
virtual void AddSolid(const G4Box &)=0

Here is the call graph for this function:

G4double G4EllipticalCone::DistanceToIn ( const G4ThreeVector p,
const G4ThreeVector v 
) const
virtual

Implements G4VSolid.

Definition at line 410 of file G4EllipticalCone.cc.

412 {
413  G4double distMin = kInfinity;
414 
415  // code from EllipticalTube
416 
417  G4double sigz = p.z()+zTopCut;
418 
419  //
420  // Check z = -dz planer surface
421  //
422 
423  if (sigz < halfCarTol)
424  {
425  //
426  // We are "behind" the shape in z, and so can
427  // potentially hit the rear face. Correct direction?
428  //
429  if (v.z() <= 0)
430  {
431  //
432  // As long as we are far enough away, we know we
433  // can't intersect
434  //
435  if (sigz < 0) return kInfinity;
436 
437  //
438  // Otherwise, we don't intersect unless we are
439  // on the surface of the ellipse
440  //
441 
442  if ( sqr(p.x()/( xSemiAxis - halfCarTol ))
443  + sqr(p.y()/( ySemiAxis - halfCarTol )) <= sqr( zheight+zTopCut ) )
444  return kInfinity;
445 
446  }
447  else
448  {
449  //
450  // How far?
451  //
452  G4double q = -sigz/v.z();
453 
454  //
455  // Where does that place us?
456  //
457  G4double xi = p.x() + q*v.x(),
458  yi = p.y() + q*v.y();
459 
460  //
461  // Is this on the surface (within ellipse)?
462  //
463  if ( sqr(xi/xSemiAxis) + sqr(yi/ySemiAxis) <= sqr( zheight + zTopCut ) )
464  {
465  //
466  // Yup. Return q, unless we are on the surface
467  //
468  return (sigz < -halfCarTol) ? q : 0;
469  }
470  else if (xi/(xSemiAxis*xSemiAxis)*v.x()
471  + yi/(ySemiAxis*ySemiAxis)*v.y() >= 0)
472  {
473  //
474  // Else, if we are traveling outwards, we know
475  // we must miss
476  //
477  // return kInfinity;
478  }
479  }
480  }
481 
482  //
483  // Check z = +dz planer surface
484  //
485  sigz = p.z() - zTopCut;
486 
487  if (sigz > -halfCarTol)
488  {
489  if (v.z() >= 0)
490  {
491 
492  if (sigz > 0) return kInfinity;
493 
494  if ( sqr(p.x()/( xSemiAxis - halfCarTol ))
495  + sqr(p.y()/( ySemiAxis - halfCarTol )) <= sqr( zheight-zTopCut ) )
496  return kInfinity;
497 
498  }
499  else {
500  G4double q = -sigz/v.z();
501 
502  G4double xi = p.x() + q*v.x(),
503  yi = p.y() + q*v.y();
504 
505  if ( sqr(xi/xSemiAxis) + sqr(yi/ySemiAxis) <= sqr( zheight - zTopCut ) )
506  {
507  return (sigz > -halfCarTol) ? q : 0;
508  }
509  else if (xi/(xSemiAxis*xSemiAxis)*v.x()
510  + yi/(ySemiAxis*ySemiAxis)*v.y() >= 0)
511  {
512  // return kInfinity;
513  }
514  }
515  }
516 
517 
518 #if 0
519 
520  // check to see if Z plane is relevant
521  //
522  if (p.z() < -zTopCut - 0.5*kCarTolerance)
523  {
524  if (v.z() <= 0.0)
525  return distMin;
526 
527  G4double lambda = (-zTopCut - p.z())/v.z();
528 
529  if ( sqr((lambda*v.x()+p.x())/xSemiAxis) +
530  sqr((lambda*v.y()+p.y())/ySemiAxis) <=
531  sqr(zTopCut + zheight + 0.5*kRadTolerance) )
532  {
533  return distMin = std::fabs(lambda);
534  }
535  }
536 
537  if (p.z() > zTopCut+0.5*kCarTolerance)
538  {
539  if (v.z() >= 0.0)
540  { return distMin; }
541 
542  G4double lambda = (zTopCut - p.z()) / v.z();
543 
544  if ( sqr((lambda*v.x() + p.x())/xSemiAxis) +
545  sqr((lambda*v.y() + p.y())/ySemiAxis) <=
546  sqr(zheight - zTopCut + 0.5*kRadTolerance) )
547  {
548  return distMin = std::fabs(lambda);
549  }
550  }
551 
552  if (p.z() > zTopCut - halfCarTol
553  && p.z() < zTopCut + halfCarTol )
554  {
555  if (v.z() > 0.)
556  { return kInfinity; }
557 
558  return distMin = 0.;
559  }
560 
561  if (p.z() < -zTopCut + halfCarTol
562  && p.z() > -zTopCut - halfCarTol)
563  {
564  if (v.z() < 0.)
565  { return distMin = kInfinity; }
566 
567  return distMin = 0.;
568  }
569 
570 #endif
571 
572  // if we are here then it either intersects or grazes the curved surface
573  // or it does not intersect at all
574  //
575  G4double A = sqr(v.x()/xSemiAxis) + sqr(v.y()/ySemiAxis) - sqr(v.z());
576  G4double B = 2*(v.x()*p.x()/sqr(xSemiAxis) +
577  v.y()*p.y()/sqr(ySemiAxis) + v.z()*(zheight-p.z()));
578  G4double C = sqr(p.x()/xSemiAxis) + sqr(p.y()/ySemiAxis) -
579  sqr(zheight - p.z());
580 
581  G4double discr = B*B - 4.*A*C;
582 
583  // if the discriminant is negative it never hits the curved object
584  //
585  if ( discr < -halfCarTol )
586  { return distMin; }
587 
588  // case below is when it hits or grazes the surface
589  //
590  if ( (discr >= - halfCarTol ) && (discr < halfCarTol ) )
591  {
592  return distMin = std::fabs(-B/(2.*A));
593  }
594 
595  G4double plus = (-B+std::sqrt(discr))/(2.*A);
596  G4double minus = (-B-std::sqrt(discr))/(2.*A);
597 
598  // Special case::Point on Surface, Check norm.dot(v)
599 
600  if ( ( std::fabs(plus) < halfCarTol )||( std::fabs(minus) < halfCarTol ) )
601  {
602  G4ThreeVector truenorm(p.x()/(xSemiAxis*xSemiAxis),
603  p.y()/(ySemiAxis*ySemiAxis),
604  -( p.z() - zheight ));
605  if ( truenorm*v >= 0) // going outside the solid from surface
606  {
607  return kInfinity;
608  }
609  else
610  {
611  return 0;
612  }
613  }
614 
615  // G4double lambda = std::fabs(plus) < std::fabs(minus) ? plus : minus;
616  G4double lambda = 0;
617 
618  if ( minus > halfCarTol && minus < distMin )
619  {
620  lambda = minus ;
621  // check normal vector n * v < 0
622  G4ThreeVector pin = p + lambda*v;
623  if(std::fabs(pin.z())<zTopCut+0.5*kCarTolerance)
624  {
625  G4ThreeVector truenorm(pin.x()/(xSemiAxis*xSemiAxis),
626  pin.y()/(ySemiAxis*ySemiAxis),
627  - ( pin.z() - zheight ));
628  if ( truenorm*v < 0)
629  { // yes, going inside the solid
630  distMin = lambda;
631  }
632  }
633  }
634  if ( plus > halfCarTol && plus < distMin )
635  {
636  lambda = plus ;
637  // check normal vector n * v < 0
638  G4ThreeVector pin = p + lambda*v;
639  if(std::fabs(pin.z())<zTopCut+0.5*kCarTolerance)
640  {
641  G4ThreeVector truenorm(pin.x()/(xSemiAxis*xSemiAxis),
642  pin.y()/(ySemiAxis*ySemiAxis),
643  - ( pin.z() - zheight ) );
644  if ( truenorm*v < 0)
645  { // yes, going inside the solid
646  distMin = lambda;
647  }
648  }
649  }
650  if (distMin < halfCarTol) distMin=0.;
651  return distMin ;
652 }
static const G4double kInfinity
Definition: geomdefs.hh:42
double x() const
double C(double temp)
double B(double temperature)
double z() const
double A(double temperature)
double y() const
G4double kCarTolerance
Definition: G4VSolid.hh:307
T sqr(const T &x)
Definition: templates.hh:145
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

G4double G4EllipticalCone::DistanceToIn ( const G4ThreeVector p) const
virtual

Implements G4VSolid.

Definition at line 659 of file G4EllipticalCone.cc.

660 {
661  G4double distR, distR2, distZ, maxDim;
662  G4double distRad;
663 
664  // check if the point lies either below z=-zTopCut in bottom elliptical
665  // region or on top within cut elliptical region
666  //
667  if( (p.z() <= -zTopCut) && (sqr(p.x()/xSemiAxis) + sqr(p.y()/ySemiAxis)
668  <= sqr(zTopCut + zheight + 0.5*kCarTolerance )) )
669  {
670  //return distZ = std::fabs(zTopCut - p.z());
671  return distZ = std::fabs(zTopCut + p.z());
672  }
673 
674  if( (p.z() >= zTopCut) && (sqr(p.x()/xSemiAxis)+sqr(p.y()/ySemiAxis)
675  <= sqr(zheight - zTopCut + kCarTolerance/2.0 )) )
676  {
677  return distZ = std::fabs(p.z() - zTopCut);
678  }
679 
680  // below we use the following approximation: we take the largest of the
681  // axes and find the shortest distance to the circular (cut) cone of that
682  // radius.
683  //
684  maxDim = xSemiAxis >= ySemiAxis ? xSemiAxis:ySemiAxis;
685  distRad = std::sqrt(p.x()*p.x()+p.y()*p.y());
686 
687  if( p.z() > maxDim*distRad + zTopCut*(1.+maxDim)-sqr(maxDim)*zheight )
688  {
689  distR2 = sqr(p.z() - zTopCut) + sqr(distRad - maxDim*(zheight - zTopCut));
690  return std::sqrt( distR2 );
691  }
692 
693  if( distRad > maxDim*( zheight - p.z() ) )
694  {
695  if( p.z() > maxDim*distRad - (zTopCut*(1.+maxDim)+sqr(maxDim)*zheight) )
696  {
697  G4double zVal = (p.z()-maxDim*(distRad-maxDim*zheight))/(1.+sqr(maxDim));
698  G4double rVal = maxDim*(zheight - zVal);
699  return distR = std::sqrt(sqr(p.z() - zVal) + sqr(distRad - rVal));
700  }
701  }
702 
703  if( distRad <= maxDim*(zheight - p.z()) )
704  {
705  distR2 = sqr(distRad - maxDim*(zheight + zTopCut)) + sqr(p.z() + zTopCut);
706  return std::sqrt( distR2 );
707  }
708 
709  return distR = 0;
710 }
double x() const
double z() const
double y() const
G4double kCarTolerance
Definition: G4VSolid.hh:307
T sqr(const T &x)
Definition: templates.hh:145
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

G4double G4EllipticalCone::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 717 of file G4EllipticalCone.cc.

722 {
723  G4double distMin, lambda;
724  enum surface_e {kPlaneSurf, kCurvedSurf, kNoSurf} surface;
725 
726  distMin = kInfinity;
727  surface = kNoSurf;
728 
729  if (v.z() < 0.0)
730  {
731  lambda = (-p.z() - zTopCut)/v.z();
732 
733  if ( (sqr((p.x() + lambda*v.x())/xSemiAxis) +
734  sqr((p.y() + lambda*v.y())/ySemiAxis)) <
735  sqr(zheight + zTopCut + 0.5*kCarTolerance) )
736  {
737  distMin = std::fabs(lambda);
738 
739  if (!calcNorm) { return distMin; }
740  }
741  distMin = std::fabs(lambda);
742  surface = kPlaneSurf;
743  }
744 
745  if (v.z() > 0.0)
746  {
747  lambda = (zTopCut - p.z()) / v.z();
748 
749  if ( (sqr((p.x() + lambda*v.x())/xSemiAxis)
750  + sqr((p.y() + lambda*v.y())/ySemiAxis) )
751  < (sqr(zheight - zTopCut + 0.5*kCarTolerance)) )
752  {
753  distMin = std::fabs(lambda);
754  if (!calcNorm) { return distMin; }
755  }
756  distMin = std::fabs(lambda);
757  surface = kPlaneSurf;
758  }
759 
760  // if we are here then it either intersects or grazes the
761  // curved surface...
762  //
763  G4double A = sqr(v.x()/xSemiAxis) + sqr(v.y()/ySemiAxis) - sqr(v.z());
764  G4double B = 2.*(v.x()*p.x()/sqr(xSemiAxis) +
765  v.y()*p.y()/sqr(ySemiAxis) + v.z()*(zheight-p.z()));
766  G4double C = sqr(p.x()/xSemiAxis) + sqr(p.y()/ySemiAxis)
767  - sqr(zheight - p.z());
768 
769  G4double discr = B*B - 4.*A*C;
770 
771  if ( discr >= - 0.5*kCarTolerance && discr < 0.5*kCarTolerance )
772  {
773  if(!calcNorm) { return distMin = std::fabs(-B/(2.*A)); }
774  }
775 
776  else if ( discr > 0.5*kCarTolerance )
777  {
778  G4double plus = (-B+std::sqrt(discr))/(2.*A);
779  G4double minus = (-B-std::sqrt(discr))/(2.*A);
780 
781  if ( plus > 0.5*kCarTolerance && minus > 0.5*kCarTolerance )
782  {
783  // take the shorter distance
784  //
785  lambda = std::fabs(plus) < std::fabs(minus) ? plus : minus;
786  }
787  else
788  {
789  // at least one solution is close to zero or negative
790  // so, take small positive solution or zero
791  //
792  lambda = plus > -0.5*kCarTolerance ? plus : 0;
793  }
794 
795  if ( std::fabs(lambda) < distMin )
796  {
797  if( std::fabs(lambda) > 0.5*kCarTolerance)
798  {
799  distMin = std::fabs(lambda);
800  surface = kCurvedSurf;
801  }
802  else // Point is On the Surface, Check Normal
803  {
804  G4ThreeVector truenorm(p.x()/(xSemiAxis*xSemiAxis),
805  p.y()/(ySemiAxis*ySemiAxis),
806  -( p.z() - zheight ));
807  if( truenorm.dot(v) > 0 )
808  {
809  distMin = 0.0;
810  surface = kCurvedSurf;
811  }
812  }
813  }
814  }
815 
816  // set normal if requested
817  //
818  if (calcNorm)
819  {
820  if (surface == kNoSurf)
821  {
822  *validNorm = false;
823  }
824  else
825  {
826  *validNorm = true;
827  switch (surface)
828  {
829  case kPlaneSurf:
830  {
831  *n = G4ThreeVector(0.,0.,(v.z() > 0.0 ? 1. : -1.));
832  }
833  break;
834 
835  case kCurvedSurf:
836  {
837  G4ThreeVector pexit = p + distMin*v;
838  G4ThreeVector truenorm( pexit.x()/(xSemiAxis*xSemiAxis),
839  pexit.y()/(ySemiAxis*ySemiAxis),
840  -( pexit.z() - zheight ) );
841  truenorm /= truenorm.mag();
842  *n= truenorm;
843  }
844  break;
845 
846  default: // Should never reach this case ...
847  DumpInfo();
848  std::ostringstream message;
849  G4int oldprc = message.precision(16);
850  message << "Undefined side for valid surface normal to solid."
851  << G4endl
852  << "Position:" << G4endl
853  << " p.x() = " << p.x()/mm << " mm" << G4endl
854  << " p.y() = " << p.y()/mm << " mm" << G4endl
855  << " p.z() = " << p.z()/mm << " mm" << G4endl
856  << "Direction:" << G4endl
857  << " v.x() = " << v.x() << G4endl
858  << " v.y() = " << v.y() << G4endl
859  << " v.z() = " << v.z() << G4endl
860  << "Proposed distance :" << G4endl
861  << " distMin = " << distMin/mm << " mm";
862  message.precision(oldprc);
863  G4Exception("G4EllipticalCone::DistanceToOut(p,v,..)",
864  "GeomSolids1002", JustWarning, message);
865  break;
866  }
867  }
868  }
869 
870  if (distMin<0.5*kCarTolerance) { distMin=0; }
871 
872  return distMin;
873 }
static constexpr double mm
Definition: G4SIunits.hh:115
static const G4double kInfinity
Definition: geomdefs.hh:42
CLHEP::Hep3Vector G4ThreeVector
double x() const
double C(double temp)
double B(double temperature)
int G4int
Definition: G4Types.hh:78
double z() const
void DumpInfo() const
double A(double temperature)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
double y() const
#define G4endl
Definition: G4ios.hh:61
G4double kCarTolerance
Definition: G4VSolid.hh:307
T sqr(const T &x)
Definition: templates.hh:145
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

G4double G4EllipticalCone::DistanceToOut ( const G4ThreeVector p) const
virtual

Implements G4VSolid.

Definition at line 879 of file G4EllipticalCone.cc.

880 {
881 #ifdef G4SPECSDEBUG
882  if( Inside(p) == kOutside )
883  {
884  DumpInfo();
885  std::ostringstream message;
886  G4int oldprc = message.precision(16);
887  message << "Point p is outside !?" << G4endl
888  << "Position:" << G4endl
889  << " p.x() = " << p.x()/mm << " mm" << G4endl
890  << " p.y() = " << p.y()/mm << " mm" << G4endl
891  << " p.z() = " << p.z()/mm << " mm";
892  message.precision(oldprc) ;
893  G4Exception("G4Ellipsoid::DistanceToOut(p)", "GeomSolids1002",
894  JustWarning, message);
895  }
896 #endif
897  // The safety is calculated in the scaled space where elliptical cone
898  // becomes a circular cone with radius equal to the smaller of the axes
899  //
900  G4double px = p.x(), py = p.y(), pz = p.z();
901  G4double axis;
902  if (xSemiAxis < ySemiAxis)
903  {
904  axis = xSemiAxis;
905  py *= xSemiAxis/ySemiAxis; // scale y
906  }
907  else
908  {
909  axis = ySemiAxis;
910  px *= ySemiAxis/xSemiAxis; // scale x
911  }
912 
913  G4double distZ = zTopCut - std::abs(pz) ;
914  if (distZ <= 0) return 0; // point is outside
915 
916  G4double rho = axis*(zheight-pz); // radius at z = p.z()
917  G4double pr = std::sqrt(px*px+py*py);
918  if (pr >= rho) return 0; // point is outside
919 
920  G4double distR = (rho-pr)/std::sqrt(1+axis*axis);
921  return std::min(distR,distZ);
922 }
static constexpr double mm
Definition: G4SIunits.hh:115
double x() const
EInside Inside(const G4ThreeVector &p) const
int G4int
Definition: G4Types.hh:78
double z() const
void DumpInfo() const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
double y() const
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

void G4EllipticalCone::Extent ( G4ThreeVector pMin,
G4ThreeVector pMax 
) const
virtual

Reimplemented from G4VSolid.

Definition at line 173 of file G4EllipticalCone.cc.

174 {
175  G4double zcut = GetZTopCut();
176  G4double height = GetZMax();
177  G4double xmax = GetSemiAxisX()*(height+zcut);
178  G4double ymax = GetSemiAxisY()*(height+zcut);
179  pMin.set(-xmax,-ymax,-zcut);
180  pMax.set( xmax, ymax, zcut);
181 
182  // Check correctness of the bounding box
183  //
184  if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
185  {
186  std::ostringstream message;
187  message << "Bad bounding box (min >= max) for solid: "
188  << GetName() << " !"
189  << "\npMin = " << pMin
190  << "\npMax = " << pMax;
191  G4Exception("G4EllipticalCone::Extent()", "GeomMgt0001",
192  JustWarning, message);
193  DumpInfo();
194  }
195 }
void set(double x, double y, double z)
G4String GetName() const
G4double GetSemiAxisX() const
double x() const
double z() const
void DumpInfo() const
G4double GetZMax() const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4double GetSemiAxisY() const
double y() const
G4double GetZTopCut() const
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

Here is the caller graph for this function:

G4double G4EllipticalCone::GetCubicVolume ( )
inlinevirtual

Reimplemented from G4VSolid.

G4GeometryType G4EllipticalCone::GetEntityType ( ) const
virtual

Implements G4VSolid.

Definition at line 928 of file G4EllipticalCone.cc.

929 {
930  return G4String("G4EllipticalCone");
931 }
G4VisExtent G4EllipticalCone::GetExtent ( ) const
virtual

Reimplemented from G4VSolid.

Definition at line 1035 of file G4EllipticalCone.cc.

1036 {
1037  // Define the sides of the box into which the solid instance would fit.
1038  //
1039  G4ThreeVector pmin,pmax;
1040  Extent(pmin,pmax);
1041  return G4VisExtent(pmin.x(),pmax.x(),
1042  pmin.y(),pmax.y(),
1043  pmin.z(),pmax.z());
1044 }
double x() const
double z() const
void Extent(G4ThreeVector &pMin, G4ThreeVector &pMax) const
double y() const

Here is the call graph for this function:

G4ThreeVector G4EllipticalCone::GetPointOnSurface ( ) const
virtual

Reimplemented from G4VSolid.

Definition at line 971 of file G4EllipticalCone.cc.

972 {
973 
974  G4double phi, sinphi, cosphi, aOne, aTwo, aThree,
975  chose, zRand, rRand1, rRand2;
976 
977  G4double rOne = std::sqrt(sqr(xSemiAxis)
978  + sqr(ySemiAxis))*(zheight - zTopCut);
979  G4double rTwo = std::sqrt(sqr(xSemiAxis)
980  + sqr(ySemiAxis))*(zheight + zTopCut);
981 
982  G4int it1=0, it2=0;
983 
984  aOne = pi*(rOne + rTwo)*std::sqrt(sqr(rOne - rTwo)+sqr(2.*zTopCut));
985  aTwo = pi*xSemiAxis*ySemiAxis*sqr(zheight+zTopCut);
986  aThree = pi*xSemiAxis*ySemiAxis*sqr(zheight-zTopCut);
987 
988  phi = G4RandFlat::shoot(0.,twopi);
989  cosphi = std::cos(phi);
990  sinphi = std::sin(phi);
991 
992  if(zTopCut >= zheight) aThree = 0.;
993 
994  chose = G4RandFlat::shoot(0.,aOne+aTwo+aThree);
995  if((chose>=0.) && (chose<aOne))
996  {
997  zRand = G4RandFlat::shoot(-zTopCut,zTopCut);
998  return G4ThreeVector(xSemiAxis*(zheight-zRand)*cosphi,
999  ySemiAxis*(zheight-zRand)*sinphi,zRand);
1000  }
1001  else if((chose>=aOne) && (chose<aOne+aTwo))
1002  {
1003  do // Loop checking, 13.08.2015, G.Cosmo
1004  {
1005  rRand1 = G4RandFlat::shoot(0.,1.) ;
1006  rRand2 = G4RandFlat::shoot(0.,1.) ;
1007  } while (( rRand2 >= rRand1 ) && (++it1 < 1000)) ;
1008 
1009  return G4ThreeVector(rRand1*xSemiAxis*(zheight+zTopCut)*cosphi,
1010  rRand1*ySemiAxis*(zheight+zTopCut)*sinphi, -zTopCut);
1011 
1012  }
1013  // else
1014  //
1015 
1016  do // Loop checking, 13.08.2015, G.Cosmo
1017  {
1018  rRand1 = G4RandFlat::shoot(0.,1.) ;
1019  rRand2 = G4RandFlat::shoot(0.,1.) ;
1020  } while (( rRand2 >= rRand1 ) && (++it2 < 1000));
1021 
1022  return G4ThreeVector(rRand1*xSemiAxis*(zheight-zTopCut)*cosphi,
1023  rRand1*ySemiAxis*(zheight-zTopCut)*sinphi, zTopCut);
1024 }
ThreeVector shoot(const G4int Ap, const G4int Af)
CLHEP::Hep3Vector G4ThreeVector
int G4int
Definition: G4Types.hh:78
static constexpr double twopi
Definition: G4SIunits.hh:76
static constexpr double pi
Definition: G4SIunits.hh:75
T sqr(const T &x)
Definition: templates.hh:145
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

G4Polyhedron * G4EllipticalCone::GetPolyhedron ( ) const
virtual

Reimplemented from G4VSolid.

Definition at line 1051 of file G4EllipticalCone.cc.

1052 {
1053  if ( (!fpPolyhedron)
1057  {
1058  G4AutoLock l(&polyhedronMutex);
1059  delete fpPolyhedron;
1061  fRebuildPolyhedron = false;
1062  l.unlock();
1063  }
1064  return fpPolyhedron;
1065 }
G4Polyhedron * CreatePolyhedron() const
G4Polyhedron * fpPolyhedron
static G4int GetNumberOfRotationSteps()
G4int GetNumberOfRotationStepsAtTimeOfCreation() const

Here is the call graph for this function:

G4double G4EllipticalCone::GetSemiAxisMax ( ) const
inline
G4double G4EllipticalCone::GetSemiAxisX ( ) const
inline

Here is the caller graph for this function:

G4double G4EllipticalCone::GetSemiAxisY ( ) const
inline

Here is the caller graph for this function:

G4double G4EllipticalCone::GetSurfaceArea ( )
inlinevirtual

Reimplemented from G4VSolid.

G4double G4EllipticalCone::GetZMax ( ) const
inline

Here is the caller graph for this function:

G4double G4EllipticalCone::GetZTopCut ( ) const
inline

Here is the caller graph for this function:

EInside G4EllipticalCone::Inside ( const G4ThreeVector p) const
virtual

Implements G4VSolid.

Definition at line 264 of file G4EllipticalCone.cc.

265 {
266  G4double rad2oo, // outside surface outer tolerance
267  rad2oi; // outside surface inner tolerance
268 
269  EInside in;
270 
271  // check this side of z cut first, because that's fast
272  //
273 
274  if ( (p.z() < -zTopCut - halfCarTol)
275  || (p.z() > zTopCut + halfCarTol ) )
276  {
277  return in = kOutside;
278  }
279 
280  rad2oo= sqr(p.x()/( xSemiAxis + halfRadTol ))
281  + sqr(p.y()/( ySemiAxis + halfRadTol ));
282 
283  if ( rad2oo > sqr( zheight-p.z() ) )
284  {
285  return in = kOutside;
286  }
287 
288  // rad2oi= sqr( p.x()*(1.0 + 0.5*kRadTolerance/(xSemiAxis*xSemiAxis)) )
289  // + sqr( p.y()*(1.0 + 0.5*kRadTolerance/(ySemiAxis*ySemiAxis)) );
290  rad2oi = sqr(p.x()/( xSemiAxis - halfRadTol ))
291  + sqr(p.y()/( ySemiAxis - halfRadTol ));
292 
293  if (rad2oi < sqr( zheight-p.z() ) )
294  {
295  in = ( ( p.z() < -zTopCut + halfRadTol )
296  || ( p.z() > zTopCut - halfRadTol ) ) ? kSurface : kInside;
297  }
298  else
299  {
300  in = kSurface;
301  }
302 
303  return in;
304 }
double x() const
double z() const
EInside
Definition: geomdefs.hh:58
double y() const
T sqr(const T &x)
Definition: templates.hh:145
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

Here is the caller graph for this function:

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

Definition at line 146 of file G4EllipticalCone.cc.

147 {
148  // Check assignment to self
149  //
150  if (this == &rhs) { return *this; }
151 
152  // Copy base class data
153  //
154  G4VSolid::operator=(rhs);
155 
156  // Copy data
157  //
158  kRadTolerance = rhs.kRadTolerance;
159  halfRadTol = rhs.halfRadTol; halfCarTol = rhs.halfCarTol;
160  fCubicVolume = rhs.fCubicVolume; fSurfaceArea = rhs.fSurfaceArea;
161  xSemiAxis = rhs.xSemiAxis; ySemiAxis = rhs.ySemiAxis;
162  zheight = rhs.zheight; semiAxisMax = rhs.semiAxisMax; zTopCut = rhs.zTopCut;
163  fRebuildPolyhedron = false;
164  delete fpPolyhedron; fpPolyhedron = 0;
165 
166  return *this;
167 }
G4Polyhedron * fpPolyhedron
G4VSolid & operator=(const G4VSolid &rhs)
Definition: G4VSolid.cc:111

Here is the call graph for this function:

void G4EllipticalCone::SetSemiAxis ( G4double  x,
G4double  y,
G4double  z 
)
inline

Here is the caller graph for this function:

void G4EllipticalCone::SetZCut ( G4double  newzTopCut)
inline

Here is the caller graph for this function:

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

Implements G4VSolid.

Definition at line 946 of file G4EllipticalCone.cc.

947 {
948  G4int oldprc = os.precision(16);
949  os << "-----------------------------------------------------------\n"
950  << " *** Dump for solid - " << GetName() << " ***\n"
951  << " ===================================================\n"
952  << " Solid type: G4EllipticalCone\n"
953  << " Parameters: \n"
954 
955  << " semi-axis x: " << xSemiAxis/mm << " mm \n"
956  << " semi-axis y: " << ySemiAxis/mm << " mm \n"
957  << " height z: " << zheight/mm << " mm \n"
958  << " half length in z: " << zTopCut/mm << " mm \n"
959  << "-----------------------------------------------------------\n";
960  os.precision(oldprc);
961 
962  return os;
963 }
G4String GetName() const
static constexpr double mm
Definition: G4SIunits.hh:115
int G4int
Definition: G4Types.hh:78

Here is the call graph for this function:

G4ThreeVector G4EllipticalCone::SurfaceNormal ( const G4ThreeVector p) const
virtual

Implements G4VSolid.

Definition at line 310 of file G4EllipticalCone.cc.

311 {
312 
313  G4double rx = sqr(p.x()/xSemiAxis),
314  ry = sqr(p.y()/ySemiAxis);
315 
316  G4double rds = std::sqrt(rx + ry);
317 
318  G4ThreeVector norm;
319 
320  if( (p.z() < -zTopCut) && ((rx+ry) < sqr(zTopCut + zheight)) )
321  {
322  return G4ThreeVector( 0., 0., -1. );
323  }
324 
325  if( (p.z() > (zheight > zTopCut ? zheight : zTopCut)) &&
326  ((rx+ry) < sqr(zheight-zTopCut)) )
327  {
328  return G4ThreeVector( 0., 0., 1. );
329  }
330 
331  if( p.z() > rds + 2.*zTopCut - zheight )
332  {
333  if ( p.z() > zTopCut )
334  {
335  if( p.x() == 0. )
336  {
337  norm = G4ThreeVector( 0., p.y() < 0. ? -1. : 1., 1. );
338  return norm /= norm.mag();
339  }
340  if( p.y() == 0. )
341  {
342  norm = G4ThreeVector( p.x() < 0. ? -1. : 1., 0., 1. );
343  return norm /= norm.mag();
344  }
345 
346  G4double k = std::fabs(p.x()/p.y());
347  G4double c2 = sqr(zheight-zTopCut)/(1./sqr(xSemiAxis)+sqr(k/ySemiAxis));
348  G4double x = std::sqrt(c2);
349  G4double y = k*x;
350 
351  x /= sqr(xSemiAxis);
352  y /= sqr(ySemiAxis);
353 
354  norm = G4ThreeVector( p.x() < 0. ? -x : x,
355  p.y() < 0. ? -y : y,
356  - ( zheight - zTopCut ) );
357  norm /= norm.mag();
358  norm += G4ThreeVector( 0., 0., 1. );
359  return norm /= norm.mag();
360  }
361 
362  return G4ThreeVector( 0., 0., 1. );
363  }
364 
365  if( p.z() < rds - 2.*zTopCut - zheight )
366  {
367  if( p.x() == 0. )
368  {
369  norm = G4ThreeVector( 0., p.y() < 0. ? -1. : 1., -1. );
370  return norm /= norm.mag();
371  }
372  if( p.y() == 0. )
373  {
374  norm = G4ThreeVector( p.x() < 0. ? -1. : 1., 0., -1. );
375  return norm /= norm.mag();
376  }
377 
378  G4double k = std::fabs(p.x()/p.y());
379  G4double c2 = sqr(zheight+zTopCut)/(1./sqr(xSemiAxis)+sqr(k/ySemiAxis));
380  G4double x = std::sqrt(c2);
381  G4double y = k*x;
382 
383  x /= sqr(xSemiAxis);
384  y /= sqr(ySemiAxis);
385 
386  norm = G4ThreeVector( p.x() < 0. ? -x : x,
387  p.y() < 0. ? -y : y,
388  - ( zheight - zTopCut ) );
389  norm /= norm.mag();
390  norm += G4ThreeVector( 0., 0., -1. );
391  return norm /= norm.mag();
392  }
393 
394  norm = G4ThreeVector(p.x()/sqr(xSemiAxis), p.y()/sqr(ySemiAxis), rds);
395 
396  G4double k = std::tan(pi/8.);
397  G4double c = -zTopCut - k*(zTopCut + zheight);
398 
399  if( p.z() < -k*rds + c )
400  return G4ThreeVector (0.,0.,-1.);
401 
402  return norm /= norm.mag();
403 }
CLHEP::Hep3Vector G4ThreeVector
double x() const
double z() const
double y() const
static constexpr double pi
Definition: G4SIunits.hh:75
T sqr(const T &x)
Definition: templates.hh:145
double G4double
Definition: G4Types.hh:76
double mag() const

Here is the call graph for this function:

Member Data Documentation

G4Polyhedron* G4EllipticalCone::fpPolyhedron
mutableprotected

Definition at line 165 of file G4EllipticalCone.hh.

G4bool G4EllipticalCone::fRebuildPolyhedron
mutableprotected

Definition at line 164 of file G4EllipticalCone.hh.


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