78 const G4int replicaNo,
99 coord = std::fabs(localPoint(axis))-width*0.5;
110 if ( localPoint.
y()||localPoint.
x() )
112 coord = std::fabs(std::atan2(localPoint.
y(),localPoint.
x()))-width*0.5;
128 rad2 = localPoint.
perp2();
129 rmax = (replicaNo+1)*width+offset;
131 tolRMax2 *= tolRMax2;
135 tolRMax2 *= tolRMax2;
136 if ( rad2<=tolRMax2 )
145 if ( replicaNo||offset )
149 tolRMin2 *= tolRMin2;
153 tolRMin2 *= tolRMin2;
154 if ( rad2>=tolRMin2 )
171 G4Exception(
"G4ReplicaNavigation::Inside()",
"GeomNav0002",
184 const G4int replicaNo,
204 coord = localPoint(axis);
205 safe1 = width*0.5-coord;
206 safe2 = width*0.5+coord;
207 safety = (safe1<=safe2) ? safe1 : safe2;
210 if ( localPoint.
y()<=0 )
212 safety = localPoint.
x()*std::sin(width*0.5)
213 + localPoint.
y()*std::cos(width*0.5);
217 safety = localPoint.
x()*std::sin(width*0.5)
218 - localPoint.
y()*std::cos(width*0.5);
222 rho = localPoint.
perp();
223 rmax = width*(replicaNo+1)+offset;
224 if ( replicaNo||offset )
229 safety = (safe1<=safe2) ? safe1 : safe2;
237 G4Exception(
"G4ReplicaNavigation::DistanceToOut()",
"GeomNav0002",
250 const G4int replicaNo,
273 coord = localPoint(axis);
274 Comp = localDirection(axis);
277 lindist = width*0.5-coord;
278 Dist = (lindist>0) ? lindist/Comp : 0;
283 lindist = width*0.5+coord;
284 Dist = (lindist>0) ? -lindist/Comp : 0;
292 candidateNormal.
exitNormal = ( signC * VecCartAxes[axis]);
296 (Comp>0) ? SideCartAxesPlus[axis] : SideCartAxesMinus[axis];
304 replicaNo,candidateNormal);
308 G4Exception(
"G4ReplicaNavigation::DistanceToOut()",
"GeomNav0002",
313 arExitNormal= candidateNormal;
331 G4double sinSPhi = -2.0, cosSPhi = -2.0;
332 G4double pDistS, pDistE, compS, compE, Dist, dist2, yi;
336 if ( (localPoint.
x()!=0.0) || (localPoint.
y()!=0.0) )
338 sinSPhi = std::sin(-width*0.5);
339 cosSPhi = std::cos(width*0.5);
343 pDistS = localPoint.
x()*sinSPhi-localPoint.
y()*cosSPhi;
345 pDistE = localPoint.
x()*sinSPhi+localPoint.
y()*cosSPhi;
350 compS = -sinSPhi*localDirection.
x()+cosSPhi*localDirection.
y();
351 compE = -sinSPhi*localDirection.
x()-cosSPhi*localDirection.
y();
359 dist2 = pDistS/compS;
360 yi = localPoint.
y()+dist2*localDirection.
y();
380 dist2 = pDistE/compE;
386 yi = localPoint.
y()+dist2*localDirection.
y();
405 Dist = ((compS>=0)&&(compE>=0)) ?
kInfinity : 0;
413 dist2 = pDistE/compE;
414 yi = localPoint.
y()+dist2*localDirection.
y();
436 dist2 = pDistS/compS;
437 yi = localPoint.
y()+dist2*localDirection.
y();
463 if( (std::fabs(localDirection.
phi())<=width*0.5) )
501 const G4int replicaNo,
524 rmin = replicaNo*width+
offset;
525 rmax = (replicaNo+1)*width+offset;
527 t1 = 1.0-localDirection.
z()*localDirection.
z();
528 t2 = localPoint.
x()*localDirection.
x()+localPoint.
y()*localDirection.
y();
529 t3 = localPoint.
x()*localPoint.
x()+localPoint.
y()*localPoint.
y();
539 deltaR = t3-rmax*
rmax;
548 srd = -b+std::sqrt(b*b-c);
566 deltaR = t3-rmin*rmin;
585 deltaR = t3-rmax*
rmax;
588 srd = (d2 < 0.) ? 0.0 : -b+std::sqrt(d2);
596 deltaR = t3-rmax*
rmax;
600 srd = (d2 < 0.) ? 0.0 : -b+std::sqrt(d2);
617 xi = localPoint.
x() + srd*localDirection.
x();
618 yi = localPoint.
y() + srd*localDirection.
y();
627 normalR *= (-1.0)/rmin;
667 val = -width*0.5*(nReplicas-1)+width*replicaNo;
669 point.
setX(point.
x()-val);
672 val = -width*0.5*(nReplicas-1)+width*replicaNo;
674 point.
setY(point.
y()-val);
677 val = -width*0.5*(nReplicas-1)+width*replicaNo;
679 point.
setZ(point.
z()-val);
682 val = -(offset+width*(replicaNo+0.5));
684 cosv = std::cos(val);
685 sinv = std::sin(val);
686 tmpx = point.
x()*cosv-point.
y()*sinv;
687 tmpy = point.
x()*sinv+point.
y()*cosv;
722 val = -width*0.5*(nReplicas-1)+width*replicaNo;
726 val = -width*0.5*(nReplicas-1)+width*replicaNo;
730 val = -width*0.5*(nReplicas-1)+width*replicaNo;
734 val = -(offset+width*(replicaNo+0.5));
753 const G4double currentProposedStepLength,
758 G4bool& calculatedExitNormal,
763 G4int& blockedReplicaNo )
770 G4double ourStep=currentProposedStepLength;
772 G4double sampleStep, sampleSafety, motherStep, motherSafety;
773 G4int localNoDaughters, sampleNo;
778 calculatedExitNormal=
false;
782 if ( exiting&&validExitNormal )
788 blockedExitedVol = *pBlockedPhysical;
808 ourSafety =
std::min( ourSafety, sampleSafety);
810 if ( sampleSafety<ourStep )
818 if ( sampleStep<ourStep )
820 ourStep = sampleStep;
822 validExitNormal = normalOutStc.validConvex;
824 exitNormalStc = normalOutStc;
827 calculatedExitNormal =
true;
830 const G4int secondDepth = topDepth;
843 if ( sampleSafety < ourSafety )
845 ourSafety = sampleSafety;
847 if ( sampleSafety < ourStep )
856 if ( sampleStep < ourStep )
858 ourStep = sampleStep;
867 exitNormalStc = normalOutStc;
869 calculatedExitNormal =
true;
878 G4bool exitConvex =
false;
882 motherPhysical = history.
GetVolume(depth);
887 motherStep = motherSolid->
DistanceToOut(repPoint,repDirection,
true,
888 &exitConvex,&exitVectorMother);
891 motherNormalStc =
G4ExitNormal( exitVectorMother,
true,
false,
893 calculatedExitNormal =
true;
897 G4bool motherDeterminedStep = (motherStep<ourStep);
899 if( (!exitConvex) && motherDeterminedStep )
902 motherNormalStc =
G4ExitNormal( exitVectorMother,
true,
false,
908 calculatedExitNormal =
true;
910 if( motherDeterminedStep )
915 exitNormalStc = motherNormalStc;
916 exitNormalStc.
exitNormal = globalExitNormalTop;
932 if ( motherSafety<ourSafety )
934 ourSafety = motherSafety;
943 message <<
"Point outside volume !" <<
G4endl
944 <<
" Point " << localPoint
945 <<
" is outside current volume " << motherPhysical->
GetName()
948 message <<
" Estimated isotropic distance to solid (distToIn)= "
949 << estDistToSolid <<
G4endl;
955 "Point is far outside Current Volume !" );
960 "Point is a little outside Current Volume.");
968 if( motherDeterminedStep )
970 ourStep = motherStep;
976 if ( calculatedExitNormal )
978 if ( motherDeterminedStep )
980 exitNormalVector = motherNormalStc.
exitNormal;
985 exitNormalVector = globalToLocalTop.
TransformAxis(exitNormalGlobal);
993 exitNormalVector *= rot->
inverse();
999 validExitNormal =
false;
1003 if ( motherSafety<=ourStep )
1005 if ( motherStep<=ourStep )
1007 ourStep = motherStep;
1009 if ( validExitNormal )
1020 validExitNormal =
false;
1027 G4bool daughterDeterminedStep =
false;
1036 for ( sampleNo=localNoDaughters-1; sampleNo>=0; sampleNo-- )
1038 samplePhysical = repLogical->
GetDaughter(sampleNo);
1039 if ( samplePhysical!=blockedExitedVol )
1048 sampleTf.TransformPoint(localPoint);
1051 const G4double sampleSafetyDistance =
1053 if ( sampleSafetyDistance<ourSafety )
1055 ourSafety = sampleSafetyDistance;
1057 if ( sampleSafetyDistance<=ourStep )
1059 sampleDirection = sampleTf.TransformAxis(localDirection);
1060 const G4double sampleStepDistance =
1061 sampleSolid->
DistanceToIn(samplePoint,sampleDirection);
1062 if ( sampleStepDistance<=ourStep )
1064 daughterDeterminedStep =
true;
1066 ourStep = sampleStepDistance;
1069 *pBlockedPhysical = samplePhysical;
1070 blockedReplicaNo = sampleNo;
1072 #ifdef DAUGHTER_NORMAL_ALSO
1075 daughtNormRepCrd = sampleTf.InverseTransformAxis(localExitNorm);
1085 intersectionPoint = samplePoint
1086 + sampleStepDistance * sampleDirection;
1087 EInside insideIntPt = sampleSolid->
Inside(intersectionPoint);
1092 message <<
"Navigator gets conflicting response from Solid."
1094 <<
" Inaccurate DistanceToIn for solid "
1096 <<
" Solid gave DistanceToIn = "
1097 << sampleStepDistance <<
" yet returns " ;
1099 message <<
"-kInside-";
1100 else if ( insideIntPt ==
kOutside )
1101 message <<
"-kOutside-";
1103 message <<
"-kSurface-";
1104 message <<
" for this point !" <<
G4endl
1105 <<
" Point = " << intersectionPoint <<
G4endl;
1107 message <<
" DistanceToIn(p) = "
1111 message <<
" DistanceToOut(p) = "
1115 G4cout.precision(oldcoutPrec);
1124 calculatedExitNormal &= (!daughterDeterminedStep);
1126 #ifdef DAUGHTER_NORMAL_ALSO
1127 if( daughterDeterminedStep )
1134 validExitNormal =
false;
1136 calculatedExitNormal =
true;
1141 newSafety = ourSafety;
1165 G4int localNoDaughters, sampleNo;
1178 if ( sampleSafety<ourSafety )
1180 ourSafety = sampleSafety;
1192 if ( sampleSafety<ourSafety )
1194 ourSafety = sampleSafety;
1202 motherPhysical = history.
GetVolume(depth);
1206 if ( sampleSafety<ourSafety )
1208 ourSafety = sampleSafety;
1214 for ( sampleNo=localNoDaughters-1; sampleNo>=0; sampleNo-- )
1216 samplePhysical = repLogical->
GetDaughter(sampleNo);
1217 if ( samplePhysical!=blockedExitedVol )
1223 sampleTf.TransformPoint(localPoint);
1226 const G4double sampleSafetyDistance =
1228 if ( sampleSafetyDistance<ourSafety )
1230 ourSafety = sampleSafetyDistance;
1246 G4bool& notKnownInside )
const
1251 G4int mdepth, depth, cdepth;
1258 for ( mdepth=cdepth-1; mdepth>=0; mdepth-- )
1267 if( pNRMother ==
nullptr )
1272 G4Exception(
"G4ReplicaNavigation::BackLocate()",
"GeomNav0002",
1279 insideCode = motherSolid->
Inside(goodPoint);
1291 notKnownInside =
false;
1296 for ( depth=mdepth+1; depth<cdepth; ++depth)
1304 localPoint = goodPoint;
1310 goodPoint = repPoint;
1323 localPoint = goodPoint;