53 #if !defined(G4GEOM_USE_UCONS)
69 using namespace CLHEP;
91 :
G4CSGSolid(pName), fRmin1(pRmin1), fRmin2(pRmin2),
92 fRmax1(pRmax1), fRmax2(pRmax2), fDz(pDz), fSPhi(0.), fDPhi(0.)
106 message <<
"Invalid Z half-length for Solid: " <<
GetName() <<
G4endl
114 if (((pRmin1>=pRmax1) || (pRmin2>=pRmax2) || (pRmin1<0)) && (pRmin2<0))
117 message <<
"Invalid values of radii for Solid: " <<
GetName() <<
G4endl
118 <<
" pRmin1 = " << pRmin1 <<
", pRmin2 = " << pRmin2
119 <<
", pRmax1 = " << pRmax1 <<
", pRmax2 = " << pRmax2;
137 :
G4CSGSolid(a), kRadTolerance(0.), kAngTolerance(0.),
138 fRmin1(0.), fRmin2(0.), fRmax1(0.), fRmax2(0.), fDz(0.),
139 fSPhi(0.), fDPhi(0.), sinCPhi(0.), cosCPhi(0.), cosHDPhiOT(0.),
140 cosHDPhiIT(0.), sinSPhi(0.), cosSPhi(0.), sinEPhi(0.), cosEPhi(0.),
141 fPhiFullCone(false), halfCarTolerance(0.), halfRadTolerance(0.),
159 :
G4CSGSolid(rhs), kRadTolerance(rhs.kRadTolerance),
160 kAngTolerance(rhs.kAngTolerance), fRmin1(rhs.fRmin1), fRmin2(rhs.fRmin2),
161 fRmax1(rhs.fRmax1), fRmax2(rhs.fRmax2), fDz(rhs.fDz), fSPhi(rhs.fSPhi),
162 fDPhi(rhs.fDPhi), sinCPhi(rhs.sinCPhi), cosCPhi(rhs.cosCPhi),
163 cosHDPhiOT(rhs.cosHDPhiOT), cosHDPhiIT(rhs.cosHDPhiIT),
164 sinSPhi(rhs.sinSPhi), cosSPhi(rhs.cosSPhi), sinEPhi(rhs.sinEPhi),
165 cosEPhi(rhs.cosEPhi), fPhiFullCone(rhs.fPhiFullCone),
166 halfCarTolerance(rhs.halfCarTolerance),
167 halfRadTolerance(rhs.halfRadTolerance),
168 halfAngTolerance(rhs.halfAngTolerance)
180 if (
this == &rhs) {
return *
this; }
211 G4double r2, rl, rh, pPhi, tolRMin, tolRMax;
218 r2 = p.
x()*p.
x() + p.
y()*p.
y() ;
225 if ( tolRMin < 0 ) { tolRMin = 0; }
228 if ( (r2<tolRMin*tolRMin) || (r2>tolRMax*tolRMax) ) {
return in =
kOutside; }
231 else { tolRMin = 0.0; }
236 if ( (r2 < tolRMin*tolRMin) || (r2 >= tolRMax*tolRMax) ) { in =
kSurface; }
240 pPhi = std::atan2(p.
y(),p.
x()) ;
290 pMin.
set(vmin.
x(),vmin.
y(),-dz);
291 pMax.
set(vmax.
x(),vmax.
y(), dz);
295 pMin.
set(-rmax,-rmax,-dz);
296 pMax.
set( rmax, rmax, dz);
301 if (pMin.
x() >= pMax.
x() || pMin.
y() >= pMax.
y() || pMin.
z() >= pMax.
z())
304 message <<
"Bad bounding box (min >= max) for solid: "
306 <<
"\npMin = " << pMin
307 <<
"\npMax = " << pMax;
308 G4Exception(
"G4Cons::BoundingLimits()",
"GeomMgt0001",
333 if (
true)
return bbox.
CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
337 return exist = (pMin < pMax) ?
true :
false;
350 const G4int NSTEPS = 24;
352 G4int ksteps = (dphi <= astep) ? 1 : (
G4int)((dphi-
deg)/astep) + 1;
355 G4double sinHalf = std::sin(0.5*ang);
356 G4double cosHalf = std::cos(0.5*ang);
357 G4double sinStep = 2.*sinHalf*cosHalf;
358 G4double cosStep = 1. - 2.*sinHalf*sinHalf;
364 if (rmin1 == 0 && rmin2 == 0 && dphi ==
twopi)
370 for (
G4int k=0; k<NSTEPS; ++k)
372 baseA[k].set(rext1*cosCur,rext1*sinCur,-dz);
373 baseB[k].set(rext2*cosCur,rext2*sinCur, dz);
376 sinCur = sinCur*cosStep + cosCur*sinStep;
377 cosCur = cosCur*cosStep - sinTmp*sinStep;
379 std::vector<const G4ThreeVectorList *> polygons(2);
380 polygons[0] = &baseA;
381 polygons[1] = &baseB;
391 G4double sinCur = sinStart*cosHalf + cosStart*sinHalf;
392 G4double cosCur = cosStart*cosHalf - sinStart*sinHalf;
396 for (
G4int k=0; k<ksteps+2; ++k) pols[k].resize(4);
397 pols[0][0].set(rmin2*cosStart,rmin2*sinStart, dz);
398 pols[0][1].set(rmin1*cosStart,rmin1*sinStart,-dz);
399 pols[0][2].set(rmax1*cosStart,rmax1*sinStart,-dz);
400 pols[0][3].set(rmax2*cosStart,rmax2*sinStart, dz);
401 for (
G4int k=1; k<ksteps+1; ++k)
403 pols[k][0].set(rmin2*cosCur,rmin2*sinCur, dz);
404 pols[k][1].set(rmin1*cosCur,rmin1*sinCur,-dz);
405 pols[k][2].set(rext1*cosCur,rext1*sinCur,-dz);
406 pols[k][3].set(rext2*cosCur,rext2*sinCur, dz);
409 sinCur = sinCur*cosStep + cosCur*sinStep;
410 cosCur = cosCur*cosStep - sinTmp*sinStep;
412 pols[ksteps+1][0].set(rmin2*cosEnd,rmin2*sinEnd, dz);
413 pols[ksteps+1][1].set(rmin1*cosEnd,rmin1*sinEnd,-dz);
414 pols[ksteps+1][2].set(rmax1*cosEnd,rmax1*sinEnd,-dz);
415 pols[ksteps+1][3].set(rmax2*cosEnd,rmax2*sinEnd, dz);
418 std::vector<const G4ThreeVectorList *> polygons;
419 polygons.resize(ksteps+2);
420 for (
G4int k=0; k<ksteps+2; ++k) polygons[k] = &pols[k];
435 G4int noSurfaces = 0;
439 G4double tanRMin, secRMin, pRMin, widRMin;
440 G4double tanRMax, secRMax, pRMax, widRMax;
445 distZ = std::fabs(std::fabs(p.
z()) -
fDz);
446 rho = std::sqrt(p.
x()*p.
x() + p.
y()*p.
y());
449 secRMin = std::sqrt(1 + tanRMin*tanRMin);
450 pRMin = rho - p.
z()*tanRMin;
452 distRMin = std::fabs(pRMin - widRMin)/secRMin;
455 secRMax = std::sqrt(1+tanRMax*tanRMax);
456 pRMax = rho - p.
z()*tanRMax;
458 distRMax = std::fabs(pRMax - widRMax)/secRMax;
464 pPhi = std::atan2(p.
y(),p.
x());
469 distSPhi = std::fabs( pPhi -
fSPhi );
482 nR =
G4ThreeVector(p.
x()/rho/secRMax, p.
y()/rho/secRMax, -tanRMax/secRMax);
485 nr =
G4ThreeVector(-p.
x()/rho/secRMin,-p.
y()/rho/secRMin,tanRMin/secRMin);
515 if ( p.
z() >= 0.) { sumnorm += nZ; }
516 else { sumnorm -= nZ; }
518 if ( noSurfaces == 0 )
521 G4Exception(
"G4Cons::SurfaceNormal(p)",
"GeomSolids1002",
526 else if ( noSurfaces == 1 ) { norm = sumnorm; }
527 else { norm = sumnorm.
unit(); }
542 G4double distZ, distRMin, distRMax, distSPhi, distEPhi, distMin ;
543 G4double tanRMin, secRMin, pRMin, widRMin ;
544 G4double tanRMax, secRMax, pRMax, widRMax ;
546 distZ = std::fabs(std::fabs(p.
z()) -
fDz) ;
547 rho = std::sqrt(p.
x()*p.
x() + p.
y()*p.
y()) ;
550 secRMin = std::sqrt(1 + tanRMin*tanRMin) ;
551 pRMin = rho - p.
z()*tanRMin ;
553 distRMin = std::fabs(pRMin - widRMin)/secRMin ;
556 secRMax = std::sqrt(1+tanRMax*tanRMax) ;
557 pRMax = rho - p.
z()*tanRMax ;
559 distRMax = std::fabs(pRMax - widRMax)/secRMax ;
561 if (distRMin < distRMax)
563 if (distZ < distRMin)
576 if (distZ < distRMax)
589 phi = std::atan2(p.
y(),p.
x()) ;
591 if (phi < 0) { phi +=
twopi; }
594 else { distSPhi = std::fabs(phi -
fSPhi)*rho; }
596 distEPhi = std::fabs(phi -
fSPhi -
fDPhi)*rho ;
600 if (distSPhi < distEPhi)
602 if (distSPhi < distMin) { side =
kNSPhi; }
606 if (distEPhi < distMin) { side =
kNEPhi; }
633 "Undefined side for valid surface normal to solid.");
669 G4double tanRMax,secRMax,rMaxAv,rMaxOAv ;
670 G4double tanRMin,secRMin,rMinAv,rMinOAv ;
673 G4double tolORMin,tolORMin2,tolIRMin,tolIRMin2 ;
674 G4double tolORMax2,tolIRMax,tolIRMax2 ;
677 G4double Dist,sd,xi,yi,zi,ri=0.,risec,rhoi2,cosPsi ;
688 secRMin = std::sqrt(1.0 + tanRMin*tanRMin) ;
700 secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
709 if (std::fabs(p.
z()) >= tolIDz)
711 if ( p.
z()*v.
z() < 0 )
713 sd = (std::fabs(p.
z()) -
fDz)/std::fabs(v.
z()) ;
715 if( sd < 0.0 ) { sd = 0.0; }
717 xi = p.
x() + sd*v.
x() ;
718 yi = p.
y() + sd*v.
y() ;
719 rhoi2 = xi*xi + yi*yi ;
726 tolORMin =
fRmin1 - halfRadTolerance*secRMin ;
727 tolIRMin =
fRmin1 + halfRadTolerance*secRMin ;
728 tolIRMax =
fRmax1 - halfRadTolerance*secRMin ;
734 tolORMin =
fRmin2 - halfRadTolerance*secRMin ;
735 tolIRMin =
fRmin2 + halfRadTolerance*secRMin ;
736 tolIRMax =
fRmax2 - halfRadTolerance*secRMin ;
743 tolIRMin2 = tolIRMin*tolIRMin ;
750 if ( tolIRMax > 0 ) { tolIRMax2 = tolIRMax*tolIRMax; }
751 else { tolIRMax2 = 0.0; }
753 if ( (tolIRMin2 <= rhoi2) && (rhoi2 <= tolIRMax2) )
794 t1 = 1.0 - v.
z()*v.
z() ;
795 t2 = p.
x()*v.
x() + p.
y()*v.
y() ;
796 t3 = p.
x()*p.
x() + p.
y()*p.
y() ;
797 rin = tanRMin*p.
z() + rMinAv ;
798 rout = tanRMax*p.
z() + rMaxAv ;
803 nt1 = t1 - (tanRMax*v.
z())*(tanRMax*v.
z()) ;
804 nt2 = t2 - tanRMax*v.
z()*rout ;
805 nt3 = t3 - rout*rout ;
821 if ((rout < 0) && (nt3 <= 0))
826 if (b>0) { sd = c/(-b-std::sqrt(d)); }
827 else { sd = -b + std::sqrt(d); }
831 if ((b <= 0) && (c >= 0))
833 sd=c/(-b+std::sqrt(d));
839 sd = -b + std::sqrt(d) ;
840 if((sd<0) & (sd>-halfRadTolerance)) sd=0;
852 G4double fTerm = sd-std::fmod(sd,dRmax);
855 zi = p.
z() + sd*v.
z() ;
857 if (std::fabs(zi) <= tolODz)
864 xi = p.
x() + sd*v.
x() ;
865 yi = p.
y() + sd*v.
y() ;
866 ri = rMaxAv + zi*tanRMax ;
880 if ( ( t3 > (rin + halfRadTolerance*secRMin)*
881 (rin + halfRadTolerance*secRMin) )
882 && (nt2 < 0) && (d >= 0) && (std::fabs(p.
z()) <= tolIDz) )
889 risec = std::sqrt(xi*xi + yi*yi)*secRMax ;
896 if ( Normal.
dot(v) <= 0 ) {
return 0.0; }
901 if ( Normal.
dot(v) <= 0 ) {
return 0.0; }
915 zi = p.
z() + sd*v.
z() ;
917 if ((std::fabs(zi) <= tolODz) && (nt2 < 0))
924 xi = p.
x() + sd*v.
x() ;
925 yi = p.
y() + sd*v.
y() ;
926 ri = rMaxAv + zi*tanRMax ;
951 nt1 = t1 - (tanRMin*v.
z())*(tanRMin*v.
z()) ;
952 nt2 = t2 - tanRMin*v.
z()*rin ;
967 if(b>0){sd = c/( -b-std::sqrt(d));}
968 else {sd = -b + std::sqrt(d) ;}
974 G4double fTerm = sd-std::fmod(sd,dRmax);
977 zi = p.
z() + sd*v.
z() ;
979 if ( std::fabs(zi) <= tolODz )
983 xi = p.
x() + sd*v.
x() ;
984 yi = p.
y() + sd*v.
y() ;
985 ri = rMinAv + zi*tanRMin ;
990 if ( sd > halfRadTolerance ) { snxt=sd; }
995 risec = std::sqrt(xi*xi + yi*yi)*secRMin ;
997 if ( Normal.
dot(v) <= 0 ) { snxt = sd; }
1003 if ( sd > halfRadTolerance ) {
return sd; }
1008 xi = p.
x() + sd*v.
x() ;
1009 yi = p.
y() + sd*v.
y() ;
1010 risec = std::sqrt(xi*xi + yi*yi)*secRMin ;
1011 Normal =
G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin) ;
1012 if ( Normal.
dot(v) <= 0 ) {
return sd; }
1032 if (b>0) { sd = c/(-b-std::sqrt(d)); }
1033 else { sd = -b + std::sqrt(d); }
1034 zi = p.
z() + sd*v.
z() ;
1035 ri = rMinAv + zi*tanRMin ;
1039 if ( (sd >= 0) && (std::fabs(zi) <= tolODz) )
1043 G4double fTerm = sd-std::fmod(sd,dRmax);
1048 xi = p.
x() + sd*v.
x() ;
1049 yi = p.
y() + sd*v.
y() ;
1054 if ( sd > halfRadTolerance ) { snxt=sd; }
1059 risec = std::sqrt(xi*xi + yi*yi)*secRMin ;
1060 Normal =
G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin);
1061 if ( Normal.
dot(v) <= 0 ) { snxt = sd; }
1067 if( sd > halfRadTolerance ) {
return sd; }
1072 xi = p.
x() + sd*v.
x() ;
1073 yi = p.
y() + sd*v.
y() ;
1074 risec = std::sqrt(xi*xi + yi*yi)*secRMin ;
1075 Normal =
G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin) ;
1076 if ( Normal.
dot(v) <= 0 ) {
return sd; }
1083 if (b>0) { sd = -b - std::sqrt(d); }
1084 else { sd = c/(-b+std::sqrt(d)); }
1085 zi = p.
z() + sd*v.
z() ;
1086 ri = rMinAv + zi*tanRMin ;
1088 if ( (sd >= 0) && (ri > 0) && (std::fabs(zi) <= tolODz) )
1092 G4double fTerm = sd-std::fmod(sd,dRmax);
1097 xi = p.
x() + sd*v.
x() ;
1098 yi = p.
y() + sd*v.
y() ;
1103 if ( sd > halfRadTolerance ) { snxt=sd; }
1108 risec = std::sqrt(xi*xi + yi*yi)*secRMin ;
1109 Normal =
G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin);
1110 if ( Normal.
dot(v) <= 0 ) { snxt = sd; }
1116 if ( sd > halfRadTolerance ) {
return sd; }
1121 xi = p.
x() + sd*v.
x() ;
1122 yi = p.
y() + sd*v.
y() ;
1123 risec = std::sqrt(xi*xi + yi*yi)*secRMin ;
1124 Normal =
G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin) ;
1125 if ( Normal.
dot(v) <= 0 ) {
return sd; }
1139 if ( std::fabs(p.
z()) <= tolODz )
1151 else {
return 0.0; }
1164 if (b>0) { sd = -b - std::sqrt(d); }
1165 else { sd = c/(-b+std::sqrt(d)); }
1166 zi = p.
z() + sd*v.
z() ;
1167 ri = rMinAv + zi*tanRMin ;
1171 if (b>0) { sd = c/(-b-std::sqrt(d)); }
1172 else { sd = -b + std::sqrt(d); }
1174 zi = p.
z() + sd*v.
z() ;
1176 if ( (sd >= 0) && (std::fabs(zi) <= tolODz) )
1180 G4double fTerm = sd-std::fmod(sd,dRmax);
1185 xi = p.
x() + sd*v.
x() ;
1186 yi = p.
y() + sd*v.
y() ;
1187 ri = rMinAv + zi*tanRMin ;
1207 if (b>0) { sd = c/(-b-std::sqrt(d)); }
1208 else { sd = -b + std::sqrt(d) ; }
1209 zi = p.
z() + sd*v.
z() ;
1211 if ( (sd >= 0) && (std::fabs(zi) <= tolODz) )
1215 G4double fTerm = sd-std::fmod(sd,dRmax);
1220 xi = p.
x() + sd*v.
x();
1221 yi = p.
y() + sd*v.
y();
1222 ri = rMinAv + zi*tanRMin ;
1254 if (Dist < halfCarTolerance)
1260 if ( sd < 0 ) { sd = 0.0; }
1262 zi = p.
z() + sd*v.
z() ;
1264 if ( std::fabs(zi) <= tolODz )
1266 xi = p.
x() + sd*v.
x() ;
1267 yi = p.
y() + sd*v.
y() ;
1268 rhoi2 = xi*xi + yi*yi ;
1269 tolORMin2 = (rMinOAv + zi*tanRMin)*(rMinOAv + zi*tanRMin) ;
1270 tolORMax2 = (rMaxOAv + zi*tanRMax)*(rMaxOAv + zi*tanRMax) ;
1272 if ( (rhoi2 >= tolORMin2) && (rhoi2 <= tolORMax2) )
1291 if (Dist < halfCarTolerance)
1297 if ( sd < 0 ) { sd = 0.0; }
1299 zi = p.
z() + sd*v.
z() ;
1301 if (std::fabs(zi) <= tolODz)
1303 xi = p.
x() + sd*v.
x() ;
1304 yi = p.
y() + sd*v.
y() ;
1305 rhoi2 = xi*xi + yi*yi ;
1306 tolORMin2 = (rMinOAv + zi*tanRMin)*(rMinOAv + zi*tanRMin) ;
1307 tolORMax2 = (rMaxOAv + zi*tanRMax)*(rMaxOAv + zi*tanRMax) ;
1309 if ( (rhoi2 >= tolORMin2) && (rhoi2 <= tolORMax2) )
1321 if (snxt < halfCarTolerance) { snxt = 0.; }
1335 G4double safe=0.0, rho, safeR1, safeR2, safeZ, safePhi, cosPsi ;
1339 rho = std::sqrt(p.
x()*p.
x() + p.
y()*p.
y()) ;
1340 safeZ = std::fabs(p.
z()) -
fDz ;
1345 secRMin = std::sqrt(1.0 + tanRMin*tanRMin) ;
1347 safeR1 = (pRMin - rho)/secRMin ;
1350 secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
1352 safeR2 = (rho - pRMax)/secRMax ;
1354 if ( safeR1 > safeR2) { safe = safeR1; }
1355 else { safe = safeR2; }
1360 secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
1362 safe = (rho - pRMax)/secRMax ;
1364 if ( safeZ > safe ) { safe = safeZ; }
1372 if ( cosPsi < std::cos(
fDPhi*0.5) )
1376 safePhi = std::fabs(p.
x()*std::sin(
fSPhi)-p.
y()*std::cos(
fSPhi));
1382 if ( safePhi > safe ) { safe = safePhi; }
1385 if ( safe < 0.0 ) { safe = 0.0; }
1405 G4double tanRMax, secRMax, rMaxAv ;
1406 G4double tanRMin, secRMin, rMinAv ;
1418 G4double pDistS, compS, pDistE, compE, sphi2, xi, yi, risec, vphi ;
1425 pdist =
fDz - p.
z() ;
1429 snxt = pdist/v.
z() ;
1442 else if ( v.
z() < 0.0 )
1444 pdist =
fDz + p.
z() ;
1448 snxt = -pdist/v.
z() ;
1485 secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
1489 t1 = 1.0 - v.
z()*v.
z() ;
1490 t2 = p.
x()*v.
x() + p.
y()*v.
y() ;
1491 t3 = p.
x()*p.
x() + p.
y()*p.
y() ;
1492 rout = tanRMax*p.
z() + rMaxAv ;
1494 nt1 = t1 - (tanRMax*v.
z())*(tanRMax*v.
z()) ;
1495 nt2 = t2 - tanRMax*v.
z()*rout ;
1496 nt3 = t3 - rout*rout ;
1500 deltaRoi2 = snxt*snxt*t1 + 2*snxt*t2 + t3
1503 else if ( v.
z() < 0.0 )
1505 deltaRoi2 = snxt*snxt*t1 + 2*snxt*t2 + t3
1513 if ( nt1 && (deltaRoi2 > 0.0) )
1530 risec = std::sqrt(t3)*secRMax ;
1539 if (b>0) { srd = -b - std::sqrt(d); }
1540 else { srd = c/(-b+std::sqrt(d)) ; }
1542 zi = p.
z() + srd*v.
z() ;
1543 ri = tanRMax*zi + rMaxAv ;
1558 if (b>0) { sr2 = c/(-b-std::sqrt(d)); }
1559 else { sr2 = -b + std::sqrt(d); }
1560 zi = p.
z() + sr2*v.
z() ;
1561 ri = tanRMax*zi + rMaxAv ;
1590 risec = std::sqrt(t3)*secRMax;
1597 else if ( nt2 && (deltaRoi2 > 0.0) )
1603 risec = std::sqrt(t3)*secRMax;
1630 xi = p.
x() + slentol*v.
x();
1631 yi = p.
y() + slentol*v.
y();
1632 risec = std::sqrt(xi*xi + yi*yi)*secRMax;
1635 if ( Normal.
dot(v) > 0 )
1639 *n = Normal.
unit() ;
1655 nt1 = t1 - (tanRMin*v.
z())*(tanRMin*v.
z()) ;
1659 secRMin = std::sqrt(1.0 + tanRMin*tanRMin) ;
1661 rin = tanRMin*p.
z() + rMinAv ;
1662 nt2 = t2 - tanRMin*v.
z()*rin ;
1663 nt3 = t3 - rin*rin ;
1680 if (calcNorm) { *validNorm =
false; }
1686 if (b>0) { sr2 = -b - std::sqrt(d); }
1687 else { sr2 = c/(-b+std::sqrt(d)); }
1688 zi = p.
z() + sr2*v.
z() ;
1689 ri = tanRMin*zi + rMinAv ;
1701 if (b>0) { sr3 = c/(-b-std::sqrt(d)); }
1702 else { sr3 = -b + std::sqrt(d) ; }
1711 zi = p.
z() + sr3*v.
z() ;
1712 ri = tanRMin*zi + rMinAv ;
1750 xi = p.
x() + slentol*v.
x() ;
1751 yi = p.
y() + slentol*v.
y() ;
1752 if( sidetol==
kRMax )
1754 risec = std::sqrt(xi*xi + yi*yi)*secRMax ;
1755 Normal =
G4ThreeVector(xi/risec,yi/risec,-tanRMax/secRMax) ;
1759 risec = std::sqrt(xi*xi + yi*yi)*secRMin ;
1760 Normal =
G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin) ;
1762 if( Normal.
dot(v) > 0 )
1768 *n = Normal.
unit() ;
1795 vphi = std::atan2(v.
y(),v.
x()) ;
1800 if ( p.
x() || p.
y() )
1822 sphi = pDistS/compS ;
1825 xi = p.
x() + sphi*v.
x() ;
1826 yi = p.
y() + sphi*v.
y() ;
1867 sphi2 = pDistE/compE ;
1873 xi = p.
x() + sphi2*v.
x() ;
1874 yi = p.
y() + sphi2*v.
y() ;
1888 else { sphi = 0.0; }
1898 else { sphi = 0.0; }
1940 xi = p.
x() + snxt*v.
x() ;
1941 yi = p.
y() + snxt*v.
y() ;
1942 risec = std::sqrt(xi*xi + yi*yi)*secRMax ;
1947 *validNorm = false ;
1957 *validNorm = false ;
1968 *validNorm = false ;
1983 G4int oldprc = message.precision(16) ;
1984 message <<
"Undefined side for valid surface normal to solid."
1986 <<
"Position:" << G4endl << G4endl
1987 <<
"p.x() = " << p.
x()/
mm <<
" mm" << G4endl
1988 <<
"p.y() = " << p.
y()/
mm <<
" mm" << G4endl
1989 <<
"p.z() = " << p.
z()/
mm <<
" mm" << G4endl << G4endl
1990 <<
"pho at z = " << std::sqrt( p.
x()*p.
x()+p.
y()*p.
y() )/
mm
1991 <<
" mm" << G4endl << G4endl ;
1992 if( p.
x() != 0. || p.
y() != 0.)
1994 message <<
"point phi = " << std::atan2(p.
y(),p.
x())/
degree
1995 <<
" degree" << G4endl << G4endl ;
1997 message <<
"Direction:" << G4endl << G4endl
1998 <<
"v.x() = " << v.
x() << G4endl
1999 <<
"v.y() = " << v.
y() << G4endl
2000 <<
"v.z() = " << v.
z() << G4endl<< G4endl
2001 <<
"Proposed distance :" << G4endl<< G4endl
2002 <<
"snxt = " << snxt/
mm <<
" mm" <<
G4endl ;
2003 message.precision(oldprc) ;
2004 G4Exception(
"G4Cons::DistanceToOut(p,v,..)",
"GeomSolids1002",
2020 G4double safe=0.0, rho, safeR1, safeR2, safeZ, safePhi;
2034 G4cout <<
"pho at z = " << std::sqrt( p.
x()*p.
x()+p.
y()*p.
y() )/
mm
2035 <<
" mm" << G4endl << G4endl ;
2036 if( (p.
x() != 0.) || (p.
x() != 0.) )
2039 <<
" degree" << G4endl << G4endl ;
2041 G4cout.precision(oldprc) ;
2042 G4Exception(
"G4Cons::DistanceToOut(p)",
"GeomSolids1002",
2047 rho = std::sqrt(p.
x()*p.
x() + p.
y()*p.
y()) ;
2048 safeZ =
fDz - std::fabs(p.
z()) ;
2053 secRMin = std::sqrt(1.0 + tanRMin*tanRMin) ;
2055 safeR1 = (rho - pRMin)/secRMin ;
2063 secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
2065 safeR2 = (pRMax - rho)/secRMax ;
2067 if (safeR1 < safeR2) { safe = safeR1; }
2068 else { safe = safeR2; }
2069 if (safeZ < safe) { safe = safeZ ; }
2085 if (safePhi < safe) { safe = safePhi; }
2087 if ( safe < 0 ) { safe = 0; }
2107 return new G4Cons(*
this);
2116 G4int oldprc = os.precision(16);
2117 os <<
"-----------------------------------------------------------\n"
2118 <<
" *** Dump for solid - " <<
GetName() <<
" ***\n"
2119 <<
" ===================================================\n"
2120 <<
" Solid type: G4Cons\n"
2121 <<
" Parameters: \n"
2122 <<
" inside -fDz radius: " <<
fRmin1/
mm <<
" mm \n"
2123 <<
" outside -fDz radius: " <<
fRmax1/
mm <<
" mm \n"
2124 <<
" inside +fDz radius: " <<
fRmin2/
mm <<
" mm \n"
2125 <<
" outside +fDz radius: " <<
fRmax2/
mm <<
" mm \n"
2126 <<
" half length in Z : " <<
fDz/
mm <<
" mm \n"
2127 <<
" starting angle of segment: " <<
fSPhi/
degree <<
" degrees \n"
2128 <<
" delta angle of segment : " <<
fDPhi/
degree <<
" degrees \n"
2129 <<
"-----------------------------------------------------------\n";
2130 os.precision(oldprc);
2145 G4double Aone, Atwo, Athree, Afour, Afive, slin, slout, phi;
2146 G4double zRand, cosu, sinu, rRand1, rRand2, chose, rone, rtwo, qone, qtwo;
2161 cosu = std::cos(phi); sinu = std::sin(phi);
2168 if( (chose >= 0.) && (chose < Aone) )
2174 rtwo*sinu*(qtwo-zRand), zRand);
2182 else if( (chose >= Aone) && (chose <= Aone + Atwo) )
2188 rone*sinu*(qone-zRand), zRand);
2196 else if( (chose >= Aone + Atwo) && (chose < Aone + Atwo + Athree) )
2200 else if( (chose >= Aone + Atwo + Athree)
2201 && (chose < Aone + Atwo + Athree + Afour) )
2205 else if( (chose >= Aone + Atwo + Athree + Afour)
2206 && (chose < Aone + Atwo + Athree + Afour + Afive) )
2212 rRand1*std::sin(
fSPhi), zRand);
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimits, const G4Transform3D &pTransform3D, G4double &pMin, G4double &pMax) const
void set(double x, double y, double z)
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=G4bool(false), G4bool *validNorm=0, G4ThreeVector *n=0) const
G4double halfRadTolerance
G4double GetInnerRadiusMinusZ() const
CLHEP::Hep3Vector G4ThreeVector
std::vector< ExP01TrackerHit * > a
virtual void AddSolid(const G4Box &)=0
static const G4double kInfinity
static constexpr double mm
G4CSGSolid & operator=(const G4CSGSolid &rhs)
virtual void ComputeDimensions(G4Box &, const G4int, const G4VPhysicalVolume *) const
G4Cons & operator=(const G4Cons &rhs)
void message(RunManager *runmanager)
G4double GetRadiusInRing(G4double rmin, G4double rmax) const
G4ThreeVector ApproxSurfaceNormal(const G4ThreeVector &p) const
G4double GetAngularTolerance() const
double dot(const Hep3Vector &) const
G4double GetRadialTolerance() const
G4Cons(const G4String &pName, G4double pRmin1, G4double pRmax1, G4double pRmin2, G4double pRmax2, G4double pDz, G4double pSPhi, G4double pDPhi)
G4double GetCosStartPhi() const
G4double GetOuterRadiusMinusZ() const
G4double GetOuterRadiusPlusZ() const
G4double GetZHalfLength() const
G4double GetCosEndPhi() const
G4double halfCarTolerance
std::ostream & StreamInfo(std::ostream &os) const
G4bool BoundingBoxVsVoxelLimits(const EAxis pAxis, const G4VoxelLimits &pVoxelLimits, const G4Transform3D &pTransform3D, G4double &pMin, G4double &pMax) const
EInside Inside(const G4ThreeVector &p) const
std::vector< G4ThreeVector > G4ThreeVectorList
static constexpr double deg
static constexpr double twopi
G4Polyhedron * CreatePolyhedron() const
G4double GetInnerRadiusPlusZ() const
ThreeVector shoot(const G4int Ap, const G4int Af)
G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
void BoundingLimits(G4ThreeVector &pMin, G4ThreeVector &pMax) const
G4double halfAngTolerance
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
G4double GetDeltaPhiAngle() const
void DescribeYourselfTo(G4VGraphicsScene &scene) const
static constexpr double degree
G4GLOB_DLL std::ostream G4cout
void CheckPhiAngles(G4double sPhi, G4double dPhi)
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pMin, G4double &pMax) const
static G4GeometryTolerance * GetInstance()
G4ThreeVector GetPointOnSurface() const
static constexpr double pi
G4double GetSinEndPhi() const
G4double GetSinStartPhi() const
void ComputeDimensions(G4VPVParameterisation *p, const G4int n, const G4VPhysicalVolume *pRep)
G4GeometryType GetEntityType() const
T min(const T t1, const T t2)
brief Return the smallest of the two arguments