59 #if !defined(G4GEOM_USE_USPHERE)
74 using namespace CLHEP;
93 :
G4CSGSolid(pName), fEpsilon(2.e-11), fSPhi(0.0),
94 fFullPhiSphere(true), fFullThetaSphere(true)
104 if ( (pRmin >= pRmax) || (pRmax < 1.1*
kRadTolerance) || (pRmin < 0) )
106 std::ostringstream message;
107 message <<
"Invalid radii for Solid: " <<
GetName() <<
G4endl
108 <<
" pRmin = " << pRmin <<
", pRmax = " << pRmax;
109 G4Exception(
"G4Sphere::G4Sphere()",
"GeomSolids0002",
128 :
G4CSGSolid(a), fRminTolerance(0.), fRmaxTolerance(0.),
129 kAngTolerance(0.), kRadTolerance(0.), fEpsilon(0.),
130 fRmin(0.), fRmax(0.), fSPhi(0.), fDPhi(0.), fSTheta(0.),
131 fDTheta(0.), sinCPhi(0.), cosCPhi(0.), cosHDPhiOT(0.), cosHDPhiIT(0.),
132 sinSPhi(0.), cosSPhi(0.), sinEPhi(0.), cosEPhi(0.), hDPhi(0.), cPhi(0.),
133 ePhi(0.), sinSTheta(0.), cosSTheta(0.), sinETheta(0.), cosETheta(0.),
134 tanSTheta(0.), tanSTheta2(0.), tanETheta(0.), tanETheta2(0.), eTheta(0.),
135 fFullPhiSphere(false), fFullThetaSphere(false), fFullSphere(true),
136 halfCarTolerance(0.), halfAngTolerance(0.)
153 :
G4CSGSolid(rhs), fRminTolerance(rhs.fRminTolerance),
154 fRmaxTolerance(rhs.fRmaxTolerance), kAngTolerance(rhs.kAngTolerance),
155 kRadTolerance(rhs.kRadTolerance), fEpsilon(rhs.fEpsilon),
156 fRmin(rhs.fRmin), fRmax(rhs.fRmax), fSPhi(rhs.fSPhi), fDPhi(rhs.fDPhi),
157 fSTheta(rhs.fSTheta), fDTheta(rhs.fDTheta),
158 sinCPhi(rhs.sinCPhi), cosCPhi(rhs.cosCPhi),
159 cosHDPhiOT(rhs.cosHDPhiOT), cosHDPhiIT(rhs.cosHDPhiIT),
160 sinSPhi(rhs.sinSPhi), cosSPhi(rhs.cosSPhi),
161 sinEPhi(rhs.sinEPhi), cosEPhi(rhs.cosEPhi),
162 hDPhi(rhs.hDPhi), cPhi(rhs.cPhi), ePhi(rhs.ePhi),
163 sinSTheta(rhs.sinSTheta), cosSTheta(rhs.cosSTheta),
164 sinETheta(rhs.sinETheta), cosETheta(rhs.cosETheta),
165 tanSTheta(rhs.tanSTheta), tanSTheta2(rhs.tanSTheta2),
166 tanETheta(rhs.tanETheta), tanETheta2(rhs.tanETheta2), eTheta(rhs.eTheta),
167 fFullPhiSphere(rhs.fFullPhiSphere), fFullThetaSphere(rhs.fFullThetaSphere),
168 fFullSphere(rhs.fFullSphere),
169 halfCarTolerance(rhs.halfCarTolerance),
170 halfAngTolerance(rhs.halfAngTolerance)
182 if (
this == &rhs) {
return *
this; }
244 G4double diff1,diff2,delta,maxDiff,newMin,newMax;
323 if ((yoff1>=0) && (yoff2>=0))
335 delta=fRmax*fRmax-yoff1*yoff1;
336 diff1=(delta>0.) ? std::sqrt(delta) : 0.;
337 delta=fRmax*fRmax-yoff2*yoff2;
338 diff2=(delta>0.) ? std::sqrt(delta) : 0.;
339 maxDiff=(diff1>diff2) ? diff1:diff2;
340 newMin=xoffset-maxDiff;
341 newMax=xoffset+maxDiff;
342 pMin=(newMin<xMin) ? xMin : newMin;
343 pMax=(newMax>xMax) ? xMax : newMax;
349 if ((xoff1>=0) && (xoff2>=0))
361 delta=fRmax*fRmax-xoff1*xoff1;
362 diff1=(delta>0.) ? std::sqrt(delta) : 0.;
363 delta=fRmax*fRmax-xoff2*xoff2;
364 diff2=(delta>0.) ? std::sqrt(delta) : 0.;
365 maxDiff=(diff1>diff2) ? diff1:diff2;
366 newMin=yoffset-maxDiff;
367 newMax=yoffset+maxDiff;
368 pMin=(newMin<yMin) ? yMin : newMin;
369 pMax=(newMax>yMax) ? yMax : newMax;
386 G4int i,j,noEntries,noBetweenSections;
387 G4bool existsAfterClip=
false;
392 G4int noPolygonVertices ;
398 noEntries=vertices->size();
399 noBetweenSections=noEntries-noPolygonVertices;
402 for (i=0;i<noEntries;i+=noPolygonVertices)
404 for(j=0;j<(noPolygonVertices/2)-1;j++)
406 ThetaPolygon.push_back((*vertices)[i+j]) ;
407 ThetaPolygon.push_back((*vertices)[i+j+1]) ;
408 ThetaPolygon.push_back((*vertices)[i+noPolygonVertices-2-j]) ;
409 ThetaPolygon.push_back((*vertices)[i+noPolygonVertices-1-j]) ;
411 ThetaPolygon.clear() ;
414 for (i=0;i<noBetweenSections;i+=noPolygonVertices)
416 for(j=0;j<noPolygonVertices-1;j++)
418 ThetaPolygon.push_back((*vertices)[i+j]) ;
419 ThetaPolygon.push_back((*vertices)[i+j+1]) ;
420 ThetaPolygon.push_back((*vertices)[i+noPolygonVertices+j+1]) ;
421 ThetaPolygon.push_back((*vertices)[i+noPolygonVertices+j]) ;
423 ThetaPolygon.clear() ;
425 ThetaPolygon.push_back((*vertices)[i+noPolygonVertices-1]) ;
426 ThetaPolygon.push_back((*vertices)[i]) ;
427 ThetaPolygon.push_back((*vertices)[i+noPolygonVertices]) ;
428 ThetaPolygon.push_back((*vertices)[i+2*noPolygonVertices-1]) ;
430 ThetaPolygon.clear() ;
435 existsAfterClip=
true;
456 existsAfterClip=
true;
462 return existsAfterClip;
474 G4double rho,rho2,rad2,tolRMin,tolRMax;
483 rho2 = p.x()*p.x() + p.y()*p.y() ;
484 rad2 = rho2 + p.z()*p.z() ;
489 tolRMax = Rmax_minus;
507 if ( (rad2 <= Rmax_minus*Rmax_minus) && (rad2 >= Rmin_plus*Rmin_plus) )
513 tolRMax =
fRmax + halfRmaxTolerance;
515 if ( (rad2 <= tolRMax*tolRMax) && (rad2 >= tolRMin*tolRMin) )
529 pPhi = std::atan2(p.y(),p.x()) ;
548 rho = std::sqrt(rho2);
549 pTheta = std::atan2(rho,p.z());
588 G4int noSurfaces = 0;
589 G4double rho, rho2, radius, pTheta, pPhi=0.;
596 rho2 = p.x()*p.x()+p.y()*p.y();
597 radius = std::sqrt(rho2+p.z()*p.z());
598 rho = std::sqrt(rho2);
601 if (
fRmin) distRMin = std::fabs(radius-
fRmin);
605 pPhi = std::atan2(p.y(),p.x());
614 distSPhi = std::fabs( pPhi-
fSPhi );
615 distEPhi = std::fabs( pPhi-
ePhi );
629 pTheta = std::atan2(rho,p.z());
630 distSTheta = std::fabs(pTheta-
fSTheta);
631 distETheta = std::fabs(pTheta-
eTheta);
655 if( radius ) { nR =
G4ThreeVector(p.x()/radius,p.y()/radius,p.z()/radius); }
686 else { sumnorm += nTs; }
692 else { sumnorm += nTe; }
693 if(sumnorm.z() == 0.) { sumnorm += nZ; }
696 if ( noSurfaces == 0 )
699 G4Exception(
"G4Sphere::SurfaceNormal(p)",
"GeomSolids1002",
704 else if ( noSurfaces == 1 ) { norm = sumnorm; }
705 else { norm = sumnorm.unit(); }
719 G4double rho,rho2,radius,pPhi,pTheta;
720 G4double distRMin,distRMax,distSPhi,distEPhi,
721 distSTheta,distETheta,distMin;
723 rho2=p.x()*p.x()+p.y()*p.y();
724 radius=std::sqrt(rho2+p.z()*p.z());
731 distRMax=std::fabs(radius-
fRmax);
734 distRMin=std::fabs(radius-
fRmin);
736 if (distRMin<distRMax)
758 pPhi = std::atan2(p.y(),p.x());
759 if (pPhi<0) { pPhi += twopi; }
765 distSPhi=std::fabs(pPhi-(
fSPhi+twopi))*rho;
769 distSPhi=std::fabs(pPhi-
fSPhi)*rho;
776 if (distSPhi<distEPhi)
778 if (distSPhi<distMin)
786 if (distEPhi<distMin)
800 pTheta=std::atan2(rho,p.z());
801 distSTheta=std::fabs(pTheta-
fSTheta)*radius;
806 if (distSTheta<distETheta)
808 if (distSTheta<distMin)
810 distMin = distSTheta ;
816 if (distETheta<distMin)
818 distMin = distETheta ;
827 norm=
G4ThreeVector(-p.x()/radius,-p.y()/radius,-p.z()/radius);
852 "Undefined side for valid surface normal to solid.");
892 G4double rho2, rad2, pDotV2d, pDotV3d, pTheta ;
893 G4double tolSTheta=0., tolETheta=0. ;
899 ? (
fRmin-halfRminTolerance)*(
fRmin-halfRminTolerance) : 0;
901 (
fRmin+halfRminTolerance)*(
fRmin+halfRminTolerance);
903 (
fRmax+halfRmaxTolerance)*(
fRmax+halfRmaxTolerance);
905 (
fRmax-halfRmaxTolerance)*(
fRmax-halfRmaxTolerance);
909 G4double xi, yi, zi, rhoi, rhoi2, radi2, iTheta ;
926 rho2 = p.x()*p.x() + p.y()*p.y() ;
927 rad2 = rho2 + p.z()*p.z() ;
928 pTheta = std::atan2(std::sqrt(rho2),p.z()) ;
930 pDotV2d = p.x()*v.x() + p.y()*v.y() ;
931 pDotV3d = pDotV2d + p.z()*v.z() ;
942 if ((rad2!=0.0) || (
fRmin!=0.0))
948 G4double vTheta = std::atan2(std::sqrt(v.x()*v.x()+v.y()*v.y()),v.z()) ;
949 if ( (vTheta < tolSTheta) || (vTheta > tolETheta) )
978 d2 = pDotV3d*pDotV3d - c ;
982 sd = -pDotV3d - std::sqrt(d2) ;
988 G4double fTerm = sd-std::fmod(sd,dRmax);
991 xi = p.x() + sd*v.x() ;
992 yi = p.y() + sd*v.y() ;
993 rhoi = std::sqrt(xi*xi + yi*yi) ;
1003 zi = p.z() + sd*v.z() ;
1008 iTheta = std::atan2(rhoi,zi) ;
1009 if ( (iTheta >= tolSTheta) && (iTheta <= tolETheta) )
1024 zi = p.z() + sd*v.z() ;
1029 iTheta = std::atan2(rhoi,zi) ;
1030 if ( (iTheta >= tolSTheta) && (iTheta <= tolETheta) )
1052 d2 = pDotV3d*pDotV3d - c ;
1054 if ( (rad2 > tolIRMax2)
1108 d2 = pDotV3d*pDotV3d - c ;
1113 if ( (c > -halfRminTolerance) && (rad2 < tolIRMin2)
1160 sd = -pDotV3d + std::sqrt(d2) ;
1161 if ( sd >= halfRminTolerance )
1163 xi = p.x() + sd*v.x() ;
1164 yi = p.y() + sd*v.y() ;
1165 rhoi = std::sqrt(xi*xi+yi*yi) ;
1175 zi = p.z() + sd*v.z() ;
1180 iTheta = std::atan2(rhoi,zi) ;
1181 if ( (iTheta >= tolSTheta) && (iTheta<=tolETheta) )
1196 zi = p.z() + sd*v.z() ;
1201 iTheta = std::atan2(rhoi,zi) ;
1202 if ( (iTheta >= tolSTheta) && (iTheta <= tolETheta) )
1245 xi = p.x() + sd*v.x() ;
1246 yi = p.y() + sd*v.y() ;
1247 zi = p.z() + sd*v.z() ;
1248 rhoi2 = xi*xi + yi*yi ;
1249 radi2 = rhoi2 + zi*zi ;
1260 if ( (radi2 <= tolORMax2)
1261 && (radi2 >= tolORMin2)
1270 iTheta = std::atan2(std::sqrt(rhoi2),zi) ;
1271 if ( (iTheta >= tolSTheta) && (iTheta <= tolETheta) )
1276 if ((yi*
cosCPhi-xi*sinCPhi) <= 0)
1307 xi = p.x() + sd*v.x() ;
1308 yi = p.y() + sd*v.y() ;
1309 zi = p.z() + sd*v.z() ;
1310 rhoi2 = xi*xi + yi*yi ;
1311 radi2 = rhoi2 + zi*zi ;
1322 if ( (radi2 <= tolORMax2)
1323 && (radi2 >= tolORMin2)
1332 iTheta = std::atan2(std::sqrt(rhoi2),zi) ;
1333 if ( (iTheta >= tolSTheta) && (iTheta <= tolETheta) )
1338 if ((yi*
cosCPhi-xi*sinCPhi) >= 0)
1380 dist2STheta = rho2 - p.z()*p.z()*
tanSTheta2 ;
1394 if ( pTheta < tolSTheta )
1404 c = dist2STheta/t1 ;
1411 zi = p.z() + sd*v.z();
1413 if ( (sd < 0) || (zi*(
fSTheta - halfpi) > 0) )
1417 if ((sd >= 0) && (sd < snxt))
1419 xi = p.x() + sd*v.x();
1420 yi = p.y() + sd*v.y();
1421 zi = p.z() + sd*v.z();
1422 rhoi2 = xi*xi + yi*yi;
1423 radi2 = rhoi2 + zi*zi;
1424 if ( (radi2 <= tolORMax2)
1425 && (radi2 >= tolORMin2)
1426 && (zi*(
fSTheta - halfpi) <= 0) )
1455 c = dist2ETheta/t1 ;
1463 if ( (sd >= 0) && (sd < snxt) )
1465 xi = p.x() + sd*v.x() ;
1466 yi = p.y() + sd*v.y() ;
1467 zi = p.z() + sd*v.z() ;
1468 rhoi2 = xi*xi + yi*yi ;
1469 radi2 = rhoi2 + zi*zi ;
1471 if ( (radi2 <= tolORMax2)
1472 && (radi2 >= tolORMin2)
1473 && (zi*(
eTheta - halfpi) <= 0) )
1493 else if ( pTheta > tolETheta )
1504 c = dist2ETheta/t1 ;
1511 zi = p.z() + sd*v.z();
1513 if ( (sd < 0) || (zi*(
eTheta - halfpi) > 0) )
1517 if ( (sd >= 0) && (sd < snxt) )
1519 xi = p.x() + sd*v.x() ;
1520 yi = p.y() + sd*v.y() ;
1521 zi = p.z() + sd*v.z() ;
1522 rhoi2 = xi*xi + yi*yi ;
1523 radi2 = rhoi2 + zi*zi ;
1525 if ( (radi2 <= tolORMax2)
1526 && (radi2 >= tolORMin2)
1527 && (zi*(
eTheta - halfpi) <= 0) )
1556 c = dist2STheta/t1 ;
1564 if ( (sd >= 0) && (sd < snxt) )
1566 xi = p.x() + sd*v.x() ;
1567 yi = p.y() + sd*v.y() ;
1568 zi = p.z() + sd*v.z() ;
1569 rhoi2 = xi*xi + yi*yi ;
1570 radi2 = rhoi2 + zi*zi ;
1572 if ( (radi2 <= tolORMax2)
1573 && (radi2 >= tolORMin2)
1574 && (zi*(
fSTheta - halfpi) <= 0) )
1602 if ( (t2>=0 && tolIRMin2<rad2 && rad2<tolIRMax2 &&
fSTheta<halfpi)
1603 || (t2<0 && tolIRMin2<rad2 && rad2<tolIRMax2 && fSTheta>halfpi)
1604 || (v.z()<0 && tolIRMin2<rad2 && rad2<tolIRMax2 &&
fSTheta==halfpi) )
1626 c = dist2STheta/t1 ;
1635 xi = p.x() + sd*v.x() ;
1636 yi = p.y() + sd*v.y() ;
1637 zi = p.z() + sd*v.z() ;
1638 rhoi2 = xi*xi + yi*yi ;
1639 radi2 = rhoi2 + zi*zi ;
1641 if ( (radi2 <= tolORMax2)
1642 && (radi2 >= tolORMin2)
1643 && (zi*(fSTheta - halfpi) <= 0) )
1671 if ( ((t2<0) && (
eTheta < halfpi)
1672 && (tolIRMin2 < rad2) && (rad2 < tolIRMax2))
1673 || ((t2>=0) && (
eTheta > halfpi)
1674 && (tolIRMin2 < rad2) && (rad2 < tolIRMax2))
1675 || ((v.z()>0) && (
eTheta == halfpi)
1676 && (tolIRMin2 < rad2) && (rad2 < tolIRMax2)) )
1698 c = dist2ETheta/t1 ;
1707 && (sd < snxt) && (
eTheta > halfpi) )
1709 xi = p.x() + sd*v.x() ;
1710 yi = p.y() + sd*v.y() ;
1711 zi = p.z() + sd*v.z() ;
1712 rhoi2 = xi*xi + yi*yi ;
1713 radi2 = rhoi2 + zi*zi ;
1715 if ( (radi2 <= tolORMax2)
1716 && (radi2 >= tolORMin2)
1717 && (zi*(
eTheta - halfpi) <= 0) )
1746 c = dist2STheta/t1 ;
1754 if ((sd >= 0) && (sd < snxt))
1756 xi = p.x() + sd*v.x() ;
1757 yi = p.y() + sd*v.y() ;
1758 zi = p.z() + sd*v.z() ;
1759 rhoi2 = xi*xi + yi*yi ;
1760 radi2 = rhoi2 + zi*zi ;
1762 if ( (radi2 <= tolORMax2)
1763 && (radi2 >= tolORMin2)
1764 && (zi*(
fSTheta - halfpi) <= 0) )
1787 c = dist2ETheta/t1 ;
1795 if ((sd >= 0) && (sd < snxt))
1797 xi = p.x() + sd*v.x() ;
1798 yi = p.y() + sd*v.y() ;
1799 zi = p.z() + sd*v.z() ;
1800 rhoi2 = xi*xi + yi*yi ;
1801 radi2 = rhoi2 + zi*zi ;
1803 if ( (radi2 <= tolORMax2)
1804 && (radi2 >= tolORMin2)
1805 && (zi*(
eTheta - halfpi) <= 0) )
1838 G4double safe=0.0,safeRMin,safeRMax,safePhi,safeTheta;
1842 rho2=p.x()*p.x()+p.y()*p.y();
1843 rds=std::sqrt(rho2+p.z()*p.z());
1844 rho=std::sqrt(rho2);
1853 if (safeRMin>safeRMax)
1875 if (cosPsi<std::cos(
hDPhi))
1887 if (safePhi>safe) { safe=safePhi; }
1895 pTheta=std::acos(p.z()/rds);
1896 if (pTheta<0) { pTheta+=
pi; }
1899 if (dTheta1>dTheta2)
1903 safeTheta=rds*std::sin(dTheta1);
1904 if (safe<=safeTheta)
1914 safeTheta=rds*std::sin(dTheta2);
1915 if (safe<=safeTheta)
1923 if (safe<0) { safe=0; }
1951 G4double pDistS,compS,pDistE,compE,sphi2,vphi;
1953 G4double rho2,rad2,pDotV2d,pDotV3d;
1960 G4double dist2STheta, dist2ETheta, distTheta;
1965 rho2 = p.x()*p.x()+p.y()*p.y();
1966 rad2 = rho2+p.z()*p.z();
1968 pDotV2d = p.x()*v.x()+p.y()*v.y();
1969 pDotV3d = pDotV2d+p.z()*v.z();
1987 if( (rad2 <= Rmax_plus*Rmax_plus) && (rad2 >= Rmin_minus*Rmin_minus) )
2002 d2 = pDotV3d*pDotV3d - c;
2005 && ((pDotV3d >=0) || (d2 < 0)) )
2017 snxt = -pDotV3d+std::sqrt(d2);
2029 d2 = pDotV3d*pDotV3d - c;
2036 if(calcNorm) { *validNorm =
false; }
2043 sd = -pDotV3d-std::sqrt(d2);
2088 if ( std::fabs( p.z() ) <= halfRmaxTolerance )
2097 stheta = -p.z()/v.z();
2107 distTheta = std::sqrt(rho2)-p.z()*
tanSTheta;
2113 if(std::fabs(distTheta) < halfRmaxTolerance)
2115 if( (
fSTheta < halfpi) && (p.z() > 0.) )
2117 if( calcNorm ) { *validNorm =
false; }
2120 else if( (
fSTheta > halfpi) && (p.z() <= 0) )
2127 rhoSecTheta = std::sqrt(rho2*(1+
tanSTheta2));
2138 stheta = -0.5*dist2STheta/t2;
2144 if( std::fabs(distTheta) < halfRmaxTolerance )
2146 if( (
fSTheta > halfpi) && (t2 >= 0.) )
2153 rhoSecTheta = std::sqrt(rho2*(1+
tanSTheta2));
2163 else if( (
fSTheta < halfpi) && (t2 < 0.) && (p.z() >=0.) )
2165 if( calcNorm ) { *validNorm =
false; }
2181 if ( ((std::fabs(
s) < halfRmaxTolerance) && (t2 < 0.))
2182 || (sd < 0.) || ( (sd > 0.) && (p.z() + sd*v.z() > 0.) ) )
2186 if( (sd > halfRmaxTolerance) && (p.z() + sd*v.z() <= 0.) )
2196 if ( ( (std::fabs(sd) < halfRmaxTolerance) && (t2 >= 0.) )
2197 || (sd < 0.) || ( (sd > 0.) && (p.z() + sd*v.z() < 0.) ) )
2201 if( (sd > halfRmaxTolerance) && (p.z() + sd*v.z() >= 0.) )
2217 if ( std::fabs( p.z() ) <= halfRmaxTolerance )
2241 distTheta = std::sqrt(rho2)-p.z()*
tanETheta;
2247 if(std::fabs(distTheta) < halfRmaxTolerance)
2249 if( (
eTheta > halfpi) && (p.z() < 0.) )
2251 if( calcNorm ) { *validNorm =
false; }
2254 else if ( (
eTheta < halfpi) && (p.z() >= 0) )
2261 rhoSecTheta = std::sqrt(rho2*(1+
tanETheta2));
2271 sd = -0.5*dist2ETheta/t2;
2282 if ( std::fabs(distTheta) < halfRmaxTolerance )
2284 if( (
eTheta < halfpi) && (t2 >= 0.) )
2291 rhoSecTheta = std::sqrt(rho2*(1+
tanETheta2));
2300 else if ( (
eTheta > halfpi)
2301 && (t2 < 0.) && (p.z() <=0.) )
2303 if( calcNorm ) { *validNorm =
false; }
2310 if ( (d2 <halfRmaxTolerance) && (d2 > -halfRmaxTolerance) )
2322 if( ((std::fabs(sd) < halfRmaxTolerance) && (t2 < 0.))
2327 if( sd > halfRmaxTolerance )
2340 if ( ((std::fabs(sd) < halfRmaxTolerance) && (t2 >= 0.))
2341 || (sd < 0.) || ( (sd > 0.) && (p.z() + sd*v.z() > halfRmaxTolerance) ) )
2345 if( (sd > halfRmaxTolerance) && (p.z() + sd*v.z() <= halfRmaxTolerance) )
2365 if ( p.x() || p.y() )
2378 if ( (pDistS <= 0) && (pDistE <= 0) )
2384 sphi = pDistS/compS ;
2385 xi = p.x()+sphi*v.x() ;
2386 yi = p.y()+sphi*v.y() ;
2392 vphi = std::atan2(v.y(),v.x());
2414 sphi2=pDistE/compE ;
2417 xi = p.x()+sphi2*v.x() ;
2418 yi = p.y()+sphi2*v.y() ;
2426 vphi = std::atan2(v.y(),v.x()) ;
2433 else { sphi = 0.0; }
2451 else if ((pDistS >= 0) && (pDistE >= 0))
2453 if ( pDistS <= pDistE )
2463 if ( (compS < 0) && (compE < 0) ) { sphi = 0; }
2471 if ( (compS >= 0) && (compE >= 0) ) { sphi =
kInfinity; }
2475 else if ( (pDistS > 0) && (pDistE < 0) )
2483 sphi = pDistE/compE ;
2484 xi = p.x() + sphi*v.x() ;
2485 yi = p.y() + sphi*v.y() ;
2492 vphi = std::atan2(v.y(),v.x());
2521 sphi = pDistE/compE ;
2522 xi = p.x() + sphi*v.x() ;
2523 yi = p.y() + sphi*v.y() ;
2530 vphi = std::atan2(v.y(),v.x());
2566 xi=p.x()+sphi*v.x();
2567 yi=p.y()+sphi*v.y();
2574 vphi = std::atan2(v.y(),v.x()) ;
2603 sphi = pDistS/compS ;
2604 xi = p.x()+sphi*v.x() ;
2605 yi = p.y()+sphi*v.y() ;
2612 vphi = std::atan2(v.y(),v.x()) ;
2647 if ( v.x() || v.y() )
2649 vphi = std::atan2(v.y(),v.x()) ;
2682 xi=p.x()+snxt*v.x();
2683 yi=p.y()+snxt*v.y();
2684 zi=p.z()+snxt*v.z();
2699 else { *validNorm=
false; }
2708 else { *validNorm=
false; }
2719 xi = p.x() + snxt*v.x();
2720 yi = p.y() + snxt*v.y();
2724 rhoSecTheta = std::sqrt(rho2*(1+
tanSTheta2));
2734 else { *validNorm=
false; }
2743 else if (
eTheta < halfpi )
2745 xi=p.x()+snxt*v.x();
2746 yi=p.y()+snxt*v.y();
2750 rhoSecTheta = std::sqrt(rho2*(1+
tanETheta2));
2760 else { *validNorm=
false; }
2766 std::ostringstream message;
2767 G4int oldprc = message.precision(16);
2768 message <<
"Undefined side for valid surface normal to solid."
2770 <<
"Position:" << G4endl << G4endl
2771 <<
"p.x() = " << p.x()/
mm <<
" mm" << G4endl
2772 <<
"p.y() = " << p.y()/
mm <<
" mm" << G4endl
2773 <<
"p.z() = " << p.z()/
mm <<
" mm" << G4endl << G4endl
2774 <<
"Direction:" << G4endl << G4endl
2775 <<
"v.x() = " << v.x() << G4endl
2776 <<
"v.y() = " << v.y() << G4endl
2777 <<
"v.z() = " << v.z() << G4endl << G4endl
2778 <<
"Proposed distance :" << G4endl << G4endl
2779 <<
"snxt = " << snxt/
mm <<
" mm" <<
G4endl;
2780 message.precision(oldprc);
2790 std::ostringstream message;
2791 G4int oldprc = message.precision(16);
2792 message <<
"Logic error: snxt = kInfinity ???" << G4endl
2793 <<
"Position:" << G4endl << G4endl
2794 <<
"p.x() = " << p.x()/
mm <<
" mm" << G4endl
2795 <<
"p.y() = " << p.y()/
mm <<
" mm" << G4endl
2796 <<
"p.z() = " << p.z()/
mm <<
" mm" << G4endl << G4endl
2797 <<
"Rp = "<< std::sqrt( p.x()*p.x()+p.y()*p.y()+p.z()*p.z() )/
mm
2798 <<
" mm" << G4endl << G4endl
2799 <<
"Direction:" << G4endl << G4endl
2800 <<
"v.x() = " << v.x() << G4endl
2801 <<
"v.y() = " << v.y() << G4endl
2802 <<
"v.z() = " << v.z() << G4endl << G4endl
2803 <<
"Proposed distance :" << G4endl << G4endl
2804 <<
"snxt = " << snxt/
mm <<
" mm" <<
G4endl;
2805 message.precision(oldprc);
2819 G4double safe=0.0,safeRMin,safeRMax,safePhi,safeTheta;
2822 rho2=p.x()*p.x()+p.y()*p.y();
2823 rds=std::sqrt(rho2+p.z()*p.z());
2824 rho=std::sqrt(rho2);
2836 G4cout.precision(old_prc) ;
2838 "GeomSolids1002",
JustWarning,
"Point p is outside !?" );
2844 safeRMax =
fRmax-rds;
2848 safeRMin = rds-
fRmin;
2849 safe =
std::min( safeRMin, safeRMax );
2883 pTheta=std::acos(p.z()/rds);
2884 if (pTheta<0) { pTheta+=
pi; }
2888 { dTheta2=
eTheta-pTheta;}
2890 safeTheta=rds*std::sin(
std::min(dTheta1, dTheta2) );
2897 safe =
std::min( safe, safeTheta );
2900 if (safe<0.0) { safe=0; }
2919 G4int& noPolygonVertices )
const
2923 G4double meshAnglePhi,meshRMax,crossAnglePhi,
2924 coscrossAnglePhi,sincrossAnglePhi,sAnglePhi;
2925 G4double meshTheta,crossTheta,startTheta;
2926 G4double rMaxX,rMaxY,rMinX,rMinY,rMinZ,rMaxZ;
2927 G4int crossSectionPhi,noPhiCrossSections,crossSectionTheta,noThetaSections;
2941 meshAnglePhi=
fDPhi/(noPhiCrossSections-1);
2948 sAnglePhi = -meshAnglePhi*0.5;
2967 meshTheta=
fDTheta/(noThetaSections-1);
2974 startTheta = -meshTheta*0.5;
2981 meshRMax = (meshAnglePhi >= meshTheta) ?
2982 fRmax/std::cos(meshAnglePhi*0.5) :
fRmax/std::cos(meshTheta*0.5);
2986 if (vertices && cosCrossTheta && sinCrossTheta)
2988 vertices->reserve(noPhiCrossSections*(noThetaSections*2));
2989 for (crossSectionPhi=0;
2990 crossSectionPhi<noPhiCrossSections; crossSectionPhi++)
2992 crossAnglePhi=sAnglePhi+crossSectionPhi*meshAnglePhi;
2993 coscrossAnglePhi=std::cos(crossAnglePhi);
2994 sincrossAnglePhi=std::sin(crossAnglePhi);
2995 for (crossSectionTheta=0;
2996 crossSectionTheta<noThetaSections;crossSectionTheta++)
3000 crossTheta=startTheta+crossSectionTheta*meshTheta;
3001 cosCrossTheta[crossSectionTheta]=std::cos(crossTheta);
3002 sinCrossTheta[crossSectionTheta]=std::sin(crossTheta);
3004 rMinX=
fRmin*sinCrossTheta[crossSectionTheta]*coscrossAnglePhi;
3005 rMinY=
fRmin*sinCrossTheta[crossSectionTheta]*sincrossAnglePhi;
3006 rMinZ=
fRmin*cosCrossTheta[crossSectionTheta];
3013 for (crossSectionTheta=noThetaSections-1;
3014 crossSectionTheta>=0; crossSectionTheta--)
3016 rMaxX=meshRMax*sinCrossTheta[crossSectionTheta]*coscrossAnglePhi;
3017 rMaxY=meshRMax*sinCrossTheta[crossSectionTheta]*sincrossAnglePhi;
3018 rMaxZ=meshRMax*cosCrossTheta[crossSectionTheta];
3025 noPolygonVertices = noThetaSections*2 ;
3032 "Error in allocation of vertices. Out of memory !");
3035 delete [] cosCrossTheta;
3036 delete [] sinCrossTheta;
3065 G4int oldprc = os.precision(16);
3066 os <<
"-----------------------------------------------------------\n"
3067 <<
" *** Dump for solid - " <<
GetName() <<
" ***\n"
3068 <<
" ===================================================\n"
3069 <<
" Solid type: G4Sphere\n"
3070 <<
" Parameters: \n"
3071 <<
" inner radius: " <<
fRmin/
mm <<
" mm \n"
3072 <<
" outer radius: " <<
fRmax/
mm <<
" mm \n"
3073 <<
" starting phi of segment : " <<
fSPhi/
degree <<
" degrees \n"
3074 <<
" delta phi of segment : " <<
fDPhi/
degree <<
" degrees \n"
3075 <<
" starting theta of segment: " <<
fSTheta/
degree <<
" degrees \n"
3076 <<
" delta theta of segment : " <<
fDTheta/
degree <<
" degrees \n"
3077 <<
"-----------------------------------------------------------\n";
3078 os.precision(oldprc);
3089 G4double zRand, aOne, aTwo, aThr, aFou, aFiv, chose, phi, sinphi, cosphi;
3090 G4double height1, height2, slant1, slant2, costheta, sintheta, rRand;
3105 cosphi = std::cos(phi);
3106 sinphi = std::sin(phi);
3108 sintheta = std::sqrt(1.-
sqr(costheta));
3117 if( (chose>=0.) && (chose<aOne) )
3122 else if( (chose>=aOne) && (chose<aOne+aTwo) )
3127 else if( (chose>=aOne+aTwo) && (chose<aOne+aTwo+aThr) )
3140 else if( (chose>=aOne+aTwo+aThr) && (chose<aOne+aTwo+aThr+aFou) )
3153 else if( (chose>=aOne+aTwo+aThr+aFou) && (chose<aOne+aTwo+aThr+aFou+aFiv) )
3156 rRand*sintheta*
sinSPhi,rRand*costheta);
3161 rRand*sintheta*
sinEPhi,rRand*costheta);
ThreeVector shoot(const G4int Ap, const G4int Af)
G4double halfCarTolerance
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
static const G4double kInfinity
G4double GetMinYExtent() const
CLHEP::Hep3Vector G4ThreeVector
G4bool IsYLimited() const
G4double GetRadiusInRing(G4double rmin, G4double rmax) const
G4Polyhedron * CreatePolyhedron() const
G4bool IsXLimited() const
const G4int kMinMeshSections
virtual void AddSolid(const G4Box &)=0
G4Sphere(const G4String &pName, G4double pRmin, G4double pRmax, G4double pSPhi, G4double pDPhi, G4double pSTheta, G4double pDTheta)
G4ThreeVector GetPointOnSurface() const
G4double GetMaxXExtent() const
G4double GetMinZExtent() const
G4GLOB_DLL std::ostream G4cout
void CalculateClippedPolygonExtent(G4ThreeVectorList &pPolygon, const G4VoxelLimits &pVoxelLimit, const EAxis pAxis, G4double &pMin, G4double &pMax) const
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pmin, G4double &pmax) const
void CheckThetaAngles(G4double sTheta, G4double dTheta)
G4double GetRadialTolerance() const
G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=G4bool(false), G4bool *validNorm=0, G4ThreeVector *n=0) const
G4ThreeVectorList * CreateRotatedVertices(const G4AffineTransform &pTransform, G4int &noPolygonVertices) const
EInside Inside(const G4ThreeVector &p) const
void CheckPhiAngles(G4double sPhi, G4double dPhi)
std::vector< G4ThreeVector > G4ThreeVectorList
G4GeometryType GetEntityType() const
G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
G4ThreeVector ApproxSurfaceNormal(const G4ThreeVector &p) const
virtual void ComputeDimensions(G4Box &, const G4int, const G4VPhysicalVolume *) const
G4Sphere & operator=(const G4Sphere &rhs)
G4double GetMinXExtent() const
G4double GetMaxZExtent() const
T max(const T t1, const T t2)
brief Return the largest of the two arguments
std::ostream & StreamInfo(std::ostream &os) const
void ComputeDimensions(G4VPVParameterisation *p, const G4int n, const G4VPhysicalVolume *pRep)
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
static const double degree
G4VisExtent GetExtent() const
G4double GetSurfaceArea()
G4double GetMaxYExtent() const
G4double GetMaxExtent(const EAxis pAxis) const
G4CSGSolid & operator=(const G4CSGSolid &rhs)
G4bool IsZLimited() const
G4double halfAngTolerance
G4double GetAngularTolerance() const
G4double GetMinExtent(const EAxis pAxis) const
const G4double kMeshAngleDefault
static G4GeometryTolerance * GetInstance()
const G4int kMaxMeshSections
void DescribeYourselfTo(G4VGraphicsScene &scene) const