58 : fCheck(false), fVerbose(0)
83 const G4int replicaNo,
95 G4double coord, rad2, rmin, tolRMax2, rmax, tolRMin2;
104 coord = std::fabs(localPoint(axis))-width*0.5;
115 if ( localPoint.
y()||localPoint.
x() )
117 coord = std::fabs(std::atan2(localPoint.
y(),localPoint.
x()))-width*0.5;
133 rad2 = localPoint.
perp2();
134 rmax = (replicaNo+1)*width+offset;
136 tolRMax2 *= tolRMax2;
140 tolRMax2 *= tolRMax2;
141 if ( rad2<=tolRMax2 )
150 if ( replicaNo||offset )
154 tolRMin2 *= tolRMin2;
158 tolRMin2 *= tolRMin2;
159 if ( rad2>=tolRMin2 )
176 G4Exception(
"G4ReplicaNavigation::Inside()",
"GeomNav0002",
189 const G4int replicaNo,
209 coord = localPoint(axis);
210 safe1 = width*0.5-coord;
211 safe2 = width*0.5+coord;
212 safety = (safe1<=safe2) ? safe1 : safe2;
215 if ( localPoint.
y()<=0 )
217 safety = localPoint.
x()*std::sin(width*0.5)
218 + localPoint.
y()*std::cos(width*0.5);
222 safety = localPoint.
x()*std::sin(width*0.5)
223 - localPoint.
y()*std::cos(width*0.5);
227 rho = localPoint.
perp();
228 rmax = width*(replicaNo+1)+offset;
229 if ( replicaNo||offset )
234 safety = (safe1<=safe2) ? safe1 : safe2;
242 G4Exception(
"G4ReplicaNavigation::DistanceToOut()",
"GeomNav0002",
255 const G4int replicaNo,
278 coord = localPoint(axis);
279 Comp = localDirection(axis);
282 lindist = width*0.5-coord;
283 Dist = (lindist>0) ? lindist/Comp : 0;
288 lindist = width*0.5+coord;
289 Dist = (lindist>0) ? -lindist/Comp : 0;
297 candidateNormal.
exitNormal = ( signC * VecCartAxes[axis]);
301 (Comp>0) ? SideCartAxesPlus[axis] : SideCartAxesMinus[axis];
309 replicaNo,candidateNormal);
313 G4Exception(
"G4ReplicaNavigation::DistanceToOut()",
"GeomNav0002",
318 arExitNormal= candidateNormal;
336 G4double sinSPhi= -2.0, cosSPhi= -2.0;
337 G4double pDistS, pDistE, compS, compE, Dist, dist2, yi;
341 if ( (localPoint.
x()!=0.0) || (localPoint.
y()!=0.0) )
343 sinSPhi = std::sin(-width*0.5);
344 cosSPhi = std::cos(width*0.5);
348 pDistS = localPoint.
x()*sinSPhi-localPoint.
y()*cosSPhi;
350 pDistE = localPoint.
x()*sinSPhi+localPoint.
y()*cosSPhi;
355 compS = -sinSPhi*localDirection.
x()+cosSPhi*localDirection.
y();
356 compE = -sinSPhi*localDirection.
x()-cosSPhi*localDirection.
y();
364 dist2 = pDistS/compS;
365 yi = localPoint.
y()+dist2*localDirection.
y();
385 dist2 = pDistE/compE;
391 yi = localPoint.
y()+dist2*localDirection.
y();
410 Dist = ((compS>=0)&&(compE>=0)) ?
kInfinity : 0;
418 dist2 = pDistE/compE;
419 yi = localPoint.
y()+dist2*localDirection.
y();
441 dist2 = pDistS/compS;
442 yi = localPoint.
y()+dist2*localDirection.
y();
468 if( (std::fabs(localDirection.
phi())<=width*0.5) )
506 const G4int replicaNo,
529 rmin = replicaNo*width+offset;
530 rmax = (replicaNo+1)*width+offset;
532 t1 = 1.0-localDirection.
z()*localDirection.
z();
533 t2 = localPoint.
x()*localDirection.
x()+localPoint.
y()*localDirection.
y();
534 t3 = localPoint.
x()*localPoint.
x()+localPoint.
y()*localPoint.
y();
544 deltaR = t3-rmax*rmax;
553 srd = -b+std::sqrt(b*b-c);
571 deltaR = t3-rmin*rmin;
590 deltaR = t3-rmax*rmax;
593 srd = (d2 < 0.) ? 0.0 : -b+std::sqrt(d2);
601 deltaR = t3-rmax*rmax;
605 srd = (d2 < 0.) ? 0.0 : -b+std::sqrt(d2);
622 xi = localPoint.
x() + srd*localDirection.
x();
623 yi = localPoint.
y() + srd*localDirection.
y();
632 normalR *= (-1.0)/rmin;
672 val = -width*0.5*(nReplicas-1)+width*replicaNo;
674 point.
setX(point.
x()-val);
677 val = -width*0.5*(nReplicas-1)+width*replicaNo;
679 point.
setY(point.
y()-val);
682 val = -width*0.5*(nReplicas-1)+width*replicaNo;
684 point.
setZ(point.
z()-val);
687 val = -(offset+width*(replicaNo+0.5));
689 cosv = std::cos(val);
690 sinv = std::sin(val);
691 tmpx = point.
x()*cosv-point.
y()*sinv;
692 tmpy = point.
x()*sinv+point.
y()*cosv;
727 val = -width*0.5*(nReplicas-1)+width*replicaNo;
731 val = -width*0.5*(nReplicas-1)+width*replicaNo;
735 val = -width*0.5*(nReplicas-1)+width*replicaNo;
739 val = -(offset+width*(replicaNo+0.5));
758 const G4double currentProposedStepLength,
763 G4bool &calculatedExitNormal,
768 G4int &blockedReplicaNo )
775 G4double ourStep=currentProposedStepLength;
777 G4double sampleStep, sampleSafety, motherStep, motherSafety;
778 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;
832 calculatedExitNormal=
true;
835 const G4int secondDepth= topDepth;
848 if ( sampleSafety < ourSafety )
850 ourSafety = sampleSafety;
852 if ( sampleSafety < ourStep )
861 if ( sampleStep < ourStep )
863 ourStep = sampleStep;
872 exitNormalStc= normalOutStc;
874 calculatedExitNormal=
true;
887 motherPhysical = history.
GetVolume(depth);
892 motherStep = motherSolid->
DistanceToOut(repPoint,repDirection,
true,
893 &exitConvex,&exitVectorMother);
896 motherNormalStc =
G4ExitNormal( exitVectorMother,
true,
false,
898 calculatedExitNormal=
true;
902 G4bool motherDeterminedStep= (motherStep<ourStep);
904 if( (!exitConvex) && motherDeterminedStep )
907 motherNormalStc=
G4ExitNormal( exitVectorMother,
true,
false,
913 calculatedExitNormal=
true;
915 if( motherDeterminedStep)
920 exitNormalStc= motherNormalStc;
921 exitNormalStc.
exitNormal= globalExitNormalTop;
937 if ( motherSafety<ourSafety )
939 ourSafety = motherSafety;
948 message <<
"Point outside volume !" <<
G4endl
949 <<
" Point " << localPoint
950 <<
" is outside current volume " << motherPhysical->
GetName()
953 message <<
" Estimated isotropic distance to solid (distToIn)= "
954 << estDistToSolid <<
G4endl;
960 "Point is far outside Current Volume !" );
965 "Point is a little outside Current Volume.");
973 if( motherDeterminedStep)
975 ourStep = motherStep;
981 if ( calculatedExitNormal )
983 if ( motherDeterminedStep )
990 exitNormalVector= globalToLocalTop.
TransformAxis(exitNormalGlobal);
998 exitNormalVector *= rot->
inverse();
1004 validExitNormal =
false;
1008 if ( motherSafety<=ourStep )
1010 if ( motherStep<=ourStep )
1012 ourStep = motherStep;
1014 if ( validExitNormal )
1025 validExitNormal =
false;
1032 G4bool daughterDeterminedStep=
false;
1041 for ( sampleNo=localNoDaughters-1; sampleNo>=0; sampleNo-- )
1043 samplePhysical = repLogical->
GetDaughter(sampleNo);
1044 if ( samplePhysical!=blockedExitedVol )
1053 sampleTf.TransformPoint(localPoint);
1056 const G4double sampleSafetyDistance =
1058 if ( sampleSafetyDistance<ourSafety )
1060 ourSafety = sampleSafetyDistance;
1062 if ( sampleSafetyDistance<=ourStep )
1064 sampleDirection = sampleTf.TransformAxis(localDirection);
1065 const G4double sampleStepDistance =
1066 sampleSolid->
DistanceToIn(samplePoint,sampleDirection);
1067 if ( sampleStepDistance<=ourStep )
1069 daughterDeterminedStep=
true;
1071 ourStep = sampleStepDistance;
1074 *pBlockedPhysical = samplePhysical;
1075 blockedReplicaNo = sampleNo;
1077 #ifdef DAUGHTER_NORMAL_ALSO
1080 daughtNormRepCrd = sampleTf.InverseTransformAxis(localExitNorm);
1090 intersectionPoint= samplePoint
1091 + sampleStepDistance * sampleDirection;
1092 EInside insideIntPt= sampleSolid->
Inside(intersectionPoint);
1097 message <<
"Navigator gets conflicting response from Solid."
1099 <<
" Inaccurate DistanceToIn for solid "
1101 <<
" Solid gave DistanceToIn = "
1102 << sampleStepDistance <<
" yet returns " ;
1104 message <<
"-kInside-";
1105 else if ( insideIntPt ==
kOutside )
1106 message <<
"-kOutside-";
1108 message <<
"-kSurface-";
1109 message <<
" for this point !" <<
G4endl
1110 <<
" Point = " << intersectionPoint <<
G4endl;
1112 message <<
" DistanceToIn(p) = "
1116 message <<
" DistanceToOut(p) = "
1120 G4cout.precision(oldcoutPrec);
1129 calculatedExitNormal &= (!daughterDeterminedStep);
1131 #ifdef DAUGHTER_NORMAL_ALSO
1132 if( daughterDeterminedStep )
1139 validExitNormal =
false;
1141 calculatedExitNormal=
true;
1146 newSafety = ourSafety;
1170 G4int localNoDaughters, sampleNo;
1183 if ( sampleSafety<ourSafety )
1185 ourSafety = sampleSafety;
1197 if ( sampleSafety<ourSafety )
1199 ourSafety = sampleSafety;
1207 motherPhysical = history.
GetVolume(depth);
1211 if ( sampleSafety<ourSafety )
1213 ourSafety = sampleSafety;
1219 for ( sampleNo=localNoDaughters-1; sampleNo>=0; sampleNo-- )
1221 samplePhysical = repLogical->
GetDaughter(sampleNo);
1222 if ( samplePhysical!=blockedExitedVol )
1228 sampleTf.TransformPoint(localPoint);
1231 const G4double sampleSafetyDistance =
1233 if ( sampleSafetyDistance<ourSafety )
1235 ourSafety = sampleSafetyDistance;
1251 G4bool ¬KnownInside )
const
1256 G4int mdepth, depth, cdepth;
1263 for ( mdepth=cdepth-1; mdepth>=0; mdepth-- )
1277 G4Exception(
"G4ReplicaNavigation::BackLocate()",
"GeomNav0002",
1284 insideCode = motherSolid->
Inside(goodPoint);
1296 notKnownInside =
false;
1301 for ( depth=mdepth+1; depth<cdepth; depth++)
1309 localPoint = goodPoint;
1315 goodPoint = repPoint;
1328 localPoint = goodPoint;
EInside BackLocate(G4NavigationHistory &history, const G4ThreeVector &globalPoint, G4ThreeVector &localPoint, const G4bool &exiting, G4bool ¬KnownInside) const
void SetPhiTransformation(const G4double ang, G4VPhysicalVolume *pVol=0) const
CLHEP::Hep3Vector G4ThreeVector
G4LogicalVolume * GetLogicalVolume() const
static const G4double kInfinity
void message(RunManager *runmanager)
G4int GetNoDaughters() const
G4double GetAngularTolerance() const
double dot(const Hep3Vector &) const
void SetTranslation(const G4ThreeVector &v)
G4double GetRadialTolerance() const
virtual G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const =0
const G4ThreeVector & GetTranslation() const
virtual G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const =0
void ComputeTransformation(const G4int replicaNo, G4VPhysicalVolume *pVol, G4ThreeVector &point) const
virtual EInside Inside(const G4ThreeVector &p) const =0
G4double halfkRadTolerance
G4double DistanceToOutPhi(const G4ThreeVector &localPoint, const G4ThreeVector &localDirection, const G4double width, G4ExitNormal &foundNormal) const
const G4AffineTransform & GetTopTransform() 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)
const G4AffineTransform & GetTransform(G4int n) const
static const double kMinExitingNormalCosine
G4VPhysicalVolume * GetVolume(G4int n) const
G4VPhysicalVolume * GetTopVolume() const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
G4int GetReplicaNo(G4int n) const
G4double ComputeSafety(const G4ThreeVector &globalPoint, const G4ThreeVector &localPoint, G4NavigationHistory &history, const G4double pProposedMaxLength=DBL_MAX)
virtual void GetReplicationData(EAxis &axis, G4int &nReplicas, G4double &width, G4double &offset, G4bool &consuming) const =0
G4double DistanceToOut(const G4VPhysicalVolume *pVol, const G4int replicaNo, const G4ThreeVector &localPoint) const
EInside Inside(const G4VPhysicalVolume *pVol, const G4int replicaNo, const G4ThreeVector &localPoint) const
G4double DistanceToOutRad(const G4ThreeVector &localPoint, const G4ThreeVector &localDirection, const G4double width, const G4double offset, const G4int replicaNo, G4ExitNormal &foundNormal) const
G4double halfkAngTolerance
G4GLOB_DLL std::ostream G4cout
EVolume GetVolumeType(G4int n) const
virtual G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=false, G4bool *validNorm=0, G4ThreeVector *n=0) const =0
G4VSolid * GetSolid() const
static G4GeometryTolerance * GetInstance()
G4double halfkCarTolerance
const G4RotationMatrix * GetRotation() const
G4int GetTopReplicaNo() const
HepRotation inverse() const
G4double GetSurfaceTolerance() const
const G4String & GetName() const
G4VPhysicalVolume * GetDaughter(const G4int i) const
T min(const T t1, const T t2)
brief Return the smallest of the two arguments