55 using namespace CLHEP;
68 :
G4VSolid(pName), fpPolyhedron(0), fCubicVolume(0.), fSurfaceArea(0.),
69 zBottomCut(0.), zTopCut(0.)
80 if ( (pxSemiAxis<=0.) || (pySemiAxis<=0.) || (pzSemiAxis<=0.) )
82 std::ostringstream message;
83 message <<
"Invalid semi-axis - " <<
GetName();
84 G4Exception(
"G4Ellipsoid::G4Ellipsoid()",
"GeomSolids0002",
89 if ( pzBottomCut == 0 && pzTopCut == 0 )
93 else if ( (pzBottomCut < pzSemiAxis) && (pzTopCut > -pzSemiAxis)
94 && (pzBottomCut < pzTopCut) )
100 std::ostringstream message;
101 message <<
"Invalid z-coordinate for cutting plane - " <<
GetName();
102 G4Exception(
"G4Ellipsoid::G4Ellipsoid()",
"GeomSolids0002",
113 :
G4VSolid(a), fpPolyhedron(0), kRadTolerance(0.),
114 halfCarTolerance(0.), halfRadTolerance(0.), fCubicVolume(0.),
115 fSurfaceArea(0.), xSemiAxis(0.), ySemiAxis(0.), zSemiAxis(0.),
116 semiAxisMax(0.), zBottomCut(0.), zTopCut(0.)
134 fpPolyhedron(0), kRadTolerance(rhs.kRadTolerance),
135 halfCarTolerance(rhs.halfCarTolerance),
136 halfRadTolerance(rhs.halfRadTolerance),
137 fCubicVolume(rhs.fCubicVolume), fSurfaceArea(rhs.fSurfaceArea),
138 xSemiAxis(rhs.xSemiAxis), ySemiAxis(rhs.ySemiAxis),
139 zSemiAxis(rhs.zSemiAxis), semiAxisMax(rhs.semiAxisMax),
140 zBottomCut(rhs.zBottomCut), zTopCut(rhs.zTopCut)
153 if (
this == &rhs) {
return *
this; }
280 xoff = (xoffset < xMin) ? (xMin-xoffset)
281 : (xoffset > xMax) ? (xoffset-xMax) : 0.0;
282 yoff = (yoffset < yMin) ? (yMin-yoffset)
283 : (yoffset > yMax) ? (yoffset-yMax) : 0.0;
305 maxDiff= 1.0-
sqr(yoff/ySemiAxis);
306 if (maxDiff < 0.0) {
return false; }
307 maxDiff= xSemiAxis * std::sqrt(maxDiff);
308 newMin=xoffset-maxDiff;
309 newMax=xoffset+maxDiff;
310 pMin=(newMin<xMin) ? xMin : newMin;
311 pMax=(newMax>xMax) ? xMax : newMax;
327 maxDiff= 1.0-
sqr(xoff/xSemiAxis);
328 if (maxDiff < 0.0) {
return false; }
329 maxDiff= ySemiAxis * std::sqrt(maxDiff);
330 newMin=yoffset-maxDiff;
331 newMax=yoffset+maxDiff;
332 pMin=(newMin<yMin) ? yMin : newMin;
333 pMax=(newMax>yMax) ? yMax : newMax;
350 G4int i,j,noEntries,noBetweenSections;
351 G4bool existsAfterClip=
false;
355 G4int noPolygonVertices=0;
362 noEntries=vertices->size();
363 noBetweenSections=noEntries-noPolygonVertices;
366 for (i=0;i<noEntries;i+=noPolygonVertices)
368 for(j=0;j<(noPolygonVertices/2)-1;j++)
370 ThetaPolygon.push_back((*vertices)[i+j]);
371 ThetaPolygon.push_back((*vertices)[i+j+1]);
372 ThetaPolygon.push_back((*vertices)[i+noPolygonVertices-2-j]);
373 ThetaPolygon.push_back((*vertices)[i+noPolygonVertices-1-j]);
375 ThetaPolygon.clear();
378 for (i=0;i<noBetweenSections;i+=noPolygonVertices)
380 for(j=0;j<noPolygonVertices-1;j++)
382 ThetaPolygon.push_back((*vertices)[i+j]);
383 ThetaPolygon.push_back((*vertices)[i+j+1]);
384 ThetaPolygon.push_back((*vertices)[i+noPolygonVertices+j+1]);
385 ThetaPolygon.push_back((*vertices)[i+noPolygonVertices+j]);
387 ThetaPolygon.clear();
389 ThetaPolygon.push_back((*vertices)[i+noPolygonVertices-1]);
390 ThetaPolygon.push_back((*vertices)[i]);
391 ThetaPolygon.push_back((*vertices)[i+noPolygonVertices]);
392 ThetaPolygon.push_back((*vertices)[i+2*noPolygonVertices-1]);
394 ThetaPolygon.clear();
398 existsAfterClip=
true;
420 existsAfterClip=
true;
426 return existsAfterClip;
451 if (rad2oo > 1.0) {
return in=
kOutside; }
480 G4double distR, distZBottom, distZTop;
492 distR = std::fabs( (p*norm - 1.0) * radius ) / 2.0;
496 distZBottom = std::fabs( p.z() -
zBottomCut );
497 distZTop = std::fabs( p.z() -
zTopCut );
499 if ( (distZBottom < distR) || (distZTop < distR) )
501 return G4ThreeVector(0.,0.,(distZBottom < distZTop) ? -1.0 : 1.0);
503 return ( norm *= radius );
522 if (v.z() <= 0.0) {
return distMin; }
529 return distMin= distZ;
534 if (v.z() >= 0.0) {
return distMin;}
540 return distMin= distZ;
557 G4double distR= (-B - std::sqrt(C)) / (2.0*A);
573 distR = (-B + std::sqrt(C) ) / (2.0*A);
574 if(distR>0.) { distMin=0.; }
578 distR= (-B + std::sqrt(C)) / (2.0*A);
579 intZ = p.z()+distR*v.z();
585 if (norm.dot(v)<0.) { distMin = distR; }
588 if ( (distMin!=
kInfinity) && (distMin>dRmax) )
591 G4double fTerm = distMin-std::fmod(distMin,dRmax);
618 distR= (p*norm - 1.0) * radius / 2.0;
632 return (distR < 0.0) ? 0.0 : distR;
634 else if (distR < 0.0)
640 return (distZ < distR) ? distZ : distR;
655 enum surface_e {kPlaneSurf, kCurvedSurf, kNoSurf} surface;
668 if (!calcNorm) {
return 0.0;}
679 if (!calcNorm) {
return 0.0;}
696 C= (p * nearnorm) - 1.0;
697 B= 2.0 * (v * nearnorm);
702 G4double distR= (-B + std::sqrt(C) ) / (2.0*A);
706 if (!calcNorm) {
return 0.0;}
711 surface= kCurvedSurf;
719 if (surface == kNoSurf)
737 truenorm *= 1.0/truenorm.mag();
742 std::ostringstream message;
743 G4int oldprc = message.precision(16);
744 message <<
"Undefined side for valid surface normal to solid."
747 <<
" p.x() = " << p.x()/
mm <<
" mm" <<
G4endl
748 <<
" p.y() = " << p.y()/
mm <<
" mm" <<
G4endl
749 <<
" p.z() = " << p.z()/
mm <<
" mm" <<
G4endl
751 <<
" v.x() = " << v.x() <<
G4endl
752 <<
" v.y() = " << v.y() <<
G4endl
753 <<
" v.z() = " << v.z() <<
G4endl
754 <<
"Proposed distance :" <<
G4endl
755 <<
" distMin = " << distMin/
mm <<
" mm";
756 message.precision(oldprc);
779 std::ostringstream message;
780 G4int oldprc = message.precision(16);
781 message <<
"Point p is outside !?" <<
G4endl
783 <<
" p.x() = " << p.x()/
mm <<
" mm" <<
G4endl
784 <<
" p.y() = " << p.y()/
mm <<
" mm" <<
G4endl
785 <<
" p.z() = " << p.z()/
mm <<
" mm";
786 message.precision(oldprc) ;
787 G4Exception(
"G4Ellipsoid::DistanceToOut(p)",
"GeomSolids1002",
802 if ( (tmp > 0.0) && (1.0 < radius*tmp) ) {radius = 1.0/tmp;}
806 distR = (1.0 - p*norm) * radius / 2.0;
811 if (distZ < 0.0) {distZ=
zTopCut - p.z();}
815 if ( (distZ < 0.0) || (distR < 0.0) )
821 return (distZ < distR) ? distZ : distR;
838 G4int& noPolygonVertices)
const
842 G4double meshAnglePhi, meshRMaxFactor,
843 crossAnglePhi, coscrossAnglePhi, sincrossAnglePhi, sAnglePhi;
844 G4double meshTheta, crossTheta, startTheta;
845 G4double rMaxX, rMaxY, rMaxZ, rMaxMax, rx, ry, rz;
846 G4int crossSectionPhi, noPhiCrossSections, crossSectionTheta, noThetaSections;
862 meshAnglePhi=twopi/(noPhiCrossSections-1);
867 sAnglePhi = -meshAnglePhi*0.5;
883 meshTheta=
pi/(noThetaSections-2);
888 startTheta = -meshTheta*0.5;
890 meshRMaxFactor = 1.0/std::cos(0.5*
891 std::sqrt(meshAnglePhi*meshAnglePhi+meshTheta*meshTheta));
894 rMaxX=
xSemiAxis + rMaxMax*(meshRMaxFactor-1.0);
895 rMaxY=
ySemiAxis + rMaxMax*(meshRMaxFactor-1.0);
896 rMaxZ=
zSemiAxis + rMaxMax*(meshRMaxFactor-1.0);
900 if (vertices && cosCrossTheta && sinCrossTheta)
902 for (crossSectionTheta=0; crossSectionTheta<noThetaSections;
907 crossTheta=startTheta+crossSectionTheta*meshTheta;
908 cosCrossTheta[crossSectionTheta]=std::cos(crossTheta);
909 sinCrossTheta[crossSectionTheta]=std::sin(crossTheta);
911 for (crossSectionPhi=0; crossSectionPhi<noPhiCrossSections;
914 crossAnglePhi=sAnglePhi+crossSectionPhi*meshAnglePhi;
915 coscrossAnglePhi=std::cos(crossAnglePhi);
916 sincrossAnglePhi=std::sin(crossAnglePhi);
917 for (crossSectionTheta=0; crossSectionTheta<noThetaSections;
922 rx= sinCrossTheta[crossSectionTheta]*coscrossAnglePhi*rMaxX;
923 ry= sinCrossTheta[crossSectionTheta]*sincrossAnglePhi*rMaxY;
924 rz= cosCrossTheta[crossSectionTheta]*rMaxZ;
933 noPolygonVertices = noThetaSections ;
938 G4Exception(
"G4Ellipsoid::CreateRotatedVertices()",
940 "Error in allocation of vertices. Out of memory !");
943 delete[] cosCrossTheta;
944 delete[] sinCrossTheta;
973 G4int oldprc = os.precision(16);
974 os <<
"-----------------------------------------------------------\n"
975 <<
" *** Dump for solid - " <<
GetName() <<
" ***\n"
976 <<
" ===================================================\n"
977 <<
" Solid type: G4Ellipsoid\n"
984 <<
" lower cut plane level z: " <<
zBottomCut/
mm <<
" mm \n"
985 <<
" upper cut plane level z: " <<
zTopCut/
mm <<
" mm \n"
986 <<
"-----------------------------------------------------------\n";
987 os.precision(oldprc);
998 G4double aTop, aBottom, aCurved, chose, xRand, yRand, zRand, phi;
999 G4double cosphi, sinphi, costheta, sintheta,
alpha, beta, max1, max2, max3;
1009 cosphi = std::cos(phi); sinphi = std::sin(phi);
1011 sintheta = std::sqrt(1.-
sqr(costheta));
1013 alpha = 1.-
sqr(max2/max1); beta = 1.-
sqr(max3/max1);
1020 aCurved = 4.*
pi*max1*max2*(1.-1./6.*(alpha+beta)-
1021 1./120.*(3.*
sqr(alpha)+2.*alpha*beta+3.*
sqr(beta)));
1028 aTop = 0; aBottom = 0;
1040 else if(chose >= aCurved && chose < aCurved + aTop)
G4Polyhedron * GetPolyhedron() const
ThreeVector shoot(const G4int Ap, const G4int Af)
G4double halfRadTolerance
G4GeometryType GetEntityType() const
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
virtual void AddSolid(const G4Box &)=0
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pmin, G4double &pmax) const
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
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
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
static const G4double A[nN]
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()