51 #if !defined(G4GEOM_USE_UCONS)
65 using namespace CLHEP;
87 :
G4CSGSolid(pName), fRmin1(pRmin1), fRmin2(pRmin2),
88 fRmax1(pRmax1), fRmax2(pRmax2), fDz(pDz), fSPhi(0.), fDPhi(0.)
101 std::ostringstream message;
102 message <<
"Invalid Z half-length for Solid: " <<
GetName() <<
G4endl
110 if (((pRmin1>=pRmax1) || (pRmin2>=pRmax2) || (pRmin1<0)) && (pRmin2<0))
112 std::ostringstream message;
113 message <<
"Invalid values of radii for Solid: " <<
GetName() <<
G4endl
114 <<
" pRmin1 = " << pRmin1 <<
", pRmin2 = " << pRmin2
115 <<
", pRmax1 = " << pRmax1 <<
", pRmax2 = " << pRmax2;
133 :
G4CSGSolid(a), kRadTolerance(0.), kAngTolerance(0.),
134 fRmin1(0.), fRmin2(0.), fRmax1(0.), fRmax2(0.), fDz(0.),
135 fSPhi(0.), fDPhi(0.), sinCPhi(0.), cosCPhi(0.), cosHDPhiOT(0.),
136 cosHDPhiIT(0.), sinSPhi(0.), cosSPhi(0.), sinEPhi(0.), cosEPhi(0.),
137 fPhiFullCone(false), halfCarTolerance(0.), halfRadTolerance(0.),
155 :
G4CSGSolid(rhs), kRadTolerance(rhs.kRadTolerance),
156 kAngTolerance(rhs.kAngTolerance), fRmin1(rhs.fRmin1), fRmin2(rhs.fRmin2),
157 fRmax1(rhs.fRmax1), fRmax2(rhs.fRmax2), fDz(rhs.fDz), fSPhi(rhs.fSPhi),
158 fDPhi(rhs.fDPhi), sinCPhi(rhs.sinCPhi), cosCPhi(rhs.cosCPhi),
159 cosHDPhiOT(rhs.cosHDPhiOT), cosHDPhiIT(rhs.cosHDPhiIT),
160 sinSPhi(rhs.sinSPhi), cosSPhi(rhs.cosSPhi), sinEPhi(rhs.sinEPhi),
161 cosEPhi(rhs.cosEPhi), fPhiFullCone(rhs.fPhiFullCone),
162 halfCarTolerance(rhs.halfCarTolerance),
163 halfRadTolerance(rhs.halfRadTolerance),
164 halfAngTolerance(rhs.halfAngTolerance)
176 if (
this == &rhs) {
return *
this; }
207 G4double r2, rl, rh, pPhi, tolRMin, tolRMax;
214 r2 = p.x()*p.x() + p.y()*p.y() ;
221 if ( tolRMin < 0 ) { tolRMin = 0; }
224 if ( (r2<tolRMin*tolRMin) || (r2>tolRMax*tolRMax) ) {
return in =
kOutside; }
227 else { tolRMin = 0.0; }
232 if ( (r2 < tolRMin*tolRMin) || (r2 >= tolRMax*tolRMax) ) { in =
kSurface; }
234 if ( !
fPhiFullCone && ((p.x() != 0.0) || (p.y() != 0.0)) )
236 pPhi = std::atan2(p.y(),p.x()) ;
290 G4double diff1, diff2, delta, maxDiff, newMin, newMax, RMax ;
291 G4double xoff1, xoff2, yoff1, yoff2 ;
294 zMin = zoffset -
fDz ;
295 zMax = zoffset +
fDz ;
320 xMin = 2*xoffset-xMax ;
344 yMin = 2*yoffset-yMax ;
345 RMax = yMax - yoffset ;
369 yoff1 = yoffset - yMin ;
370 yoff2 = yMax - yoffset ;
372 if ((yoff1 >= 0) && (yoff2 >= 0))
381 delta=RMax*RMax-yoff1*yoff1;
382 diff1=(delta>0.) ? std::sqrt(delta) : 0.;
383 delta=RMax*RMax-yoff2*yoff2;
384 diff2=(delta>0.) ? std::sqrt(delta) : 0.;
385 maxDiff = (diff1>diff2) ? diff1:diff2 ;
386 newMin = xoffset - maxDiff ;
387 newMax = xoffset + maxDiff ;
388 pMin = ( newMin < xMin ) ? xMin : newMin ;
389 pMax = ( newMax > xMax) ? xMax : newMax ;
394 xoff1 = xoffset - xMin ;
395 xoff2 = xMax - xoffset ;
397 if ((xoff1 >= 0) && (xoff2 >= 0) )
406 delta=RMax*RMax-xoff1*xoff1;
407 diff1=(delta>0.) ? std::sqrt(delta) : 0.;
408 delta=RMax*RMax-xoff2*xoff2;
409 diff2=(delta>0.) ? std::sqrt(delta) : 0.;
410 maxDiff = (diff1 > diff2) ? diff1:diff2 ;
411 newMin = yoffset - maxDiff ;
412 newMax = yoffset + maxDiff ;
413 pMin = (newMin < yMin) ? yMin : newMin ;
414 pMax = (newMax > yMax) ? yMax : newMax ;
433 G4int i, noEntries, noBetweenSections4 ;
434 G4bool existsAfterClip = false ;
440 noEntries = vertices->size() ;
441 noBetweenSections4 = noEntries-4 ;
443 for ( i = 0 ; i < noEntries ; i += 4 )
447 for ( i = 0 ; i < noBetweenSections4 ; i += 4 )
453 existsAfterClip = true ;
474 existsAfterClip = true ;
480 return existsAfterClip ;
492 G4int noSurfaces = 0;
496 G4double tanRMin, secRMin, pRMin, widRMin;
497 G4double tanRMax, secRMax, pRMax, widRMax;
502 distZ = std::fabs(std::fabs(p.z()) -
fDz);
503 rho = std::sqrt(p.x()*p.x() + p.y()*p.y());
506 secRMin = std::sqrt(1 + tanRMin*tanRMin);
507 pRMin = rho - p.z()*tanRMin;
509 distRMin = std::fabs(pRMin - widRMin)/secRMin;
512 secRMax = std::sqrt(1+tanRMax*tanRMax);
513 pRMax = rho - p.z()*tanRMax;
515 distRMax = std::fabs(pRMax - widRMax)/secRMax;
521 pPhi = std::atan2(p.y(),p.x());
526 distSPhi = std::fabs( pPhi -
fSPhi );
539 nR =
G4ThreeVector(p.x()/rho/secRMax, p.y()/rho/secRMax, -tanRMax/secRMax);
542 nr =
G4ThreeVector(-p.x()/rho/secRMin,-p.y()/rho/secRMin,tanRMin/secRMin);
572 if ( p.z() >= 0.) { sumnorm += nZ; }
573 else { sumnorm -= nZ; }
575 if ( noSurfaces == 0 )
578 G4Exception(
"G4Cons::SurfaceNormal(p)",
"GeomSolids1002",
583 else if ( noSurfaces == 1 ) { norm = sumnorm; }
584 else { norm = sumnorm.unit(); }
599 G4double distZ, distRMin, distRMax, distSPhi, distEPhi, distMin ;
600 G4double tanRMin, secRMin, pRMin, widRMin ;
601 G4double tanRMax, secRMax, pRMax, widRMax ;
603 distZ = std::fabs(std::fabs(p.z()) -
fDz) ;
604 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
607 secRMin = std::sqrt(1 + tanRMin*tanRMin) ;
608 pRMin = rho - p.z()*tanRMin ;
610 distRMin = std::fabs(pRMin - widRMin)/secRMin ;
613 secRMax = std::sqrt(1+tanRMax*tanRMax) ;
614 pRMax = rho - p.z()*tanRMax ;
616 distRMax = std::fabs(pRMax - widRMax)/secRMax ;
618 if (distRMin < distRMax)
620 if (distZ < distRMin)
633 if (distZ < distRMax)
646 phi = std::atan2(p.y(),p.x()) ;
648 if (phi < 0) { phi +=
twopi; }
651 else { distSPhi = std::fabs(phi -
fSPhi)*rho; }
653 distEPhi = std::fabs(phi -
fSPhi -
fDPhi)*rho ;
657 if (distSPhi < distEPhi)
659 if (distSPhi < distMin) { side =
kNSPhi; }
663 if (distEPhi < distMin) { side =
kNEPhi; }
670 norm =
G4ThreeVector(-p.x()/rho, -p.y()/rho, tanRMin/secRMin) ;
674 norm =
G4ThreeVector(p.x()/rho, p.y()/rho, -tanRMax/secRMax) ;
690 "Undefined side for valid surface normal to solid.");
726 G4double tanRMax,secRMax,rMaxAv,rMaxOAv ;
727 G4double tanRMin,secRMin,rMinAv,rMinOAv ;
730 G4double tolORMin,tolORMin2,tolIRMin,tolIRMin2 ;
731 G4double tolORMax2,tolIRMax,tolIRMax2 ;
734 G4double Dist,sd,xi,yi,zi,ri=0.,risec,rhoi2,cosPsi ;
745 secRMin = std::sqrt(1.0 + tanRMin*tanRMin) ;
757 secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
766 if (std::fabs(p.z()) >= tolIDz)
768 if ( p.z()*v.z() < 0 )
770 sd = (std::fabs(p.z()) -
fDz)/std::fabs(v.z()) ;
772 if( sd < 0.0 ) { sd = 0.0; }
774 xi = p.x() + sd*v.x() ;
775 yi = p.y() + sd*v.y() ;
776 rhoi2 = xi*xi + yi*yi ;
783 tolORMin =
fRmin1 - halfRadTolerance*secRMin ;
784 tolIRMin =
fRmin1 + halfRadTolerance*secRMin ;
785 tolIRMax =
fRmax1 - halfRadTolerance*secRMin ;
791 tolORMin =
fRmin2 - halfRadTolerance*secRMin ;
792 tolIRMin =
fRmin2 + halfRadTolerance*secRMin ;
793 tolIRMax =
fRmax2 - halfRadTolerance*secRMin ;
800 tolIRMin2 = tolIRMin*tolIRMin ;
807 if ( tolIRMax > 0 ) { tolIRMax2 = tolIRMax*tolIRMax; }
808 else { tolIRMax2 = 0.0; }
810 if ( (tolIRMin2 <= rhoi2) && (rhoi2 <= tolIRMax2) )
851 t1 = 1.0 - v.z()*v.z() ;
852 t2 = p.x()*v.x() + p.y()*v.y() ;
853 t3 = p.x()*p.x() + p.y()*p.y() ;
854 rin = tanRMin*p.z() + rMinAv ;
855 rout = tanRMax*p.z() + rMaxAv ;
860 nt1 = t1 - (tanRMax*v.z())*(tanRMax*v.z()) ;
861 nt2 = t2 - tanRMax*v.z()*rout ;
862 nt3 = t3 - rout*rout ;
878 if ((rout < 0) && (nt3 <= 0))
883 if (b>0) { sd = c/(-b-std::sqrt(d)); }
884 else { sd = -b + std::sqrt(d); }
888 if ((b <= 0) && (c >= 0))
890 sd=c/(-b+std::sqrt(d));
896 sd = -b + std::sqrt(d) ;
897 if((sd<0) & (sd>-halfRadTolerance)) sd=0;
909 G4double fTerm = sd-std::fmod(sd,dRmax);
912 zi = p.z() + sd*v.z() ;
914 if (std::fabs(zi) <= tolODz)
921 xi = p.x() + sd*v.x() ;
922 yi = p.y() + sd*v.y() ;
923 ri = rMaxAv + zi*tanRMax ;
937 if ( ( t3 > (rin + halfRadTolerance*secRMin)*
938 (rin + halfRadTolerance*secRMin) )
939 && (nt2 < 0) && (d >= 0) && (std::fabs(p.z()) <= tolIDz) )
946 risec = std::sqrt(xi*xi + yi*yi)*secRMax ;
953 if ( Normal.dot(v) <= 0 ) {
return 0.0; }
958 if ( Normal.dot(v) <= 0 ) {
return 0.0; }
972 zi = p.z() + sd*v.z() ;
974 if ((std::fabs(zi) <= tolODz) && (nt2 < 0))
981 xi = p.x() + sd*v.x() ;
982 yi = p.y() + sd*v.y() ;
983 ri = rMaxAv + zi*tanRMax ;
1008 nt1 = t1 - (tanRMin*v.z())*(tanRMin*v.z()) ;
1009 nt2 = t2 - tanRMin*v.z()*rin ;
1010 nt3 = t3 - rin*rin ;
1024 if(b>0){sd = c/( -b-std::sqrt(d));}
1025 else {sd = -b + std::sqrt(d) ;}
1031 G4double fTerm = sd-std::fmod(sd,dRmax);
1034 zi = p.z() + sd*v.z() ;
1036 if ( std::fabs(zi) <= tolODz )
1040 xi = p.x() + sd*v.x() ;
1041 yi = p.y() + sd*v.y() ;
1042 ri = rMinAv + zi*tanRMin ;
1047 if ( sd > halfRadTolerance ) { snxt=sd; }
1052 risec = std::sqrt(xi*xi + yi*yi)*secRMin ;
1053 Normal =
G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin);
1054 if ( Normal.dot(v) <= 0 ) { snxt = sd; }
1060 if ( sd > halfRadTolerance ) {
return sd; }
1065 xi = p.x() + sd*v.x() ;
1066 yi = p.y() + sd*v.y() ;
1067 risec = std::sqrt(xi*xi + yi*yi)*secRMin ;
1068 Normal =
G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin) ;
1069 if ( Normal.dot(v) <= 0 ) {
return sd; }
1089 if (b>0) { sd = c/(-b-std::sqrt(d)); }
1090 else { sd = -b + std::sqrt(d); }
1091 zi = p.z() + sd*v.z() ;
1092 ri = rMinAv + zi*tanRMin ;
1096 if ( (sd >= 0) && (std::fabs(zi) <= tolODz) )
1100 G4double fTerm = sd-std::fmod(sd,dRmax);
1105 xi = p.x() + sd*v.x() ;
1106 yi = p.y() + sd*v.y() ;
1111 if ( sd > halfRadTolerance ) { snxt=sd; }
1116 risec = std::sqrt(xi*xi + yi*yi)*secRMin ;
1117 Normal =
G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin);
1118 if ( Normal.dot(v) <= 0 ) { snxt = sd; }
1124 if( sd > halfRadTolerance ) {
return sd; }
1129 xi = p.x() + sd*v.x() ;
1130 yi = p.y() + sd*v.y() ;
1131 risec = std::sqrt(xi*xi + yi*yi)*secRMin ;
1132 Normal =
G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin) ;
1133 if ( Normal.dot(v) <= 0 ) {
return sd; }
1140 if (b>0) { sd = -b - std::sqrt(d); }
1141 else { sd = c/(-b+std::sqrt(d)); }
1142 zi = p.z() + sd*v.z() ;
1143 ri = rMinAv + zi*tanRMin ;
1145 if ( (sd >= 0) && (ri > 0) && (std::fabs(zi) <= tolODz) )
1149 G4double fTerm = sd-std::fmod(sd,dRmax);
1154 xi = p.x() + sd*v.x() ;
1155 yi = p.y() + sd*v.y() ;
1160 if ( sd > halfRadTolerance ) { snxt=sd; }
1165 risec = std::sqrt(xi*xi + yi*yi)*secRMin ;
1166 Normal =
G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin);
1167 if ( Normal.dot(v) <= 0 ) { snxt = sd; }
1173 if ( sd > halfRadTolerance ) {
return sd; }
1178 xi = p.x() + sd*v.x() ;
1179 yi = p.y() + sd*v.y() ;
1180 risec = std::sqrt(xi*xi + yi*yi)*secRMin ;
1181 Normal =
G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin) ;
1182 if ( Normal.dot(v) <= 0 ) {
return sd; }
1196 if ( std::fabs(p.z()) <= tolODz )
1208 else {
return 0.0; }
1221 if (b>0) { sd = -b - std::sqrt(d); }
1222 else { sd = c/(-b+std::sqrt(d)); }
1223 zi = p.z() + sd*v.z() ;
1224 ri = rMinAv + zi*tanRMin ;
1228 if (b>0) { sd = c/(-b-std::sqrt(d)); }
1229 else { sd = -b + std::sqrt(d); }
1231 zi = p.z() + sd*v.z() ;
1233 if ( (sd >= 0) && (std::fabs(zi) <= tolODz) )
1237 G4double fTerm = sd-std::fmod(sd,dRmax);
1242 xi = p.x() + sd*v.x() ;
1243 yi = p.y() + sd*v.y() ;
1244 ri = rMinAv + zi*tanRMin ;
1264 if (b>0) { sd = c/(-b-std::sqrt(d)); }
1265 else { sd = -b + std::sqrt(d) ; }
1266 zi = p.z() + sd*v.z() ;
1268 if ( (sd >= 0) && (std::fabs(zi) <= tolODz) )
1272 G4double fTerm = sd-std::fmod(sd,dRmax);
1277 xi = p.x() + sd*v.x();
1278 yi = p.y() + sd*v.y();
1279 ri = rMinAv + zi*tanRMin ;
1311 if (Dist < halfCarTolerance)
1317 if ( sd < 0 ) { sd = 0.0; }
1319 zi = p.z() + sd*v.z() ;
1321 if ( std::fabs(zi) <= tolODz )
1323 xi = p.x() + sd*v.x() ;
1324 yi = p.y() + sd*v.y() ;
1325 rhoi2 = xi*xi + yi*yi ;
1326 tolORMin2 = (rMinOAv + zi*tanRMin)*(rMinOAv + zi*tanRMin) ;
1327 tolORMax2 = (rMaxOAv + zi*tanRMax)*(rMaxOAv + zi*tanRMax) ;
1329 if ( (rhoi2 >= tolORMin2) && (rhoi2 <= tolORMax2) )
1348 if (Dist < halfCarTolerance)
1354 if ( sd < 0 ) { sd = 0.0; }
1356 zi = p.z() + sd*v.z() ;
1358 if (std::fabs(zi) <= tolODz)
1360 xi = p.x() + sd*v.x() ;
1361 yi = p.y() + sd*v.y() ;
1362 rhoi2 = xi*xi + yi*yi ;
1363 tolORMin2 = (rMinOAv + zi*tanRMin)*(rMinOAv + zi*tanRMin) ;
1364 tolORMax2 = (rMaxOAv + zi*tanRMax)*(rMaxOAv + zi*tanRMax) ;
1366 if ( (rhoi2 >= tolORMin2) && (rhoi2 <= tolORMax2) )
1378 if (snxt < halfCarTolerance) { snxt = 0.; }
1392 G4double safe=0.0, rho, safeR1, safeR2, safeZ, safePhi, cosPsi ;
1396 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
1397 safeZ = std::fabs(p.z()) -
fDz ;
1402 secRMin = std::sqrt(1.0 + tanRMin*tanRMin) ;
1404 safeR1 = (pRMin - rho)/secRMin ;
1407 secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
1409 safeR2 = (rho - pRMax)/secRMax ;
1411 if ( safeR1 > safeR2) { safe = safeR1; }
1412 else { safe = safeR2; }
1417 secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
1419 safe = (rho - pRMax)/secRMax ;
1421 if ( safeZ > safe ) { safe = safeZ; }
1429 if ( cosPsi < std::cos(
fDPhi*0.5) )
1433 safePhi = std::fabs(p.x()*std::sin(
fSPhi)-p.y()*std::cos(
fSPhi));
1439 if ( safePhi > safe ) { safe = safePhi; }
1442 if ( safe < 0.0 ) { safe = 0.0; }
1462 G4double tanRMax, secRMax, rMaxAv ;
1463 G4double tanRMin, secRMin, rMinAv ;
1465 G4double t1, t2, t3, rout, rin, nt1, nt2, nt3 ;
1475 G4double pDistS, compS, pDistE, compE, sphi2, xi, yi, risec, vphi ;
1482 pdist =
fDz - p.z() ;
1486 snxt = pdist/v.z() ;
1499 else if ( v.z() < 0.0 )
1501 pdist =
fDz + p.z() ;
1505 snxt = -pdist/v.z() ;
1542 secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
1546 t1 = 1.0 - v.z()*v.z() ;
1547 t2 = p.x()*v.x() + p.y()*v.y() ;
1548 t3 = p.x()*p.x() + p.y()*p.y() ;
1549 rout = tanRMax*p.z() + rMaxAv ;
1551 nt1 = t1 - (tanRMax*v.z())*(tanRMax*v.z()) ;
1552 nt2 = t2 - tanRMax*v.z()*rout ;
1553 nt3 = t3 - rout*rout ;
1557 deltaRoi2 = snxt*snxt*t1 + 2*snxt*t2 + t3
1560 else if ( v.z() < 0.0 )
1562 deltaRoi2 = snxt*snxt*t1 + 2*snxt*t2 + t3
1570 if ( nt1 && (deltaRoi2 > 0.0) )
1587 risec = std::sqrt(t3)*secRMax ;
1589 *n =
G4ThreeVector(p.x()/risec,p.y()/risec,-tanRMax/secRMax);
1596 if (b>0) { srd = -b - std::sqrt(d); }
1597 else { srd = c/(-b+std::sqrt(d)) ; }
1599 zi = p.z() + srd*v.z() ;
1600 ri = tanRMax*zi + rMaxAv ;
1615 if (b>0) { sr2 = c/(-b-std::sqrt(d)); }
1616 else { sr2 = -b + std::sqrt(d); }
1617 zi = p.z() + sr2*v.z() ;
1618 ri = tanRMax*zi + rMaxAv ;
1647 risec = std::sqrt(t3)*secRMax;
1649 *n =
G4ThreeVector(p.x()/risec,p.y()/risec,-tanRMax/secRMax);
1654 else if ( nt2 && (deltaRoi2 > 0.0) )
1660 risec = std::sqrt(t3)*secRMax;
1662 *n =
G4ThreeVector(p.x()/risec,p.y()/risec,-tanRMax/secRMax);
1687 xi = p.x() + slentol*v.x();
1688 yi = p.y() + slentol*v.y();
1689 risec = std::sqrt(xi*xi + yi*yi)*secRMax;
1692 if ( Normal.dot(v) > 0 )
1696 *n = Normal.unit() ;
1712 nt1 = t1 - (tanRMin*v.z())*(tanRMin*v.z()) ;
1716 secRMin = std::sqrt(1.0 + tanRMin*tanRMin) ;
1718 rin = tanRMin*p.z() + rMinAv ;
1719 nt2 = t2 - tanRMin*v.z()*rin ;
1720 nt3 = t3 - rin*rin ;
1737 if (calcNorm) { *validNorm =
false; }
1743 if (b>0) { sr2 = -b - std::sqrt(d); }
1744 else { sr2 = c/(-b+std::sqrt(d)); }
1745 zi = p.z() + sr2*v.z() ;
1746 ri = tanRMin*zi + rMinAv ;
1758 if (b>0) { sr3 = c/(-b-std::sqrt(d)); }
1759 else { sr3 = -b + std::sqrt(d) ; }
1768 zi = p.z() + sr3*v.z() ;
1769 ri = tanRMin*zi + rMinAv ;
1807 xi = p.x() + slentol*v.x() ;
1808 yi = p.y() + slentol*v.y() ;
1809 if( sidetol==
kRMax )
1811 risec = std::sqrt(xi*xi + yi*yi)*secRMax ;
1812 Normal =
G4ThreeVector(xi/risec,yi/risec,-tanRMax/secRMax) ;
1816 risec = std::sqrt(xi*xi + yi*yi)*secRMin ;
1817 Normal =
G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin) ;
1819 if( Normal.dot(v) > 0 )
1825 *n = Normal.unit() ;
1852 vphi = std::atan2(v.y(),v.x()) ;
1857 if ( p.x() || p.y() )
1879 sphi = pDistS/compS ;
1882 xi = p.x() + sphi*v.x() ;
1883 yi = p.y() + sphi*v.y() ;
1924 sphi2 = pDistE/compE ;
1930 xi = p.x() + sphi2*v.x() ;
1931 yi = p.y() + sphi2*v.y() ;
1945 else { sphi = 0.0; }
1955 else { sphi = 0.0; }
1997 xi = p.x() + snxt*v.x() ;
1998 yi = p.y() + snxt*v.y() ;
1999 risec = std::sqrt(xi*xi + yi*yi)*secRMax ;
2004 *validNorm = false ;
2014 *validNorm = false ;
2025 *validNorm = false ;
2039 std::ostringstream message;
2040 G4int oldprc = message.precision(16) ;
2041 message <<
"Undefined side for valid surface normal to solid."
2043 <<
"Position:" << G4endl << G4endl
2044 <<
"p.x() = " << p.x()/
mm <<
" mm" << G4endl
2045 <<
"p.y() = " << p.y()/
mm <<
" mm" << G4endl
2046 <<
"p.z() = " << p.z()/
mm <<
" mm" << G4endl << G4endl
2047 <<
"pho at z = " << std::sqrt( p.x()*p.x()+p.y()*p.y() )/
mm
2048 <<
" mm" << G4endl << G4endl ;
2049 if( p.x() != 0. || p.y() != 0.)
2051 message <<
"point phi = " << std::atan2(p.y(),p.x())/
degree
2052 <<
" degree" << G4endl << G4endl ;
2054 message <<
"Direction:" << G4endl << G4endl
2055 <<
"v.x() = " << v.x() << G4endl
2056 <<
"v.y() = " << v.y() << G4endl
2057 <<
"v.z() = " << v.z() << G4endl<< G4endl
2058 <<
"Proposed distance :" << G4endl<< G4endl
2059 <<
"snxt = " << snxt/
mm <<
" mm" <<
G4endl ;
2060 message.precision(oldprc) ;
2061 G4Exception(
"G4Cons::DistanceToOut(p,v,..)",
"GeomSolids1002",
2077 G4double safe=0.0, rho, safeR1, safeR2, safeZ, safePhi;
2091 G4cout <<
"pho at z = " << std::sqrt( p.x()*p.x()+p.y()*p.y() )/
mm
2092 <<
" mm" << G4endl << G4endl ;
2093 if( (p.x() != 0.) || (p.x() != 0.) )
2095 G4cout <<
"point phi = " << std::atan2(p.y(),p.x())/
degree
2096 <<
" degree" << G4endl << G4endl ;
2098 G4cout.precision(oldprc) ;
2099 G4Exception(
"G4Cons::DistanceToOut(p)",
"GeomSolids1002",
2104 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
2105 safeZ =
fDz - std::fabs(p.z()) ;
2110 secRMin = std::sqrt(1.0 + tanRMin*tanRMin) ;
2112 safeR1 = (rho - pRMin)/secRMin ;
2120 secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
2122 safeR2 = (pRMax - rho)/secRMax ;
2124 if (safeR1 < safeR2) { safe = safeR1; }
2125 else { safe = safeR2; }
2126 if (safeZ < safe) { safe = safeZ ; }
2142 if (safePhi < safe) { safe = safePhi; }
2144 if ( safe < 0 ) { safe = 0; }
2165 G4double meshAngle, meshRMax1, meshRMax2, crossAngle;
2166 G4double cosCrossAngle, sinCrossAngle, sAngle ;
2167 G4double rMaxX1, rMaxX2, rMaxY1, rMaxY2, rMinX1, rMinX2, rMinY1, rMinY2 ;
2168 G4int crossSection, noCrossSections ;
2182 meshAngle =
fDPhi/(noCrossSections - 1) ;
2184 meshRMax1 =
fRmax1/std::cos(meshAngle*0.5) ;
2185 meshRMax2 =
fRmax2/std::cos(meshAngle*0.5) ;
2192 sAngle = -meshAngle*0.5 ;
2202 vertices->reserve(noCrossSections*4) ;
2203 for (crossSection = 0 ; crossSection < noCrossSections ; crossSection++)
2207 crossAngle = sAngle + crossSection*meshAngle ;
2208 cosCrossAngle = std::cos(crossAngle) ;
2209 sinCrossAngle = std::sin(crossAngle) ;
2211 rMaxX1 = meshRMax1*cosCrossAngle ;
2212 rMaxY1 = meshRMax1*sinCrossAngle ;
2213 rMaxX2 = meshRMax2*cosCrossAngle ;
2214 rMaxY2 = meshRMax2*sinCrossAngle ;
2216 rMinX1 =
fRmin1*cosCrossAngle ;
2217 rMinY1 =
fRmin1*sinCrossAngle ;
2218 rMinX2 =
fRmin2*cosCrossAngle ;
2219 rMinY2 =
fRmin2*sinCrossAngle ;
2237 "Error in allocation of vertices. Out of memory !");
2258 return new G4Cons(*
this);
2267 G4int oldprc = os.precision(16);
2268 os <<
"-----------------------------------------------------------\n"
2269 <<
" *** Dump for solid - " <<
GetName() <<
" ***\n"
2270 <<
" ===================================================\n"
2271 <<
" Solid type: G4Cons\n"
2272 <<
" Parameters: \n"
2273 <<
" inside -fDz radius: " <<
fRmin1/
mm <<
" mm \n"
2274 <<
" outside -fDz radius: " <<
fRmax1/
mm <<
" mm \n"
2275 <<
" inside +fDz radius: " <<
fRmin2/
mm <<
" mm \n"
2276 <<
" outside +fDz radius: " <<
fRmax2/
mm <<
" mm \n"
2277 <<
" half length in Z : " <<
fDz/
mm <<
" mm \n"
2278 <<
" starting angle of segment: " <<
fSPhi/
degree <<
" degrees \n"
2279 <<
" delta angle of segment : " <<
fDPhi/
degree <<
" degrees \n"
2280 <<
"-----------------------------------------------------------\n";
2281 os.precision(oldprc);
2296 G4double Aone, Atwo, Athree, Afour, Afive, slin, slout, phi;
2297 G4double zRand, cosu, sinu, rRand1, rRand2, chose, rone, rtwo, qone, qtwo;
2312 cosu = std::cos(phi); sinu = std::sin(phi);
2319 if( (chose >= 0.) && (chose < Aone) )
2325 rtwo*sinu*(qtwo-zRand), zRand);
2333 else if( (chose >= Aone) && (chose <= Aone + Atwo) )
2339 rone*sinu*(qone-zRand), zRand);
2347 else if( (chose >= Aone + Atwo) && (chose < Aone + Atwo + Athree) )
2351 else if( (chose >= Aone + Atwo + Athree)
2352 && (chose < Aone + Atwo + Athree + Afour) )
2356 else if( (chose >= Aone + Atwo + Athree + Afour)
2357 && (chose < Aone + Atwo + Athree + Afour + Afive) )
2363 rRand1*std::sin(
fSPhi), zRand);
G4Cons(const G4String &pName, G4double pRmin1, G4double pRmax1, G4double pRmin2, G4double pRmax2, G4double pDz, G4double pSPhi, G4double pDPhi)
void ClipCrossSection(G4ThreeVectorList *pVertices, const G4int pSectionIndex, const G4VoxelLimits &pVoxelLimit, const EAxis pAxis, G4double &pMin, G4double &pMax) const
ThreeVector shoot(const G4int Ap, const G4int Af)
G4double halfRadTolerance
static const G4double kInfinity
G4double GetMinYExtent() const
CLHEP::Hep3Vector G4ThreeVector
G4ThreeVectorList * CreateRotatedVertices(const G4AffineTransform &pTransform) const
G4Cons & operator=(const G4Cons &rhs)
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
G4bool IsYLimited() const
G4double GetRadiusInRing(G4double rmin, G4double rmax) const
EInside Inside(const G4ThreeVector &p) const
G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=G4bool(false), G4bool *validNorm=0, G4ThreeVector *n=0) const
G4ThreeVector ApproxSurfaceNormal(const G4ThreeVector &p) const
std::ostream & StreamInfo(std::ostream &os) const
G4bool IsXLimited() const
const G4int kMinMeshSections
virtual void AddSolid(const G4Box &)=0
G4GeometryType GetEntityType() const
void ComputeDimensions(G4VPVParameterisation *p, const G4int n, const G4VPhysicalVolume *pRep)
void CheckPhiAngles(G4double sPhi, G4double dPhi)
G4double GetMaxXExtent() const
G4double GetMinZExtent() const
G4GLOB_DLL std::ostream G4cout
G4double GetRadialTolerance() const
G4double halfCarTolerance
static const double twopi
std::vector< G4ThreeVector > G4ThreeVectorList
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
G4double GetMaxZExtent() const
G4double halfAngTolerance
G4Polyhedron * CreatePolyhedron() const
static const double degree
void DescribeYourselfTo(G4VGraphicsScene &scene) const
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pmin, G4double &pmax) const
G4double GetMaxYExtent() const
G4ThreeVector GetPointOnSurface() const
G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const
G4double GetMaxExtent(const EAxis pAxis) const
G4CSGSolid & operator=(const G4CSGSolid &rhs)
G4bool IsZLimited() const
G4double GetAngularTolerance() const
void ClipBetweenSections(G4ThreeVectorList *pVertices, const G4int pSectionIndex, const G4VoxelLimits &pVoxelLimit, const EAxis pAxis, G4double &pMin, G4double &pMax) const
G4double GetMinExtent(const EAxis pAxis) const
const G4double kMeshAngleDefault
static G4GeometryTolerance * GetInstance()
const G4int kMaxMeshSections