64 using namespace CLHEP;
75 :
G4VSolid(pName), fRebuildPolyhedron(false), fpPolyhedron(0),
76 fCubicVolume(0.), fSurfaceArea(0.), zTopCut(0.)
86 if ( (pxSemiAxis <= 0.) || (pySemiAxis <= 0.) || (pzMax <= 0.) )
88 std::ostringstream message;
89 message <<
"Invalid semi-axis or height - " <<
GetName();
90 G4Exception(
"G4EllipticalCone::G4EllipticalCone()",
"GeomSolids0002",
95 std::ostringstream message;
96 message <<
"Invalid z-coordinate for cutting plane - " <<
GetName();
97 G4Exception(
"G4EllipticalCone::G4EllipticalCone()",
"InvalidSetup",
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.)
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)
150 if (
this == &rhs) {
return *
this; }
198 G4double rFudge = 1.0/std::cos(0.5*sigPhi);
212 cosPhi = std::cos(phi),
213 sinPhi = std::sin(phi);
215 v1( dxFudgeBot*cosPhi, dyFudgeBot*sinPhi, -
zTopCut ),
222 if (numPhi == 1) phi = 0;
223 cosPhi = std::cos(phi),
224 sinPhi = std::sin(phi);
252 phiPoly.
SetNormal( (v1-v0).cross(w0-v0).unit() );
262 }
while( --numPhi > 0 );
367 return norm /= norm.mag();
372 return norm /= norm.mag();
375 G4double k = std::fabs(p.x()/p.y());
388 return norm /= norm.mag();
399 return norm /= norm.mag();
404 return norm /= norm.mag();
407 G4double k = std::fabs(p.x()/p.y());
420 return norm /= norm.mag();
428 if( p.z() < -k*rds + c )
431 return norm /= norm.mag();
487 yi = p.y() + q*v.y();
532 yi = p.y() + q*v.y();
562 return distMin = std::fabs(lambda);
577 return distMin = std::fabs(lambda);
621 return distMin = std::fabs(-B/(2.*A));
624 G4double plus = (-B+std::sqrt(discr))/(2.*A);
625 G4double minus = (-B-std::sqrt(discr))/(2.*A);
634 if ( truenorm*v >= 0)
690 G4double distR, distR2, distZ, maxDim;
700 return distZ = std::fabs(
zTopCut + p.z());
706 return distZ = std::fabs(p.z() -
zTopCut);
714 distRad = std::sqrt(p.x()*p.x()+p.y()*p.y());
719 return std::sqrt( distR2 );
722 if( distRad > maxDim*(
zheight - p.z() ) )
728 return distR = std::sqrt(
sqr(p.z() - zVal) +
sqr(distRad - rVal));
732 if( distRad <= maxDim*(
zheight - p.z()) )
735 return std::sqrt( distR2 );
753 enum surface_e {kPlaneSurf, kCurvedSurf, kNoSurf} surface;
760 lambda = (-p.z() -
zTopCut)/v.z();
766 distMin = std::fabs(lambda);
768 if (!calcNorm) {
return distMin; }
770 distMin = std::fabs(lambda);
771 surface = kPlaneSurf;
776 lambda = (
zTopCut - p.z()) / v.z();
782 distMin = std::fabs(lambda);
783 if (!calcNorm) {
return distMin; }
785 distMin = std::fabs(lambda);
786 surface = kPlaneSurf;
802 if(!calcNorm) {
return distMin = std::fabs(-B/(2.*A)); }
807 G4double plus = (-B+std::sqrt(discr))/(2.*A);
808 G4double minus = (-B-std::sqrt(discr))/(2.*A);
814 lambda = std::fabs(plus) < std::fabs(minus) ? plus : minus;
824 if ( std::fabs(lambda) < distMin )
828 distMin = std::fabs(lambda);
829 surface = kCurvedSurf;
836 if( truenorm.dot(v) > 0 )
839 surface = kCurvedSurf;
849 if (surface == kNoSurf)
870 truenorm /= truenorm.mag();
877 std::ostringstream message;
878 G4int oldprc = message.precision(16);
879 message <<
"Undefined side for valid surface normal to solid."
882 <<
" p.x() = " << p.x()/
mm <<
" mm" <<
G4endl
883 <<
" p.y() = " << p.y()/
mm <<
" mm" <<
G4endl
884 <<
" p.z() = " << p.z()/
mm <<
" mm" <<
G4endl
886 <<
" v.x() = " << v.x() <<
G4endl
887 <<
" v.y() = " << v.y() <<
G4endl
888 <<
" v.z() = " << v.z() <<
G4endl
889 <<
"Proposed distance :" <<
G4endl
890 <<
" distMin = " << distMin/
mm <<
" mm";
891 message.precision(oldprc);
892 G4Exception(
"G4EllipticalCone::DistanceToOut(p,v,..)",
910 G4double rds,roo,roo1, distR, distZ, distMin=0.;
917 std::ostringstream message;
918 G4int oldprc = message.precision(16);
919 message <<
"Point p is outside !?" <<
G4endl
921 <<
" p.x() = " << p.x()/
mm <<
" mm" <<
G4endl
922 <<
" p.y() = " << p.y()/
mm <<
" mm" <<
G4endl
923 <<
" p.z() = " << p.z()/
mm <<
" mm";
924 message.precision(oldprc) ;
925 G4Exception(
"G4Ellipsoid::DistanceToOut(p)",
"GeomSolids1002",
936 rds = std::sqrt(
sqr(p.x()) +
sqr(p.y()));
940 distZ=
zTopCut - std::fabs(p.z()) ;
941 distR=(roo-rds)/(std::sqrt(1+
sqr(minAxis)));
945 distMin=(
zTopCut-p.z())*(roo-rds)/(roo-roo1);
960 return G4String(
"G4EllipticalCone");
978 G4int oldprc = os.precision(16);
979 os <<
"-----------------------------------------------------------\n"
980 <<
" *** Dump for solid - " <<
GetName() <<
" ***\n"
981 <<
" ===================================================\n"
982 <<
" Solid type: G4EllipticalCone\n"
987 <<
" height z: " <<
zheight/
mm <<
" mm \n"
988 <<
" half length in z: " <<
zTopCut/
mm <<
" mm \n"
989 <<
"-----------------------------------------------------------\n";
990 os.precision(oldprc);
1004 G4double phi, sinphi, cosphi, aOne, aTwo, aThree,
1005 chose, zRand, rRand1, rRand2;
1019 cosphi = std::cos(phi);
1020 sinphi = std::sin(phi);
1025 if((chose>=0.) && (chose<aOne))
1031 else if((chose>=aOne) && (chose<aOne+aTwo))
1037 }
while (( rRand2 >= rRand1 ) && (++it1 < 1000)) ;
1051 }
while (( rRand2 >= rRand1 ) && (++it2 < 1000));
ThreeVector shoot(const G4int Ap, const G4int Af)
static c2_factory< G4double > c2
G4EllipticalCone & operator=(const G4EllipticalCone &rhs)
static const G4double kInfinity
CLHEP::Hep3Vector G4ThreeVector
G4bool GetExtent(G4double &min, G4double &max) const
EInside Inside(const G4ThreeVector &p) const
void SetNormal(const G4ThreeVector &newNormal)
virtual G4bool PartialClip(const G4VoxelLimits &voxelLimit, const EAxis IgnoreMe)
virtual void AddVertexInOrder(const G4ThreeVector vertex)
virtual void ClearAllVertices()
double B(double temperature)
virtual void AddSolid(const G4Box &)=0
G4Polyhedron * CreatePolyhedron() const
#define G4MUTEX_INITIALIZER
static double normal(HepRandomEngine *eptr)
std::ostream & StreamInfo(std::ostream &os) const
G4bool fRebuildPolyhedron
G4ThreeVector GetPointOnSurface() const
double A(double temperature)
G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const
G4Polyhedron * fpPolyhedron
G4double GetRadialTolerance() const
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
G4EllipticalCone(const G4String &pName, G4double pxSemiAxis, G4double pySemiAxis, G4double zMax, G4double pzTopCut)
static const double twopi
void AddSurface(const G4ClippablePolygon &surface)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
void DescribeYourselfTo(G4VGraphicsScene &scene) const
void SetSemiAxis(G4double x, G4double y, G4double z)
T max(const T t1, const T t2)
brief Return the largest of the two arguments
const G4double x[NPOINTSGL]
void SetZCut(G4double newzTopCut)
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pmin, G4double &pmax) const
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
G4VSolid & operator=(const G4VSolid &rhs)
G4int GetNumberOfRotationStepsAtTimeOfCreation() const
G4Polyhedron * GetPolyhedron() const
G4GeometryType GetEntityType() const
G4VisExtent GetExtent() const
virtual ~G4EllipticalCone()
G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=G4bool(false), G4bool *validNorm=0, G4ThreeVector *n=0) const
static G4GeometryTolerance * GetInstance()
const G4int kMaxMeshSections