49 : fCheck(false), fVerbose(0)
70 const G4int replicaNo,
82 G4double coord, rad2, rmin, tolRMax2, rmax, tolRMin2;
91 coord = std::fabs(localPoint(axis))-width*0.5;
102 if ( localPoint.y()||localPoint.x() )
104 coord = std::fabs(std::atan2(localPoint.y(),localPoint.x()))-width*0.5;
120 rad2 = localPoint.perp2();
121 rmax = (replicaNo+1)*width+offset;
123 tolRMax2 *= tolRMax2;
127 tolRMax2 *= tolRMax2;
128 if ( rad2<=tolRMax2 )
137 if ( replicaNo||offset )
141 tolRMin2 *= tolRMin2;
145 tolRMin2 *= tolRMin2;
146 if ( rad2>=tolRMin2 )
163 G4Exception(
"G4ReplicaNavigation::Inside()",
"GeomNav0002",
176 const G4int replicaNo,
196 coord = localPoint(axis);
197 safe1 = width*0.5-coord;
198 safe2 = width*0.5+coord;
199 safety = (safe1<=safe2) ? safe1 : safe2;
202 if ( localPoint.y()<=0 )
204 safety = localPoint.x()*std::sin(width*0.5)
205 + localPoint.y()*std::cos(width*0.5);
209 safety = localPoint.x()*std::sin(width*0.5)
210 - localPoint.y()*std::cos(width*0.5);
214 rho = localPoint.perp();
215 rmax = width*(replicaNo+1)+offset;
216 if ( replicaNo||offset )
221 safety = (safe1<=safe2) ? safe1 : safe2;
229 G4Exception(
"G4ReplicaNavigation::DistanceToOut()",
"GeomNav0002",
242 const G4int replicaNo,
272 coord = localPoint(axis);
273 Comp = localDirection(axis);
276 lindist = width*0.5-coord;
277 Dist = (lindist>0) ? lindist/Comp : 0;
282 lindist = width*0.5+coord;
283 Dist = (lindist>0) ? -lindist/Comp : 0;
291 candidateNormal.
exitNormal = ( signC * VecCartAxes[axis]);
295 (Comp>0) ? SideCartAxesPlus[axis] : SideCartAxesMinus[axis];
303 replicaNo,candidateNormal);
307 G4Exception(
"G4ReplicaNavigation::DistanceToOut()",
"GeomNav0002",
312 arExitNormal= candidateNormal;
330 G4double sinSPhi= -2.0, cosSPhi= -2.0;
331 G4double pDistS, pDistE, compS, compE, Dist, dist2, yi;
335 if ( (localPoint.x()!=0.0) || (localPoint.y()!=0.0) )
337 sinSPhi = std::sin(-width*0.5);
338 cosSPhi = std::cos(width*0.5);
342 pDistS = localPoint.x()*sinSPhi-localPoint.y()*cosSPhi;
344 pDistE = localPoint.x()*sinSPhi+localPoint.y()*cosSPhi;
349 compS = -sinSPhi*localDirection.x()+cosSPhi*localDirection.y();
350 compE = -sinSPhi*localDirection.x()-cosSPhi*localDirection.y();
352 if ( (pDistS<=0)&&(pDistE<=0) )
358 dist2 = pDistS/compS;
359 yi = localPoint.y()+dist2*localDirection.y();
379 dist2 = pDistE/compE;
385 yi = localPoint.y()+dist2*localDirection.y();
399 else if ( (pDistS>=0)&&(pDistE>=0) )
404 Dist = ((compS>=0)&&(compE>=0)) ?
kInfinity : 0;
406 else if ( (pDistS>0)&&(pDistE<0) )
412 dist2 = pDistE/compE;
413 yi = localPoint.y()+dist2*localDirection.y();
435 dist2 = pDistS/compS;
436 yi = localPoint.y()+dist2*localDirection.y();
462 if( (std::fabs(localDirection.phi())<=width*0.5) )
500 const G4int replicaNo,
503 G4double rmin, rmax, t1, t2, t3, deltaR;
523 rmin = replicaNo*width+offset;
524 rmax = (replicaNo+1)*width+offset;
526 t1 = 1.0-localDirection.z()*localDirection.z();
527 t2 = localPoint.x()*localDirection.x()+localPoint.y()*localDirection.y();
528 t3 = localPoint.x()*localPoint.x()+localPoint.y()*localPoint.y();
538 deltaR = t3-rmax*rmax;
547 srd = -b+std::sqrt(b*b-c);
565 deltaR = t3-rmin*rmin;
584 deltaR = t3-rmax*rmax;
586 srd = -b+std::sqrt(b*b-c);
594 deltaR = t3-rmax*rmax;
597 srd = -b+std::sqrt(b*b-c);
613 xi = localPoint.x() + srd*localDirection.x() ;
614 yi = localPoint.y() + srd*localDirection.y() ;
620 normalR *= (-1.0)/rmin;
658 val = -width*0.5*(nReplicas-1)+width*replicaNo;
660 point.setX(point.x()-val);
663 val = -width*0.5*(nReplicas-1)+width*replicaNo;
665 point.setY(point.y()-val);
668 val = -width*0.5*(nReplicas-1)+width*replicaNo;
670 point.setZ(point.z()-val);
673 val = -(offset+width*(replicaNo+0.5));
675 cosv = std::cos(val);
676 sinv = std::sin(val);
677 tmpx = point.x()*cosv-point.y()*sinv;
678 tmpy = point.x()*sinv+point.y()*cosv;
713 val = -width*0.5*(nReplicas-1)+width*replicaNo;
717 val = -width*0.5*(nReplicas-1)+width*replicaNo;
721 val = -width*0.5*(nReplicas-1)+width*replicaNo;
725 val = -(offset+width*(replicaNo+0.5));
744 const G4double currentProposedStepLength,
749 G4bool &calculatedExitNormal,
754 G4int &blockedReplicaNo )
761 G4double ourStep=currentProposedStepLength;
763 G4double sampleStep, sampleSafety, motherStep, motherSafety;
764 G4int localNoDaughters, sampleNo;
769 calculatedExitNormal=
false;
773 if ( exiting&&validExitNormal )
779 blockedExitedVol = *pBlockedPhysical;
799 if ( sampleSafety<ourSafety )
801 ourSafety = sampleSafety;
803 if ( sampleSafety<ourStep )
811 if ( sampleStep<ourStep )
813 ourStep = sampleStep;
815 validExitNormal = normalOutStc.validConvex;
817 exitNormalStc= normalOutStc;
819 TransformAxis(normalOutStc.exitNormal);
820 calculatedExitNormal=
true;
823 const G4int secondDepth= topDepth;
835 if ( sampleSafety < ourSafety )
837 ourSafety = sampleSafety;
839 if ( sampleSafety < ourStep )
847 if ( sampleStep < ourStep )
849 ourStep = sampleStep;
858 exitNormalStc= normalOutStc;
860 calculatedExitNormal=
true;
873 motherPhysical = history.
GetVolume(depth);
878 motherStep = motherSolid->
DistanceToOut(repPoint,repDirection,
true,
879 &exitConvex,&exitVectorMother);
882 motherNormalStc =
G4ExitNormal( exitVectorMother,
true,
false,
884 calculatedExitNormal=
true;
888 G4bool motherDeterminedStep= (motherStep<ourStep);
890 if( (!exitConvex) && motherDeterminedStep )
893 motherNormalStc=
G4ExitNormal( exitVectorMother,
true,
false,
899 calculatedExitNormal=
true;
901 if( motherDeterminedStep)
906 exitNormalStc= motherNormalStc;
907 exitNormalStc.
exitNormal= globalExitNormalTop;
923 if ( motherSafety<ourSafety )
925 ourSafety = motherSafety;
933 std::ostringstream message;
934 message <<
"Point outside volume !" <<
G4endl
935 <<
" Point " << localPoint
936 <<
" is outside current volume " << motherPhysical->
GetName()
939 message <<
" Estimated isotropic distance to solid (distToIn)= "
940 << estDistToSolid <<
G4endl;
946 "Point is far outside Current Volume !" );
951 "Point is a little outside Current Volume.");
959 if( motherDeterminedStep)
961 ourStep = motherStep;
967 if ( calculatedExitNormal )
969 if ( motherDeterminedStep )
974 exitNormalVector= globalToLocalTop.
TransformAxis(exitNormalGlobal);
982 exitNormalVector *= rot->inverse();
988 validExitNormal =
false;
992 if ( motherSafety<=ourStep )
994 if ( motherStep<=ourStep )
996 ourStep = motherStep;
998 if ( validExitNormal )
1003 exitNormal *= rot->inverse();
1009 validExitNormal =
false;
1016 G4bool daughterDeterminedStep=
false;
1025 for ( sampleNo=localNoDaughters-1; sampleNo>=0; sampleNo-- )
1027 samplePhysical = repLogical->
GetDaughter(sampleNo);
1028 if ( samplePhysical!=blockedExitedVol )
1040 const G4double sampleSafetyDistance =
1042 if ( sampleSafetyDistance<ourSafety )
1044 ourSafety = sampleSafetyDistance;
1046 if ( sampleSafetyDistance<=ourStep )
1048 sampleDirection = sampleTf.TransformAxis(localDirection);
1049 const G4double sampleStepDistance =
1050 sampleSolid->
DistanceToIn(samplePoint,sampleDirection);
1051 if ( sampleStepDistance<=ourStep )
1053 daughterDeterminedStep=
true;
1055 ourStep = sampleStepDistance;
1058 *pBlockedPhysical = samplePhysical;
1059 blockedReplicaNo = -1;
1061 #ifdef DAUGHTER_NORMAL_ALSO
1064 daughtNormRepCrd = sampleTf.Inverse().TransformAxis(localExitNorm);
1074 intersectionPoint= samplePoint
1075 + sampleStepDistance * sampleDirection;
1076 EInside insideIntPt= sampleSolid->
Inside(intersectionPoint);
1080 std::ostringstream message;
1081 message <<
"Navigator gets conflicting response from Solid."
1083 <<
" Inaccurate DistanceToIn for solid "
1085 <<
" Solid gave DistanceToIn = "
1086 << sampleStepDistance <<
" yet returns " ;
1088 message <<
"-kInside-";
1089 else if ( insideIntPt ==
kOutside )
1090 message <<
"-kOutside-";
1092 message <<
"-kSurface-";
1093 message <<
" for this point !" <<
G4endl
1094 <<
" Point = " << intersectionPoint <<
G4endl;
1096 message <<
" DistanceToIn(p) = "
1100 message <<
" DistanceToOut(p) = "
1104 G4cout.precision(oldcoutPrec);
1113 calculatedExitNormal &= (!daughterDeterminedStep);
1115 #ifdef DAUGHTER_NORMAL_ALSO
1116 if( daughterDeterminedStep )
1123 validExitNormal=
false;
1125 calculatedExitNormal=
true;
1130 newSafety = ourSafety;
1154 G4int localNoDaughters, sampleNo;
1167 if ( sampleSafety<ourSafety )
1169 ourSafety = sampleSafety;
1179 if ( sampleSafety<ourSafety )
1181 ourSafety = sampleSafety;
1189 motherPhysical = history.
GetVolume(depth);
1193 if ( sampleSafety<ourSafety )
1195 ourSafety = sampleSafety;
1201 for ( sampleNo=localNoDaughters-1; sampleNo>=0; sampleNo-- )
1203 samplePhysical = repLogical->
GetDaughter(sampleNo);
1204 if ( samplePhysical!=blockedExitedVol )
1213 const G4double sampleSafetyDistance =
1215 if ( sampleSafetyDistance<ourSafety )
1217 ourSafety = sampleSafetyDistance;
1233 G4bool ¬KnownInside )
const
1238 G4int mdepth, depth, cdepth;
1245 for ( mdepth=cdepth-1; mdepth>=0; mdepth-- )
1259 G4Exception(
"G4ReplicaNavigation::BackLocate()",
"GeomNav0002",
1266 insideCode = motherSolid->
Inside(goodPoint);
1278 notKnownInside =
false;
1283 for ( depth=mdepth+1; depth<cdepth; depth++)
1291 localPoint = goodPoint;
1297 goodPoint = repPoint;
1310 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
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