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.; }
G4double halfRadTolerance
static const G4double kInfinity
CLHEP::Hep3Vector G4ThreeVector
G4double halfCarTolerance
double dot(const Hep3Vector &) const
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const