58 : fCheck(false), fVerbose(0)
79 const G4int replicaNo,
95 G4double coord, rad2, rmin, tolRMax2, rmax, tolRMin2;
104 coord = std::fabs(localPoint(axis))-width*0.5;
105 if ( coord<=-halfkCarTolerance )
109 else if ( coord<=halfkCarTolerance )
115 if ( localPoint.y()||localPoint.x() )
117 coord = std::fabs(std::atan2(localPoint.y(),localPoint.x()))-width*0.5;
118 if ( coord<=-halfkAngTolerance )
122 else if ( coord<=halfkAngTolerance )
133 rad2 = localPoint.perp2();
134 rmax = (replicaNo+1)*width+offset;
135 tolRMax2 = rmax-halfkRadTolerance;
136 tolRMax2 *= tolRMax2;
139 tolRMax2 = rmax+halfkRadTolerance;
140 tolRMax2 *= tolRMax2;
141 if ( rad2<=tolRMax2 )
150 if ( replicaNo||offset )
153 tolRMin2 = rmin-halfkRadTolerance;
154 tolRMin2 *= tolRMin2;
157 tolRMin2 = rmin+halfkRadTolerance;
158 tolRMin2 *= tolRMin2;
159 if ( rad2>=tolRMin2 )
176 G4Exception(
"G4ReplicaNavigation::Inside()",
"GeomNav0002",
189 const G4int replicaNo,
210 coord = localPoint(axis);
211 safe1 = width*0.5-coord;
212 safe2 = width*0.5+coord;
213 safety = (safe1<=safe2) ? safe1 : safe2;
216 if ( localPoint.y()<=0 )
218 safety = localPoint.x()*std::sin(width*0.5)
219 + localPoint.y()*std::cos(width*0.5);
223 safety = localPoint.x()*std::sin(width*0.5)
224 - localPoint.y()*std::cos(width*0.5);
228 rho = localPoint.perp();
229 rmax = width*(replicaNo+1)+offset;
230 if ( replicaNo||offset )
235 safety = (safe1<=safe2) ? safe1 : safe2;
243 G4Exception(
"G4ReplicaNavigation::DistanceToOut()",
"GeomNav0002",
247 return (safety >= halfkCarTolerance) ? safety : 0;
256 const G4int replicaNo,
279 coord = localPoint(axis);
280 Comp = localDirection(axis);
283 lindist = width*0.5-coord;
284 Dist = (lindist>0) ? lindist/Comp : 0;
289 lindist = width*0.5+coord;
290 Dist = (lindist>0) ? -lindist/Comp : 0;
298 candidateNormal.
exitNormal = ( signC * VecCartAxes[axis]);
302 (Comp>0) ? SideCartAxesPlus[axis] : SideCartAxesMinus[axis];
310 replicaNo,candidateNormal);
314 G4Exception(
"G4ReplicaNavigation::DistanceToOut()",
"GeomNav0002",
319 arExitNormal= candidateNormal;
337 G4double sinSPhi= -2.0, cosSPhi= -2.0;
338 G4double pDistS, pDistE, compS, compE, Dist, dist2, yi;
343 if ( (localPoint.x()!=0.0) || (localPoint.y()!=0.0) )
345 sinSPhi = std::sin(-width*0.5);
346 cosSPhi = std::cos(width*0.5);
350 pDistS = localPoint.x()*sinSPhi-localPoint.y()*cosSPhi;
352 pDistE = localPoint.x()*sinSPhi+localPoint.y()*cosSPhi;
357 compS = -sinSPhi*localDirection.x()+cosSPhi*localDirection.y();
358 compE = -sinSPhi*localDirection.x()-cosSPhi*localDirection.y();
360 if ( (pDistS<=0)&&(pDistE<=0) )
366 dist2 = pDistS/compS;
367 yi = localPoint.y()+dist2*localDirection.y();
373 Dist = (pDistS<=-halfkCarTolerance) ? dist2 : 0;
387 dist2 = pDistE/compE;
393 yi = localPoint.y()+dist2*localDirection.y();
401 Dist = (pDistE<=-halfkCarTolerance) ? dist2 : 0;
407 else if ( (pDistS>=0)&&(pDistE>=0) )
412 Dist = ((compS>=0)&&(compE>=0)) ?
kInfinity : 0;
414 else if ( (pDistS>0)&&(pDistE<0) )
420 dist2 = pDistE/compE;
421 yi = localPoint.y()+dist2*localDirection.y();
443 dist2 = pDistS/compS;
444 yi = localPoint.y()+dist2*localDirection.y();
470 if( (std::fabs(localDirection.phi())<=width*0.5) )
508 const G4int replicaNo,
511 G4double rmin, rmax, t1, t2, t3, deltaR;
532 rmin = replicaNo*width+offset;
533 rmax = (replicaNo+1)*width+offset;
535 t1 = 1.0-localDirection.z()*localDirection.z();
536 t2 = localPoint.x()*localDirection.x()+localPoint.y()*localDirection.y();
537 t3 = localPoint.x()*localPoint.x()+localPoint.y()*localPoint.y();
547 deltaR = t3-rmax*rmax;
552 if ( deltaR<-halfkRadTolerance )
556 srd = -b+std::sqrt(b*b-c);
574 deltaR = t3-rmin*rmin;
584 srd = (deltaR>halfkRadTolerance) ? -b-std::sqrt(d2) : 0.0;
593 deltaR = t3-rmax*rmax;
596 srd = (d2 < 0.) ? 0.0 : -b+std::sqrt(d2);
604 deltaR = t3-rmax*rmax;
608 srd = (d2 < 0.) ? 0.0 : -b+std::sqrt(d2);
624 xi = localPoint.x() + srd*localDirection.x() ;
625 yi = localPoint.y() + srd*localDirection.y() ;
631 normalR *= (-1.0)/rmin;
669 val = -width*0.5*(nReplicas-1)+width*replicaNo;
671 point.setX(point.x()-val);
674 val = -width*0.5*(nReplicas-1)+width*replicaNo;
676 point.setY(point.y()-val);
679 val = -width*0.5*(nReplicas-1)+width*replicaNo;
681 point.setZ(point.z()-val);
684 val = -(offset+width*(replicaNo+0.5));
686 cosv = std::cos(val);
687 sinv = std::sin(val);
688 tmpx = point.x()*cosv-point.y()*sinv;
689 tmpy = point.x()*sinv+point.y()*cosv;
724 val = -width*0.5*(nReplicas-1)+width*replicaNo;
728 val = -width*0.5*(nReplicas-1)+width*replicaNo;
732 val = -width*0.5*(nReplicas-1)+width*replicaNo;
736 val = -(offset+width*(replicaNo+0.5));
755 const G4double currentProposedStepLength,
760 G4bool &calculatedExitNormal,
765 G4int &blockedReplicaNo )
772 G4double ourStep=currentProposedStepLength;
774 G4double sampleStep, sampleSafety, motherStep, motherSafety;
775 G4int localNoDaughters, sampleNo;
783 calculatedExitNormal=
false;
787 if ( exiting&&validExitNormal )
793 blockedExitedVol = *pBlockedPhysical;
813 ourSafety =
std::min( ourSafety, sampleSafety);
815 if ( sampleSafety<ourStep )
823 if ( sampleStep<ourStep )
825 ourStep = sampleStep;
827 validExitNormal = normalOutStc.validConvex;
829 exitNormalStc= normalOutStc;
831 TransformAxis(normalOutStc.exitNormal);
832 calculatedExitNormal=
true;
835 const G4int secondDepth= topDepth;
847 if ( sampleSafety < ourSafety )
849 ourSafety = sampleSafety;
851 if ( sampleSafety < ourStep )
860 if ( sampleStep < ourStep )
862 ourStep = sampleStep;
871 exitNormalStc= normalOutStc;
873 calculatedExitNormal=
true;
886 motherPhysical = history.
GetVolume(depth);
891 motherStep = motherSolid->
DistanceToOut(repPoint,repDirection,
true,
892 &exitConvex,&exitVectorMother);
895 motherNormalStc =
G4ExitNormal( exitVectorMother,
true,
false,
897 calculatedExitNormal=
true;
901 G4bool motherDeterminedStep= (motherStep<ourStep);
903 if( (!exitConvex) && motherDeterminedStep )
906 motherNormalStc=
G4ExitNormal( exitVectorMother,
true,
false,
912 calculatedExitNormal=
true;
914 if( motherDeterminedStep)
919 exitNormalStc= motherNormalStc;
920 exitNormalStc.
exitNormal= globalExitNormalTop;
930 if ( ( (ourStep<minStep) && (sampleSafety<halfkCarTolerance) )
936 if ( motherSafety<ourSafety )
938 ourSafety = motherSafety;
946 std::ostringstream message;
947 message <<
"Point outside volume !" <<
G4endl
948 <<
" Point " << localPoint
949 <<
" is outside current volume " << motherPhysical->
GetName()
952 message <<
" Estimated isotropic distance to solid (distToIn)= "
953 << estDistToSolid <<
G4endl;
959 "Point is far outside Current Volume !" );
964 "Point is a little outside Current Volume.");
972 if( motherDeterminedStep)
974 ourStep = motherStep;
980 if ( calculatedExitNormal )
982 if ( motherDeterminedStep )
989 exitNormalVector= globalToLocalTop.
TransformAxis(exitNormalGlobal);
997 exitNormalVector *= rot->inverse();
1003 validExitNormal =
false;
1007 if ( motherSafety<=ourStep )
1009 if ( motherStep<=ourStep )
1011 ourStep = motherStep;
1013 if ( validExitNormal )
1018 exitNormal *= rot->inverse();
1024 validExitNormal =
false;
1031 G4bool daughterDeterminedStep=
false;
1040 for ( sampleNo=localNoDaughters-1; sampleNo>=0; sampleNo-- )
1042 samplePhysical = repLogical->
GetDaughter(sampleNo);
1043 if ( samplePhysical!=blockedExitedVol )
1055 const G4double sampleSafetyDistance =
1057 if ( sampleSafetyDistance<ourSafety )
1059 ourSafety = sampleSafetyDistance;
1061 if ( sampleSafetyDistance<=ourStep )
1063 sampleDirection = sampleTf.TransformAxis(localDirection);
1064 const G4double sampleStepDistance =
1065 sampleSolid->
DistanceToIn(samplePoint,sampleDirection);
1066 if ( sampleStepDistance<=ourStep )
1068 daughterDeterminedStep=
true;
1070 ourStep = sampleStepDistance;
1073 *pBlockedPhysical = samplePhysical;
1074 blockedReplicaNo = sampleNo;
1076 #ifdef DAUGHTER_NORMAL_ALSO
1079 daughtNormRepCrd = sampleTf.Inverse().TransformAxis(localExitNorm);
1089 intersectionPoint= samplePoint
1090 + sampleStepDistance * sampleDirection;
1091 EInside insideIntPt= sampleSolid->
Inside(intersectionPoint);
1095 std::ostringstream message;
1096 message <<
"Navigator gets conflicting response from Solid."
1098 <<
" Inaccurate DistanceToIn for solid "
1100 <<
" Solid gave DistanceToIn = "
1101 << sampleStepDistance <<
" yet returns " ;
1103 message <<
"-kInside-";
1104 else if ( insideIntPt ==
kOutside )
1105 message <<
"-kOutside-";
1107 message <<
"-kSurface-";
1108 message <<
" for this point !" <<
G4endl
1109 <<
" Point = " << intersectionPoint <<
G4endl;
1111 message <<
" DistanceToIn(p) = "
1115 message <<
" DistanceToOut(p) = "
1119 G4cout.precision(oldcoutPrec);
1128 calculatedExitNormal &= (!daughterDeterminedStep);
1130 #ifdef DAUGHTER_NORMAL_ALSO
1131 if( daughterDeterminedStep )
1138 validExitNormal=
false;
1140 calculatedExitNormal=
true;
1145 newSafety = ourSafety;
1169 G4int localNoDaughters, sampleNo;
1182 if ( sampleSafety<ourSafety )
1184 ourSafety = sampleSafety;
1194 if ( sampleSafety<ourSafety )
1196 ourSafety = sampleSafety;
1204 motherPhysical = history.
GetVolume(depth);
1208 if ( sampleSafety<ourSafety )
1210 ourSafety = sampleSafety;
1216 for ( sampleNo=localNoDaughters-1; sampleNo>=0; sampleNo-- )
1218 samplePhysical = repLogical->
GetDaughter(sampleNo);
1219 if ( samplePhysical!=blockedExitedVol )
1228 const G4double sampleSafetyDistance =
1230 if ( sampleSafetyDistance<ourSafety )
1232 ourSafety = sampleSafetyDistance;
1248 G4bool ¬KnownInside )
const
1253 G4int mdepth, depth, cdepth;
1260 for ( mdepth=cdepth-1; mdepth>=0; mdepth-- )
1274 G4Exception(
"G4ReplicaNavigation::BackLocate()",
"GeomNav0002",
1281 insideCode = motherSolid->
Inside(goodPoint);
1293 notKnownInside =
false;
1298 for ( depth=mdepth+1; depth<cdepth; depth++)
1306 localPoint = goodPoint;
1312 goodPoint = repPoint;
1325 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