66 using namespace CLHEP;
79 :
G4VSolid(pName), fRebuildPolyhedron(false), fpPolyhedron(0),
80 fCubicVolume(0.), fSurfaceArea(0.), zBottomCut(0.), zTopCut(0.)
91 if ( (pxSemiAxis<=0.) || (pySemiAxis<=0.) || (pzSemiAxis<=0.) )
93 std::ostringstream message;
94 message <<
"Invalid semi-axis - " <<
GetName();
95 G4Exception(
"G4Ellipsoid::G4Ellipsoid()",
"GeomSolids0002",
100 if ( pzBottomCut == 0 && pzTopCut == 0 )
104 else if ( (pzBottomCut < pzSemiAxis) && (pzTopCut > -pzSemiAxis)
105 && (pzBottomCut < pzTopCut) )
111 std::ostringstream message;
112 message <<
"Invalid z-coordinate for cutting plane - " <<
GetName();
113 G4Exception(
"G4Ellipsoid::G4Ellipsoid()",
"GeomSolids0002",
124 :
G4VSolid(a), fRebuildPolyhedron(false), fpPolyhedron(0), kRadTolerance(0.),
125 halfCarTolerance(0.), halfRadTolerance(0.), fCubicVolume(0.),
126 fSurfaceArea(0.), xSemiAxis(0.), ySemiAxis(0.), zSemiAxis(0.),
127 semiAxisMax(0.), zBottomCut(0.), zTopCut(0.)
146 fRebuildPolyhedron(false), fpPolyhedron(0),
147 kRadTolerance(rhs.kRadTolerance),
148 halfCarTolerance(rhs.halfCarTolerance),
149 halfRadTolerance(rhs.halfRadTolerance),
150 fCubicVolume(rhs.fCubicVolume), fSurfaceArea(rhs.fSurfaceArea),
151 xSemiAxis(rhs.xSemiAxis), ySemiAxis(rhs.ySemiAxis),
152 zSemiAxis(rhs.zSemiAxis), semiAxisMax(rhs.semiAxisMax),
153 zBottomCut(rhs.zBottomCut), zTopCut(rhs.zTopCut)
165 if (
this == &rhs) {
return *
this; }
209 pMin.set(-dx,-dy,zmin);
210 pMax.set( dx, dy,zmax);
214 if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
216 std::ostringstream message;
217 message <<
"Bad bounding box (min >= max) for solid: "
219 <<
"\npMin = " << pMin
220 <<
"\npMax = " << pMax;
267 if (rad2oo > 1.0) {
return in=
kOutside; }
296 G4double distR, distZBottom, distZTop;
308 distR = std::fabs( (p*norm - 1.0) * radius ) / 2.0;
312 distZBottom = std::fabs( p.z() -
zBottomCut );
313 distZTop = std::fabs( p.z() -
zTopCut );
315 if ( (distZBottom < distR) || (distZTop < distR) )
317 return G4ThreeVector(0.,0.,(distZBottom < distZTop) ? -1.0 : 1.0);
319 return ( norm *= radius );
338 if (v.z() <= 0.0) {
return distMin; }
345 return distMin= distZ;
350 if (v.z() >= 0.0) {
return distMin;}
356 return distMin= distZ;
373 G4double distR= (-B - std::sqrt(C)) / (2.0*A);
389 distR = (-B + std::sqrt(C) ) / (2.0*A);
390 if(distR>0.) { distMin=0.; }
394 distR= (-B + std::sqrt(C)) / (2.0*A);
395 intZ = p.z()+distR*v.z();
401 if (norm.dot(v)<0.) { distMin = distR; }
404 if ( (distMin!=
kInfinity) && (distMin>dRmax) )
407 G4double fTerm = distMin-std::fmod(distMin,dRmax);
434 distR= (p*norm - 1.0) * radius / 2.0;
448 return (distR < 0.0) ? 0.0 : distR;
450 else if (distR < 0.0)
456 return (distZ < distR) ? distZ : distR;
471 enum surface_e {kPlaneSurf, kCurvedSurf, kNoSurf} surface;
484 if (!calcNorm) {
return 0.0;}
495 if (!calcNorm) {
return 0.0;}
512 C= (p * nearnorm) - 1.0;
513 B= 2.0 * (v * nearnorm);
518 G4double distR= (-B + std::sqrt(C) ) / (2.0*A);
522 if (!calcNorm) {
return 0.0;}
527 surface= kCurvedSurf;
535 if (surface == kNoSurf)
553 truenorm *= 1.0/truenorm.mag();
558 std::ostringstream message;
559 G4int oldprc = message.precision(16);
560 message <<
"Undefined side for valid surface normal to solid."
563 <<
" p.x() = " << p.x()/
mm <<
" mm" <<
G4endl
564 <<
" p.y() = " << p.y()/
mm <<
" mm" <<
G4endl
565 <<
" p.z() = " << p.z()/
mm <<
" mm" <<
G4endl
567 <<
" v.x() = " << v.x() <<
G4endl
568 <<
" v.y() = " << v.y() <<
G4endl
569 <<
" v.z() = " << v.z() <<
G4endl
570 <<
"Proposed distance :" <<
G4endl
571 <<
" distMin = " << distMin/
mm <<
" mm";
572 message.precision(oldprc);
595 std::ostringstream message;
596 G4int oldprc = message.precision(16);
597 message <<
"Point p is outside !?" <<
G4endl
599 <<
" p.x() = " << p.x()/
mm <<
" mm" <<
G4endl
600 <<
" p.y() = " << p.y()/
mm <<
" mm" <<
G4endl
601 <<
" p.z() = " << p.z()/
mm <<
" mm";
602 message.precision(oldprc) ;
603 G4Exception(
"G4Ellipsoid::DistanceToOut(p)",
"GeomSolids1002",
618 if ( (tmp > 0.0) && (1.0 < radius*tmp) ) {radius = 1.0/tmp;}
622 distR = (1.0 - p*norm) * radius / 2.0;
627 if (distZ < 0.0) {distZ=
zTopCut - p.z();}
631 if ( (distZ < 0.0) || (distR < 0.0) )
637 return (distZ < distR) ? distZ : distR;
665 G4int oldprc = os.precision(16);
666 os <<
"-----------------------------------------------------------\n"
667 <<
" *** Dump for solid - " <<
GetName() <<
" ***\n"
668 <<
" ===================================================\n"
669 <<
" Solid type: G4Ellipsoid\n"
676 <<
" lower cut plane level z: " <<
zBottomCut/
mm <<
" mm \n"
677 <<
" upper cut plane level z: " <<
zTopCut/
mm <<
" mm \n"
678 <<
"-----------------------------------------------------------\n";
679 os.precision(oldprc);
690 G4double aTop, aBottom, aCurved, chose, xRand, yRand, zRand, phi;
701 cosphi = std::cos(phi); sinphi = std::sin(phi);
703 sintheta = std::sqrt(1.-
sqr(costheta));
705 alpha = 1.-
sqr(max2/max1); beta = 1.-
sqr(max3/max1);
712 aCurved = 4.*
pi*max1*max2*(1.-1./6.*(alpha+beta)-
713 1./120.*(3.*
sqr(alpha)+2.*alpha*beta+3.*
sqr(beta)));
720 aTop = 0; aBottom = 0;
732 else if(chose >= aCurved && chose < aCurved + aTop)
G4Polyhedron * GetPolyhedron() const
ThreeVector shoot(const G4int Ap, const G4int Af)
G4double halfRadTolerance
G4GeometryType GetEntityType() const
G4bool fRebuildPolyhedron
static constexpr double mm
static const G4double kInfinity
CLHEP::Hep3Vector G4ThreeVector
std::vector< ExP01TrackerHit * > a
void Extent(G4ThreeVector &pMin, G4ThreeVector &pMax) const
G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=G4bool(false), G4bool *validNorm=0, G4ThreeVector *n=0) const
G4double GetZTopCut() const
EInside Inside(const G4ThreeVector &p) const
G4double GetZBottomCut() const
double B(double temperature)
virtual void AddSolid(const G4Box &)=0
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pmin, G4double &pmax) const
G4double GetSemiAxisMax(G4int i) const
#define G4MUTEX_INITIALIZER
static constexpr double twopi
void SetZCuts(G4double newzBottomCut, G4double newzTopCut)
void DescribeYourselfTo(G4VGraphicsScene &scene) const
G4Ellipsoid(const G4String &pName, G4double pxSemiAxis, G4double pySemiAxis, G4double pzSemiAxis, G4double pzBottomCut=0, G4double pzTopCut=0)
double A(double temperature)
G4VisExtent GetExtent() const
G4Polyhedron * fpPolyhedron
G4double GetRadialTolerance() const
G4Ellipsoid & operator=(const G4Ellipsoid &rhs)
G4Polyhedron * CreatePolyhedron() const
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimits, const G4Transform3D &pTransform3D, G4double &pMin, G4double &pMax) const
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const
G4double halfCarTolerance
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
virtual void ComputeDimensions(G4Box &, const G4int, const G4VPhysicalVolume *) const
std::ostream & StreamInfo(std::ostream &os) const
T max(const T t1, const T t2)
brief Return the largest of the two arguments
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
void SetSemiAxis(G4double x, G4double y, G4double z)
G4VSolid & operator=(const G4VSolid &rhs)
G4int GetNumberOfRotationStepsAtTimeOfCreation() const
static constexpr double pi
G4ThreeVector GetPointOnSurface() const
static const G4double alpha
void ComputeDimensions(G4VPVParameterisation *p, const G4int n, const G4VPhysicalVolume *pRep)
static G4GeometryTolerance * GetInstance()
static int max3(int a, int b, int c)