37 #if !defined(G4GEOM_USE_UTUBS)
53 using namespace CLHEP;
64 :
G4CSGSolid(pName), fRMin(pRMin), fRMax(pRMax), fDz(pDz),
66 fInvRmax( pRMax > 0.0 ? 1.0/pRMax : 0.0 ),
67 fInvRmin( pRMin > 0.0 ? 1.0/pRMin : 0.0 )
79 message <<
"Negative Z half-length (" << pDz <<
") in solid: " <<
GetName();
82 if ( (pRMin >= pRMax) || (pRMin < 0) )
85 message <<
"Invalid values for radii in solid: " <<
GetName()
87 <<
" pRMin = " << pRMin <<
", pRMax = " << pRMax;
102 :
G4CSGSolid(a), kRadTolerance(0.), kAngTolerance(0.),
103 fRMin(0.), fRMax(0.), fDz(0.), fSPhi(0.), fDPhi(0.),
104 sinCPhi(0.), cosCPhi(0.), cosHDPhi(0.), cosHDPhiOT(0.), cosHDPhiIT(0.),
105 sinSPhi(0.), cosSPhi(0.), sinEPhi(0.), cosEPhi(0.),
106 fPhiFullTube(
false), fInvRmax(0.), fInvRmin(0.),
107 halfCarTolerance(0.), halfRadTolerance(0.),
126 kRadTolerance(rhs.kRadTolerance), kAngTolerance(rhs.kAngTolerance),
127 fRMin(rhs.fRMin), fRMax(rhs.fRMax), fDz(rhs.fDz),
128 fSPhi(rhs.fSPhi), fDPhi(rhs.fDPhi),
129 sinCPhi(rhs.sinCPhi), cosCPhi(rhs.cosCPhi), cosHDPhi(rhs.cosHDPhi),
130 cosHDPhiOT(rhs.cosHDPhiOT), cosHDPhiIT(rhs.cosHDPhiIT),
131 sinSPhi(rhs.sinSPhi), cosSPhi(rhs.cosSPhi),
132 sinEPhi(rhs.sinEPhi), cosEPhi(rhs.cosEPhi), fPhiFullTube(rhs.fPhiFullTube),
133 fInvRmax(rhs.fInvRmax), fInvRmin(rhs.fInvRmin),
134 halfCarTolerance(rhs.halfCarTolerance),
135 halfRadTolerance(rhs.halfRadTolerance),
136 halfAngTolerance(rhs.halfAngTolerance)
148 if (
this == &rhs) {
return *
this; }
204 pMin.
set(vmin.
x(),vmin.
y(),-
dz);
205 pMax.
set(vmax.
x(),vmax.
y(),
dz);
209 pMin.
set(-rmax,-rmax,-dz);
210 pMax.
set( rmax, rmax, dz);
215 if (pMin.
x() >= pMax.
x() || pMin.
y() >= pMax.
y() || pMin.
z() >= pMax.
z())
218 message <<
"Bad bounding box (min >= max) for solid: "
220 <<
"\npMin = " << pMin
221 <<
"\npMax = " <<
pMax;
222 G4Exception(
"G4Tubs::BoundingLimits()",
"GeomMgt0001",
251 return exist = (pMin <
pMax) ?
true :
false;
262 const G4int NSTEPS = 24;
264 G4int ksteps = (dphi <= astep) ? 1 : (
G4int)((dphi-
deg)/astep) + 1;
267 G4double sinHalf = std::sin(0.5*ang);
268 G4double cosHalf = std::cos(0.5*ang);
269 G4double sinStep = 2.*sinHalf*cosHalf;
270 G4double cosStep = 1. - 2.*sinHalf*sinHalf;
275 if (rmin == 0 && dphi ==
twopi)
283 baseA[
k].set(rext*cosCur,rext*sinCur,-dz);
284 baseB[
k].set(rext*cosCur,rext*sinCur, dz);
287 sinCur = sinCur*cosStep + cosCur*sinStep;
288 cosCur = cosCur*cosStep - sinTmp*sinStep;
290 std::vector<const G4ThreeVectorList *> polygons(2);
291 polygons[0] = &baseA;
292 polygons[1] = &baseB;
302 G4double sinCur = sinStart*cosHalf + cosStart*sinHalf;
303 G4double cosCur = cosStart*cosHalf - sinStart*sinHalf;
307 for (
G4int k=0;
k<ksteps+2; ++
k) pols[
k].resize(4);
308 pols[0][0].set(rmin*cosStart,rmin*sinStart, dz);
309 pols[0][1].set(rmin*cosStart,rmin*sinStart,-dz);
310 pols[0][2].set(rmax*cosStart,rmax*sinStart,-dz);
311 pols[0][3].set(rmax*cosStart,rmax*sinStart, dz);
314 pols[
k][0].set(rmin*cosCur,rmin*sinCur, dz);
315 pols[
k][1].set(rmin*cosCur,rmin*sinCur,-dz);
316 pols[
k][2].set(rext*cosCur,rext*sinCur,-dz);
317 pols[
k][3].set(rext*cosCur,rext*sinCur, dz);
320 sinCur = sinCur*cosStep + cosCur*sinStep;
321 cosCur = cosCur*cosStep - sinTmp*sinStep;
323 pols[ksteps+1][0].set(rmin*cosEnd,rmin*sinEnd, dz);
324 pols[ksteps+1][1].set(rmin*cosEnd,rmin*sinEnd,-dz);
325 pols[ksteps+1][2].set(rmax*cosEnd,rmax*sinEnd,-dz);
326 pols[ksteps+1][3].set(rmax*cosEnd,rmax*sinEnd, dz);
329 std::vector<const G4ThreeVectorList *> polygons;
330 polygons.resize(ksteps+2);
331 for (
G4int k=0;
k<ksteps+2; ++
k) polygons[
k] = &pols[
k];
349 r2 = p.
x()*p.
x() + p.
y()*p.
y() ;
352 else { tolRMin = 0 ; }
356 if ((r2 >= tolRMin*tolRMin) && (r2 <= tolRMax*tolRMax))
374 pPhi = std::atan2(p.
y(),p.
x()) ;
417 if ( tolRMin < 0 ) { tolRMin = 0; }
419 if ( (r2 >= tolRMin*tolRMin) && (r2 <= tolRMax*tolRMax) )
421 if (
fPhiFullTube || (r2 <=halfRadTolerance*halfRadTolerance) )
427 pPhi = std::atan2(p.
y(),p.
x()) ;
458 r2 = p.
x()*p.
x() + p.
y()*p.
y() ;
462 if ( tolRMin < 0 ) { tolRMin = 0; }
464 if ( (r2 >= tolRMin*tolRMin) && (r2 <= tolRMax*tolRMax) )
466 if (
fPhiFullTube || (r2 <=halfRadTolerance*halfRadTolerance))
472 pPhi = std::atan2(p.
y(),p.
x()) ;
511 G4int noSurfaces = 0;
520 rho = std::sqrt(p.
x()*p.
x() + p.
y()*p.
y());
522 distRMin = std::fabs(rho -
fRMin);
523 distRMax = std::fabs(rho -
fRMax);
524 distZ = std::fabs(std::fabs(p.
z()) -
fDz);
530 pPhi = std::atan2(p.
y(),p.
x());
535 distSPhi = std::fabs( pPhi -
fSPhi );
574 if ( p.
z() >= 0.) { sumnorm += nZ; }
575 else { sumnorm -= nZ; }
577 if ( noSurfaces == 0 )
580 G4Exception(
"G4Tubs::SurfaceNormal(p)",
"GeomSolids1002",
583 G4cout<<
"G4Tubs::SN ( "<<p.
x()<<
", "<<p.
y()<<
", "<<p.
z()<<
" ); "
585 G4cout.precision(oldprc) ;
589 else if ( noSurfaces == 1 ) { norm = sumnorm; }
590 else { norm = sumnorm.
unit(); }
605 G4double distZ, distRMin, distRMax, distSPhi, distEPhi, distMin ;
607 rho = std::sqrt(p.
x()*p.
x() + p.
y()*p.
y()) ;
609 distRMin = std::fabs(rho -
fRMin) ;
610 distRMax = std::fabs(rho -
fRMax) ;
611 distZ = std::fabs(std::fabs(p.
z()) -
fDz) ;
613 if (distRMin < distRMax)
615 if ( distZ < distRMin )
628 if ( distZ < distRMax )
641 phi = std::atan2(p.
y(),p.
x()) ;
643 if ( phi < 0 ) { phi +=
twopi; }
647 distSPhi = std::fabs(phi - (
fSPhi +
twopi))*rho ;
651 distSPhi = std::fabs(phi -
fSPhi)*rho ;
653 distEPhi = std::fabs(phi -
fSPhi -
fDPhi)*rho ;
655 if (distSPhi < distEPhi)
657 if ( distSPhi < distMin )
664 if ( distEPhi < distMin )
703 "Undefined side for valid surface normal to solid.");
737 G4double tolORMax2, tolIRMin2, tolODz, tolIDz ;
742 G4double Dist, sd, xi, yi, zi, rho2, inum, iden, cosPsi, Comp ;
765 if (std::fabs(p.
z()) >= tolIDz)
767 if ( p.
z()*v.
z() < 0 )
769 sd = (std::fabs(p.
z()) -
fDz)/std::fabs(v.
z()) ;
771 if(sd < 0.0) { sd = 0.0; }
773 xi = p.
x() + sd*v.
x() ;
774 yi = p.
y() + sd*v.
y() ;
775 rho2 = xi*xi + yi*yi ;
779 if ((tolIRMin2 <= rho2) && (rho2 <= tolIRMax2))
786 iden = std::sqrt(rho2) ;
798 if ( snxt<halfCarTolerance ) { snxt=0; }
815 t1 = 1.0 - v.
z()*v.
z() ;
816 t2 = p.
x()*v.
x() + p.
y()*v.
y() ;
817 t3 = p.
x()*p.
x() + p.
y()*p.
y() ;
823 if ((t3 >= tolORMax2) && (t2<0))
833 sd = c/(-b+std::sqrt(d));
838 G4double fTerm = sd-std::fmod(sd,dRmax);
843 zi = p.
z() + sd*v.
z() ;
844 if (std::fabs(zi)<=tolODz)
854 xi = p.
x() + sd*v.
x() ;
855 yi = p.
y() + sd*v.
y() ;
868 if ((t3 > tolIRMin2) && (t2 < 0) && (std::fabs(p.
z()) <= tolIDz))
875 iden = std::sqrt(t3) ;
896 snxt = c/(-b+std::sqrt(d));
898 if ( snxt < halfCarTolerance ) { snxt=0; }
916 c = t3 - fRMax*
fRMax;
927 snxt= c/(-b+std::sqrt(d));
929 if ( snxt < halfCarTolerance ) { snxt=0; }
949 sd =( b > 0. )? c/(-b - std::sqrt(d)) : (-b + std::sqrt(d));
950 if (sd >= -halfCarTolerance)
954 if(sd < 0.0) { sd = 0.0; }
957 G4double fTerm = sd-std::fmod(sd,dRmax);
960 zi = p.
z() + sd*v.
z() ;
961 if (std::fabs(zi) <= tolODz)
971 xi = p.
x() + sd*v.
x() ;
972 yi = p.
y() + sd*v.
y() ;
1007 if ( Dist < halfCarTolerance )
1013 if ( sd < 0 ) { sd = 0.0; }
1014 zi = p.
z() + sd*v.
z() ;
1015 if ( std::fabs(zi) <= tolODz )
1017 xi = p.
x() + sd*v.
x() ;
1018 yi = p.
y() + sd*v.
y() ;
1019 rho2 = xi*xi + yi*yi ;
1021 if ( ( (rho2 >= tolIRMin2) && (rho2 <= tolIRMax2) )
1022 || ( (rho2 > tolORMin2) && (rho2 < tolIRMin2)
1025 || ( (rho2 > tolIRMax2) && (rho2 < tolORMax2)
1047 if ( Dist < halfCarTolerance )
1053 if ( sd < 0 ) { sd = 0; }
1054 zi = p.
z() + sd*v.
z() ;
1055 if ( std::fabs(zi) <= tolODz )
1057 xi = p.
x() + sd*v.
x() ;
1058 yi = p.
y() + sd*v.
y() ;
1059 rho2 = xi*xi + yi*yi ;
1060 if ( ( (rho2 >= tolIRMin2) && (rho2 <= tolIRMax2) )
1061 || ( (rho2 > tolORMin2) && (rho2 < tolIRMin2)
1064 || ( (rho2 > tolIRMax2) && (rho2 < tolORMax2)
1078 if ( snxt<halfCarTolerance ) { snxt=0; }
1110 G4double safe=0.0, rho, safe1, safe2, safe3 ;
1113 rho = std::sqrt(p.
x()*p.
x() + p.
y()*p.
y()) ;
1114 safe1 =
fRMin - rho ;
1115 safe2 = rho -
fRMax ;
1116 safe3 = std::fabs(p.
z()) -
fDz ;
1118 if ( safe1 > safe2 ) { safe = safe1; }
1119 else { safe = safe2; }
1120 if ( safe3 > safe ) { safe = safe3; }
1140 if ( safePhi > safe ) { safe = safePhi; }
1143 if ( safe < 0 ) { safe = 0; }
1164 G4double pDistS, compS, pDistE, compE, sphi2, xi, yi, vphi, roi2 ;
1170 pdist =
fDz - p.
z() ;
1173 snxt = pdist/v.
z() ;
1186 else if ( v.
z() < 0 )
1188 pdist =
fDz + p.
z() ;
1192 snxt = -pdist/v.
z() ;
1222 t1 = 1.0 - v.
z()*v.
z() ;
1223 t2 = p.
x()*v.
x() + p.
y()*v.
y() ;
1224 t3 = p.
x()*p.
x() + p.
y()*p.
y() ;
1227 else { roi2 = snxt*snxt*t1 + 2*snxt*t2 + t3; }
1247 if( d2 >= 0 ) { srd = c/( -b - std::sqrt(d2)); }
1267 roMin2 = t3 - t2*t2/
t1 ;
1283 srd = c/(-b+std::sqrt(d2));
1301 srd = -b + std::sqrt(d2) ;
1326 srd = -b + std::sqrt(d2) ;
1350 vphi = std::atan2(v.
y(),v.
x()) ;
1356 if ( p.
x() || p.
y() )
1379 sphi = pDistS/compS ;
1383 xi = p.
x() + sphi*v.
x() ;
1384 yi = p.
y() + sphi*v.
y() ;
1423 sphi2 = pDistE/compE ;
1429 xi = p.
x() + sphi2*v.
x() ;
1430 yi = p.
y() + sphi2*v.
y() ;
1441 else { sphi = 0.0 ; }
1452 else { sphi = 0.0 ; }
1498 xi = p.
x() + snxt*v.
x() ;
1499 yi = p.
y() + snxt*v.
y() ;
1505 *validNorm =
false ;
1516 *validNorm =
false ;
1528 *validNorm =
false ;
1546 G4int oldprc = message.precision(16);
1547 message <<
"Undefined side for valid surface normal to solid."
1549 <<
"Position:" << G4endl << G4endl
1550 <<
"p.x() = " << p.
x()/
mm <<
" mm" << G4endl
1551 <<
"p.y() = " << p.
y()/
mm <<
" mm" << G4endl
1552 <<
"p.z() = " << p.
z()/
mm <<
" mm" << G4endl << G4endl
1553 <<
"Direction:" << G4endl << G4endl
1554 <<
"v.x() = " << v.
x() << G4endl
1555 <<
"v.y() = " << v.
y() << G4endl
1556 <<
"v.z() = " << v.
z() << G4endl << G4endl
1557 <<
"Proposed distance :" << G4endl << G4endl
1558 <<
"snxt = " << snxt/
mm <<
" mm" <<
G4endl ;
1559 message.precision(oldprc) ;
1560 G4Exception(
"G4Tubs::DistanceToOut(p,v,..)",
"GeomSolids1002",
1576 G4double safe=0.0, rho, safeR1, safeR2, safeZ, safePhi ;
1577 rho = std::sqrt(p.
x()*p.
x() + p.
y()*p.
y()) ;
1589 G4cout.precision(oldprc) ;
1590 G4Exception(
"G4Tubs::DistanceToOut(p)",
"GeomSolids1002",
1597 safeR1 = rho -
fRMin ;
1598 safeR2 =
fRMax - rho ;
1600 if ( safeR1 < safeR2 ) { safe = safeR1 ; }
1601 else { safe = safeR2 ; }
1605 safe =
fRMax - rho ;
1607 safeZ =
fDz - std::fabs(p.
z()) ;
1609 if ( safeZ < safe ) { safe = safeZ ; }
1623 if (safePhi < safe) { safe = safePhi ; }
1625 if ( safe < 0 ) { safe = 0 ; }
1645 return new G4Tubs(*
this);
1654 G4int oldprc = os.precision(16);
1655 os <<
"-----------------------------------------------------------\n"
1656 <<
" *** Dump for solid - " <<
GetName() <<
" ***\n"
1657 <<
" ===================================================\n"
1658 <<
" Solid type: G4Tubs\n"
1659 <<
" Parameters: \n"
1660 <<
" inner radius : " <<
fRMin/
mm <<
" mm \n"
1661 <<
" outer radius : " <<
fRMax/
mm <<
" mm \n"
1662 <<
" half length Z: " <<
fDz/
mm <<
" mm \n"
1663 <<
" starting phi : " <<
fSPhi/
degree <<
" degrees \n"
1664 <<
" delta phi : " <<
fDPhi/
degree <<
" degrees \n"
1665 <<
"-----------------------------------------------------------\n";
1666 os.precision(oldprc);
1677 G4double xRand, yRand, zRand,
phi, cosphi, sinphi, chose,
1678 aOne, aTwo, aThr, aFou;
1687 cosphi = std::cos(phi);
1688 sinphi = std::sin(phi);
1696 if( (chose >=0) && (chose < aOne) )
1698 xRand = fRMax*cosphi;
1699 yRand = fRMax*sinphi;
1703 else if( (chose >= aOne) && (chose < aOne + aTwo) )
1705 xRand = fRMin*cosphi;
1706 yRand = fRMin*sinphi;
1710 else if( (chose >= aOne + aTwo) && (chose < aOne + aTwo + aThr) )
1712 xRand = rRand*cosphi;
1713 yRand = rRand*sinphi;
1717 else if( (chose >= aOne + aTwo + aThr) && (chose < aOne + aTwo + 2.*aThr) )
1719 xRand = rRand*cosphi;
1720 yRand = rRand*sinphi;
1724 else if( (chose >= aOne + aTwo + 2.*aThr)
1725 && (chose < aOne + aTwo + 2.*aThr + aFou) )