876 G4double rho2, rad2, pDotV2d, pDotV3d, pTheta ;
877 G4double tolSTheta=0., tolETheta=0. ;
880 const G4double halfRmaxTolerance = fRmaxTolerance*0.5;
881 const G4double halfRminTolerance = fRminTolerance*0.5;
882 const G4double tolORMin2 = (fRmin>halfRminTolerance)
883 ? (fRmin-halfRminTolerance)*(fRmin-halfRminTolerance) : 0;
885 (fRmin+halfRminTolerance)*(fRmin+halfRminTolerance);
887 (fRmax+halfRmaxTolerance)*(fRmax+halfRmaxTolerance);
889 (fRmax-halfRmaxTolerance)*(fRmax-halfRmaxTolerance);
893 G4double xi, yi, zi, rhoi, rhoi2, radi2, iTheta ;
910 rho2 = p.
x()*p.
x() + p.
y()*p.
y() ;
911 rad2 = rho2 + p.
z()*p.
z() ;
912 pTheta = std::atan2(std::sqrt(rho2),p.
z()) ;
914 pDotV2d = p.
x()*v.
x() + p.
y()*v.
y() ;
915 pDotV3d = pDotV2d + p.
z()*v.
z() ;
919 if (!fFullThetaSphere)
921 tolSTheta = fSTheta - halfAngTolerance ;
922 tolETheta = eTheta + halfAngTolerance ;
939 c = rad2 - fRmax*fRmax ;
941 if (c > fRmaxTolerance*fRmax)
946 d2 = pDotV3d*pDotV3d -
c ;
950 sd = -pDotV3d - std::sqrt(d2) ;
956 G4double fTerm = sd-std::fmod(sd,dRmax);
959 xi = p.
x() + sd*v.
x() ;
960 yi = p.
y() + sd*v.
y() ;
961 rhoi = std::sqrt(xi*xi + yi*yi) ;
963 if (!fFullPhiSphere && rhoi)
965 cosPsi = (xi*cosCPhi + yi*sinCPhi)/rhoi ;
967 if (cosPsi >= cosHDPhiOT)
969 if (!fFullThetaSphere)
971 zi = p.
z() + sd*v.
z() ;
976 iTheta = std::atan2(rhoi,zi) ;
977 if ( (iTheta >= tolSTheta) && (iTheta <= tolETheta) )
990 if (!fFullThetaSphere)
992 zi = p.
z() + sd*v.
z() ;
997 iTheta = std::atan2(rhoi,zi) ;
998 if ( (iTheta >= tolSTheta) && (iTheta <= tolETheta) )
1012 return snxt=kInfinity;
1020 d2 = pDotV3d*pDotV3d -
c ;
1022 if ( (rad2 > tolIRMax2)
1023 && ( (d2 >= fRmaxTolerance*fRmax) && (pDotV3d < 0) ) )
1025 if (!fFullPhiSphere)
1030 cosPsi = (p.
x()*cosCPhi + p.
y()*sinCPhi)/std::sqrt(rho2) ;
1032 if (cosPsi>=cosHDPhiIT)
1036 if ( !fFullThetaSphere )
1038 if ( (pTheta >= tolSTheta + kAngTolerance)
1039 && (pTheta <= tolETheta - kAngTolerance) )
1052 if ( !fFullThetaSphere )
1054 if ( (pTheta >= tolSTheta + kAngTolerance)
1055 && (pTheta <= tolETheta - kAngTolerance) )
1075 c = rad2 - fRmin*fRmin ;
1076 d2 = pDotV3d*pDotV3d -
c ;
1081 if ( (c > -halfRminTolerance) && (rad2 < tolIRMin2)
1084 if ( !fFullPhiSphere )
1089 cosPsi = (p.
x()*cosCPhi+p.
y()*sinCPhi)/std::sqrt(rho2) ;
1090 if (cosPsi >= cosHDPhiIT)
1094 if ( !fFullThetaSphere )
1096 if ( (pTheta >= tolSTheta + kAngTolerance)
1097 && (pTheta <= tolETheta - kAngTolerance) )
1110 if ( !fFullThetaSphere )
1112 if ( (pTheta >= tolSTheta + kAngTolerance)
1113 && (pTheta <= tolETheta - kAngTolerance) )
1128 sd = -pDotV3d + std::sqrt(d2) ;
1129 if ( sd >= halfRminTolerance )
1131 xi = p.
x() + sd*v.
x() ;
1132 yi = p.
y() + sd*v.
y() ;
1133 rhoi = std::sqrt(xi*xi+yi*yi) ;
1135 if ( !fFullPhiSphere && rhoi )
1137 cosPsi = (xi*cosCPhi + yi*sinCPhi)/rhoi ;
1139 if (cosPsi >= cosHDPhiOT)
1141 if ( !fFullThetaSphere )
1143 zi = p.
z() + sd*v.
z() ;
1148 iTheta = std::atan2(rhoi,zi) ;
1149 if ( (iTheta >= tolSTheta) && (iTheta<=tolETheta) )
1162 if ( !fFullThetaSphere )
1164 zi = p.
z() + sd*v.
z() ;
1169 iTheta = std::atan2(rhoi,zi) ;
1170 if ( (iTheta >= tolSTheta) && (iTheta <= tolETheta) )
1194 if ( !fFullPhiSphere )
1199 Comp = v.
x()*sinSPhi - v.
y()*cosSPhi ;
1203 Dist = p.
y()*cosSPhi - p.
x()*sinSPhi ;
1205 if (Dist < halfCarTolerance)
1213 xi = p.
x() + sd*v.
x() ;
1214 yi = p.
y() + sd*v.
y() ;
1215 zi = p.
z() + sd*v.
z() ;
1216 rhoi2 = xi*xi + yi*yi ;
1217 radi2 = rhoi2 + zi*zi ;
1228 if ( (radi2 <= tolORMax2)
1229 && (radi2 >= tolORMin2)
1230 && ((yi*cosCPhi-xi*sinCPhi) <= 0) )
1236 if ( !fFullThetaSphere )
1238 iTheta = std::atan2(std::sqrt(rhoi2),zi) ;
1239 if ( (iTheta >= tolSTheta) && (iTheta <= tolETheta) )
1244 if ((yi*cosCPhi-xi*sinCPhi) <= 0)
1262 Comp = -( v.
x()*sinEPhi-v.
y()*cosEPhi ) ;
1266 Dist = -(p.
y()*cosEPhi-p.
x()*sinEPhi) ;
1267 if ( Dist < halfCarTolerance )
1275 xi = p.
x() + sd*v.
x() ;
1276 yi = p.
y() + sd*v.
y() ;
1277 zi = p.
z() + sd*v.
z() ;
1278 rhoi2 = xi*xi + yi*yi ;
1279 radi2 = rhoi2 + zi*zi ;
1290 if ( (radi2 <= tolORMax2)
1291 && (radi2 >= tolORMin2)
1292 && ((yi*cosCPhi-xi*sinCPhi) >= 0) )
1298 if ( !fFullThetaSphere )
1300 iTheta = std::atan2(std::sqrt(rhoi2),zi) ;
1301 if ( (iTheta >= tolSTheta) && (iTheta <= tolETheta) )
1306 if ((yi*cosCPhi-xi*sinCPhi) >= 0)
1324 if ( !fFullThetaSphere )
1348 dist2STheta = rho2 - p.
z()*p.
z()*tanSTheta2 ;
1352 dist2STheta = kInfinity ;
1356 dist2ETheta=rho2-p.
z()*p.
z()*tanETheta2;
1360 dist2ETheta=kInfinity;
1362 if ( pTheta < tolSTheta )
1367 t1 = 1 - v.
z()*v.
z()*(1 + tanSTheta2) ;
1368 t2 = pDotV2d - p.
z()*v.
z()*tanSTheta2 ;
1372 c = dist2STheta/
t1 ;
1379 zi = p.
z() + sd*v.
z();
1381 if ( (sd < 0) || (zi*(fSTheta -
halfpi) > 0) )
1385 if ((sd >= 0) && (sd < snxt))
1387 xi = p.
x() + sd*v.
x();
1388 yi = p.
y() + sd*v.
y();
1389 zi = p.
z() + sd*v.
z();
1390 rhoi2 = xi*xi + yi*yi;
1391 radi2 = rhoi2 + zi*zi;
1392 if ( (radi2 <= tolORMax2)
1393 && (radi2 >= tolORMin2)
1394 && (zi*(fSTheta -
halfpi) <= 0) )
1396 if ( !fFullPhiSphere && rhoi2 )
1398 cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
1399 if (cosPsi >= cosHDPhiOT)
1418 t1 = 1 - v.
z()*v.
z()*(1 + tanETheta2) ;
1419 t2 = pDotV2d - p.
z()*v.
z()*tanETheta2 ;
1423 c = dist2ETheta/
t1 ;
1431 if ( (sd >= 0) && (sd < snxt) )
1433 xi = p.
x() + sd*v.
x() ;
1434 yi = p.
y() + sd*v.
y() ;
1435 zi = p.
z() + sd*v.
z() ;
1436 rhoi2 = xi*xi + yi*yi ;
1437 radi2 = rhoi2 + zi*zi ;
1439 if ( (radi2 <= tolORMax2)
1440 && (radi2 >= tolORMin2)
1441 && (zi*(eTheta -
halfpi) <= 0) )
1443 if (!fFullPhiSphere && rhoi2)
1445 cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
1446 if (cosPsi >= cosHDPhiOT)
1461 else if ( pTheta > tolETheta )
1467 t1 = 1 - v.
z()*v.
z()*(1 + tanETheta2) ;
1468 t2 = pDotV2d - p.
z()*v.
z()*tanETheta2 ;
1472 c = dist2ETheta/
t1 ;
1479 zi = p.
z() + sd*v.
z();
1481 if ( (sd < 0) || (zi*(eTheta -
halfpi) > 0) )
1485 if ( (sd >= 0) && (sd < snxt) )
1487 xi = p.
x() + sd*v.
x() ;
1488 yi = p.
y() + sd*v.
y() ;
1489 zi = p.
z() + sd*v.
z() ;
1490 rhoi2 = xi*xi + yi*yi ;
1491 radi2 = rhoi2 + zi*zi ;
1493 if ( (radi2 <= tolORMax2)
1494 && (radi2 >= tolORMin2)
1495 && (zi*(eTheta -
halfpi) <= 0) )
1497 if (!fFullPhiSphere && rhoi2)
1499 cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
1500 if (cosPsi >= cosHDPhiOT)
1519 t1 = 1 - v.
z()*v.
z()*(1 + tanSTheta2) ;
1520 t2 = pDotV2d - p.
z()*v.
z()*tanSTheta2 ;
1524 c = dist2STheta/
t1 ;
1532 if ( (sd >= 0) && (sd < snxt) )
1534 xi = p.
x() + sd*v.
x() ;
1535 yi = p.
y() + sd*v.
y() ;
1536 zi = p.
z() + sd*v.
z() ;
1537 rhoi2 = xi*xi + yi*yi ;
1538 radi2 = rhoi2 + zi*zi ;
1540 if ( (radi2 <= tolORMax2)
1541 && (radi2 >= tolORMin2)
1542 && (zi*(fSTheta -
halfpi) <= 0) )
1544 if (!fFullPhiSphere && rhoi2)
1546 cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
1547 if (cosPsi >= cosHDPhiOT)
1562 else if ( (pTheta < tolSTheta + kAngTolerance)
1563 && (fSTheta > halfAngTolerance) )
1569 t2 = pDotV2d - p.
z()*v.
z()*tanSTheta2 ;
1570 if ( (t2>=0 && tolIRMin2<rad2 && rad2<tolIRMax2 && fSTheta<
halfpi)
1571 || (t2<0 && tolIRMin2<rad2 && rad2<tolIRMax2 && fSTheta>
halfpi)
1572 || (v.
z()<0 && tolIRMin2<rad2 && rad2<tolIRMax2 && fSTheta==
halfpi) )
1574 if (!fFullPhiSphere && rho2)
1576 cosPsi = (p.
x()*cosCPhi + p.
y()*sinCPhi)/std::sqrt(rho2) ;
1577 if (cosPsi >= cosHDPhiIT)
1590 t1 = 1 - v.
z()*v.
z()*(1 + tanSTheta2) ;
1594 c = dist2STheta/
t1 ;
1601 if ( (sd >= halfCarTolerance) && (sd < snxt) && (fSTheta < halfpi) )
1603 xi = p.
x() + sd*v.
x() ;
1604 yi = p.
y() + sd*v.
y() ;
1605 zi = p.
z() + sd*v.
z() ;
1606 rhoi2 = xi*xi + yi*yi ;
1607 radi2 = rhoi2 + zi*zi ;
1609 if ( (radi2 <= tolORMax2)
1610 && (radi2 >= tolORMin2)
1611 && (zi*(fSTheta - halfpi) <= 0) )
1613 if ( !fFullPhiSphere && rhoi2 )
1615 cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
1616 if ( cosPsi >= cosHDPhiOT )
1630 else if ((pTheta > tolETheta-kAngTolerance) && (eTheta <
pi-kAngTolerance))
1637 t2 = pDotV2d - p.
z()*v.
z()*tanETheta2 ;
1639 if ( ((t2<0) && (eTheta < halfpi)
1640 && (tolIRMin2 < rad2) && (rad2 < tolIRMax2))
1641 || ((t2>=0) && (eTheta > halfpi)
1642 && (tolIRMin2 < rad2) && (rad2 < tolIRMax2))
1643 || ((v.
z()>0) && (eTheta == halfpi)
1644 && (tolIRMin2 < rad2) && (rad2 < tolIRMax2)) )
1646 if (!fFullPhiSphere && rho2)
1648 cosPsi = (p.
x()*cosCPhi + p.
y()*sinCPhi)/std::sqrt(rho2) ;
1649 if (cosPsi >= cosHDPhiIT)
1662 t1 = 1 - v.
z()*v.
z()*(1 + tanETheta2) ;
1666 c = dist2ETheta/
t1 ;
1674 if ( (sd >= halfCarTolerance)
1675 && (sd < snxt) && (eTheta > halfpi) )
1677 xi = p.
x() + sd*v.
x() ;
1678 yi = p.
y() + sd*v.
y() ;
1679 zi = p.
z() + sd*v.
z() ;
1680 rhoi2 = xi*xi + yi*yi ;
1681 radi2 = rhoi2 + zi*zi ;
1683 if ( (radi2 <= tolORMax2)
1684 && (radi2 >= tolORMin2)
1685 && (zi*(eTheta - halfpi) <= 0) )
1687 if (!fFullPhiSphere && rhoi2)
1689 cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
1690 if (cosPsi >= cosHDPhiOT)
1709 t1 = 1 - v.
z()*v.
z()*(1 + tanSTheta2) ;
1710 t2 = pDotV2d - p.
z()*v.
z()*tanSTheta2 ;
1714 c = dist2STheta/
t1 ;
1722 if ((sd >= 0) && (sd < snxt))
1724 xi = p.
x() + sd*v.
x() ;
1725 yi = p.
y() + sd*v.
y() ;
1726 zi = p.
z() + sd*v.
z() ;
1727 rhoi2 = xi*xi + yi*yi ;
1728 radi2 = rhoi2 + zi*zi ;
1730 if ( (radi2 <= tolORMax2)
1731 && (radi2 >= tolORMin2)
1732 && (zi*(fSTheta - halfpi) <= 0) )
1734 if (!fFullPhiSphere && rhoi2)
1736 cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
1737 if (cosPsi >= cosHDPhiOT)
1750 t1 = 1 - v.
z()*v.
z()*(1 + tanETheta2) ;
1751 t2 = pDotV2d - p.
z()*v.
z()*tanETheta2 ;
1755 c = dist2ETheta/
t1 ;
1763 if ((sd >= 0) && (sd < snxt))
1765 xi = p.
x() + sd*v.
x() ;
1766 yi = p.
y() + sd*v.
y() ;
1767 zi = p.
z() + sd*v.
z() ;
1768 rhoi2 = xi*xi + yi*yi ;
1769 radi2 = rhoi2 + zi*zi ;
1771 if ( (radi2 <= tolORMax2)
1772 && (radi2 >= tolORMin2)
1773 && (zi*(eTheta - halfpi) <= 0) )
1775 if (!fFullPhiSphere && rhoi2)
1777 cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
1778 if ( cosPsi >= cosHDPhiOT )
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const