65 using namespace CLHEP;
78 :
G4VSolid(pName), fRebuildPolyhedron(false), fpPolyhedron(0),
79 fCubicVolume(0.), fSurfaceArea(0.), zBottomCut(0.), zTopCut(0.)
90 if ( (pxSemiAxis<=0.) || (pySemiAxis<=0.) || (pzSemiAxis<=0.) )
93 message <<
"Invalid semi-axis - " <<
GetName();
94 G4Exception(
"G4Ellipsoid::G4Ellipsoid()",
"GeomSolids0002",
99 if ( pzBottomCut == 0 && pzTopCut == 0 )
103 else if ( (pzBottomCut < pzSemiAxis) && (pzTopCut > -pzSemiAxis)
104 && (pzBottomCut < pzTopCut) )
111 message <<
"Invalid z-coordinate for cutting plane - " <<
GetName();
112 G4Exception(
"G4Ellipsoid::G4Ellipsoid()",
"GeomSolids0002",
123 :
G4VSolid(a), fRebuildPolyhedron(false), fpPolyhedron(0), kRadTolerance(0.),
124 halfCarTolerance(0.), halfRadTolerance(0.), fCubicVolume(0.),
125 fSurfaceArea(0.), xSemiAxis(0.), ySemiAxis(0.), zSemiAxis(0.),
126 semiAxisMax(0.), zBottomCut(0.), zTopCut(0.)
145 fRebuildPolyhedron(false), fpPolyhedron(0),
146 kRadTolerance(rhs.kRadTolerance),
147 halfCarTolerance(rhs.halfCarTolerance),
148 halfRadTolerance(rhs.halfRadTolerance),
149 fCubicVolume(rhs.fCubicVolume), fSurfaceArea(rhs.fSurfaceArea),
150 xSemiAxis(rhs.xSemiAxis), ySemiAxis(rhs.ySemiAxis),
151 zSemiAxis(rhs.zSemiAxis), semiAxisMax(rhs.semiAxisMax),
152 zBottomCut(rhs.zBottomCut), zTopCut(rhs.zTopCut)
164 if (
this == &rhs) {
return *
this; }
208 pMin.
set(-dx,-dy,zmin);
209 pMax.
set( dx, dy,zmax);
213 if (pMin.
x() >= pMax.
x() || pMin.
y() >= pMax.
y() || pMin.
z() >= pMax.
z())
216 message <<
"Bad bounding box (min >= max) for solid: "
218 <<
"\npMin = " << pMin
219 <<
"\npMax = " << pMax;
220 G4Exception(
"G4Ellipsoid::BoundingLimits()",
"GeomMgt0001",
267 if (rad2oo > 1.0) {
return in=
kOutside; }
296 G4double distR, distZBottom, distZTop;
308 distR = std::fabs( (p*norm - 1.0) * radius ) / 2.0;
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();
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);
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)
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimits, const G4Transform3D &pTransform3D, G4double &pMin, G4double &pMax) const
void set(double x, double y, double z)
T max(const T t1, const T t2)
brief Return the largest of the two arguments
CLHEP::Hep3Vector G4ThreeVector
G4double GetSemiAxisMax(G4int i) const
std::vector< ExP01TrackerHit * > a
virtual void AddSolid(const G4Box &)=0
static const G4double kInfinity
void BoundingLimits(G4ThreeVector &pMin, G4ThreeVector &pMax) const
G4bool fRebuildPolyhedron
G4GeometryType GetEntityType() const
static constexpr double mm
virtual void ComputeDimensions(G4Box &, const G4int, const G4VPhysicalVolume *) const
void message(RunManager *runmanager)
double dot(const Hep3Vector &) const
G4double GetRadialTolerance() const
G4VisExtent GetExtent() const
std::ostream & StreamInfo(std::ostream &os) const
void SetSemiAxis(G4double x, G4double y, G4double z)
G4VSolid & operator=(const G4VSolid &rhs)
#define G4MUTEX_INITIALIZER
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pmin, G4double &pmax) const
static int max3(int a, int b, int c)
G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const
G4double halfCarTolerance
G4ThreeVector GetPointOnSurface() const
G4double halfRadTolerance
G4Polyhedron * GetPolyhedron() const
static const G4double alpha
static constexpr double twopi
double A(double temperature)
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
G4double GetZBottomCut() const
void SetZCuts(G4double newzBottomCut, G4double newzTopCut)
ThreeVector shoot(const G4int Ap, const G4int Af)
static G4int GetNumberOfRotationSteps()
G4Ellipsoid(const G4String &pName, G4double pxSemiAxis, G4double pySemiAxis, G4double pzSemiAxis, G4double pzBottomCut=0, G4double pzTopCut=0)
EInside Inside(const G4ThreeVector &p) const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
G4int GetNumberOfRotationStepsAtTimeOfCreation() const
G4Polyhedron * fpPolyhedron
G4Polyhedron * CreatePolyhedron() const
void DescribeYourselfTo(G4VGraphicsScene &scene) const
G4double GetZTopCut() const
static G4GeometryTolerance * GetInstance()
void ComputeDimensions(G4VPVParameterisation *p, const G4int n, const G4VPhysicalVolume *pRep)
static constexpr double pi
G4Ellipsoid & operator=(const G4Ellipsoid &rhs)
G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=G4bool(false), G4bool *validNorm=0, G4ThreeVector *n=0) const
double B(double temperature)
T min(const T t1, const T t2)
brief Return the smallest of the two arguments