67 using namespace CLHEP;
78 :
G4VSolid(pName), fRebuildPolyhedron(false), fpPolyhedron(0),
79 fCubicVolume(0.), fSurfaceArea(0.), zTopCut(0.)
85 if ( (pxSemiAxis <= 0.) || (pySemiAxis <= 0.) || (pzMax <= 0.) )
88 message <<
"Invalid semi-axis or height for solid: " <<
GetName()
89 <<
"\n X semi-axis, Y semi-axis, height = "
90 << pxSemiAxis <<
", " << pySemiAxis <<
", " << pzMax;
91 G4Exception(
"G4EllipticalCone::G4EllipticalCone()",
"GeomSolids0002",
98 message <<
"Invalid z-coordinate for cutting plane for solid: " <<
GetName()
99 <<
"\n Z top cut = " << pzTopCut;
100 G4Exception(
"G4EllipticalCone::G4EllipticalCone()",
"GeomSolids0002",
114 :
G4VSolid(a), fRebuildPolyhedron(false), fpPolyhedron(0),
115 halfCarTol(0.), fCubicVolume(0.), fSurfaceArea(0.),
116 xSemiAxis(0.), ySemiAxis(0.), zheight(0.), zTopCut(0.),
117 cosAxisMin(0.), invXX(0.), invYY(0.)
135 :
G4VSolid(rhs), fRebuildPolyhedron(false), fpPolyhedron(0),
136 halfCarTol(rhs.halfCarTol),
137 fCubicVolume(rhs.fCubicVolume), fSurfaceArea(rhs.fSurfaceArea),
138 xSemiAxis(rhs.xSemiAxis), ySemiAxis(rhs.ySemiAxis),
139 zheight(rhs.zheight), zTopCut(rhs.zTopCut),
140 cosAxisMin(rhs.cosAxisMin), invXX(rhs.invXX), invYY(rhs.invYY)
152 if (
this == &rhs) {
return *
this; }
183 pMin.
set(-xmax,-ymax,-zcut);
184 pMax.
set( xmax, ymax, zcut);
188 if (pMin.
x() >= pMax.
x() || pMin.
y() >= pMax.
y() || pMin.
z() >= pMax.
z())
191 message <<
"Bad bounding box (min >= max) for solid: "
193 <<
"\npMin = " << pMin
194 <<
"\npMax = " << pMax;
195 G4Exception(
"G4EllipticalCone::BoundingLimits()",
"GeomMgt0001",
219 if (
true)
return bbox.
CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
223 return exist = (pMin < pMax) ?
true :
false;
228 static const G4int NSTEPS = 48;
230 static const G4double sinHalf = std::sin(0.5*ang);
231 static const G4double cosHalf = std::cos(0.5*ang);
232 static const G4double sinStep = 2.*sinHalf*cosHalf;
233 static const G4double cosStep = 1. - 2.*sinHalf*sinHalf;
244 for (
G4int k=0; k<NSTEPS; ++k)
246 baseA[k].set(sxmax*cosCur,symax*sinCur,-zcut);
247 baseB[k].set(sxmin*cosCur,symin*sinCur, zcut);
250 sinCur = sinCur*cosStep + cosCur*sinStep;
251 cosCur = cosCur*cosStep - sinTmp*sinStep;
254 std::vector<const G4ThreeVectorList *> polygons(2);
255 polygons[0] = &baseA;
256 polygons[1] = &baseB;
303 if (nsurf == 1)
return norm;
304 else if (nsurf > 1)
return norm.
unit();
311 G4int oldprc = message.precision(16);
312 message <<
"Point p is not on surface (!?) of solid: "
314 message <<
"Position:\n";
315 message <<
" p.x() = " << p.
x()/
mm <<
" mm\n";
316 message <<
" p.y() = " << p.
y()/
mm <<
" mm\n";
317 message <<
" p.z() = " << p.
z()/
mm <<
" mm";
319 G4Exception(
"G4EllipticalCone::SurfaceNormal(p)",
"GeomSolids1002",
398 yi = p.
y() + q*v.
y();
443 yi = p.
y() + q*v.
y();
473 return distMin = std::fabs(lambda);
488 return distMin = std::fabs(lambda);
532 return distMin = std::fabs(-B/(2.*A));
535 G4double plus = (-B+std::sqrt(discr))/(2.*A);
536 G4double minus = (-B-std::sqrt(discr))/(2.*A);
545 if ( truenorm*v >= 0)
605 return (dist > 0) ? dist : 0.;
620 enum surface_e {kPlaneSurf, kCurvedSurf, kNoSurf} surface;
633 distMin = std::fabs(lambda);
635 if (!calcNorm) {
return distMin; }
637 distMin = std::fabs(lambda);
638 surface = kPlaneSurf;
649 distMin = std::fabs(lambda);
650 if (!calcNorm) {
return distMin; }
652 distMin = std::fabs(lambda);
653 surface = kPlaneSurf;
669 if(!calcNorm) {
return distMin = std::fabs(-B/(2.*A)); }
674 G4double plus = (-B+std::sqrt(discr))/(2.*A);
675 G4double minus = (-B-std::sqrt(discr))/(2.*A);
681 lambda = std::fabs(plus) < std::fabs(minus) ? plus : minus;
691 if ( std::fabs(lambda) < distMin )
695 distMin = std::fabs(lambda);
696 surface = kCurvedSurf;
703 if( truenorm.
dot(v) > 0 )
706 surface = kCurvedSurf;
716 if (surface == kNoSurf)
737 truenorm /= truenorm.
mag();
745 G4int oldprc = message.precision(16);
746 message <<
"Undefined side for valid surface normal to solid."
749 <<
" p.x() = " << p.
x()/
mm <<
" mm" <<
G4endl
750 <<
" p.y() = " << p.
y()/
mm <<
" mm" <<
G4endl
751 <<
" p.z() = " << p.
z()/
mm <<
" mm" <<
G4endl
753 <<
" v.x() = " << v.
x() <<
G4endl
754 <<
" v.y() = " << v.
y() <<
G4endl
755 <<
" v.z() = " << v.
z() <<
G4endl
756 <<
"Proposed distance :" <<
G4endl
757 <<
" distMin = " << distMin/
mm <<
" mm";
758 message.precision(oldprc);
759 G4Exception(
"G4EllipticalCone::DistanceToOut(p,v,..)",
781 G4int oldprc = message.precision(16);
782 message <<
"Point p is outside (!?) of solid: " <<
GetName() <<
"\n"
784 <<
" p.x() = " << p.
x()/
mm <<
" mm\n"
785 <<
" p.y() = " << p.
y()/
mm <<
" mm\n"
786 <<
" p.z() = " << p.
z()/
mm <<
" mm";
787 message.precision(oldprc) ;
788 G4Exception(
"G4Ellipsoid::DistanceToOut(p)",
"GeomSolids1002",
797 return (dist > 0) ? dist : 0.;
806 return G4String(
"G4EllipticalCone");
824 G4int oldprc = os.precision(16);
825 os <<
"-----------------------------------------------------------\n"
826 <<
" *** Dump for solid - " <<
GetName() <<
" ***\n"
827 <<
" ===================================================\n"
828 <<
" Solid type: G4EllipticalCone\n"
833 <<
" height z: " <<
zheight/
mm <<
" mm \n"
834 <<
" half length in z: " <<
zTopCut/
mm <<
" mm \n"
835 <<
"-----------------------------------------------------------\n";
836 os.precision(oldprc);
857 G4double sside = s0*(kmax*kmax - kmin*kmin);
858 G4double ssurf[3] = { szmin, sside, szmax };
859 for (
G4int i=1; i<3; ++i) { ssurf[i] += ssurf[i-1]; }
865 if (select <= ssurf[1]) k = 1;
866 if (select <= ssurf[0]) k = 0;
890 G4double mu_max = R*std::sqrt(hh + R*R);
893 for (
G4int i=0; i<1000; ++i)
933 fCubicVolume = (kmax - kmin)*(kmax*kmax + kmax*kmin + kmin*kmin)*v0;
952 +
CLHEP::pi*x0*y0*(kmin*kmin + kmax*kmax);
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
virtual ~G4EllipticalCone()
CLHEP::Hep3Vector G4ThreeVector
std::vector< ExP01TrackerHit * > a
virtual void AddSolid(const G4Box &)=0
static const G4double kInfinity
static constexpr double mm
void message(RunManager *runmanager)
G4EllipticalCone & operator=(const G4EllipticalCone &rhs)
double dot(const Hep3Vector &) const
G4double GetSurfaceArea()
G4VisExtent GetExtent() const
G4double GetSemiAxisY() const
G4VSolid & operator=(const G4VSolid &rhs)
G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=G4bool(false), G4bool *validNorm=0, G4ThreeVector *n=0) const
G4bool BoundingBoxVsVoxelLimits(const EAxis pAxis, const G4VoxelLimits &pVoxelLimits, const G4Transform3D &pTransform3D, G4double &pMin, G4double &pMax) const
std::vector< G4ThreeVector > G4ThreeVectorList
#define G4MUTEX_INITIALIZER
void DescribeYourselfTo(G4VGraphicsScene &scene) const
void SetZCut(G4double newzTopCut)
EInside Inside(const G4ThreeVector &p) const
G4ThreeVector GetPointOnSurface() const
G4double GetZTopCut() const
static constexpr double twopi
double A(double temperature)
G4bool fRebuildPolyhedron
G4Polyhedron * fpPolyhedron
G4EllipticalCone(const G4String &pName, G4double pxSemiAxis, G4double pySemiAxis, G4double zMax, G4double pzTopCut)
std::ostream & StreamInfo(std::ostream &os) const
G4Polyhedron * CreatePolyhedron() const
static G4int GetNumberOfRotationSteps()
static constexpr double twopi
G4GeometryType GetEntityType() const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
G4double GetSemiAxisX() const
G4int GetNumberOfRotationStepsAtTimeOfCreation() const
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const
G4GLOB_DLL std::ostream G4cout
G4double GetCubicVolume()
void SetSemiAxis(G4double x, G4double y, G4double z)
G4ThreeVector ApproxSurfaceNormal(const G4ThreeVector &p) const
static constexpr double pi
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pMin, G4double &pMax) const
double B(double temperature)
void BoundingLimits(G4ThreeVector &pMin, G4ThreeVector &pMax) const
static constexpr double pi
G4Polyhedron * GetPolyhedron() const
T min(const T t1, const T t2)
brief Return the smallest of the two arguments