57 #define fHistory fpNavigatorState->fHistory
58 #define fLastTriedStepComputation fpNavigatorState->fLastTriedStepComputation
59 #define fEnteredDaughter fpNavigatorState->fEnteredDaughter
60 #define fExitedMother fpNavigatorState->fExitedMother
61 #define fWasLimitedByGeometry fpNavigatorState->fWasLimitedByGeometry
62 #define fStepEndPoint fpNavigatorState->fStepEndPoint
63 #define fLastStepEndPointLocal fpNavigatorState->fLastStepEndPointLocal
64 #define fPushed fpNavigatorState->fPushed
65 #define fLastTriedStepComputation fpNavigatorState->fLastTriedStepComputation
66 #define fEntering fpNavigatorState->fEntering
67 #define fExiting fpNavigatorState->fExiting
68 #define fBlockedPhysicalVolume fpNavigatorState->fBlockedPhysicalVolume
69 #define fBlockedReplicaNo fpNavigatorState->fBlockedReplicaNo
70 #define fLastLocatedPointLocal fpNavigatorState->fLastLocatedPointLocal
71 #define fLocatedOutsideWorld fpNavigatorState->fLocatedOutsideWorld
72 #define fValidExitNormal fpNavigatorState->fValidExitNormal
73 #define fExitNormal fpNavigatorState->fExitNormal
74 #define fGrandMotherExitNormal fpNavigatorState->fGrandMotherExitNormal
75 #define fChangedGrandMotherRefFrame fpNavigatorState->fChangedGrandMotherRefFrame
76 #define fExitNormalGlobalFrame fpNavigatorState->fExitNormalGlobalFrame
77 #define fCalculatedExitNormal fpNavigatorState->fCalculatedExitNormal
78 #define fLastStepWasZero fpNavigatorState->fLastStepWasZero
79 #define fLocatedOnEdge fpNavigatorState->fLocatedOnEdge
80 #define fNumberZeroSteps fpNavigatorState->fNumberZeroSteps
81 #define fPreviousSftOrigin fpNavigatorState->fPreviousSftOrigin
82 #define fPreviousSafety fpNavigatorState->fPreviousSafety
101 fTopPhysical(0), fCheck(
false),
160 const G4bool relativeSearch,
161 const G4bool ignoreDirection )
164 G4bool notKnownContained=
true, noResult;
176 if( considerDirection && pGlobalDirection != 0 )
178 globalDirection=*pGlobalDirection;
185 static int nCalls = 0;
187 G4cout <<
"*** G4ITNavigator2::LocateGlobalPointAndSetup: ***" <<
G4endl;
188 G4cout <<
" Called with arguments: " << G4endl
189 <<
" Globalpoint = " << globalPoint << G4endl
190 <<
" RelativeSearch = " << relativeSearch <<
G4endl;
196 G4cout.precision(oldcoutPrec);
202 G4int noLevelsExited=0 ;
203 G4int noLevelsEntered= 0;
205 if ( !relativeSearch )
302 G4Exception(
"G4ITNavigator2::LocateGlobalPointAndSetup()",
304 "Not applicable for external volumes.");
309 localPoint =
fHistory.GetTopTransform().TransformPoint(globalPoint);
310 notKnownContained =
false;
332 while (notKnownContained)
336 targetSolid =
fHistory.GetTopVolume()->GetLogicalVolume()->GetSolid();
337 localPoint =
fHistory.GetTopTransform().TransformPoint(globalPoint);
338 insideCode = targetSolid->
Inside(localPoint);
342 G4String solidResponse =
"-kInside-";
344 solidResponse =
"-kOutside-";
346 solidResponse =
"-kSurface-";
347 G4cout <<
"*** G4ITNavigator2::LocateGlobalPointAndSetup(): ***" <<
G4endl
348 <<
" Invoked Inside() for solid: " << targetSolid->
GetName()
349 <<
". Solid replied: " << solidResponse <<
G4endl
350 <<
" For local point p: " << localPoint <<
G4endl;
375 if( noLevelsExited > 1 )
399 if( (!
fExiting)&&considerDirection )
404 G4bool directionExiting =
false;
406 fHistory.GetTopTransform().TransformAxis(globalDirection);
411 localPoint=
fHistory.GetTopTransform().TransformPoint(globalPoint);
415 directionExiting = normal.
dot(localDirection) > 0.0;
416 isExiting = isExiting || directionExiting;
432 if( noLevelsExited > 1 )
455 notKnownContained=
false;
460 notKnownContained=
false;
480 targetPhysical =
fHistory.GetTopVolume();
481 if (!targetPhysical) {
break; }
539 G4Exception(
"G4ITNavigator2::LocateGlobalPointAndSetup()",
541 "Not applicable for external volumes.");
579 #ifdef G4DEBUG_NAVIGATION
583 G4cout <<
"*** G4ITNavigator2::LocateGlobalPointAndSetup() ***" <<
G4endl;
598 if (targetPhysical) { curPhysVol_Name = targetPhysical->
GetName(); }
599 G4cout <<
" Return value = new volume = " << curPhysVol_Name <<
G4endl;
604 G4cout <<
"Upon exiting LocateGlobalPointAndSetup():" <<
G4endl;
607 G4cout.precision(oldcoutPrec);
613 return targetPhysical;
634 #ifdef G4DEBUG_NAVIGATION
644 #ifdef G4DEBUG_NAVIGATION
647 G4cout <<
"Entering LocateGlobalWithinVolume(): History = " <<
G4endl;
681 G4Exception(
"G4ITNavigator2::LocateGlobalPointWithinVolume()",
683 "Not applicable for replicated volumes.");
686 G4Exception(
"G4ITNavigator2::LocateGlobalPointWithinVolume()",
688 "Not applicable for external volumes.");
713 exceptionDescription <<
"The navigator state is NULL. ";
714 exceptionDescription <<
"Either NewNavigatorStateAndLocate was not called ";
715 exceptionDescription <<
"or the provided navigator state was already NULL.";
717 G4Exception(
"G4ITNavigator::CheckNavigatorStateIsValid",
747 exceptionDescription <<
"No World Volume";
764 exceptionDescription <<
"No World Volume";
784 exceptionDescription <<
"No World Volume";
786 G4Exception(
"G4ITNavigator::NewNavigatorStateAndLocate",
926 const G4double pCurrentProposedStepLength,
953 G4cout <<
"*** G4ITNavigator2::ComputeStep: ***" <<
G4endl;
955 <<
" - Proposed step length = " << pCurrentProposedStepLength
957 #ifdef G4DEBUG_NAVIGATION
960 G4cout <<
" Called with the arguments: " << G4endl
961 <<
" Globalpoint = " << std::setw(25) << pGlobalpoint << G4endl
962 <<
" Direction = " << std::setw(25) << pDirection <<
G4endl;
979 G4double moveLenSq = (newLocalPoint-oldLocalPoint).mag2();
1002 pCurrentProposedStepLength,
1019 pCurrentProposedStepLength,
1043 if(
fHistory.GetTopVolume()->GetRegularStructureId() == 0 )
1047 "Point is relocated in voxels, while it should be outside!");
1050 pCurrentProposedStepLength,
1065 pCurrentProposedStepLength,
1084 pCurrentProposedStepLength,
1098 pCurrentProposedStepLength,
1110 G4Exception(
"G4ITNavigator2::ComputeStep()",
"GeomNav0001",
1114 G4Exception(
"G4ITNavigator2::ComputeStep()",
"GeomNav0001",
1125 G4bool calculatedExitNormal;
1130 pCurrentProposedStepLength,
1134 calculatedExitNormal,
1171 #ifdef G4DEBUG_NAVIGATION
1174 G4cout <<
"G4ITNavigator2::ComputeStep(): another zero step, # "
1176 <<
" at " << pGlobalpoint
1177 <<
" in volume " << motherPhysical->
GetName()
1178 <<
" nav-comp-step calls # " << sNavCScalls
1191 message <<
"Track stuck or not moving." <<
G4endl
1192 <<
" Track stuck, not moving for "
1194 <<
" in volume -" << motherPhysical->
GetName()
1195 <<
"- at point " << pGlobalpoint <<
G4endl
1196 <<
" direction: " << pDirection <<
"." <<
G4endl
1197 <<
" Potential geometry or navigation problem !"
1199 <<
" Trying pushing it of " << Step <<
" mm ...";
1200 G4Exception(
"G4ITNavigator2::ComputeStep()",
"GeomNav1002",
1201 JustWarning, message,
"Potential overlap in geometry!");
1211 message <<
"Stuck Track: potential geometry or navigation problem."
1213 <<
" Track stuck, not moving for "
1215 <<
" in volume -" << motherPhysical->
GetName()
1216 <<
"- at point " << pGlobalpoint <<
G4endl
1217 <<
" direction: " << pDirection <<
".";
1219 G4Exception(
"G4ITNavigator2::ComputeStep()",
"GeomNav0003",
1232 +
std::min(Step,pCurrentProposedStepLength) * pDirection;
1237 #ifdef G4DEBUG_NAVIGATION
1240 G4cout <<
" At G4Nav CompStep End - if(exiting) - fExiting= " <<
fExiting
1301 #ifdef G4DEBUG_NAVIGATION
1304 desc <<
"Problem in ComputeStep: Replica Navigation did not provide"
1305 <<
" valid exit Normal. " <<
G4endl;
1306 desc <<
" Do not know how calculate it in this case." <<
G4endl;
1307 desc <<
" Location = " << finalLocalPoint <<
G4endl;
1308 desc <<
" Volume name = " << motherPhysical->
GetName()
1310 G4Exception(
"G4ITNavigator2::ComputeStep()",
"GeomNav0003",
1311 JustWarning, desc,
"Normal not available for exiting.");
1324 fHistory.GetTransform(depth-1).Inverse();
1355 G4cout <<
" Returned step= " << Step;
1359 G4cout <<
" Requested step= " << pCurrentProposedStepLength ;
1377 const G4double pCurrentProposedStepLength,
1391 pCurrentProposedStepLength,
1460 for ( i=1; i<=cdepth; i++ )
1463 switch (
fHistory.GetVolumeType(i) )
1474 replicaNo =
fHistory.GetReplicaNo(i);
1499 ComputeMaterial(replicaNo, current, pTouchable) );
1507 "Not applicable for external volumes.");
1536 if( candidateLogical )
1558 currentSolid= candidateLogical->
GetSolid();
1559 inSideIt = currentSolid->
Inside(daughterPointOwnLocal);
1560 onSurface = (inSideIt ==
kSurface);
1565 safety = (currentSolid->
DistanceToIn(daughterPointOwnLocal));
1568 else if (inSideIt ==
kInside )
1570 safety = (currentSolid->
DistanceToOut(daughterPointOwnLocal));
1577 nextSolidExitNormal =
1582 ExitNormal = -nextSolidExitNormal;
1591 message <<
"Point not on surface ! " <<
G4endl
1593 << daughterPointOwnLocal <<
G4endl
1594 <<
" Physical volume = "
1596 <<
" Logical volume = "
1598 <<
" Solid = " << currentSolid->
GetName()
1601 << *currentSolid <<
G4endl;
1604 message <<
"Point is Outside. " << G4endl
1605 <<
" Safety (from outside) = " << safety <<
G4endl;
1609 message <<
"Point is Inside. " << G4endl
1610 <<
" Safety (from inside) = " << safety <<
G4endl;
1612 G4Exception(
"G4ITNavigator2::GetLocalExitNormal()",
"GeomNav1001",
1630 G4Exception(
"G4ITNavigator2::GetLocalExitNormal()",
1632 "Incorrect call to GetLocalSurfaceNormal." );
1645 desc <<
" Parameters of solid: " << *daughterSolid
1647 G4Exception(
"G4ITNavigator2::GetLocalExitNormal()",
1649 "Surface Normal returned by Solid is not a Unit Vector." );
1667 message <<
"Function called when *NOT* at a Boundary." <<
G4endl;
1668 G4Exception(
"G4ITNavigator2::GetLocalExitNormal()",
1684 G4int enteringReplicaNo,
1688 switch (enteringVolumeType)
1693 G4Exception(
"G4ITNavigator2::GetMotherToDaughterTransform()",
1695 "Method NOT Implemented yet for replica volumes.");
1703 pParam->
ComputeSolid(enteringReplicaNo, pEnteringPhysVol);
1717 G4Exception(
"G4ITNavigator2::GetMotherToDaughterTransform()",
1719 "Not applicable for external volumes.");
1744 #ifdef G4DEBUG_NAVIGATION
1749 G4ThreeVector ExpectedBoundaryPointLocal;
1752 ExpectedBoundaryPointLocal =
1772 G4bool* pNormalCalculated)
1796 if( std::fabs ( normMag2 - 1.0 ) <
perMillion )
1798 *pNormalCalculated =
true;
1805 message <<
" ERROR> Expected normal-global-frame to valid (unit vector) "
1806 <<
" - but |normal| = " << std::sqrt(normMag2)
1807 <<
" - and |normal|^ = " << normMag2
1808 <<
" which differs from 1.0 by " << normMag2 - 1.0 <<
G4endl
1810 message <<
"============================================================"
1814 message <<
" State of Navigator: " <<
G4endl;
1815 message << *
this <<
G4endl;
1817 message <<
"============================================================"
1820 G4Exception(
"G4ITNavigator2::GetGlobalExitNormal()",
1822 "Value obtained from stored global-normal is not a unit vector.");
1839 #ifdef G4DEBUG_NAVIGATION
1849 edN <<
" State of Navigator: " <<
G4endl;
1853 G4Exception(
"G4ITNavigator2::GetGlobalExitNormal()",
1855 "LocalExitNormalAndCheck() did not calculate Normal.");
1864 edN <<
"G4ITNavigator2::GetGlobalExitNormal: "
1865 <<
" Using Local Normal - from call to GetLocalExitNormalAndCheck. "
1867 <<
" Local Exit Normal : " <<
" || = " << std::sqrt(localMag2)
1868 <<
" vec = " << localNormal <<
G4endl
1869 <<
" Global Exit Normal : " <<
" || = " << globalNormal.
mag()
1870 <<
" vec = " << globalNormal <<
G4endl;
1873 G4Exception(
"G4ITNavigator2::GetGlobalExitNormal()",
1875 "Value obtained from new local *solid* is incorrect.");
1876 localNormal = localNormal.
unit();
1882 #ifdef G4DEBUG_NAVIGATION
1897 edDfn <<
"Found difference in normals in case of exiting mother "
1898 <<
"- when Get is called after ComputingStep " <<
G4endl;
1899 edDfn <<
" Magnitude of diff = " << diffNorm.
mag() <<
G4endl;
1900 edDfn <<
" Normal stored (Global) = " << fExitNormalGlobalFrame
1902 edDfn <<
" Global Computed from Local = " << globalNormAgn <<
G4endl;
1903 G4Exception(
"G4ITNavigator::GetGlobalExitNormal()",
"GeomNav0003",
1909 return globalNormal;
1913 #define G4NEW_SAFETY 1
1931 #ifdef G4DEBUG_NAVIGATION
1935 G4cout <<
"*** G4ITNavigator2::ComputeSafety: ***" <<
G4endl
1936 <<
" Called at point: " << pGlobalpoint <<
G4endl;
1940 <<
" - Maximum length = " << pMaxLength <<
G4endl;
1943 G4cout <<
" ----- Upon entering Compute Safety:" <<
G4endl;
1955 if( endpointOnSurface && stayedOnEndpoint )
1957 #ifdef G4DEBUG_NAVIGATION
1960 G4cout <<
" G4Navigator::ComputeSafety() finds that point - "
1961 << pGlobalpoint <<
" - is on surface " <<
G4endl;
1991 #ifdef G4DEBUG_NAVIGATION
1994 G4cout <<
" G4ITNavigator2::ComputeSafety() relocates-in-volume to point: "
1995 << pGlobalpoint <<
G4endl;
2012 *motherPhysical, pMaxLength);
2013 newSafety= safetyTwo;
2019 newSafety= safetyOldVoxel;
2038 G4Exception(
"G4ITNavigator2::ComputeSafety()",
"GeomNav0001",
2044 "Not applicable for external volumes.");
2069 #ifdef G4DEBUG_NAVIGATION
2074 G4cout <<
" Returned value of Safety = " << newSafety <<
G4endl;
2076 G4cout.precision(oldcoutPrec);
2117 #ifdef CHECK_ORDER_OF_METHODS
2120 G4Exception(
"G4Navigator::RecheckDistanceToCurrentBoundary()",
2122 "Method must be called after ComputeStep(), before call to LocateMethod.");
2151 G4ThreeVector dgPosition= nextLevelTrf.TransformPoint(localPosition);
2152 G4ThreeVector dgDirection= nextLevelTrf.TransformAxis(localDirection);
2153 locatedDaug = candSolid->
Inside(dgPosition);
2162 true, &validExitNormal, &exitNormal);
2163 daughterStep= - distanceBackOut;
2190 daughterSafety= 0.0;
2199 *prDistance= daughterStep;
2200 if( prNewSafety ) { *prNewSafety= daughterSafety; }
2220 if( ProposedMove >= motherSafety )
2224 true, &validExitNormal, &exitNormal);
2228 motherStep= ProposedMove;
2233 motherSafety= motherSolid->
DistanceToIn(localPosition);
2234 if( ProposedMove >= motherSafety )
2244 if( prNewSafety ) { *prNewSafety= motherSafety; }
2253 *prDistance=
std::min( motherStep, daughterStep );
2256 *prNewSafety=
std::min( motherSafety, daughterSafety );
2282 G4cout <<
"The current state of G4Navigator is: " <<
G4endl;
2287 <<
" BlockedPhysicalVolume= " ;
2300 G4cout << std::setw(30) <<
" ExitNormal " <<
" "
2301 << std::setw( 5) <<
" Valid " <<
" "
2302 << std::setw( 9) <<
" Exiting " <<
" "
2303 << std::setw( 9) <<
" Entering" <<
" "
2304 << std::setw(15) <<
" Blocked:Volume " <<
" "
2305 << std::setw( 9) <<
" ReplicaNo" <<
" "
2306 << std::setw( 8) <<
" LastStepZero " <<
" "
2310 <<
", " << std::setw(7) <<
fExitNormal.z() <<
" ) "
2312 << std::setw( 9) <<
fExiting <<
" "
2316 G4cout << std::setw(15) <<
"None";
2333 G4cout.precision(oldcoutPrec);
2363 G4double shiftOrigin = std::sqrt(shiftOriginSafSq);
2366 if( diffShiftSaf > fAccuracyForWarning )
2370 std::ostringstream
message, suggestion;
2371 message <<
"Accuracy error or slightly inaccurate position shift."
2373 <<
" The Step's starting point has moved "
2374 << std::sqrt(moveLenSq)/
mm <<
" mm " <<
G4endl
2375 <<
" since the last call to a Locate method." <<
G4endl
2376 <<
" This has resulted in moving "
2377 << shiftOrigin/
mm <<
" mm "
2378 <<
" from the last point at which the safety "
2379 <<
" was calculated " <<
G4endl
2380 <<
" which is more than the computed safety= "
2381 << fPreviousSafety/
mm <<
" mm at that point." <<
G4endl
2382 <<
" This difference is "
2383 << diffShiftSaf/
mm <<
" mm." <<
G4endl
2384 <<
" The tolerated accuracy is "
2385 << fAccuracyForException/
mm <<
" mm.";
2389 if( ((++warnNow % 100) == 1) )
2392 <<
" This problem can be due to either " <<
G4endl
2393 <<
" - a process that has proposed a displacement"
2394 <<
" larger than the current safety , or" <<
G4endl
2395 <<
" - inaccuracy in the computation of the safety";
2396 suggestion <<
"We suggest that you " << G4endl
2397 <<
" - find i) what particle is being tracked, and "
2398 <<
" ii) through what part of your geometry " << G4endl
2399 <<
" for example by re-running this event with "
2401 <<
" /tracking/verbose 1 " << G4endl
2402 <<
" - check which processes you declare for"
2403 <<
" this particle (and look at non-standard ones)"
2405 <<
" - in case, create a detailed logfile"
2406 <<
" of this event using:" << G4endl
2407 <<
" /tracking/verbose 6 ";
2411 message,
G4String(suggestion.str()));
2412 G4cout.precision(oldcoutPrec);
2413 G4cerr.precision(oldcerrPrec);
2415 #ifdef G4DEBUG_NAVIGATION
2418 G4cerr <<
"WARNING - G4ITNavigator2::ComputeStep()" <<
G4endl
2419 <<
" The Step's starting point has moved "
2420 << std::sqrt(moveLenSq) <<
"," <<
G4endl
2421 <<
" which has taken it to the limit of"
2422 <<
" the current safety. " <<
G4endl;
2427 if ( shiftOriginSafSq >
sqr(safetyPlus) )
2430 message <<
"May lead to a crash or unreliable results." <<
G4endl
2431 <<
" Position has shifted considerably without"
2432 <<
" notifying the navigator !" <<
G4endl
2433 <<
" Tolerated safety: " << safetyPlus <<
G4endl
2434 <<
" Computed shift : " << shiftOriginSafSq;
2435 G4Exception(
"G4ITNavigator2::ComputeStep()",
"GeomNav1002",
2452 G4int oldcoutPrec = os.precision(4);
2455 os <<
"The current state of G4ITNavigator2 is: " <<
G4endl;
2456 os <<
" ValidExitNormal= " << n.fValidExitNormal << G4endl
2457 <<
" ExitNormal = " << n.fExitNormal << G4endl
2458 <<
" Exiting = " << n.fExiting << G4endl
2459 <<
" Entering = " << n.fEntering << G4endl
2460 <<
" BlockedPhysicalVolume= " ;
2462 if (n.fBlockedPhysicalVolume==0)
2468 os << n.fBlockedPhysicalVolume->GetName();
2472 <<
" BlockedReplicaNo = " << n.fBlockedReplicaNo << G4endl
2473 <<
" LastStepWasZero = " << n.fLastStepWasZero << G4endl
2479 os << std::setw(30) <<
" ExitNormal " <<
" "
2480 << std::setw( 5) <<
" Valid " <<
" "
2481 << std::setw( 9) <<
" Exiting " <<
" "
2482 << std::setw( 9) <<
" Entering" <<
" "
2483 << std::setw(15) <<
" Blocked:Volume " <<
" "
2484 << std::setw( 9) <<
" ReplicaNo" <<
" "
2485 << std::setw( 8) <<
" LastStepZero " <<
" "
2487 os <<
"( " << std::setw(7) << n.fExitNormal.x()
2488 <<
", " << std::setw(7) << n.fExitNormal.y()
2489 <<
", " << std::setw(7) << n.fExitNormal.z() <<
" ) "
2490 << std::setw( 5) << n.fValidExitNormal <<
" "
2491 << std::setw( 9) << n.fExiting <<
" "
2492 << std::setw( 9) << n.fEntering <<
" ";
2494 if ( n.fBlockedPhysicalVolume==0 )
2495 { os << std::setw(15) <<
"None"; }
2497 { os << std::setw(15)<< n.fBlockedPhysicalVolume->GetName(); }
2499 os << std::setw( 9) << n.fBlockedReplicaNo <<
" "
2500 << std::setw( 8) << n.fLastStepWasZero <<
" "
2506 os <<
" Current Localpoint = " << n.fLastLocatedPointLocal <<
G4endl;
2507 os <<
" PreviousSftOrigin = " << n.fPreviousSftOrigin <<
G4endl;
2508 os <<
" PreviousSafety = " << n.fPreviousSafety <<
G4endl;
2512 os <<
"Current History: " <<
G4endl << n.fHistory;
2515 os.precision(oldcoutPrec);
2529 return solid->Inside(localPoint);
2542 std::vector<std::vector<double> > fExtend;
2560 for(
size_t i = 0 ; i < 3 ; ++i)