59 : fCheck(false), fVerbose(0)
80 const G4int replicaNo,
92 G4double coord, rad2, rmin, tolRMax2, rmax, tolRMin2;
101 coord = std::fabs(localPoint(axis))-width*0.5;
112 if ( localPoint.y()||localPoint.x() )
114 coord = std::fabs(std::atan2(localPoint.y(),localPoint.x()))-width*0.5;
130 rad2 = localPoint.perp2();
131 rmax = (replicaNo+1)*width+offset;
133 tolRMax2 *= tolRMax2;
137 tolRMax2 *= tolRMax2;
138 if ( rad2<=tolRMax2 )
147 if ( replicaNo||offset )
151 tolRMin2 *= tolRMin2;
155 tolRMin2 *= tolRMin2;
156 if ( rad2>=tolRMin2 )
173 G4Exception(
"G4ReplicaNavigation::Inside()",
"GeomNav0002",
186 const G4int replicaNo,
206 coord = localPoint(axis);
207 safe1 = width*0.5-coord;
208 safe2 = width*0.5+coord;
209 safety = (safe1<=safe2) ? safe1 : safe2;
212 if ( localPoint.y()<=0 )
214 safety = localPoint.x()*std::sin(width*0.5)
215 + localPoint.y()*std::cos(width*0.5);
219 safety = localPoint.x()*std::sin(width*0.5)
220 - localPoint.y()*std::cos(width*0.5);
224 rho = localPoint.perp();
225 rmax = width*(replicaNo+1)+offset;
226 if ( replicaNo||offset )
231 safety = (safe1<=safe2) ? safe1 : safe2;
239 G4Exception(
"G4ReplicaNavigation::DistanceToOut()",
"GeomNav0002",
252 const G4int replicaNo,
275 coord = localPoint(axis);
276 Comp = localDirection(axis);
279 lindist = width*0.5-coord;
280 Dist = (lindist>0) ? lindist/Comp : 0;
285 lindist = width*0.5+coord;
286 Dist = (lindist>0) ? -lindist/Comp : 0;
294 candidateNormal.
exitNormal = ( signC * VecCartAxes[axis]);
298 (Comp>0) ? SideCartAxesPlus[axis] : SideCartAxesMinus[axis];
306 replicaNo,candidateNormal);
310 G4Exception(
"G4ReplicaNavigation::DistanceToOut()",
"GeomNav0002",
315 arExitNormal= candidateNormal;
333 G4double sinSPhi= -2.0, cosSPhi= -2.0;
334 G4double pDistS, pDistE, compS, compE, Dist, dist2, yi;
338 if ( (localPoint.x()!=0.0) || (localPoint.y()!=0.0) )
340 sinSPhi = std::sin(-width*0.5);
341 cosSPhi = std::cos(width*0.5);
345 pDistS = localPoint.x()*sinSPhi-localPoint.y()*cosSPhi;
347 pDistE = localPoint.x()*sinSPhi+localPoint.y()*cosSPhi;
352 compS = -sinSPhi*localDirection.x()+cosSPhi*localDirection.y();
353 compE = -sinSPhi*localDirection.x()-cosSPhi*localDirection.y();
355 if ( (pDistS<=0)&&(pDistE<=0) )
361 dist2 = pDistS/compS;
362 yi = localPoint.y()+dist2*localDirection.y();
382 dist2 = pDistE/compE;
388 yi = localPoint.y()+dist2*localDirection.y();
402 else if ( (pDistS>=0)&&(pDistE>=0) )
407 Dist = ((compS>=0)&&(compE>=0)) ?
kInfinity : 0;
409 else if ( (pDistS>0)&&(pDistE<0) )
415 dist2 = pDistE/compE;
416 yi = localPoint.y()+dist2*localDirection.y();
438 dist2 = pDistS/compS;
439 yi = localPoint.y()+dist2*localDirection.y();
465 if( (std::fabs(localDirection.phi())<=width*0.5) )
503 const G4int replicaNo,
506 G4double rmin, rmax, t1, t2, t3, deltaR;
526 rmin = replicaNo*width+offset;
527 rmax = (replicaNo+1)*width+offset;
529 t1 = 1.0-localDirection.z()*localDirection.z();
530 t2 = localPoint.x()*localDirection.x()+localPoint.y()*localDirection.y();
531 t3 = localPoint.x()*localPoint.x()+localPoint.y()*localPoint.y();
541 deltaR = t3-rmax*rmax;
550 srd = -b+std::sqrt(b*b-c);
568 deltaR = t3-rmin*rmin;
587 deltaR = t3-rmax*rmax;
590 srd = (d2 < 0.) ? 0.0 : -b+std::sqrt(d2);
598 deltaR = t3-rmax*rmax;
602 srd = (d2 < 0.) ? 0.0 : -b+std::sqrt(d2);
618 xi = localPoint.x() + srd*localDirection.x() ;
619 yi = localPoint.y() + srd*localDirection.y() ;
625 normalR *= (-1.0)/rmin;
663 val = -width*0.5*(nReplicas-1)+width*replicaNo;
665 point.setX(point.x()-val);
668 val = -width*0.5*(nReplicas-1)+width*replicaNo;
670 point.setY(point.y()-val);
673 val = -width*0.5*(nReplicas-1)+width*replicaNo;
675 point.setZ(point.z()-val);
678 val = -(offset+width*(replicaNo+0.5));
680 cosv = std::cos(val);
681 sinv = std::sin(val);
682 tmpx = point.x()*cosv-point.y()*sinv;
683 tmpy = point.x()*sinv+point.y()*cosv;
718 val = -width*0.5*(nReplicas-1)+width*replicaNo;
722 val = -width*0.5*(nReplicas-1)+width*replicaNo;
726 val = -width*0.5*(nReplicas-1)+width*replicaNo;
730 val = -(offset+width*(replicaNo+0.5));
749 const G4double currentProposedStepLength,
754 G4bool &calculatedExitNormal,
759 G4int &blockedReplicaNo )
766 G4double ourStep=currentProposedStepLength;
768 G4double sampleStep, sampleSafety, motherStep, motherSafety;
769 G4int localNoDaughters, sampleNo;
774 calculatedExitNormal=
false;
778 if ( exiting&&validExitNormal )
784 blockedExitedVol = *pBlockedPhysical;
804 ourSafety =
std::min( ourSafety, sampleSafety);
806 if ( sampleSafety<ourStep )
814 if ( sampleStep<ourStep )
816 ourStep = sampleStep;
818 validExitNormal = normalOutStc.validConvex;
820 exitNormalStc= normalOutStc;
822 TransformAxis(normalOutStc.exitNormal);
823 calculatedExitNormal=
true;
826 const G4int secondDepth= topDepth;
838 if ( sampleSafety < ourSafety )
840 ourSafety = sampleSafety;
842 if ( sampleSafety < ourStep )
850 if ( sampleStep < ourStep )
852 ourStep = sampleStep;
861 exitNormalStc= normalOutStc;
863 calculatedExitNormal=
true;
876 motherPhysical = history.
GetVolume(depth);
881 motherStep = motherSolid->
DistanceToOut(repPoint,repDirection,
true,
882 &exitConvex,&exitVectorMother);
885 motherNormalStc =
G4ExitNormal( exitVectorMother,
true,
false,
887 calculatedExitNormal=
true;
891 G4bool motherDeterminedStep= (motherStep<ourStep);
893 if( (!exitConvex) && motherDeterminedStep )
896 motherNormalStc=
G4ExitNormal( exitVectorMother,
true,
false,
902 calculatedExitNormal=
true;
904 if( motherDeterminedStep)
909 exitNormalStc= motherNormalStc;
910 exitNormalStc.
exitNormal= globalExitNormalTop;
926 if ( motherSafety<ourSafety )
928 ourSafety = motherSafety;
936 std::ostringstream message;
937 message <<
"Point outside volume !" <<
G4endl
938 <<
" Point " << localPoint
939 <<
" is outside current volume " << motherPhysical->
GetName()
942 message <<
" Estimated isotropic distance to solid (distToIn)= "
943 << estDistToSolid <<
G4endl;
949 "Point is far outside Current Volume !" );
954 "Point is a little outside Current Volume.");
962 if( motherDeterminedStep)
964 ourStep = motherStep;
970 if ( calculatedExitNormal )
972 if ( motherDeterminedStep )
977 exitNormalVector= globalToLocalTop.
TransformAxis(exitNormalGlobal);
985 exitNormalVector *= rot->inverse();
991 validExitNormal =
false;
995 if ( motherSafety<=ourStep )
997 if ( motherStep<=ourStep )
999 ourStep = motherStep;
1001 if ( validExitNormal )
1006 exitNormal *= rot->inverse();
1012 validExitNormal =
false;
1019 G4bool daughterDeterminedStep=
false;
1028 for ( sampleNo=localNoDaughters-1; sampleNo>=0; sampleNo-- )
1030 samplePhysical = repLogical->
GetDaughter(sampleNo);
1031 if ( samplePhysical!=blockedExitedVol )
1043 const G4double sampleSafetyDistance =
1045 if ( sampleSafetyDistance<ourSafety )
1047 ourSafety = sampleSafetyDistance;
1049 if ( sampleSafetyDistance<=ourStep )
1051 sampleDirection = sampleTf.TransformAxis(localDirection);
1052 const G4double sampleStepDistance =
1053 sampleSolid->
DistanceToIn(samplePoint,sampleDirection);
1054 if ( sampleStepDistance<=ourStep )
1056 daughterDeterminedStep=
true;
1058 ourStep = sampleStepDistance;
1061 *pBlockedPhysical = samplePhysical;
1062 blockedReplicaNo = -1;
1064 #ifdef DAUGHTER_NORMAL_ALSO
1067 daughtNormRepCrd = sampleTf.Inverse().TransformAxis(localExitNorm);
1077 intersectionPoint= samplePoint
1078 + sampleStepDistance * sampleDirection;
1079 EInside insideIntPt= sampleSolid->
Inside(intersectionPoint);
1083 std::ostringstream message;
1084 message <<
"Navigator gets conflicting response from Solid."
1086 <<
" Inaccurate DistanceToIn for solid "
1088 <<
" Solid gave DistanceToIn = "
1089 << sampleStepDistance <<
" yet returns " ;
1091 message <<
"-kInside-";
1092 else if ( insideIntPt ==
kOutside )
1093 message <<
"-kOutside-";
1095 message <<
"-kSurface-";
1096 message <<
" for this point !" <<
G4endl
1097 <<
" Point = " << intersectionPoint <<
G4endl;
1099 message <<
" DistanceToIn(p) = "
1103 message <<
" DistanceToOut(p) = "
1107 G4cout.precision(oldcoutPrec);
1116 calculatedExitNormal &= (!daughterDeterminedStep);
1118 #ifdef DAUGHTER_NORMAL_ALSO
1119 if( daughterDeterminedStep )
1126 validExitNormal=
false;
1128 calculatedExitNormal=
true;
1133 newSafety = ourSafety;
1157 G4int localNoDaughters, sampleNo;
1170 if ( sampleSafety<ourSafety )
1172 ourSafety = sampleSafety;
1182 if ( sampleSafety<ourSafety )
1184 ourSafety = sampleSafety;
1192 motherPhysical = history.
GetVolume(depth);
1196 if ( sampleSafety<ourSafety )
1198 ourSafety = sampleSafety;
1204 for ( sampleNo=localNoDaughters-1; sampleNo>=0; sampleNo-- )
1206 samplePhysical = repLogical->
GetDaughter(sampleNo);
1207 if ( samplePhysical!=blockedExitedVol )
1216 const G4double sampleSafetyDistance =
1218 if ( sampleSafetyDistance<ourSafety )
1220 ourSafety = sampleSafetyDistance;
1236 G4bool ¬KnownInside )
const
1241 G4int mdepth, depth, cdepth;
1248 for ( mdepth=cdepth-1; mdepth>=0; mdepth-- )
1262 G4Exception(
"G4ReplicaNavigation::BackLocate()",
"GeomNav0002",
1269 insideCode = motherSolid->
Inside(goodPoint);
1281 notKnownInside =
false;
1286 for ( depth=mdepth+1; depth<cdepth; depth++)
1294 localPoint = goodPoint;
1300 goodPoint = repPoint;
1313 localPoint = goodPoint;
G4VPhysicalVolume * GetTopVolume() const
const G4ThreeVector & GetTranslation() const
static const G4double kInfinity
CLHEP::Hep3Vector G4ThreeVector
CLHEP::HepRotation G4RotationMatrix
EInside Inside(const G4VPhysicalVolume *pVol, const G4int replicaNo, const G4ThreeVector &localPoint) const
G4double GetSurfaceTolerance() const
G4VPhysicalVolume * GetDaughter(const G4int i) const
EVolume GetVolumeType(G4int n) const
G4double ComputeStep(const G4ThreeVector &globalPoint, const G4ThreeVector &globalDirection, const G4ThreeVector &localPoint, const G4ThreeVector &localDirection, const G4double currentProposedStepLength, G4double &newSafety, G4NavigationHistory &history, G4bool &validExitNormal, G4bool &calculatedExitNormal, G4ThreeVector &exitNormal, G4bool &exiting, G4bool &entering, G4VPhysicalVolume *(*pBlockedPhysical), G4int &blockedReplicaNo)
G4double DistanceToOut(const G4VPhysicalVolume *pVol, const G4int replicaNo, const G4ThreeVector &localPoint) const
void ComputeTransformation(const G4int replicaNo, G4VPhysicalVolume *pVol, G4ThreeVector &point) const
G4int GetTopReplicaNo() const
G4GLOB_DLL std::ostream G4cout
const G4String & GetName() const
virtual EInside Inside(const G4ThreeVector &p) const =0
G4double GetRadialTolerance() const
virtual G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const =0
G4int GetReplicaNo(G4int n) const
void SetTranslation(const G4ThreeVector &v)
virtual G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const =0
static const double kMinExitingNormalCosine
G4int GetNoDaughters() const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
G4LogicalVolume * GetLogicalVolume() const
const G4AffineTransform & GetTransform(G4int n) const
G4double DistanceToOutRad(const G4ThreeVector &localPoint, const G4ThreeVector &localDirection, const G4double width, const G4double offset, const G4int replicaNo, G4ExitNormal &foundNormal) const
const G4RotationMatrix * GetRotation() const
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
G4double DistanceToOutPhi(const G4ThreeVector &localPoint, const G4ThreeVector &localDirection, const G4double width, G4ExitNormal &foundNormal) const
virtual void GetReplicationData(EAxis &axis, G4int &nReplicas, G4double &width, G4double &offset, G4bool &consuming) const =0
const G4AffineTransform & GetTopTransform() const
G4double ComputeSafety(const G4ThreeVector &globalPoint, const G4ThreeVector &localPoint, G4NavigationHistory &history, const G4double pProposedMaxLength=DBL_MAX)
G4VPhysicalVolume * GetVolume(G4int n) const
void SetPhiTransformation(const G4double ang, G4VPhysicalVolume *pVol=0) const
virtual G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=false, G4bool *validNorm=0, G4ThreeVector *n=0) const =0
EInside BackLocate(G4NavigationHistory &history, const G4ThreeVector &globalPoint, G4ThreeVector &localPoint, const G4bool &exiting, G4bool ¬KnownInside) const
G4double GetAngularTolerance() const
static G4GeometryTolerance * GetInstance()
G4VSolid * GetSolid() const