56 using namespace CLHEP;
69 :
G4VSolid(pName), fpPolyhedron(0), fCubicVolume(0.), fSurfaceArea(0.),
70 zBottomCut(0.), zTopCut(0.)
78 halfRadTolerance = kRadTolerance*0.5;
81 if ( (pxSemiAxis<=0.) || (pySemiAxis<=0.) || (pzSemiAxis<=0.) )
83 std::ostringstream message;
84 message <<
"Invalid semi-axis - " <<
GetName();
85 G4Exception(
"G4Ellipsoid::G4Ellipsoid()",
"GeomSolids0002",
90 if ( pzBottomCut == 0 && pzTopCut == 0 )
94 else if ( (pzBottomCut < pzSemiAxis) && (pzTopCut > -pzSemiAxis)
95 && (pzBottomCut < pzTopCut) )
101 std::ostringstream message;
102 message <<
"Invalid z-coordinate for cutting plane - " <<
GetName();
103 G4Exception(
"G4Ellipsoid::G4Ellipsoid()",
"GeomSolids0002",
114 :
G4VSolid(a), fpPolyhedron(0), kRadTolerance(0.),
115 halfCarTolerance(0.), halfRadTolerance(0.), fCubicVolume(0.),
116 fSurfaceArea(0.), xSemiAxis(0.), ySemiAxis(0.), zSemiAxis(0.),
117 semiAxisMax(0.), zBottomCut(0.), zTopCut(0.)
135 fpPolyhedron(0), kRadTolerance(rhs.kRadTolerance),
136 halfCarTolerance(rhs.halfCarTolerance),
137 halfRadTolerance(rhs.halfRadTolerance),
138 fCubicVolume(rhs.fCubicVolume), fSurfaceArea(rhs.fSurfaceArea),
139 xSemiAxis(rhs.xSemiAxis), ySemiAxis(rhs.ySemiAxis),
140 zSemiAxis(rhs.zSemiAxis), semiAxisMax(rhs.semiAxisMax),
141 zBottomCut(rhs.zBottomCut), zTopCut(rhs.zTopCut)
153 if (
this == &rhs) {
return *
this; }
162 halfCarTolerance = rhs.halfCarTolerance;
163 halfRadTolerance = rhs.halfRadTolerance;
164 fCubicVolume = rhs.fCubicVolume; fSurfaceArea = rhs.fSurfaceArea;
165 xSemiAxis = rhs.xSemiAxis; ySemiAxis = rhs.ySemiAxis;
166 zSemiAxis = rhs.zSemiAxis; semiAxisMax = rhs.semiAxisMax;
167 zBottomCut = rhs.zBottomCut; zTopCut = rhs.zTopCut;
209 xMin=xoffset - xSemiAxis;
210 xMax=xoffset + xSemiAxis;
232 yMin=yoffset - ySemiAxis;
233 yMax=yoffset + ySemiAxis;
255 zMin=zoffset + (-zSemiAxis > zBottomCut ? -zSemiAxis : zBottomCut);
256 zMax=zoffset + ( zSemiAxis < zTopCut ? zSemiAxis : zTopCut);
279 xoff = (xoffset < xMin) ? (xMin-xoffset)
280 : (xoffset > xMax) ? (xoffset-xMax) : 0.0;
281 yoff = (yoffset < yMin) ? (yMin-yoffset)
282 : (yoffset > yMax) ? (yoffset-yMax) : 0.0;
304 maxDiff= 1.0-
sqr(yoff/ySemiAxis);
305 if (maxDiff < 0.0) {
return false; }
306 maxDiff= xSemiAxis * std::sqrt(maxDiff);
307 newMin=xoffset-maxDiff;
308 newMax=xoffset+maxDiff;
309 pMin=(newMin<xMin) ? xMin : newMin;
310 pMax=(newMax>xMax) ? xMax : newMax;
326 maxDiff= 1.0-
sqr(xoff/xSemiAxis);
327 if (maxDiff < 0.0) {
return false; }
328 maxDiff= ySemiAxis * std::sqrt(maxDiff);
329 newMin=yoffset-maxDiff;
330 newMax=yoffset+maxDiff;
331 pMin=(newMin<yMin) ? yMin : newMin;
332 pMax=(newMax>yMax) ? yMax : newMax;
349 G4int i,j,noEntries,noBetweenSections;
350 G4bool existsAfterClip=
false;
354 G4int noPolygonVertices=0;
361 noEntries=vertices->size();
362 noBetweenSections=noEntries-noPolygonVertices;
365 for (i=0;i<noEntries;i+=noPolygonVertices)
367 for(j=0;j<(noPolygonVertices/2)-1;j++)
369 ThetaPolygon.push_back((*vertices)[i+j]);
370 ThetaPolygon.push_back((*vertices)[i+j+1]);
371 ThetaPolygon.push_back((*vertices)[i+noPolygonVertices-2-j]);
372 ThetaPolygon.push_back((*vertices)[i+noPolygonVertices-1-j]);
374 ThetaPolygon.clear();
377 for (i=0;i<noBetweenSections;i+=noPolygonVertices)
379 for(j=0;j<noPolygonVertices-1;j++)
381 ThetaPolygon.push_back((*vertices)[i+j]);
382 ThetaPolygon.push_back((*vertices)[i+j+1]);
383 ThetaPolygon.push_back((*vertices)[i+noPolygonVertices+j+1]);
384 ThetaPolygon.push_back((*vertices)[i+noPolygonVertices+j]);
386 ThetaPolygon.clear();
388 ThetaPolygon.push_back((*vertices)[i+noPolygonVertices-1]);
389 ThetaPolygon.push_back((*vertices)[i]);
390 ThetaPolygon.push_back((*vertices)[i+noPolygonVertices]);
391 ThetaPolygon.push_back((*vertices)[i+2*noPolygonVertices-1]);
393 ThetaPolygon.clear();
395 if ( (pMin!=kInfinity) || (pMax!=-kInfinity) )
397 existsAfterClip=
true;
419 existsAfterClip=
true;
425 return existsAfterClip;
443 if (p.
z() < zBottomCut-halfRadTolerance) {
return in=
kOutside; }
444 if (p.
z() > zTopCut+halfRadTolerance) {
return in=
kOutside; }
446 rad2oo=
sqr(p.
x()/(xSemiAxis+halfRadTolerance))
447 +
sqr(p.
y()/(ySemiAxis+halfRadTolerance))
448 +
sqr(p.
z()/(zSemiAxis+halfRadTolerance));
450 if (rad2oo > 1.0) {
return in=
kOutside; }
452 rad2oi=
sqr(p.
x()*(1.0+halfRadTolerance/xSemiAxis)/xSemiAxis)
453 +
sqr(p.
y()*(1.0+halfRadTolerance/ySemiAxis)/ySemiAxis)
454 +
sqr(p.
z()*(1.0+halfRadTolerance/zSemiAxis)/zSemiAxis);
461 in = ( (p.
z() < zBottomCut+halfRadTolerance)
463 if ( rad2oi > 1.0-halfRadTolerance ) { in=
kSurface; }
479 G4double distR, distZBottom, distZTop;
485 p.
y()/(ySemiAxis*ySemiAxis),
486 p.
z()/(zSemiAxis*zSemiAxis));
491 distR = std::fabs( (p*norm - 1.0) * radius ) / 2.0;
495 distZBottom = std::fabs( p.
z() - zBottomCut );
496 distZTop = std::fabs( p.
z() - zTopCut );
498 if ( (distZBottom < distR) || (distZTop < distR) )
500 return G4ThreeVector(0.,0.,(distZBottom < distZTop) ? -1.0 : 1.0);
502 return ( norm *= radius );
519 if (p.
z() <= zBottomCut+halfCarTolerance)
521 if (v.
z() <= 0.0) {
return distMin; }
522 G4double distZ = (zBottomCut - p.
z()) / v.
z();
524 if ( (distZ > -halfRadTolerance) && (
Inside(p+distZ*v) !=
kOutside) )
527 if ( std::fabs(distZ) < halfRadTolerance ) { distZ=0.; }
528 return distMin= distZ;
531 if (p.
z() >= zTopCut-halfCarTolerance)
533 if (v.
z() >= 0.0) {
return distMin;}
535 if ( (distZ > -halfRadTolerance) && (
Inside(p+distZ*v) !=
kOutside) )
538 if ( std::fabs(distZ) < halfRadTolerance ) { distZ=0.; }
539 return distMin= distZ;
547 A=
sqr(v.
x()/xSemiAxis) +
sqr(v.
y()/ySemiAxis) +
sqr(v.
z()/zSemiAxis);
548 C=
sqr(p.
x()/xSemiAxis) +
sqr(p.
y()/ySemiAxis) +
sqr(p.
z()/zSemiAxis) - 1.0;
549 B= 2.0 * ( p.
x()*v.
x()/(xSemiAxis*xSemiAxis)
550 + p.
y()*v.
y()/(ySemiAxis*ySemiAxis)
551 + p.
z()*v.
z()/(zSemiAxis*zSemiAxis) );
556 G4double distR= (-B - std::sqrt(C)) / (2.0*A);
558 if ( (distR > halfRadTolerance)
559 && (intZ >= zBottomCut-halfRadTolerance)
560 && (intZ <= zTopCut+halfRadTolerance) )
564 else if( (distR >- halfRadTolerance)
565 && (intZ >= zBottomCut-halfRadTolerance)
566 && (intZ <= zTopCut+halfRadTolerance) )
572 distR = (-B + std::sqrt(C) ) / (2.0*A);
573 if(distR>0.) { distMin=0.; }
577 distR= (-B + std::sqrt(C)) / (2.0*A);
578 intZ = p.
z()+distR*v.
z();
579 if ( (distR > halfRadTolerance)
580 && (intZ >= zBottomCut-halfRadTolerance)
581 && (intZ <= zTopCut+halfRadTolerance) )
584 if (norm.
dot(v)<0.) { distMin = distR; }
587 if ( (distMin!=kInfinity) && (distMin>dRmax) )
590 G4double fTerm = distMin-std::fmod(distMin,dRmax);
595 if (std::fabs(distMin)<halfRadTolerance) { distMin=0.; }
611 p.
y()/(ySemiAxis*ySemiAxis),
612 p.
z()/(zSemiAxis*zSemiAxis));
617 distR= (p*norm - 1.0) * radius / 2.0;
621 distZ= zBottomCut - p.
z();
624 distZ = p.
z() - zTopCut;
631 return (distR < 0.0) ? 0.0 : distR;
633 else if (distR < 0.0)
639 return (distZ < distR) ? distZ : distR;
654 enum surface_e {kPlaneSurf, kCurvedSurf, kNoSurf} surface;
663 G4double distZ = (zBottomCut - p.
z()) / v.
z();
667 if (!calcNorm) {
return 0.0;}
678 if (!calcNorm) {
return 0.0;}
687 p.
y()/(ySemiAxis*ySemiAxis),
688 p.
z()/(zSemiAxis*zSemiAxis));
694 A=
sqr(v.
x()/xSemiAxis) +
sqr(v.
y()/ySemiAxis) +
sqr(v.
z()/zSemiAxis);
695 C= (p * nearnorm) - 1.0;
696 B= 2.0 * (v * nearnorm);
701 G4double distR= (-B + std::sqrt(C) ) / (2.0*A);
705 if (!calcNorm) {
return 0.0;}
710 surface= kCurvedSurf;
718 if (surface == kNoSurf)
734 pexit.
y()/(ySemiAxis*ySemiAxis),
735 pexit.
z()/(zSemiAxis*zSemiAxis));
736 truenorm *= 1.0/truenorm.
mag();
741 std::ostringstream message;
742 G4int oldprc = message.precision(16);
743 message <<
"Undefined side for valid surface normal to solid."
746 <<
" p.x() = " << p.
x()/
mm <<
" mm" <<
G4endl
747 <<
" p.y() = " << p.
y()/
mm <<
" mm" <<
G4endl
748 <<
" p.z() = " << p.
z()/
mm <<
" mm" <<
G4endl
750 <<
" v.x() = " << v.
x() <<
G4endl
751 <<
" v.y() = " << v.
y() <<
G4endl
752 <<
" v.z() = " << v.
z() <<
G4endl
753 <<
"Proposed distance :" <<
G4endl
754 <<
" distMin = " << distMin/
mm <<
" mm";
755 message.precision(oldprc);
778 std::ostringstream message;
779 G4int oldprc = message.precision(16);
780 message <<
"Point p is outside !?" <<
G4endl
782 <<
" p.x() = " << p.
x()/
mm <<
" mm" <<
G4endl
783 <<
" p.y() = " << p.
y()/
mm <<
" mm" <<
G4endl
784 <<
" p.z() = " << p.
z()/
mm <<
" mm";
785 message.precision(oldprc) ;
786 G4Exception(
"G4Ellipsoid::DistanceToOut(p)",
"GeomSolids1002",
794 p.
y()/(ySemiAxis*ySemiAxis),
795 p.
z()/(zSemiAxis*zSemiAxis));
801 if ( (tmp > 0.0) && (1.0 < radius*tmp) ) {radius = 1.0/
tmp;}
805 distR = (1.0 - p*
norm) * radius / 2.0;
809 distZ = p.
z() - zBottomCut;
810 if (distZ < 0.0) {distZ= zTopCut - p.
z();}
814 if ( (distZ < 0.0) || (distR < 0.0) )
820 return (distZ < distR) ? distZ : distR;
837 G4int& noPolygonVertices)
const
841 G4double meshAnglePhi, meshRMaxFactor,
842 crossAnglePhi, coscrossAnglePhi, sincrossAnglePhi, sAnglePhi;
843 G4double meshTheta, crossTheta, startTheta;
844 G4double rMaxX, rMaxY, rMaxZ, rMaxMax, rx, ry, rz;
845 G4int crossSectionPhi, noPhiCrossSections, crossSectionTheta, noThetaSections;
861 meshAnglePhi=
twopi/(noPhiCrossSections-1);
866 sAnglePhi = -meshAnglePhi*0.5;
882 meshTheta=
pi/(noThetaSections-2);
887 startTheta = -meshTheta*0.5;
889 meshRMaxFactor = 1.0/std::cos(0.5*
890 std::sqrt(meshAnglePhi*meshAnglePhi+meshTheta*meshTheta));
891 rMaxMax= (xSemiAxis > ySemiAxis ? xSemiAxis : ySemiAxis);
892 if (zSemiAxis > rMaxMax) rMaxMax= zSemiAxis;
893 rMaxX= xSemiAxis + rMaxMax*(meshRMaxFactor-1.0);
894 rMaxY= ySemiAxis + rMaxMax*(meshRMaxFactor-1.0);
895 rMaxZ= zSemiAxis + rMaxMax*(meshRMaxFactor-1.0);
899 if (vertices && cosCrossTheta && sinCrossTheta)
901 for (crossSectionTheta=0; crossSectionTheta<noThetaSections;
906 crossTheta=startTheta+crossSectionTheta*meshTheta;
907 cosCrossTheta[crossSectionTheta]=std::cos(crossTheta);
908 sinCrossTheta[crossSectionTheta]=std::sin(crossTheta);
910 for (crossSectionPhi=0; crossSectionPhi<noPhiCrossSections;
913 crossAnglePhi=sAnglePhi+crossSectionPhi*meshAnglePhi;
914 coscrossAnglePhi=std::cos(crossAnglePhi);
915 sincrossAnglePhi=std::sin(crossAnglePhi);
916 for (crossSectionTheta=0; crossSectionTheta<noThetaSections;
921 rx= sinCrossTheta[crossSectionTheta]*coscrossAnglePhi*rMaxX;
922 ry= sinCrossTheta[crossSectionTheta]*sincrossAnglePhi*rMaxY;
923 rz= cosCrossTheta[crossSectionTheta]*rMaxZ;
932 noPolygonVertices = noThetaSections ;
937 G4Exception(
"G4Ellipsoid::CreateRotatedVertices()",
939 "Error in allocation of vertices. Out of memory !");
942 delete[] cosCrossTheta;
943 delete[] sinCrossTheta;
972 G4int oldprc = os.precision(16);
973 os <<
"-----------------------------------------------------------\n"
974 <<
" *** Dump for solid - " <<
GetName() <<
" ***\n"
975 <<
" ===================================================\n"
976 <<
" Solid type: G4Ellipsoid\n"
979 <<
" semi-axis x: " << xSemiAxis/
mm <<
" mm \n"
980 <<
" semi-axis y: " << ySemiAxis/
mm <<
" mm \n"
981 <<
" semi-axis z: " << zSemiAxis/
mm <<
" mm \n"
982 <<
" max semi-axis: " << semiAxisMax/
mm <<
" mm \n"
983 <<
" lower cut plane level z: " << zBottomCut/
mm <<
" mm \n"
984 <<
" upper cut plane level z: " << zTopCut/
mm <<
" mm \n"
985 <<
"-----------------------------------------------------------\n";
986 os.precision(oldprc);
997 G4double aTop, aBottom, aCurved, chose, xRand, yRand, zRand, phi;
998 G4double cosphi, sinphi, costheta, sintheta, alpha, beta, max1, max2, max3;
1000 max1 = xSemiAxis > ySemiAxis ? xSemiAxis : ySemiAxis;
1001 max1 = max1 > zSemiAxis ? max1 : zSemiAxis;
1002 if (max1 == xSemiAxis) { max2 = ySemiAxis; max3 = zSemiAxis; }
1003 else if (max1 == ySemiAxis) { max2 = xSemiAxis; max3 = zSemiAxis; }
1004 else { max2 = xSemiAxis; max3 = ySemiAxis; }
1008 cosphi = std::cos(phi); sinphi = std::sin(phi);
1010 sintheta = std::sqrt(1.-
sqr(costheta));
1012 alpha = 1.-
sqr(max2/max1); beta = 1.-
sqr(max3/max1);
1014 aTop =
pi*xSemiAxis*ySemiAxis*(1 -
sqr(zTopCut/zSemiAxis));
1015 aBottom =
pi*xSemiAxis*ySemiAxis*(1 -
sqr(zBottomCut/zSemiAxis));
1019 aCurved = 4.*
pi*max1*max2*(1.-1./6.*(alpha+beta)-
1020 1./120.*(3.*
sqr(alpha)+2.*alpha*beta+3.*
sqr(beta)));
1022 aCurved *= 0.5*(1.2*zTopCut/zSemiAxis - 1.2*zBottomCut/zSemiAxis);
1024 if( ( zTopCut >= zSemiAxis && zBottomCut <= -1.*zSemiAxis )
1025 || ( zTopCut == 0 && zBottomCut ==0 ) )
1027 aTop = 0; aBottom = 0;
1034 xRand = xSemiAxis*sintheta*cosphi;
1035 yRand = ySemiAxis*sintheta*sinphi;
1036 zRand = zSemiAxis*costheta;
1039 else if(chose >= aCurved && chose < aCurved + aTop)
1042 * std::sqrt(1-
sqr(zTopCut/zSemiAxis));
1044 * std::sqrt(1.-
sqr(zTopCut/zSemiAxis)-
sqr(xRand/xSemiAxis));
1051 * std::sqrt(1-
sqr(zBottomCut/zSemiAxis));
1053 * std::sqrt(1.-
sqr(zBottomCut/zSemiAxis)-
sqr(xRand/xSemiAxis));
1073 -semiAxisMax, semiAxisMax,
1074 -semiAxisMax, semiAxisMax);
1080 zBottomCut, zTopCut);
G4Polyhedron * GetPolyhedron() const
ThreeVector shoot(const G4int Ap, const G4int Af)
G4GeometryType GetEntityType() const
G4double GetMinYExtent() const
CLHEP::Hep3Vector G4ThreeVector
double dot(const Hep3Vector &) const
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
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
static G4int GetNumberOfRotationSteps()
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
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()