| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // G4Navigator class Implementation
27 //
28 // Original author: Paul Kent, July 95/96
29 // Responsible 1996-present: John Apostolakis, Gabriele Cosmo
30 // Additional revisions by: Pedro Arce, Vladimir Grichine
31 // --------------------------------------------------------------------
32
33 #include <iomanip>
34
35 #include "G4Navigator.hh"
36 #include "G4ios.hh"
37 #include "G4SystemOfUnits.hh"
38 #include "G4GeometryTolerance.hh"
39 #include "G4VPhysicalVolume.hh"
40
41 #include "G4VoxelSafety.hh"
42 #include "G4SafetyCalculator.hh"
43
44 // Constant determining how precise normals should be (how close to unit
45 // vectors). If exceeded, warnings will be issued.
46 // Can be CLHEP::perMillion (its old default) for geometry checking.
47 //
48 static const G4double kToleranceNormalCheck = CLHEP::perThousand;
49
50 // ********************************************************************
51 // Constructor
52 // ********************************************************************
53 //
54 G4Navigator::G4Navigator()
55 {
56 ResetStackAndState();
57 // Initialises also all
58 // - exit / entry flags
59 // - flags & variables for exit normals
60 // - zero step counters
61 // - blocked volume
62
63 if( fVerbose > 2 )
64 {
65 G4cout << " G4Navigator parameters: Action Threshold (No Zero Steps) = "
66 << fActionThreshold_NoZeroSteps
67 << " Abandon Threshold (No Zero Steps) = "
68 << fAbandonThreshold_NoZeroSteps << G4endl;
69 }
70 kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
71 fMinStep = 0.05*kCarTolerance;
72 fSqTol = sqr(kCarTolerance);
73
74 fregularNav.SetNormalNavigation( &fnormalNav );
75
76 fStepEndPoint = G4ThreeVector( kInfinity, kInfinity, kInfinity );
77 fLastStepEndPointLocal = G4ThreeVector( kInfinity, kInfinity, kInfinity );
78
79 fpVoxelSafety = new G4VoxelSafety();
80 fpvoxelNav = new G4VoxelNavigation();
81 fpSafetyCalculator = new G4SafetyCalculator( *this, fHistory );
82 fpSafetyCalculator->SetExternalNavigation(fpExternalNav);
83 }
84
85 // ********************************************************************
86 // Destructor
87 // ********************************************************************
88 //
89 G4Navigator::~G4Navigator()
90 {
91 delete fpVoxelSafety;
92 delete fpExternalNav;
93 delete fpvoxelNav;
94 delete fpSafetyCalculator;
95 }
96
97 // ********************************************************************
98 // ResetHierarchyAndLocate
99 // ********************************************************************
100 //
101 G4VPhysicalVolume*
102 G4Navigator::ResetHierarchyAndLocate(const G4ThreeVector& p,
103 const G4ThreeVector& direction,
104 const G4TouchableHistory& h)
105 {
106 ResetState();
107 fHistory = *h.GetHistory();
108 SetupHierarchy();
109 fLastTriedStepComputation = false; // Redundant, but best
110 return LocateGlobalPointAndSetup(p, &direction, true, false);
111 }
112
113 // ********************************************************************
114 // LocateGlobalPointAndSetup
115 //
116 // Locate the point in the hierarchy return 0 if outside
117 // The direction is required
118 // - if on an edge shared by more than two surfaces
119 // (to resolve likely looping in tracking)
120 // - at initial location of a particle
121 // (to resolve potential ambiguity at boundary)
122 //
123 // Flags on exit: (comments to be completed)
124 // fEntering - True if entering `daughter' volume (or replica)
125 // whether daughter of last mother directly
126 // or daughter of that volume's ancestor.
127 // fExiting - True if exited 'mother' volume
128 // (always ? - how about if going back down ? - tbc)
129 // ********************************************************************
130 //
131 G4VPhysicalVolume*
132 G4Navigator::LocateGlobalPointAndSetup( const G4ThreeVector& globalPoint,
133 const G4ThreeVector* pGlobalDirection,
134 const G4bool relativeSearch,
135 const G4bool ignoreDirection )
136 {
137 G4bool notKnownContained = true, noResult;
138 G4VPhysicalVolume *targetPhysical;
139 G4LogicalVolume *targetLogical;
140 G4VSolid *targetSolid = nullptr;
141 G4ThreeVector localPoint, globalDirection;
142 EInside insideCode;
143
144 G4bool considerDirection = (pGlobalDirection != nullptr) && ((!ignoreDirection) || fLocatedOnEdge);
145
146 fLastTriedStepComputation = false;
147 fChangedGrandMotherRefFrame = false; // For local exit normal
148
149 if( considerDirection )
150 {
151 globalDirection=*pGlobalDirection;
152 }
153
154 #ifdef G4VERBOSE
155 if( fVerbose > 2 )
156 {
157 G4long oldcoutPrec = G4cout.precision(8);
158 G4cout << "*** G4Navigator::LocateGlobalPointAndSetup: ***" << G4endl;
159 G4cout << " Called with arguments: " << G4endl
160 << " Globalpoint = " << globalPoint << G4endl
161 << " RelativeSearch = " << relativeSearch << G4endl;
162 if( fVerbose >= 4 )
163 {
164 G4cout << " ----- Upon entering:" << G4endl;
165 PrintState();
166 }
167 G4cout.precision(oldcoutPrec);
168 }
169 #endif
170
171 G4int noLevelsExited = 0;
172
173 if ( !relativeSearch )
174 {
175 ResetStackAndState();
176 }
177 else
178 {
179 if ( fWasLimitedByGeometry )
180 {
181 fWasLimitedByGeometry = false;
182 fEnteredDaughter = fEntering; // Remember
183 fExitedMother = fExiting; // Remember
184 if ( fExiting )
185 {
186 ++noLevelsExited; // count this first level entered too
187
188 if ( fHistory.GetDepth() != 0 )
189 {
190 fBlockedPhysicalVolume = fHistory.GetTopVolume();
191 fBlockedReplicaNo = fHistory.GetTopReplicaNo();
192 fHistory.BackLevel();
193 }
194 else
195 {
196 fLastLocatedPointLocal = localPoint;
197 fLocatedOutsideWorld = true;
198 fBlockedPhysicalVolume = nullptr; // to be sure
199 fBlockedReplicaNo = -1;
200 fEntering = false; // No longer
201 fEnteredDaughter = false;
202 fExitedMother = true; // ??
203
204 return nullptr; // Have exited world volume
205 }
206 // A fix for the case where a volume is "entered" at an edge
207 // and a coincident surface exists outside it.
208 // - This stops it from exiting further volumes and cycling
209 // - However ReplicaNavigator treats this case itself
210 //
211 // assert( fBlockedPhysicalVolume!=0 );
212
213 // Expect to be on edge => on surface
214 //
215 if ( fLocatedOnEdge && (VolumeType(fBlockedPhysicalVolume)!=kReplica ))
216 {
217 fExiting = false;
218 // Consider effect on Exit Normal !?
219 }
220 }
221 else
222 if ( fEntering )
223 {
224 switch (VolumeType(fBlockedPhysicalVolume))
225 {
226 case kNormal:
227 fHistory.NewLevel(fBlockedPhysicalVolume, kNormal,
228 fBlockedPhysicalVolume->GetCopyNo());
229 break;
230 case kReplica:
231 freplicaNav.ComputeTransformation(fBlockedReplicaNo,
232 fBlockedPhysicalVolume);
233 fHistory.NewLevel(fBlockedPhysicalVolume, kReplica,
234 fBlockedReplicaNo);
235 fBlockedPhysicalVolume->SetCopyNo(fBlockedReplicaNo);
236 break;
237 case kParameterised:
238 if( fBlockedPhysicalVolume->GetRegularStructureId() == 0 )
239 {
240 G4VSolid *pSolid;
241 G4VPVParameterisation *pParam;
242 G4TouchableHistory parentTouchable( fHistory );
243 pParam = fBlockedPhysicalVolume->GetParameterisation();
244 pSolid = pParam->ComputeSolid(fBlockedReplicaNo,
245 fBlockedPhysicalVolume);
246 pSolid->ComputeDimensions(pParam, fBlockedReplicaNo,
247 fBlockedPhysicalVolume);
248 pParam->ComputeTransformation(fBlockedReplicaNo,
249 fBlockedPhysicalVolume);
250 fHistory.NewLevel(fBlockedPhysicalVolume, kParameterised,
251 fBlockedReplicaNo);
252 fBlockedPhysicalVolume->SetCopyNo(fBlockedReplicaNo);
253 //
254 // Set the correct solid and material in Logical Volume
255 //
256 G4LogicalVolume *pLogical;
257 pLogical = fBlockedPhysicalVolume->GetLogicalVolume();
258 pLogical->SetSolid( pSolid );
259 pLogical->UpdateMaterial(pParam ->
260 ComputeMaterial(fBlockedReplicaNo,
261 fBlockedPhysicalVolume,
262 &parentTouchable));
263 }
264 break;
265 case kExternal:
266 G4Exception("G4Navigator::LocateGlobalPointAndSetup()",
267 "GeomNav0001", FatalException,
268 "Extra levels not applicable for external volumes.");
269 break;
270 }
271 fEntering = false;
272 fBlockedPhysicalVolume = nullptr;
273 localPoint = fHistory.GetTopTransform().TransformPoint(globalPoint);
274 notKnownContained = false;
275 }
276 }
277 else
278 {
279 fBlockedPhysicalVolume = nullptr;
280 fEntering = false;
281 fEnteredDaughter = false; // Full Step was not taken, did not enter
282 fExiting = false;
283 fExitedMother = false; // Full Step was not taken, did not exit
284 }
285 }
286 //
287 // Search from top of history up through geometry until
288 // containing volume found:
289 // If on
290 // o OUTSIDE - Back up level, not/no longer exiting volumes
291 // o SURFACE and EXITING - Back up level, setting new blocking no.s
292 // else
293 // o containing volume found
294 //
295
296 while (notKnownContained) // Loop checking, 07.10.2016, J.Apostolakis
297 {
298 EVolume topVolumeType = fHistory.GetTopVolumeType();
299 if (topVolumeType!=kReplica && topVolumeType!=kExternal)
300 {
301 targetSolid = fHistory.GetTopVolume()->GetLogicalVolume()->GetSolid();
302 localPoint = fHistory.GetTopTransform().TransformPoint(globalPoint);
303 insideCode = targetSolid->Inside(localPoint);
304 #ifdef G4VERBOSE
305 if(( fVerbose == 1 ) && ( fCheck ))
306 {
307 G4String solidResponse = "-kInside-";
308 if (insideCode == kOutside)
309 {
310 solidResponse = "-kOutside-";
311 }
312 else if (insideCode == kSurface)
313 {
314 solidResponse = "-kSurface-";
315 }
316 G4cout << "*** G4Navigator::LocateGlobalPointAndSetup(): ***" << G4endl
317 << " Invoked Inside() for solid: " << targetSolid->GetName()
318 << ". Solid replied: " << solidResponse << G4endl
319 << " For local point p: " << localPoint << G4endl;
320 }
321 #endif
322 }
323 else
324 {
325 if( topVolumeType == kReplica )
326 {
327 insideCode = freplicaNav.BackLocate(fHistory, globalPoint, localPoint,
328 fExiting, notKnownContained);
329 // !CARE! if notKnownContained returns false then the point is within
330 // the containing placement volume of the replica(s). If insidecode
331 // will result in the history being backed up one level, then the
332 // local point returned is the point in the system of this new level
333 }
334 else
335 {
336 targetSolid = fHistory.GetTopVolume()->GetLogicalVolume()->GetSolid();
337 localPoint = fHistory.GetTopTransform().TransformPoint(globalPoint);
338 G4ThreeVector localDirection =
339 fHistory.GetTopTransform().TransformAxis(globalDirection);
340 insideCode = fpExternalNav->Inside(targetSolid, localPoint, localDirection);
341 }
342 }
343
344 // Point is inside current volume, break out of the loop
345 if ( insideCode == kInside ) { break; }
346
347 // Point is outside current volume, move up a level in the hierarchy
348 if ( insideCode == kOutside )
349 {
350 ++noLevelsExited;
351
352 // Exiting world volume
353 if ( fHistory.GetDepth() == 0 )
354 {
355 fLocatedOutsideWorld = true;
356 fLastLocatedPointLocal = localPoint;
357 return nullptr;
358 }
359
360 fBlockedPhysicalVolume = fHistory.GetTopVolume();
361 fBlockedReplicaNo = fHistory.GetTopReplicaNo();
362 fHistory.BackLevel();
363 fExiting = false;
364
365 if( noLevelsExited > 1 )
366 {
367 // The first transformation was done by the sub-navigator
368 //
369 if(const auto *mRot = fBlockedPhysicalVolume->GetRotation())
370 {
371 fGrandMotherExitNormal *= (*mRot).inverse();
372 fChangedGrandMotherRefFrame = true;
373 }
374 }
375 continue;
376 }
377
378 // Point is on the surface of a volume
379 G4bool isExiting = fExiting;
380 if( (!fExiting) && considerDirection )
381 {
382 // Figure out whether we are exiting this level's volume
383 // by using the direction
384 //
385 G4bool directionExiting = false;
386 G4ThreeVector localDirection =
387 fHistory.GetTopTransform().TransformAxis(globalDirection);
388
389 // Make sure localPoint in correct reference frame
390 // ( Was it already correct ? How ? )
391 //
392 localPoint= fHistory.GetTopTransform().TransformPoint(globalPoint);
393 if ( fHistory.GetTopVolumeType() != kReplica )
394 {
395 G4ThreeVector normal = targetSolid->SurfaceNormal(localPoint);
396 directionExiting = normal.dot(localDirection) > 0.0;
397 isExiting = isExiting || directionExiting;
398 }
399 }
400
401 // Point is on a surface, but no longer exiting, break out of the loop
402 if ( !isExiting ) { break; }
403
404 ++noLevelsExited;
405
406 // Point is on the outer surface, leaving world volume
407 if ( fHistory.GetDepth() == 0 )
408 {
409 fLocatedOutsideWorld = true;
410 fLastLocatedPointLocal = localPoint;
411 return nullptr;
412 }
413
414 // Point is still on a surface, but exited a volume not necessarily convex
415 fValidExitNormal = false;
416 fBlockedPhysicalVolume = fHistory.GetTopVolume();
417 fBlockedReplicaNo = fHistory.GetTopReplicaNo();
418 fHistory.BackLevel();
419
420 if( noLevelsExited > 1 )
421 {
422 // The first transformation was done by the sub-navigator
423 //
424 const G4RotationMatrix* mRot =
425 fBlockedPhysicalVolume->GetRotation();
426 if( mRot != nullptr )
427 {
428 fGrandMotherExitNormal *= (*mRot).inverse();
429 fChangedGrandMotherRefFrame = true;
430 }
431 }
432 } // END while (notKnownContained)
433 //
434 // Search downwards until deepest containing volume found,
435 // blocking fBlockedPhysicalVolume/BlockedReplicaNum
436 //
437 // 3 Cases:
438 //
439 // o Parameterised daughters
440 // =>Must be one G4PVParameterised daughter & voxels
441 // o Positioned daughters & voxels
442 // o Positioned daughters & no voxels
443
444 noResult = true; // noResult should be renamed to
445 // something like enteredLevel, as that is its meaning.
446 do
447 {
448 // Determine `type' of current mother volume
449 //
450 targetPhysical = fHistory.GetTopVolume();
451 if (targetPhysical == nullptr) { break; }
452 targetLogical = targetPhysical->GetLogicalVolume();
453 switch( CharacteriseDaughters(targetLogical) )
454 {
455 case kNormal:
456 if ( targetLogical->GetVoxelHeader() != nullptr ) // use optimised navigation
457 {
458 noResult = GetVoxelNavigator().LevelLocate(fHistory,
459 fBlockedPhysicalVolume,
460 fBlockedReplicaNo,
461 globalPoint,
462 pGlobalDirection,
463 considerDirection,
464 localPoint);
465 }
466 else // do not use optimised navigation
467 {
468 noResult = fnormalNav.LevelLocate(fHistory,
469 fBlockedPhysicalVolume,
470 fBlockedReplicaNo,
471 globalPoint,
472 pGlobalDirection,
473 considerDirection,
474 localPoint);
475 }
476 break;
477 case kReplica:
478 noResult = freplicaNav.LevelLocate(fHistory,
479 fBlockedPhysicalVolume,
480 fBlockedReplicaNo,
481 globalPoint,
482 pGlobalDirection,
483 considerDirection,
484 localPoint);
485 break;
486 case kParameterised:
487 if( GetDaughtersRegularStructureId(targetLogical) != 1 )
488 {
489 noResult = fparamNav.LevelLocate(fHistory,
490 fBlockedPhysicalVolume,
491 fBlockedReplicaNo,
492 globalPoint,
493 pGlobalDirection,
494 considerDirection,
495 localPoint);
496 }
497 else // Regular structure
498 {
499 noResult = fregularNav.LevelLocate(fHistory,
500 fBlockedPhysicalVolume,
501 fBlockedReplicaNo,
502 globalPoint,
503 pGlobalDirection,
504 considerDirection,
505 localPoint);
506 }
507 break;
508 case kExternal:
509 noResult = fpExternalNav->LevelLocate(fHistory,
510 fBlockedPhysicalVolume,
511 fBlockedReplicaNo,
512 globalPoint,
513 pGlobalDirection,
514 considerDirection,
515 localPoint);
516 break;
517 }
518
519 // LevelLocate returns true if it finds a daughter volume
520 // in which globalPoint is inside (or on the surface).
521
522 if ( noResult )
523 {
524 // Entering a daughter after ascending
525 //
526 // The blocked volume is no longer valid - it was for another level
527 //
528 fBlockedPhysicalVolume = nullptr;
529 fBlockedReplicaNo = -1;
530
531 // fEntering should be false -- else blockedVolume is assumed good.
532 // fEnteredDaughter is used for ExitNormal
533 //
534 fEntering = false;
535 fEnteredDaughter = true;
536
537 if( fExitedMother )
538 {
539 G4VPhysicalVolume* enteredPhysical = fHistory.GetTopVolume();
540 const G4RotationMatrix* mRot = enteredPhysical->GetRotation();
541 if( mRot != nullptr )
542 {
543 // Go deeper, i.e. move 'down' in the hierarchy
544 // Apply direct rotation, not inverse
545 //
546 fGrandMotherExitNormal *= (*mRot);
547 fChangedGrandMotherRefFrame= true;
548 }
549 }
550
551 #ifdef G4DEBUG_NAVIGATION
552 if( fVerbose > 2 )
553 {
554 G4VPhysicalVolume* enteredPhysical = fHistory.GetTopVolume();
555 G4cout << "*** G4Navigator::LocateGlobalPointAndSetup() ***" << G4endl;
556 G4cout << " Entering volume: " << enteredPhysical->GetName()
557 << G4endl;
558 }
559 #endif
560 }
561 } while (noResult); // Loop checking, 07.10.2016, J.Apostolakis
562
563 fLastLocatedPointLocal = localPoint;
564
565 #ifdef G4VERBOSE
566 if( fVerbose >= 4 )
567 {
568 G4long oldcoutPrec = G4cout.precision(8);
569 G4String curPhysVol_Name("None");
570 if (targetPhysical != nullptr) { curPhysVol_Name = targetPhysical->GetName(); }
571 G4cout << " Return value = new volume = " << curPhysVol_Name << G4endl;
572 G4cout << " ----- Upon exiting:" << G4endl;
573 PrintState();
574 if( fVerbose >= 5 )
575 {
576 G4cout << "Upon exiting LocateGlobalPointAndSetup():" << G4endl;
577 G4cout << " History = " << G4endl << fHistory << G4endl << G4endl;
578 }
579 G4cout.precision(oldcoutPrec);
580 }
581 #endif
582
583 fLocatedOutsideWorld = false;
584
585 return targetPhysical;
586 }
587
588 // ********************************************************************
589 // LocateGlobalPointWithinVolume
590 //
591 // -> the state information of this Navigator and its subNavigators
592 // is updated in order to start the next step at pGlobalpoint
593 // -> no check is performed whether pGlobalpoint is inside the
594 // original volume (this must be the case).
595 //
596 // Note: a direction could be added to the arguments, to aid in future
597 // optional checking (via the old code below, flagged by OLD_LOCATE).
598 // [ This would be done only in verbose mode ]
599 // ********************************************************************
600 //
601 void
602 G4Navigator::LocateGlobalPointWithinVolume(const G4ThreeVector& pGlobalpoint)
603 {
604 #ifdef G4DEBUG_NAVIGATION
605 assert( !fWasLimitedByGeometry );
606 // Check: Either step was not limited by a boundary or
607 // else the full step is no longer being taken
608 #endif
609
610 fLastLocatedPointLocal = ComputeLocalPoint(pGlobalpoint);
611 fLastTriedStepComputation = false;
612 fChangedGrandMotherRefFrame = false; // Frame for Exit Normal
613
614 // For the case of Voxel (or Parameterised) volume the respective
615 // Navigator must be messaged to update its voxel information etc
616
617 // Update the state of the Sub Navigators
618 // - in particular any voxel information they store/cache
619 //
620 G4VPhysicalVolume* motherPhysical = fHistory.GetTopVolume();
621 G4LogicalVolume* motherLogical = motherPhysical->GetLogicalVolume();
622
623 switch( CharacteriseDaughters(motherLogical) )
624 {
625 case kNormal:
626 GetVoxelNavigator().RelocateWithinVolume( motherPhysical, fLastLocatedPointLocal );
627 break;
628 case kParameterised:
629 fparamNav.RelocateWithinVolume( motherPhysical, fLastLocatedPointLocal );
630 break;
631 case kReplica:
632 // Nothing to do
633 break;
634 case kExternal:
635 fpExternalNav->RelocateWithinVolume( motherPhysical,
636 fLastLocatedPointLocal );
637 break;
638 }
639
640 // Reset the state variables
641 // - which would have been affected
642 // by the 'equivalent' call to LocateGlobalPointAndSetup
643 // - who's values have been invalidated by the 'move'.
644 //
645 fBlockedPhysicalVolume = nullptr;
646 fBlockedReplicaNo = -1;
647 fEntering = false;
648 fEnteredDaughter = false; // Boundary not encountered, did not enter
649 fExiting = false;
650 fExitedMother = false; // Boundary not encountered, did not exit
651 }
652
653 // ********************************************************************
654 // SetSavedState
655 //
656 // Save the state, in case this is a parasitic call
657 // Save fValidExitNormal, fExitNormal, fExiting, fEntering,
658 // fBlockedPhysicalVolume, fBlockedReplicaNo, fLastStepWasZero;
659 // ********************************************************************
660 //
661 void G4Navigator::SetSavedState()
662 {
663 // Note: the state of dependent objects is not currently saved.
664 // ( This means that the full state is changed by calls between
665 // SetSavedState() and RestoreSavedState();
666
667 fSaveState.sExitNormal = fExitNormal;
668 fSaveState.sValidExitNormal = fValidExitNormal;
669 fSaveState.sExiting = fExiting;
670 fSaveState.sEntering = fEntering;
671
672 fSaveState.spBlockedPhysicalVolume = fBlockedPhysicalVolume;
673 fSaveState.sBlockedReplicaNo = fBlockedReplicaNo;
674
675 fSaveState.sLastStepWasZero = static_cast<G4int>(fLastStepWasZero);
676
677 fSaveState.sLocatedOutsideWorld = fLocatedOutsideWorld;
678 fSaveState.sLastLocatedPointLocal = fLastLocatedPointLocal;
679 fSaveState.sEnteredDaughter = fEnteredDaughter;
680 fSaveState.sExitedMother = fExitedMother;
681 fSaveState.sWasLimitedByGeometry = fWasLimitedByGeometry;
682
683 // Even the safety sphere - if you want to change it do it explicitly!
684 //
685 fSaveState.sPreviousSftOrigin = fPreviousSftOrigin;
686 fSaveState.sPreviousSafety = fPreviousSafety;
687 }
688
689 // ********************************************************************
690 // RestoreSavedState
691 //
692 // Restore the state (in Compute Step), in case this is a parasitic call
693 // ********************************************************************
694 //
695 void G4Navigator::RestoreSavedState()
696 {
697 fExitNormal = fSaveState.sExitNormal;
698 fValidExitNormal = fSaveState.sValidExitNormal;
699 fExiting = fSaveState.sExiting;
700 fEntering = fSaveState.sEntering;
701
702 fBlockedPhysicalVolume = fSaveState.spBlockedPhysicalVolume;
703 fBlockedReplicaNo = fSaveState.sBlockedReplicaNo;
704
705 fLastStepWasZero = (fSaveState.sLastStepWasZero != 0);
706
707 fLocatedOutsideWorld = fSaveState.sLocatedOutsideWorld;
708 fLastLocatedPointLocal = fSaveState.sLastLocatedPointLocal;
709 fEnteredDaughter = fSaveState.sEnteredDaughter;
710 fExitedMother = fSaveState.sExitedMother;
711 fWasLimitedByGeometry = fSaveState.sWasLimitedByGeometry;
712
713 // The 'expected' behaviour is to restore these too (fix 2014.05.26)
714 fPreviousSftOrigin = fSaveState.sPreviousSftOrigin;
715 fPreviousSafety = fSaveState.sPreviousSafety;
716 }
717
718 // ********************************************************************
719 // ComputeStep
720 //
721 // Computes the next geometric Step: intersections with current
722 // mother and `daughter' volumes.
723 //
724 // NOTE:
725 //
726 // Flags on entry:
727 // --------------
728 // fValidExitNormal - Normal of exited volume is valid (convex, not a
729 // coincident boundary)
730 // fExitNormal - Surface normal of exited volume
731 // fExiting - True if have exited solid
732 //
733 // fBlockedPhysicalVolume - Ptr to exited volume (or 0)
734 // fBlockedReplicaNo - Replication no of exited volume
735 // fLastStepWasZero - True if last Step size was almost zero.
736 //
737 // Flags on exit:
738 // -------------
739 // fValidExitNormal - True if surface normal of exited volume is valid
740 // fExitNormal - Surface normal of exited volume rotated to mothers
741 // reference system
742 // fExiting - True if exiting mother
743 // fEntering - True if entering `daughter' volume (or replica)
744 // fBlockedPhysicalVolume - Ptr to candidate (entered) volume
745 // fBlockedReplicaNo - Replication no of candidate (entered) volume
746 // fLastStepWasZero - True if this Step size was almost zero.
747 // ********************************************************************
748 //
749 G4double G4Navigator::ComputeStep( const G4ThreeVector& pGlobalpoint,
750 const G4ThreeVector& pDirection,
751 const G4double pCurrentProposedStepLength,
752 G4double& pNewSafety)
753 {
754 #ifdef G4DEBUG_NAVIGATION
755 static G4ThreadLocal G4int sNavCScalls = 0;
756 ++sNavCScalls;
757 #endif
758
759 G4ThreeVector localDirection = ComputeLocalAxis(pDirection);
760 G4double Step = kInfinity;
761 G4VPhysicalVolume *motherPhysical = fHistory.GetTopVolume();
762 G4LogicalVolume *motherLogical = motherPhysical->GetLogicalVolume();
763
764 // All state relating to exiting normals must be reset
765 //
766 fExitNormalGlobalFrame = G4ThreeVector( 0., 0., 0.);
767 // Reset value - to erase its memory
768 fChangedGrandMotherRefFrame = false;
769 // Reset - used for local exit normal
770 fGrandMotherExitNormal = G4ThreeVector( 0., 0., 0.);
771 fCalculatedExitNormal = false;
772 // Reset for new step
773
774 #ifdef G4VERBOSE
775 if( fVerbose > 0 )
776 {
777 G4cout << "*** G4Navigator::ComputeStep: ***" << G4endl;
778 G4cout << " Volume = " << motherPhysical->GetName()
779 << " - Proposed step length = " << pCurrentProposedStepLength
780 << G4endl;
781 #ifdef G4DEBUG_NAVIGATION
782 if( fVerbose >= 2 )
783 {
784 G4cout << " Called with the arguments: " << G4endl
785 << " Globalpoint = " << std::setw(25) << pGlobalpoint << G4endl
786 << " Direction = " << std::setw(25) << pDirection << G4endl;
787 if( fVerbose >= 4 )
788 {
789 G4cout << " ---- Upon entering : State" << G4endl;
790 PrintState();
791 }
792 }
793 #endif
794 }
795 #endif
796
797 G4ThreeVector newLocalPoint = ComputeLocalPoint(pGlobalpoint);
798
799 if( newLocalPoint != fLastLocatedPointLocal )
800 {
801 // Check whether the relocation is within safety
802 //
803 G4ThreeVector oldLocalPoint = fLastLocatedPointLocal;
804 G4double moveLenSq = (newLocalPoint-oldLocalPoint).mag2();
805
806 if ( moveLenSq >= fSqTol )
807 {
808 #ifdef G4VERBOSE
809 ComputeStepLog(pGlobalpoint, moveLenSq);
810 #endif
811 // Relocate the point within the same volume
812 //
813 LocateGlobalPointWithinVolume( pGlobalpoint );
814 }
815 }
816 if ( fHistory.GetTopVolumeType()!=kReplica )
817 {
818 switch( CharacteriseDaughters(motherLogical) )
819 {
820 case kNormal:
821 if ( motherLogical->GetVoxelHeader() != nullptr )
822 {
823 Step = GetVoxelNavigator().ComputeStep(fLastLocatedPointLocal,
824 localDirection,
825 pCurrentProposedStepLength,
826 pNewSafety,
827 fHistory,
828 fValidExitNormal,
829 fExitNormal,
830 fExiting,
831 fEntering,
832 &fBlockedPhysicalVolume,
833 fBlockedReplicaNo);
834
835 }
836 else
837 {
838 if( motherPhysical->GetRegularStructureId() == 0 )
839 {
840 Step = fnormalNav.ComputeStep(fLastLocatedPointLocal,
841 localDirection,
842 pCurrentProposedStepLength,
843 pNewSafety,
844 fHistory,
845 fValidExitNormal,
846 fExitNormal,
847 fExiting,
848 fEntering,
849 &fBlockedPhysicalVolume,
850 fBlockedReplicaNo);
851 }
852 else // Regular (non-voxelised) structure
853 {
854 LocateGlobalPointAndSetup( pGlobalpoint, &pDirection, true, true );
855 //
856 // if physical process limits the step, the voxel will not be the
857 // one given by ComputeStepSkippingEqualMaterials() and the local
858 // point will be wrongly calculated.
859
860 // There is a problem: when msc limits the step and the point is
861 // assigned wrongly to phantom in previous step (while it is out
862 // of the container volume). Then LocateGlobalPointAndSetup() has
863 // reset the history topvolume to world.
864 //
865 if(fHistory.GetTopVolume()->GetRegularStructureId() == 0 )
866 {
867 G4Exception("G4Navigator::ComputeStep()",
868 "GeomNav1001", JustWarning,
869 "Point is relocated in voxels, while it should be outside!");
870 Step = fnormalNav.ComputeStep(fLastLocatedPointLocal,
871 localDirection,
872 pCurrentProposedStepLength,
873 pNewSafety,
874 fHistory,
875 fValidExitNormal,
876 fExitNormal,
877 fExiting,
878 fEntering,
879 &fBlockedPhysicalVolume,
880 fBlockedReplicaNo);
881 }
882 else
883 {
884 Step = fregularNav.
885 ComputeStepSkippingEqualMaterials(fLastLocatedPointLocal,
886 localDirection,
887 pCurrentProposedStepLength,
888 pNewSafety,
889 fHistory,
890 fValidExitNormal,
891 fExitNormal,
892 fExiting,
893 fEntering,
894 &fBlockedPhysicalVolume,
895 fBlockedReplicaNo,
896 motherPhysical);
897 }
898 }
899 }
900 break;
901 case kParameterised:
902 if( GetDaughtersRegularStructureId(motherLogical) != 1 )
903 {
904 Step = fparamNav.ComputeStep(fLastLocatedPointLocal,
905 localDirection,
906 pCurrentProposedStepLength,
907 pNewSafety,
908 fHistory,
909 fValidExitNormal,
910 fExitNormal,
911 fExiting,
912 fEntering,
913 &fBlockedPhysicalVolume,
914 fBlockedReplicaNo);
915 }
916 else // Regular structure
917 {
918 Step = fregularNav.ComputeStep(fLastLocatedPointLocal,
919 localDirection,
920 pCurrentProposedStepLength,
921 pNewSafety,
922 fHistory,
923 fValidExitNormal,
924 fExitNormal,
925 fExiting,
926 fEntering,
927 &fBlockedPhysicalVolume,
928 fBlockedReplicaNo);
929 }
930 break;
931 case kReplica:
932 G4Exception("G4Navigator::ComputeStep()", "GeomNav0001",
933 FatalException, "Not applicable for replicated volumes.");
934 break;
935 case kExternal:
936 Step = fpExternalNav->ComputeStep(fLastLocatedPointLocal,
937 localDirection,
938 pCurrentProposedStepLength,
939 pNewSafety,
940 fHistory,
941 fValidExitNormal,
942 fExitNormal,
943 fExiting,
944 fEntering,
945 &fBlockedPhysicalVolume,
946 fBlockedReplicaNo);
947 break;
948 }
949 }
950 else
951 {
952 // In the case of a replica, it must handle the exiting
953 // edge/corner problem by itself
954 //
955 fExiting = fExitedMother;
956 Step = freplicaNav.ComputeStep(pGlobalpoint,
957 pDirection,
958 fLastLocatedPointLocal,
959 localDirection,
960 pCurrentProposedStepLength,
961 pNewSafety,
962 fHistory,
963 fValidExitNormal,
964 fCalculatedExitNormal,
965 fExitNormal,
966 fExiting,
967 fEntering,
968 &fBlockedPhysicalVolume,
969 fBlockedReplicaNo);
970 }
971
972 // Remember last safety origin & value.
973 //
974 fPreviousSftOrigin = pGlobalpoint;
975 fPreviousSafety = pNewSafety;
976
977 // Count zero steps - one can occur due to changing momentum at a boundary
978 // - one, two (or a few) can occur at common edges between
979 // volumes
980 // - more than two is likely a problem in the geometry
981 // description or the Navigation
982
983 // Rule of thumb: likely at an Edge if two consecutive steps are zero,
984 // because at least two candidate volumes must have been
985 // checked
986 //
987 fLocatedOnEdge = fLastStepWasZero && (Step==0.0);
988 fLastStepWasZero = (Step<fMinStep);
989 if (fPushed) { fPushed = fLastStepWasZero; }
990
991 // Handle large number of consecutive zero steps
992 //
993 if ( fLastStepWasZero )
994 {
995 ++fNumberZeroSteps;
996
997 G4bool act = fNumberZeroSteps >= fActionThreshold_NoZeroSteps;
998 G4bool actAndReport = false;
999 G4bool abandon = fNumberZeroSteps >= fAbandonThreshold_NoZeroSteps;
1000 G4bool inform = false;
1001 #ifdef G4VERBOSE
1002 actAndReport = act && (!fPushed) && fWarnPush;
1003 #endif
1004 #ifdef G4DEBUG_NAVIGATION
1005 inform = fNumberZeroSteps > 1;
1006 #endif
1007
1008 if ( act || inform )
1009 {
1010 if( act && !abandon )
1011 {
1012 // Act to recover this stuck track. Pushing it along original direction
1013 //
1014 Step += 100*kCarTolerance;
1015 fPushed = true;
1016 }
1017
1018 if( actAndReport || abandon || inform )
1019 {
1020 std::ostringstream message;
1021
1022 message.precision(16);
1023 message << "Stuck Track: potential geometry or navigation problem."
1024 << G4endl;
1025 message << " Track stuck, not moving for "
1026 << fNumberZeroSteps << " steps." << G4endl
1027 << " Current phys volume: '" << motherPhysical->GetName()
1028 << "'" << G4endl
1029 << " - at position : " << pGlobalpoint << G4endl
1030 << " in direction: " << pDirection << G4endl
1031 << " (local position: " << newLocalPoint << ")" << G4endl
1032 << " (local direction: " << localDirection << ")." << G4endl
1033 << " Previous phys volume: '"
1034 << ( fLastMotherPhys != nullptr ? fLastMotherPhys->GetName() : G4String("") )
1035 << "'" << G4endl << G4endl;
1036 if( actAndReport || abandon )
1037 {
1038 message << " Likely geometry overlap - else navigation problem !"
1039 << G4endl;
1040 }
1041 if( abandon ) // i.e. fNumberZeroSteps >= fAbandonThreshold_NoZeroSteps
1042 {
1043 // Must kill this stuck track
1044 #ifdef G4VERBOSE
1045 if ( fWarnPush ) { CheckOverlapsIterative(motherPhysical); }
1046 #endif
1047 message << " Track *abandoned* due to excessive number of Zero steps."
1048 << " Event aborted. " << G4endl << G4endl;
1049 G4Exception("G4Navigator::ComputeStep()", "GeomNav0003",
1050 EventMustBeAborted, message);
1051 }
1052 else
1053 {
1054 #ifdef G4VERBOSE
1055 if ( actAndReport ) // (!fPushed => !wasPushed) && (fWarnPush))
1056 {
1057 message << " *** Trying to get *unstuck* using a push"
1058 << " - expanding step to " << Step << " (mm) ..."
1059 << " Potential overlap in geometry !" << G4endl;
1060 G4Exception("G4Navigator::ComputeStep()", "GeomNav1002",
1061 JustWarning, message);
1062 }
1063 #endif
1064 #ifdef G4DEBUG_NAVIGATION
1065 else
1066 {
1067 if( fNumberZeroSteps > 1 )
1068 {
1069 message << ", nav-comp-step calls # " << sNavCScalls
1070 << ", Step= " << Step << G4endl;
1071 G4cout << message.str();
1072 }
1073 }
1074 #endif
1075 } // end of else if ( abandon )
1076 } // end of if( actAndReport || abandon || inform )
1077 } // end of if ( act || inform )
1078 }
1079 else
1080 {
1081 if (!fPushed) { fNumberZeroSteps = 0; }
1082 }
1083 fLastMotherPhys = motherPhysical;
1084
1085 fEnteredDaughter = fEntering; // I expect to enter a volume in this Step
1086 fExitedMother = fExiting;
1087
1088 fStepEndPoint = pGlobalpoint
1089 + std::min(Step,pCurrentProposedStepLength) * pDirection;
1090 fLastStepEndPointLocal = fLastLocatedPointLocal + Step * localDirection;
1091
1092 if( fExiting )
1093 {
1094 #ifdef G4DEBUG_NAVIGATION
1095 if( fVerbose > 2 )
1096 {
1097 G4cout << " At G4Nav CompStep End - if(exiting) - fExiting= " << fExiting
1098 << " fValidExitNormal = " << fValidExitNormal << G4endl;
1099 G4cout << " fExitNormal= " << fExitNormal << G4endl;
1100 }
1101 #endif
1102
1103 if ( fValidExitNormal || fCalculatedExitNormal )
1104 {
1105 // Convention: fExitNormal is in the 'grand-mother' coordinate system
1106 fGrandMotherExitNormal = fExitNormal;
1107 }
1108 else
1109 {
1110 // We must calculate the normal anyway (in order to have it if requested)
1111 //
1112 G4ThreeVector finalLocalPoint = fLastLocatedPointLocal
1113 + localDirection*Step;
1114
1115 if ( fHistory.GetTopVolumeType() != kReplica )
1116 {
1117 // Find normal in the 'mother' coordinate system
1118 //
1119 G4ThreeVector exitNormalMotherFrame=
1120 motherLogical->GetSolid()->SurfaceNormal(finalLocalPoint);
1121
1122 // Transform it to the 'grand-mother' coordinate system
1123 //
1124 const G4RotationMatrix* mRot = motherPhysical->GetRotation();
1125 if( mRot != nullptr )
1126 {
1127 fChangedGrandMotherRefFrame = true;
1128 fGrandMotherExitNormal = (*mRot).inverse() * exitNormalMotherFrame;
1129 }
1130 else
1131 {
1132 fGrandMotherExitNormal = exitNormalMotherFrame;
1133 }
1134
1135 // Do not set fValidExitNormal -- this signifies
1136 // that the solid is convex!
1137 }
1138 else
1139 {
1140 fCalculatedExitNormal = false;
1141 //
1142 // Nothing can be done at this stage currently - to solve this
1143 // Replica Navigation must have calculated the normal for this case
1144 // already.
1145 // Cases: mother is not convex, and exit is at previous replica level
1146
1147 #ifdef G4DEBUG_NAVIGATION
1148 G4ExceptionDescription desc;
1149
1150 desc << "Problem in ComputeStep: Replica Navigation did not provide"
1151 << " valid exit Normal. " << G4endl;
1152 desc << " Do not know how calculate it in this case." << G4endl;
1153 desc << " Location = " << finalLocalPoint << G4endl;
1154 desc << " Volume name = " << motherPhysical->GetName()
1155 << " copy/replica No = " << motherPhysical->GetCopyNo() << G4endl;
1156 G4Exception("G4Navigator::ComputeStep()", "GeomNav0003",
1157 JustWarning, desc, "Normal not available for exiting.");
1158 #endif
1159 }
1160 }
1161
1162 if ( fHistory.GetTopVolumeType() != kReplica )
1163 {
1164 fCalculatedExitNormal = true;
1165 }
1166
1167 // Now transform it to the global reference frame !!
1168 //
1169 if( fValidExitNormal || fCalculatedExitNormal )
1170 {
1171 auto depth = (G4int)fHistory.GetDepth();
1172 if( depth > 0 )
1173 {
1174 fExitNormalGlobalFrame = fHistory.GetTransform(depth-1)
1175 .InverseTransformAxis( fGrandMotherExitNormal );
1176 }
1177 else
1178 {
1179 fExitNormalGlobalFrame = fGrandMotherExitNormal;
1180 }
1181 }
1182 else
1183 {
1184 fExitNormalGlobalFrame = G4ThreeVector( 0., 0., 0.);
1185 }
1186 }
1187
1188 if( (Step == pCurrentProposedStepLength) && (!fExiting) && (!fEntering) )
1189 {
1190 // This if Step is not really limited by the geometry.
1191 // The Navigator is obliged to return "infinity"
1192 //
1193 Step = kInfinity;
1194 }
1195
1196 #ifdef G4VERBOSE
1197 if( fVerbose > 1 )
1198 {
1199 if( fVerbose >= 4 )
1200 {
1201 G4cout << " ----- Upon exiting :" << G4endl;
1202 PrintState();
1203 }
1204 G4cout << " Returned step= " << Step;
1205 if( fVerbose > 5 ) { G4cout << G4endl; }
1206 if( Step == kInfinity )
1207 {
1208 G4cout << " Requested step= " << pCurrentProposedStepLength ;
1209 if( fVerbose > 5) { G4cout << G4endl; }
1210 }
1211 G4cout << " Safety = " << pNewSafety << G4endl;
1212 }
1213 #endif
1214
1215 fLastTriedStepComputation = true;
1216
1217 return Step;
1218 }
1219
1220 // ********************************************************************
1221 // CheckNextStep
1222 //
1223 // Compute the step without altering the navigator state
1224 // ********************************************************************
1225 //
1226 G4double G4Navigator::CheckNextStep( const G4ThreeVector& pGlobalpoint,
1227 const G4ThreeVector& pDirection,
1228 const G4double pCurrentProposedStepLength,
1229 G4double& pNewSafety)
1230 {
1231 G4double step;
1232
1233 // Save the state, for this parasitic call
1234 //
1235 SetSavedState();
1236
1237 step = ComputeStep ( pGlobalpoint,
1238 pDirection,
1239 pCurrentProposedStepLength,
1240 pNewSafety );
1241
1242 // It is a parasitic call, so attempt to restore the key parts of the state
1243 //
1244 RestoreSavedState();
1245 // NOTE: the state of the current subnavigator is NOT restored.
1246 // ***> TODO: restore subnavigator state
1247 // if( last_located) Need Position of last location
1248 // if( last_computed step) Need Endposition of last step
1249
1250 return step;
1251 }
1252
1253 // ********************************************************************
1254 // ResetState
1255 //
1256 // Resets stack and minimum of navigator state `machine'
1257 // ********************************************************************
1258 //
1259 void G4Navigator::ResetState()
1260 {
1261 fWasLimitedByGeometry = false;
1262 fEntering = false;
1263 fExiting = false;
1264 fLocatedOnEdge = false;
1265 fLastStepWasZero = false;
1266 fEnteredDaughter = false;
1267 fExitedMother = false;
1268 fPushed = false;
1269
1270 fValidExitNormal = false;
1271 fChangedGrandMotherRefFrame = false;
1272 fCalculatedExitNormal = false;
1273
1274 fExitNormal = G4ThreeVector(0,0,0);
1275 fGrandMotherExitNormal = G4ThreeVector(0,0,0);
1276 fExitNormalGlobalFrame = G4ThreeVector(0,0,0);
1277
1278 fPreviousSftOrigin = G4ThreeVector(0,0,0);
1279 fPreviousSafety = 0.0;
1280
1281 fNumberZeroSteps = 0;
1282
1283 fBlockedPhysicalVolume = nullptr;
1284 fBlockedReplicaNo = -1;
1285
1286 fLastLocatedPointLocal = G4ThreeVector( kInfinity, -kInfinity, 0.0 );
1287 fLocatedOutsideWorld = false;
1288
1289 fLastMotherPhys = nullptr;
1290 }
1291
1292 // ********************************************************************
1293 // SetupHierarchy
1294 //
1295 // Renavigates & resets hierarchy described by current history
1296 // o Reset volumes
1297 // o Recompute transforms and/or solids of replicated/parameterised volumes
1298 // ********************************************************************
1299 //
1300 void G4Navigator::SetupHierarchy()
1301 {
1302 const auto depth = (G4int)fHistory.GetDepth();
1303 for ( auto i = 1; i <= depth; ++i )
1304 {
1305 switch ( fHistory.GetVolumeType(i) )
1306 {
1307 case kNormal:
1308 case kExternal:
1309 break;
1310 case kReplica:
1311 freplicaNav.ComputeTransformation(fHistory.GetReplicaNo(i), fHistory.GetVolume(i));
1312 break;
1313 case kParameterised:
1314 G4VPhysicalVolume* current = fHistory.GetVolume(i);
1315 G4int replicaNo = fHistory.GetReplicaNo(i);
1316 G4VPVParameterisation* pParam = current->GetParameterisation();
1317 G4VSolid* pSolid = pParam->ComputeSolid(replicaNo, current);
1318
1319 // Set up dimensions & transform in solid/physical volume
1320 //
1321 pSolid->ComputeDimensions(pParam, replicaNo, current);
1322 pParam->ComputeTransformation(replicaNo, current);
1323
1324 G4TouchableHistory* pTouchable = nullptr;
1325 if( pParam->IsNested() )
1326 {
1327 pTouchable= new G4TouchableHistory( fHistory );
1328 pTouchable->MoveUpHistory(); // Move up to the parent level
1329 // Adequate only if Nested at the Branch level (last)
1330 // To extend to other cases:
1331 // pTouchable->MoveUpHistory(cdepth-i-1);
1332 // Move to the parent level of *Current* level
1333 // Could replace this line and constructor with a revised
1334 // c-tor for History(levels to drop)
1335 }
1336 // Set up the correct solid and material in Logical Volume
1337 //
1338 G4LogicalVolume* pLogical = current->GetLogicalVolume();
1339 pLogical->SetSolid( pSolid );
1340 pLogical->UpdateMaterial( pParam ->
1341 ComputeMaterial(replicaNo, current, pTouchable) );
1342 delete pTouchable;
1343 break;
1344 }
1345 }
1346 }
1347
1348 // ********************************************************************
1349 // GetLocalExitNormal
1350 //
1351 // Obtains the Normal vector to a surface (in local coordinates)
1352 // pointing out of previous volume and into current volume
1353 // ********************************************************************
1354 //
1355 G4ThreeVector G4Navigator::GetLocalExitNormal( G4bool* valid )
1356 {
1357 G4ThreeVector ExitNormal(0.,0.,0.);
1358 G4VSolid* currentSolid = nullptr;
1359 G4LogicalVolume* candidateLogical;
1360
1361 if ( fLastTriedStepComputation )
1362 {
1363 // use fLastLocatedPointLocal and next candidate volume
1364 //
1365 G4ThreeVector nextSolidExitNormal(0.,0.,0.);
1366
1367 if( fEntering && (fBlockedPhysicalVolume!=nullptr) )
1368 {
1369 candidateLogical = fBlockedPhysicalVolume->GetLogicalVolume();
1370 if( candidateLogical != nullptr )
1371 {
1372 // fLastStepEndPointLocal is in the coordinates of the mother
1373 // we need it in the daughter's coordinate system.
1374
1375 // The following code should also work in case of Replica
1376 {
1377 // First transform fLastLocatedPointLocal to the new daughter
1378 // coordinates
1379 //
1380 G4AffineTransform MotherToDaughterTransform=
1381 GetMotherToDaughterTransform( fBlockedPhysicalVolume,
1382 fBlockedReplicaNo,
1383 VolumeType(fBlockedPhysicalVolume) );
1384 G4ThreeVector daughterPointOwnLocal =
1385 MotherToDaughterTransform.TransformPoint( fLastStepEndPointLocal );
1386
1387 // OK if it is a parameterised volume
1388 //
1389 EInside inSideIt;
1390 G4bool onSurface;
1391 G4double safety = -1.0;
1392 currentSolid = candidateLogical->GetSolid();
1393 inSideIt = currentSolid->Inside(daughterPointOwnLocal);
1394 onSurface = (inSideIt == kSurface);
1395 if( !onSurface )
1396 {
1397 if( inSideIt == kOutside )
1398 {
1399 safety = (currentSolid->DistanceToIn(daughterPointOwnLocal));
1400 onSurface = safety < 100.0 * kCarTolerance;
1401 }
1402 else if (inSideIt == kInside )
1403 {
1404 safety = (currentSolid->DistanceToOut(daughterPointOwnLocal));
1405 onSurface = safety < 100.0 * kCarTolerance;
1406 }
1407 }
1408
1409 if( onSurface )
1410 {
1411 nextSolidExitNormal =
1412 currentSolid->SurfaceNormal(daughterPointOwnLocal);
1413
1414 // Entering the solid ==> opposite
1415 //
1416 // First flip ( ExitNormal = -nextSolidExitNormal; )
1417 // and then rotate the the normal to the frame of the mother (current volume)
1418 ExitNormal = MotherToDaughterTransform
1419 .InverseTransformAxis( -nextSolidExitNormal );
1420 fCalculatedExitNormal = true;
1421 }
1422 else
1423 {
1424 #ifdef G4VERBOSE
1425 if(( fVerbose == 1 ) && ( fCheck ))
1426 {
1427 std::ostringstream message;
1428 message << "Point not on surface ! " << G4endl
1429 << " Point = "
1430 << daughterPointOwnLocal << G4endl
1431 << " Physical volume = "
1432 << fBlockedPhysicalVolume->GetName() << G4endl
1433 << " Logical volume = "
1434 << candidateLogical->GetName() << G4endl
1435 << " Solid = " << currentSolid->GetName()
1436 << " Type = "
1437 << currentSolid->GetEntityType() << G4endl
1438 << *currentSolid << G4endl;
1439 if( inSideIt == kOutside )
1440 {
1441 message << "Point is Outside. " << G4endl
1442 << " Safety (from outside) = " << safety << G4endl;
1443 }
1444 else // if( inSideIt == kInside )
1445 {
1446 message << "Point is Inside. " << G4endl
1447 << " Safety (from inside) = " << safety << G4endl;
1448 }
1449 G4Exception("G4Navigator::GetLocalExitNormal()", "GeomNav1001",
1450 JustWarning, message);
1451 }
1452 #endif
1453 }
1454 *valid = onSurface; // was =true;
1455 }
1456 }
1457 }
1458 else if ( fExiting )
1459 {
1460 ExitNormal = fGrandMotherExitNormal;
1461 *valid = true;
1462 fCalculatedExitNormal = true; // Should be true already
1463 }
1464 else // i.e. ( fBlockedPhysicalVolume == 0 )
1465 {
1466 *valid = false;
1467 G4Exception("G4Navigator::GetLocalExitNormal()",
1468 "GeomNav0003", JustWarning,
1469 "Incorrect call to GetLocalSurfaceNormal." );
1470 }
1471 }
1472 else // ( ! fLastTriedStepComputation ) i.e. last call was to Locate
1473 {
1474 if ( EnteredDaughterVolume() )
1475 {
1476 G4VSolid* daughterSolid = fHistory.GetTopVolume()->GetLogicalVolume()
1477 ->GetSolid();
1478 ExitNormal = -(daughterSolid->SurfaceNormal(fLastLocatedPointLocal));
1479 if( std::fabs(ExitNormal.mag2()-1.0 ) > kToleranceNormalCheck )
1480 {
1481 G4ExceptionDescription desc;
1482 desc << " Parameters of solid: " << *daughterSolid
1483 << " Point for surface = " << fLastLocatedPointLocal << std::endl;
1484 G4Exception("G4Navigator::GetLocalExitNormal()",
1485 "GeomNav0003", FatalException, desc,
1486 "Surface Normal returned by Solid is not a Unit Vector." );
1487 }
1488 fCalculatedExitNormal = true;
1489 *valid = true;
1490 }
1491 else
1492 {
1493 if( fExitedMother )
1494 {
1495 ExitNormal = fGrandMotherExitNormal;
1496 *valid = true;
1497 fCalculatedExitNormal = true;
1498 }
1499 else // We are not at a boundary. ExitNormal remains (0,0,0)
1500 {
1501 *valid = false;
1502 fCalculatedExitNormal = false;
1503 G4ExceptionDescription message;
1504 message << "Function called when *NOT* at a Boundary." << G4endl;
1505 message << "Exit Normal not calculated." << G4endl;
1506 G4Exception("G4Navigator::GetLocalExitNormal()",
1507 "GeomNav0003", JustWarning, message);
1508 }
1509 }
1510 }
1511 return ExitNormal;
1512 }
1513
1514 // ********************************************************************
1515 // GetMotherToDaughterTransform
1516 //
1517 // Obtains the mother to daughter affine transformation
1518 // ********************************************************************
1519 //
1520 G4AffineTransform
1521 G4Navigator::GetMotherToDaughterTransform( G4VPhysicalVolume* pEnteringPhysVol,
1522 G4int enteringReplicaNo,
1523 EVolume enteringVolumeType )
1524 {
1525 switch (enteringVolumeType)
1526 {
1527 case kNormal: // Nothing is needed to prepare the transformation
1528 break; // It is stored already in the physical volume (placement)
1529 case kReplica: // Sets the transform in the Replica - tbc
1530 G4Exception("G4Navigator::GetMotherToDaughterTransform()",
1531 "GeomNav0001", FatalException,
1532 "Method NOT Implemented yet for replica volumes.");
1533 break;
1534 case kParameterised:
1535 if( pEnteringPhysVol->GetRegularStructureId() == 0 )
1536 {
1537 G4VPVParameterisation *pParam =
1538 pEnteringPhysVol->GetParameterisation();
1539 G4VSolid* pSolid =
1540 pParam->ComputeSolid(enteringReplicaNo, pEnteringPhysVol);
1541 pSolid->ComputeDimensions(pParam, enteringReplicaNo, pEnteringPhysVol);
1542
1543 // Sets the transform in the Parameterisation
1544 //
1545 pParam->ComputeTransformation(enteringReplicaNo, pEnteringPhysVol);
1546
1547 // Set the correct solid and material in Logical Volume
1548 //
1549 G4LogicalVolume* pLogical = pEnteringPhysVol->GetLogicalVolume();
1550 pLogical->SetSolid( pSolid );
1551 }
1552 break;
1553 case kExternal:
1554 // Expect that nothing is needed to prepare the transformation.
1555 // It is stored already in the physical volume (placement)
1556 break;
1557 }
1558 return G4AffineTransform(pEnteringPhysVol->GetRotation(),
1559 pEnteringPhysVol->GetTranslation()).Invert();
1560 }
1561
1562
1563 // ********************************************************************
1564 // GetLocalExitNormalAndCheck
1565 //
1566 // Obtains the Normal vector to a surface (in local coordinates)
1567 // pointing out of previous volume and into current volume, and
1568 // checks the current point against expected 'local' value.
1569 // ********************************************************************
1570 //
1571 G4ThreeVector
1572 G4Navigator::GetLocalExitNormalAndCheck(
1573 #ifdef G4DEBUG_NAVIGATION
1574 const G4ThreeVector& ExpectedBoundaryPointGlobal,
1575 #else
1576 const G4ThreeVector&,
1577 #endif
1578 G4bool* pValid)
1579 {
1580 #ifdef G4DEBUG_NAVIGATION
1581 // Check Current point against expected 'local' value
1582 //
1583 if ( fLastTriedStepComputation )
1584 {
1585 G4ThreeVector ExpectedBoundaryPointLocal;
1586
1587 const G4AffineTransform& GlobalToLocal = GetGlobalToLocalTransform();
1588 ExpectedBoundaryPointLocal =
1589 GlobalToLocal.TransformPoint( ExpectedBoundaryPointGlobal );
1590
1591 // Add here: Comparison against expected position,
1592 // i.e. the endpoint of ComputeStep
1593 }
1594 #endif
1595
1596 return GetLocalExitNormal( pValid );
1597 }
1598
1599 // ********************************************************************
1600 // GetGlobalExitNormal
1601 //
1602 // Obtains the Normal vector to a surface (in global coordinates)
1603 // pointing out of previous volume and into current volume
1604 // ********************************************************************
1605 //
1606 G4ThreeVector
1607 G4Navigator::GetGlobalExitNormal(const G4ThreeVector& IntersectPointGlobal,
1608 G4bool* pNormalCalculated)
1609 {
1610 G4bool validNormal;
1611 G4ThreeVector localNormal, globalNormal;
1612
1613 G4bool usingStored = fCalculatedExitNormal && (
1614 ( fLastTriedStepComputation && fExiting ) // Just calculated it
1615 || // No locate in between
1616 ( !fLastTriedStepComputation
1617 && (IntersectPointGlobal-fStepEndPoint).mag2() < 10.0*fSqTol ) );
1618 // Calculated it 'just' before & then called locate
1619 // but it did not move position
1620
1621 if( usingStored )
1622 {
1623 // This was computed in last call to ComputeStep
1624 // and only if it arrived at boundary
1625 //
1626 globalNormal = fExitNormalGlobalFrame;
1627 G4double normMag2 = globalNormal.mag2();
1628 if( std::fabs ( normMag2 - 1.0 ) < perThousand ) // was perMillion
1629 {
1630 *pNormalCalculated = true; // ComputeStep always computes it if Exiting
1631 // (fExiting==true)
1632 }
1633 else
1634 {
1635 G4ExceptionDescription message;
1636 message.precision(10);
1637 message << " WARNING> Expected normal-global-frame to be valid, "
1638 << " i.e. a unit vector!" << G4endl
1639 << " - but |normal| = " << std::sqrt(normMag2)
1640 << " - and |normal|^2 = " << normMag2 << G4endl
1641 << " which differs from 1.0 by " << normMag2 - 1.0 << G4endl
1642 << " n = " << fExitNormalGlobalFrame << G4endl
1643 << " Global point: " << IntersectPointGlobal << G4endl
1644 << " Volume: " << fHistory.GetTopVolume()->GetName() << G4endl;
1645 #ifdef G4VERBOSE
1646 G4LogicalVolume* candLog = fHistory.GetTopVolume()->GetLogicalVolume();
1647 if ( candLog != nullptr )
1648 {
1649 message << " Solid: " << candLog->GetSolid()->GetName()
1650 << ", Type: " << candLog->GetSolid()->GetEntityType() << G4endl
1651 << *candLog->GetSolid() << G4endl;
1652 }
1653 #endif
1654 message << "============================================================"
1655 << G4endl;
1656 G4int oldVerbose = fVerbose;
1657 fVerbose = 4;
1658 message << " State of Navigator: " << G4endl;
1659 message << *this << G4endl;
1660 fVerbose = oldVerbose;
1661 message << "============================================================"
1662 << G4endl;
1663
1664 G4Exception("G4Navigator::GetGlobalExitNormal()",
1665 "GeomNav0003",JustWarning, message,
1666 "Value obtained from stored global-normal is not a unit vector.");
1667
1668 // (Re)Compute it now -- as either it was not computed, or it is wrong.
1669 //
1670 localNormal = GetLocalExitNormalAndCheck(IntersectPointGlobal,
1671 &validNormal);
1672 *pNormalCalculated = fCalculatedExitNormal;
1673 globalNormal = fHistory.GetTopTransform()
1674 .InverseTransformAxis(localNormal);
1675 }
1676 }
1677 else
1678 {
1679 localNormal = GetLocalExitNormalAndCheck(IntersectPointGlobal,&validNormal);
1680 *pNormalCalculated = fCalculatedExitNormal;
1681
1682 #ifdef G4DEBUG_NAVIGATION
1683 usingStored = false;
1684
1685 if( (!validNormal) && !fCalculatedExitNormal )
1686 {
1687 G4ExceptionDescription edN;
1688 edN << " Calculated = " << fCalculatedExitNormal << G4endl;
1689 edN << " Entering= " << fEntering << G4endl;
1690 G4int oldVerbose = this->GetVerboseLevel();
1691 this->SetVerboseLevel(4);
1692 edN << " State of Navigator: " << G4endl;
1693 edN << *this << G4endl;
1694 this->SetVerboseLevel( oldVerbose );
1695
1696 G4Exception("G4Navigator::GetGlobalExitNormal()",
1697 "GeomNav0003", JustWarning, edN,
1698 "LocalExitNormalAndCheck() did not calculate Normal.");
1699 }
1700 #endif
1701
1702 G4double localMag2 = localNormal.mag2();
1703 if( validNormal && (std::fabs(localMag2-1.0)) > kToleranceNormalCheck )
1704 {
1705 G4ExceptionDescription edN;
1706 edN.precision(10);
1707 edN << "G4Navigator::GetGlobalExitNormal: "
1708 << " Using Local Normal - from call to GetLocalExitNormalAndCheck. "
1709 << G4endl
1710 << " Local Exit Normal : " << " || = " << std::sqrt(localMag2)
1711 << " vec = " << localNormal << G4endl
1712 << " Global Exit Normal : " << " || = " << globalNormal.mag()
1713 << " vec = " << globalNormal << G4endl
1714 << " Global point: " << IntersectPointGlobal << G4endl;
1715 edN << " Calculated It = " << fCalculatedExitNormal << G4endl
1716 << " Volume: " << fHistory.GetTopVolume()->GetName() << G4endl;
1717 #ifdef G4VERBOSE
1718 G4LogicalVolume* candLog = fHistory.GetTopVolume()->GetLogicalVolume();
1719 if ( candLog != nullptr )
1720 {
1721 edN << " Solid: " << candLog->GetSolid()->GetName()
1722 << ", Type: " << candLog->GetSolid()->GetEntityType() << G4endl
1723 << *candLog->GetSolid();
1724 }
1725 #endif
1726 G4Exception("G4Navigator::GetGlobalExitNormal()",
1727 "GeomNav0003",JustWarning, edN,
1728 "Value obtained from new local *solid* is incorrect.");
1729 localNormal = localNormal.unit(); // Should we correct it ??
1730 }
1731 globalNormal = fHistory.GetTopTransform()
1732 .InverseTransformAxis(localNormal);
1733 }
1734
1735 #ifdef G4DEBUG_NAVIGATION
1736 if( usingStored )
1737 {
1738 G4ThreeVector globalNormAgn;
1739
1740 localNormal = GetLocalExitNormalAndCheck(IntersectPointGlobal,&validNormal);
1741
1742 globalNormAgn = fHistory.GetTopTransform()
1743 .InverseTransformAxis(localNormal);
1744
1745 // Check the value computed against fExitNormalGlobalFrame
1746 G4ThreeVector diffNorm = globalNormAgn - fExitNormalGlobalFrame;
1747 if( diffNorm.mag2() > kToleranceNormalCheck )
1748 {
1749 G4ExceptionDescription edDfn;
1750 edDfn << "Found difference in normals in case of exiting mother "
1751 << "- when Get is called after ComputingStep " << G4endl;
1752 edDfn << " Magnitude of diff = " << diffNorm.mag() << G4endl;
1753 edDfn << " Normal stored (Global) = " << fExitNormalGlobalFrame
1754 << G4endl;
1755 edDfn << " Global Computed from Local = " << globalNormAgn << G4endl;
1756 G4Exception("G4Navigator::GetGlobalExitNormal()", "GeomNav0003",
1757 JustWarning, edDfn);
1758 }
1759 }
1760 #endif
1761
1762 // Synchronise stored global exit normal as possibly re-computed here
1763 //
1764 fExitNormalGlobalFrame = globalNormal;
1765
1766 return globalNormal;
1767 }
1768
1769 // ********************************************************************
1770 // ComputeSafety
1771 //
1772 // It assumes that it will be
1773 // i) called at the Point in the same volume as the EndPoint of the
1774 // ComputeStep.
1775 // ii) after (or at the end of) ComputeStep OR after the relocation.
1776 // ********************************************************************
1777 //
1778 G4double G4Navigator::ComputeSafety( const G4ThreeVector& pGlobalpoint,
1779 const G4double pMaxLength,
1780 const G4bool )
1781 {
1782 G4VPhysicalVolume *motherPhysical = fHistory.GetTopVolume();
1783 G4double safety = 0.0;
1784
1785 G4double distEndpointSq = (pGlobalpoint-fStepEndPoint).mag2();
1786 G4bool stayedOnEndpoint = distEndpointSq < sqr(kCarTolerance);
1787 G4bool endpointOnSurface = fEnteredDaughter || fExitedMother;
1788
1789 G4bool onSurface = endpointOnSurface && stayedOnEndpoint;
1790 if( ! onSurface )
1791 {
1792 safety= fpSafetyCalculator->SafetyInCurrentVolume(pGlobalpoint, motherPhysical, pMaxLength);
1793 // offload to G4SafetyCalculator - avoids need to save / reload state
1794
1795 // Remember last safety origin & value
1796 //
1797 fPreviousSftOrigin = pGlobalpoint;
1798 fPreviousSafety = safety;
1799 // We overwrite the Safety 'sphere' - keeping old behaviour
1800 }
1801
1802 return safety;
1803 }
1804
1805 // ********************************************************************
1806 // CreateTouchableHistoryHandle
1807 // ********************************************************************
1808 //
1809 G4TouchableHandle G4Navigator::CreateTouchableHistoryHandle() const
1810 {
1811 return G4TouchableHandle( CreateTouchableHistory() );
1812 }
1813
1814 // ********************************************************************
1815 // PrintState
1816 // ********************************************************************
1817 //
1818 void G4Navigator::PrintState() const
1819 {
1820 G4long oldcoutPrec = G4cout.precision(4);
1821 if( fVerbose >= 4 )
1822 {
1823 G4cout << "The current state of G4Navigator is: " << G4endl;
1824 G4cout << " ValidExitNormal= " << fValidExitNormal // << G4endl
1825 << " ExitNormal = " << fExitNormal // << G4endl
1826 << " Exiting = " << fExiting // << G4endl
1827 << " Entering = " << fEntering // << G4endl
1828 << " BlockedPhysicalVolume= " ;
1829 if (fBlockedPhysicalVolume==nullptr)
1830 {
1831 G4cout << "None";
1832 }
1833 else
1834 {
1835 G4cout << fBlockedPhysicalVolume->GetName();
1836 }
1837 G4cout << G4endl
1838 << " BlockedReplicaNo = " << fBlockedReplicaNo // << G4endl
1839 << " LastStepWasZero = " << fLastStepWasZero // << G4endl
1840 << G4endl;
1841 }
1842 if( ( 1 < fVerbose) && (fVerbose < 4) )
1843 {
1844 G4cout << G4endl; // Make sure to line up
1845 G4cout << std::setw(30) << " ExitNormal " << " "
1846 << std::setw( 5) << " Valid " << " "
1847 << std::setw( 9) << " Exiting " << " "
1848 << std::setw( 9) << " Entering" << " "
1849 << std::setw(15) << " Blocked:Volume " << " "
1850 << std::setw( 9) << " ReplicaNo" << " "
1851 << std::setw( 8) << " LastStepZero " << " "
1852 << G4endl;
1853 G4cout << "( " << std::setw(7) << fExitNormal.x()
1854 << ", " << std::setw(7) << fExitNormal.y()
1855 << ", " << std::setw(7) << fExitNormal.z() << " ) "
1856 << std::setw( 5) << fValidExitNormal << " "
1857 << std::setw( 9) << fExiting << " "
1858 << std::setw( 9) << fEntering << " ";
1859 if ( fBlockedPhysicalVolume == nullptr )
1860 { G4cout << std::setw(15) << "None"; }
1861 else
1862 { G4cout << std::setw(15)<< fBlockedPhysicalVolume->GetName(); }
1863 G4cout << std::setw( 9) << fBlockedReplicaNo << " "
1864 << std::setw( 8) << fLastStepWasZero << " "
1865 << G4endl;
1866 }
1867 if( fVerbose > 2 )
1868 {
1869 G4cout.precision(8);
1870 G4cout << " Current Localpoint = " << fLastLocatedPointLocal << G4endl;
1871 G4cout << " PreviousSftOrigin = " << fPreviousSftOrigin << G4endl;
1872 G4cout << " PreviousSafety = " << fPreviousSafety << G4endl;
1873 }
1874 G4cout.precision(oldcoutPrec);
1875 }
1876
1877 // ********************************************************************
1878 // ComputeStepLog
1879 // ********************************************************************
1880 //
1881 void G4Navigator::ComputeStepLog(const G4ThreeVector& pGlobalpoint,
1882 G4double moveLenSq) const
1883 {
1884 // The following checks only make sense if the move is larger
1885 // than the tolerance.
1886
1887 const G4double fAccuracyForWarning = kCarTolerance,
1888 fAccuracyForException = 1000*kCarTolerance;
1889
1890 G4ThreeVector OriginalGlobalpoint = fHistory.GetTopTransform().
1891 InverseTransformPoint(fLastLocatedPointLocal);
1892
1893 G4double shiftOriginSafSq = (fPreviousSftOrigin-pGlobalpoint).mag2();
1894
1895 // Check that the starting point of this step is
1896 // within the isotropic safety sphere of the last point
1897 // to a accuracy/precision given by fAccuracyForWarning.
1898 // If so give warning.
1899 // If it fails by more than fAccuracyForException exit with error.
1900 //
1901 if( shiftOriginSafSq >= sqr(fPreviousSafety) )
1902 {
1903 G4double shiftOrigin = std::sqrt(shiftOriginSafSq);
1904 G4double diffShiftSaf = shiftOrigin - fPreviousSafety;
1905
1906 if( diffShiftSaf > fAccuracyForWarning )
1907 {
1908 G4long oldcoutPrec = G4cout.precision(8);
1909 G4long oldcerrPrec = G4cerr.precision(10);
1910 std::ostringstream message, suggestion;
1911 message << "Accuracy error or slightly inaccurate position shift."
1912 << G4endl
1913 << " The Step's starting point has moved "
1914 << std::sqrt(moveLenSq)/mm << " mm " << G4endl
1915 << " since the last call to a Locate method." << G4endl
1916 << " This has resulted in moving "
1917 << shiftOrigin/mm << " mm "
1918 << " from the last point at which the safety "
1919 << " was calculated " << G4endl
1920 << " which is more than the computed safety= "
1921 << fPreviousSafety/mm << " mm at that point." << G4endl
1922 << " This difference is "
1923 << diffShiftSaf/mm << " mm." << G4endl
1924 << " The tolerated accuracy is "
1925 << fAccuracyForException/mm << " mm.";
1926
1927 suggestion << " ";
1928 static G4ThreadLocal G4int warnNow = 0;
1929 if( ((++warnNow % 100) == 1) )
1930 {
1931 message << G4endl
1932 << " This problem can be due to either " << G4endl
1933 << " - a process that has proposed a displacement"
1934 << " larger than the current safety , or" << G4endl
1935 << " - inaccuracy in the computation of the safety";
1936 suggestion << "We suggest that you " << G4endl
1937 << " - find i) what particle is being tracked, and "
1938 << " ii) through what part of your geometry " << G4endl
1939 << " for example by re-running this event with "
1940 << G4endl
1941 << " /tracking/verbose 1 " << G4endl
1942 << " - check which processes you declare for"
1943 << " this particle (and look at non-standard ones)"
1944 << G4endl
1945 << " - in case, create a detailed logfile"
1946 << " of this event using:" << G4endl
1947 << " /tracking/verbose 6 ";
1948 }
1949 G4Exception("G4Navigator::ComputeStep()",
1950 "GeomNav1002", JustWarning,
1951 message, G4String(suggestion.str()));
1952 G4cout.precision(oldcoutPrec);
1953 G4cerr.precision(oldcerrPrec);
1954 }
1955 #ifdef G4DEBUG_NAVIGATION
1956 else
1957 {
1958 G4cerr << "WARNING - G4Navigator::ComputeStep()" << G4endl
1959 << " The Step's starting point has moved "
1960 << std::sqrt(moveLenSq) << "," << G4endl
1961 << " which has taken it to the limit of"
1962 << " the current safety. " << G4endl;
1963 }
1964 #endif
1965 }
1966 G4double safetyPlus = fPreviousSafety + fAccuracyForException;
1967 if ( shiftOriginSafSq > sqr(safetyPlus) )
1968 {
1969 std::ostringstream message;
1970 message << "May lead to a crash or unreliable results." << G4endl
1971 << " Position has shifted considerably without"
1972 << " notifying the navigator !" << G4endl
1973 << " Tolerated safety: " << safetyPlus << G4endl
1974 << " Computed shift : " << shiftOriginSafSq;
1975 G4Exception("G4Navigator::ComputeStep()", "GeomNav1002",
1976 JustWarning, message);
1977 }
1978 }
1979
1980 // ********************************************************************
1981 // CheckOverlapsIterative
1982 // ********************************************************************
1983 //
1984 G4bool G4Navigator::CheckOverlapsIterative(G4VPhysicalVolume* vol)
1985 {
1986 // Check and report overlaps
1987 //
1988 G4bool foundOverlap = false;
1989 G4int nPoints = 300000, ntrials = 9, numOverlaps = 5;
1990 G4double trialLength = 1.0 * CLHEP::centimeter;
1991 while ( ntrials-- > 0 && !foundOverlap )
1992 {
1993 if ( fVerbose > 1 )
1994 {
1995 G4cout << " ** Running overlap checks in volume "
1996 << vol->GetName()
1997 << " with length = " << trialLength << G4endl;
1998 }
1999 foundOverlap = vol->CheckOverlaps(nPoints, trialLength,
2000 fVerbose != 0, numOverlaps);
2001 trialLength *= 0.1;
2002 if ( trialLength <= 1.0e-5 ) { numOverlaps= 1;}
2003 }
2004 return foundOverlap;
2005 }
2006
2007 // ********************************************************************
2008 // Operator <<
2009 // ********************************************************************
2010 //
2011 std::ostream& operator << (std::ostream &os,const G4Navigator &n)
2012 {
2013 // Old version did only the following:
2014 // os << "Current History: " << G4endl << n.fHistory;
2015 // Old behaviour is recovered for fVerbose = 0
2016
2017 // Adapted from G4Navigator::PrintState() const
2018
2019 G4long oldcoutPrec = os.precision(4);
2020 if( n.fVerbose >= 4 )
2021 {
2022 os << "The current state of G4Navigator is: " << G4endl;
2023 os << " ValidExitNormal= " << n.fValidExitNormal << G4endl
2024 << " ExitNormal = " << n.fExitNormal << G4endl
2025 << " Exiting = " << n.fExiting << G4endl
2026 << " Entering = " << n.fEntering << G4endl
2027 << " BlockedPhysicalVolume= " ;
2028 if (n.fBlockedPhysicalVolume==nullptr)
2029 {
2030 os << "None";
2031 }
2032 else
2033 {
2034 os << n.fBlockedPhysicalVolume->GetName();
2035 }
2036 os << G4endl
2037 << " BlockedReplicaNo = " << n.fBlockedReplicaNo << G4endl
2038 << " LastStepWasZero = " << n.fLastStepWasZero << G4endl
2039 << G4endl;
2040 }
2041 if( ( 1 < n.fVerbose) && (n.fVerbose < 4) )
2042 {
2043 os << G4endl; // Make sure to line up
2044 os << std::setw(30) << " ExitNormal " << " "
2045 << std::setw( 5) << " Valid " << " "
2046 << std::setw( 9) << " Exiting " << " "
2047 << std::setw( 9) << " Entering" << " "
2048 << std::setw(15) << " Blocked:Volume " << " "
2049 << std::setw( 9) << " ReplicaNo" << " "
2050 << std::setw( 8) << " LastStepZero " << " "
2051 << G4endl;
2052 os << "( " << std::setw(7) << n.fExitNormal.x()
2053 << ", " << std::setw(7) << n.fExitNormal.y()
2054 << ", " << std::setw(7) << n.fExitNormal.z() << " ) "
2055 << std::setw( 5) << n.fValidExitNormal << " "
2056 << std::setw( 9) << n.fExiting << " "
2057 << std::setw( 9) << n.fEntering << " ";
2058 if ( n.fBlockedPhysicalVolume==nullptr )
2059 { os << std::setw(15) << "None"; }
2060 else
2061 { os << std::setw(15)<< n.fBlockedPhysicalVolume->GetName(); }
2062 os << std::setw( 9) << n.fBlockedReplicaNo << " "
2063 << std::setw( 8) << n.fLastStepWasZero << " "
2064 << G4endl;
2065 }
2066 if( n.fVerbose > 2 )
2067 {
2068 os.precision(8);
2069 os << " Current Localpoint = " << n.fLastLocatedPointLocal << G4endl;
2070 os << " PreviousSftOrigin = " << n.fPreviousSftOrigin << G4endl;
2071 os << " PreviousSafety = " << n.fPreviousSafety << G4endl;
2072 }
2073 if( n.fVerbose > 3 || n.fVerbose == 0 )
2074 {
2075 os << "Current History: " << G4endl << n.fHistory;
2076 }
2077
2078 os.precision(oldcoutPrec);
2079 return os;
2080 }
2081
2082 // ********************************************************************
2083 // SetVoxelNavigation: alternative navigator for voxelised geometry
2084 // ********************************************************************
2085 //
2086 void G4Navigator::SetVoxelNavigation(G4VoxelNavigation* voxelNav)
2087 {
2088 delete fpvoxelNav;
2089 fpvoxelNav = voxelNav;
2090 }
2091
2092 // ********************************************************************
2093 // InformLastStep: derived navigators can inform of its step
2094 // used to update fLastStepWasZero
2095 // ********************************************************************
2096 void G4Navigator::InformLastStep(G4double lastStep, G4bool entersDaughtVol,
2097 G4bool exitsMotherVol)
2098 {
2099 G4bool zeroStep = ( lastStep == 0.0 );
2100 fLocatedOnEdge = fLastStepWasZero && zeroStep;
2101 fLastStepWasZero = zeroStep;
2102
2103 fExiting = exitsMotherVol;
2104 fEntering = entersDaughtVol;
2105 }
2106
2107 // ********************************************************************
2108 // SetExternalNavigation
2109 // ********************************************************************
2110 //
2111 void G4Navigator::SetExternalNavigation(G4VExternalNavigation* externalNav)
2112 {
2113 fpExternalNav = externalNav;
2114 fpSafetyCalculator->SetExternalNavigation(externalNav);
2115 }
2116