62 using namespace CLHEP;
75 :
G4VSolid(pName), fRebuildPolyhedron(false), fpPolyhedron(0),
76 fCubicVolume(0.), fSurfaceArea(0.), zBottomCut(0.), zTopCut(0.)
87 if ( (pxSemiAxis<=0.) || (pySemiAxis<=0.) || (pzSemiAxis<=0.) )
89 std::ostringstream message;
90 message <<
"Invalid semi-axis - " <<
GetName();
91 G4Exception(
"G4Ellipsoid::G4Ellipsoid()",
"GeomSolids0002",
96 if ( pzBottomCut == 0 && pzTopCut == 0 )
100 else if ( (pzBottomCut < pzSemiAxis) && (pzTopCut > -pzSemiAxis)
101 && (pzBottomCut < pzTopCut) )
107 std::ostringstream message;
108 message <<
"Invalid z-coordinate for cutting plane - " <<
GetName();
109 G4Exception(
"G4Ellipsoid::G4Ellipsoid()",
"GeomSolids0002",
120 :
G4VSolid(a), fRebuildPolyhedron(false), fpPolyhedron(0), kRadTolerance(0.),
121 halfCarTolerance(0.), halfRadTolerance(0.), fCubicVolume(0.),
122 fSurfaceArea(0.), xSemiAxis(0.), ySemiAxis(0.), zSemiAxis(0.),
123 semiAxisMax(0.), zBottomCut(0.), zTopCut(0.)
142 fRebuildPolyhedron(false), fpPolyhedron(0),
143 kRadTolerance(rhs.kRadTolerance),
144 halfCarTolerance(rhs.halfCarTolerance),
145 halfRadTolerance(rhs.halfRadTolerance),
146 fCubicVolume(rhs.fCubicVolume), fSurfaceArea(rhs.fSurfaceArea),
147 xSemiAxis(rhs.xSemiAxis), ySemiAxis(rhs.ySemiAxis),
148 zSemiAxis(rhs.zSemiAxis), semiAxisMax(rhs.semiAxisMax),
149 zBottomCut(rhs.zBottomCut), zTopCut(rhs.zTopCut)
161 if (
this == &rhs) {
return *
this; }
289 xoff = (xoffset < xMin) ? (xMin-xoffset)
290 : (xoffset > xMax) ? (xoffset-xMax) : 0.0;
291 yoff = (yoffset < yMin) ? (yMin-yoffset)
292 : (yoffset > yMax) ? (yoffset-yMax) : 0.0;
314 maxDiff= 1.0-
sqr(yoff/ySemiAxis);
315 if (maxDiff < 0.0) {
return false; }
316 maxDiff= xSemiAxis * std::sqrt(maxDiff);
317 newMin=xoffset-maxDiff;
318 newMax=xoffset+maxDiff;
319 pMin=(newMin<xMin) ? xMin : newMin;
320 pMax=(newMax>xMax) ? xMax : newMax;
336 maxDiff= 1.0-
sqr(xoff/xSemiAxis);
337 if (maxDiff < 0.0) {
return false; }
338 maxDiff= ySemiAxis * std::sqrt(maxDiff);
339 newMin=yoffset-maxDiff;
340 newMax=yoffset+maxDiff;
341 pMin=(newMin<yMin) ? yMin : newMin;
342 pMax=(newMax>yMax) ? yMax : newMax;
359 G4int i,j,noEntries,noBetweenSections;
360 G4bool existsAfterClip=
false;
364 G4int noPolygonVertices=0;
371 noEntries=vertices->size();
372 noBetweenSections=noEntries-noPolygonVertices;
375 for (i=0;i<noEntries;i+=noPolygonVertices)
377 for(j=0;j<(noPolygonVertices/2)-1;j++)
379 ThetaPolygon.push_back((*vertices)[i+j]);
380 ThetaPolygon.push_back((*vertices)[i+j+1]);
381 ThetaPolygon.push_back((*vertices)[i+noPolygonVertices-2-j]);
382 ThetaPolygon.push_back((*vertices)[i+noPolygonVertices-1-j]);
384 ThetaPolygon.clear();
387 for (i=0;i<noBetweenSections;i+=noPolygonVertices)
389 for(j=0;j<noPolygonVertices-1;j++)
391 ThetaPolygon.push_back((*vertices)[i+j]);
392 ThetaPolygon.push_back((*vertices)[i+j+1]);
393 ThetaPolygon.push_back((*vertices)[i+noPolygonVertices+j+1]);
394 ThetaPolygon.push_back((*vertices)[i+noPolygonVertices+j]);
396 ThetaPolygon.clear();
398 ThetaPolygon.push_back((*vertices)[i+noPolygonVertices-1]);
399 ThetaPolygon.push_back((*vertices)[i]);
400 ThetaPolygon.push_back((*vertices)[i+noPolygonVertices]);
401 ThetaPolygon.push_back((*vertices)[i+2*noPolygonVertices-1]);
403 ThetaPolygon.clear();
407 existsAfterClip=
true;
429 existsAfterClip=
true;
435 return existsAfterClip;
460 if (rad2oo > 1.0) {
return in=
kOutside; }
489 G4double distR, distZBottom, distZTop;
501 distR = std::fabs( (p*norm - 1.0) * radius ) / 2.0;
505 distZBottom = std::fabs( p.z() -
zBottomCut );
506 distZTop = std::fabs( p.z() -
zTopCut );
508 if ( (distZBottom < distR) || (distZTop < distR) )
510 return G4ThreeVector(0.,0.,(distZBottom < distZTop) ? -1.0 : 1.0);
512 return ( norm *= radius );
531 if (v.z() <= 0.0) {
return distMin; }
538 return distMin= distZ;
543 if (v.z() >= 0.0) {
return distMin;}
549 return distMin= distZ;
566 G4double distR= (-B - std::sqrt(C)) / (2.0*A);
582 distR = (-B + std::sqrt(C) ) / (2.0*A);
583 if(distR>0.) { distMin=0.; }
587 distR= (-B + std::sqrt(C)) / (2.0*A);
588 intZ = p.z()+distR*v.z();
594 if (norm.dot(v)<0.) { distMin = distR; }
597 if ( (distMin!=
kInfinity) && (distMin>dRmax) )
600 G4double fTerm = distMin-std::fmod(distMin,dRmax);
627 distR= (p*norm - 1.0) * radius / 2.0;
641 return (distR < 0.0) ? 0.0 : distR;
643 else if (distR < 0.0)
649 return (distZ < distR) ? distZ : distR;
664 enum surface_e {kPlaneSurf, kCurvedSurf, kNoSurf} surface;
677 if (!calcNorm) {
return 0.0;}
688 if (!calcNorm) {
return 0.0;}
705 C= (p * nearnorm) - 1.0;
706 B= 2.0 * (v * nearnorm);
711 G4double distR= (-B + std::sqrt(C) ) / (2.0*A);
715 if (!calcNorm) {
return 0.0;}
720 surface= kCurvedSurf;
728 if (surface == kNoSurf)
746 truenorm *= 1.0/truenorm.mag();
751 std::ostringstream message;
752 G4int oldprc = message.precision(16);
753 message <<
"Undefined side for valid surface normal to solid."
756 <<
" p.x() = " << p.x()/
mm <<
" mm" <<
G4endl
757 <<
" p.y() = " << p.y()/
mm <<
" mm" <<
G4endl
758 <<
" p.z() = " << p.z()/
mm <<
" mm" <<
G4endl
760 <<
" v.x() = " << v.x() <<
G4endl
761 <<
" v.y() = " << v.y() <<
G4endl
762 <<
" v.z() = " << v.z() <<
G4endl
763 <<
"Proposed distance :" <<
G4endl
764 <<
" distMin = " << distMin/
mm <<
" mm";
765 message.precision(oldprc);
788 std::ostringstream message;
789 G4int oldprc = message.precision(16);
790 message <<
"Point p is outside !?" <<
G4endl
792 <<
" p.x() = " << p.x()/
mm <<
" mm" <<
G4endl
793 <<
" p.y() = " << p.y()/
mm <<
" mm" <<
G4endl
794 <<
" p.z() = " << p.z()/
mm <<
" mm";
795 message.precision(oldprc) ;
796 G4Exception(
"G4Ellipsoid::DistanceToOut(p)",
"GeomSolids1002",
811 if ( (tmp > 0.0) && (1.0 < radius*tmp) ) {radius = 1.0/tmp;}
815 distR = (1.0 - p*norm) * radius / 2.0;
820 if (distZ < 0.0) {distZ=
zTopCut - p.z();}
824 if ( (distZ < 0.0) || (distR < 0.0) )
830 return (distZ < distR) ? distZ : distR;
847 G4int& noPolygonVertices)
const
851 G4double meshAnglePhi, meshRMaxFactor,
852 crossAnglePhi, coscrossAnglePhi, sincrossAnglePhi, sAnglePhi;
853 G4double meshTheta, crossTheta, startTheta;
854 G4double rMaxX, rMaxY, rMaxZ, rMaxMax, rx, ry, rz;
855 G4int crossSectionPhi, noPhiCrossSections, crossSectionTheta, noThetaSections;
871 meshAnglePhi=
twopi/(noPhiCrossSections-1);
876 sAnglePhi = -meshAnglePhi*0.5;
892 meshTheta=
pi/(noThetaSections-2);
897 startTheta = -meshTheta*0.5;
899 meshRMaxFactor = 1.0/std::cos(0.5*
900 std::sqrt(meshAnglePhi*meshAnglePhi+meshTheta*meshTheta));
903 rMaxX=
xSemiAxis + rMaxMax*(meshRMaxFactor-1.0);
904 rMaxY=
ySemiAxis + rMaxMax*(meshRMaxFactor-1.0);
905 rMaxZ=
zSemiAxis + rMaxMax*(meshRMaxFactor-1.0);
909 if (vertices && cosCrossTheta && sinCrossTheta)
911 for (crossSectionTheta=0; crossSectionTheta<noThetaSections;
916 crossTheta=startTheta+crossSectionTheta*meshTheta;
917 cosCrossTheta[crossSectionTheta]=std::cos(crossTheta);
918 sinCrossTheta[crossSectionTheta]=std::sin(crossTheta);
920 for (crossSectionPhi=0; crossSectionPhi<noPhiCrossSections;
923 crossAnglePhi=sAnglePhi+crossSectionPhi*meshAnglePhi;
924 coscrossAnglePhi=std::cos(crossAnglePhi);
925 sincrossAnglePhi=std::sin(crossAnglePhi);
926 for (crossSectionTheta=0; crossSectionTheta<noThetaSections;
931 rx= sinCrossTheta[crossSectionTheta]*coscrossAnglePhi*rMaxX;
932 ry= sinCrossTheta[crossSectionTheta]*sincrossAnglePhi*rMaxY;
933 rz= cosCrossTheta[crossSectionTheta]*rMaxZ;
942 noPolygonVertices = noThetaSections ;
947 G4Exception(
"G4Ellipsoid::CreateRotatedVertices()",
949 "Error in allocation of vertices. Out of memory !");
952 delete[] cosCrossTheta;
953 delete[] sinCrossTheta;
982 G4int oldprc = os.precision(16);
983 os <<
"-----------------------------------------------------------\n"
984 <<
" *** Dump for solid - " <<
GetName() <<
" ***\n"
985 <<
" ===================================================\n"
986 <<
" Solid type: G4Ellipsoid\n"
993 <<
" lower cut plane level z: " <<
zBottomCut/
mm <<
" mm \n"
994 <<
" upper cut plane level z: " <<
zTopCut/
mm <<
" mm \n"
995 <<
"-----------------------------------------------------------\n";
996 os.precision(oldprc);
1007 G4double aTop, aBottom, aCurved, chose, xRand, yRand, zRand, phi;
1018 cosphi = std::cos(phi); sinphi = std::sin(phi);
1020 sintheta = std::sqrt(1.-
sqr(costheta));
1022 alpha = 1.-
sqr(max2/max1); beta = 1.-
sqr(max3/max1);
1029 aCurved = 4.*
pi*max1*max2*(1.-1./6.*(alpha+beta)-
1030 1./120.*(3.*
sqr(alpha)+2.*alpha*beta+3.*
sqr(beta)));
1037 aTop = 0; aBottom = 0;
1049 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 const G4double kInfinity
G4double GetMinYExtent() const
CLHEP::Hep3Vector G4ThreeVector
G4bool IsYLimited() const
G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=G4bool(false), G4bool *validNorm=0, G4ThreeVector *n=0) const
EInside Inside(const G4ThreeVector &p) const
G4bool IsXLimited() 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
#define G4MUTEX_INITIALIZER
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)
G4double GetMaxXExtent() const
G4double GetMinZExtent() const
double A(double temperature)
void CalculateClippedPolygonExtent(G4ThreeVectorList &pPolygon, const G4VoxelLimits &pVoxelLimit, const EAxis pAxis, G4double &pMin, G4double &pMax) const
G4VisExtent GetExtent() const
G4Polyhedron * fpPolyhedron
G4ThreeVectorList * CreateRotatedVertices(const G4AffineTransform &pT, G4int &noPV) const
G4double GetRadialTolerance() const
static const double twopi
G4Ellipsoid & operator=(const G4Ellipsoid &rhs)
G4Polyhedron * CreatePolyhedron() const
std::vector< G4ThreeVector > G4ThreeVectorList
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
G4double GetMinXExtent() const
std::ostream & StreamInfo(std::ostream &os) const
G4double GetMaxZExtent() const
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
void SetSemiAxis(G4double x, G4double y, G4double z)
G4double GetMaxYExtent() const
G4VSolid & operator=(const G4VSolid &rhs)
G4int GetNumberOfRotationStepsAtTimeOfCreation() const
G4ThreeVector GetPointOnSurface() const
G4double GetMaxExtent(const EAxis pAxis) const
static const G4double alpha
G4bool IsZLimited() const
G4double GetMinExtent(const EAxis pAxis) const
void ComputeDimensions(G4VPVParameterisation *p, const G4int n, const G4VPhysicalVolume *pRep)
const G4double kMeshAngleDefault
static G4GeometryTolerance * GetInstance()
static int max3(int a, int b, int c)