| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // G4Navigator class Implementation
27 //
28 // Original author: Paul Kent, July 95/96
29 // Responsible 1996-present: John Apostolakis,
30 // Additional revisions by: Pedro Arce, Vladim
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 sh
45 // vectors). If exceeded, warnings will be iss
46 // Can be CLHEP::perMillion (its old default)
47 //
48 static const G4double kToleranceNormalCheck =
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
66 << fActionThreshold_NoZeroSteps
67 << " Abandon Threshold (No Zero St
68 << fAbandonThreshold_NoZeroSteps <<
69 }
70 kCarTolerance = G4GeometryTolerance::GetInst
71 fMinStep = 0.05*kCarTolerance;
72 fSqTol = sqr(kCarTolerance);
73
74 fregularNav.SetNormalNavigation( &fnormalNav
75
76 fStepEndPoint = G4ThreeVector( kInfinity, kI
77 fLastStepEndPointLocal = G4ThreeVector( kInf
78
79 fpVoxelSafety = new G4VoxelSafety();
80 fpvoxelNav = new G4VoxelNavigation();
81 fpSafetyCalculator = new G4SafetyCalculator(
82 fpSafetyCalculator->SetExternalNavigation(fp
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 G4T
103 const G4T
104 const G4T
105 {
106 ResetState();
107 fHistory = *h.GetHistory();
108 SetupHierarchy();
109 fLastTriedStepComputation = false; // Redun
110 return LocateGlobalPointAndSetup(p, &directi
111 }
112
113 // *******************************************
114 // LocateGlobalPointAndSetup
115 //
116 // Locate the point in the hierarchy return 0
117 // The direction is required
118 // - if on an edge shared by more than two
119 // (to resolve likely looping in tracking
120 // - at initial location of a particle
121 // (to resolve potential ambiguity at bou
122 //
123 // Flags on exit: (comments to be completed)
124 // fEntering - True if entering `daugh
125 // whether daughter of las
126 // or daughter of that vol
127 // fExiting - True if exited 'mother'
128 // (always ? - how about i
129 // *******************************************
130 //
131 G4VPhysicalVolume*
132 G4Navigator::LocateGlobalPointAndSetup( const
133 const
134 const
135 const
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
145
146 fLastTriedStepComputation = false;
147 fChangedGrandMotherRefFrame = false; // For
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::LocateGlobalPo
159 G4cout << " Called with arguments: " <<
160 << " Globalpoint = " << glob
161 << " RelativeSearch = " << r
162 if( fVerbose >= 4 )
163 {
164 G4cout << " ----- Upon entering:" <<
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; // Remem
183 fExitedMother = fExiting; // Remem
184 if ( fExiting )
185 {
186 ++noLevelsExited; // count this first
187
188 if ( fHistory.GetDepth() != 0 )
189 {
190 fBlockedPhysicalVolume = fHistory.Ge
191 fBlockedReplicaNo = fHistory.GetTopR
192 fHistory.BackLevel();
193 }
194 else
195 {
196 fLastLocatedPointLocal = localPoint;
197 fLocatedOutsideWorld = true;
198 fBlockedPhysicalVolume = nullptr;
199 fBlockedReplicaNo = -1;
200 fEntering = false; // No
201 fEnteredDaughter = false;
202 fExitedMother = true; // ??
203
204 return nullptr; // Have ex
205 }
206 // A fix for the case where a volume i
207 // and a coincident surface exists out
208 // - This stops it from exiting furth
209 // - However ReplicaNavigator treats
210 //
211 // assert( fBlockedPhysicalVolume!=0 )
212
213 // Expect to be on edge => on surface
214 //
215 if ( fLocatedOnEdge && (VolumeType(fBl
216 {
217 fExiting = false;
218 // Consider effect on Exit Normal !?
219 }
220 }
221 else
222 if ( fEntering )
223 {
224 switch (VolumeType(fBlockedPhysicalV
225 {
226 case kNormal:
227 fHistory.NewLevel(fBlockedPhysic
228 fBlockedPhysic
229 break;
230 case kReplica:
231 freplicaNav.ComputeTransformatio
232
233 fHistory.NewLevel(fBlockedPhysic
234 fBlockedReplic
235 fBlockedPhysicalVolume->SetCopyN
236 break;
237 case kParameterised:
238 if( fBlockedPhysicalVolume->GetR
239 {
240 G4VSolid *pSolid;
241 G4VPVParameterisation *pParam;
242 G4TouchableHistory parentTouch
243 pParam = fBlockedPhysicalVolum
244 pSolid = pParam->ComputeSolid(
245
246 pSolid->ComputeDimensions(pPar
247 fBlo
248 pParam->ComputeTransformation(
249
250 fHistory.NewLevel(fBlockedPhys
251 fBlockedRepl
252 fBlockedPhysicalVolume->SetCop
253 //
254 // Set the correct solid and m
255 //
256 G4LogicalVolume *pLogical;
257 pLogical = fBlockedPhysicalVol
258 pLogical->SetSolid( pSolid );
259 pLogical->UpdateMaterial(pPara
260 ComputeMaterial(fBlockedRepl
261 fBlockedPhys
262 &parentTouch
263 }
264 break;
265 case kExternal:
266 G4Exception("G4Navigator::Locate
267 "GeomNav0001", Fatal
268 "Extra levels not ap
269 break;
270 }
271 fEntering = false;
272 fBlockedPhysicalVolume = nullptr;
273 localPoint = fHistory.GetTopTransfor
274 notKnownContained = false;
275 }
276 }
277 else
278 {
279 fBlockedPhysicalVolume = nullptr;
280 fEntering = false;
281 fEnteredDaughter = false; // Full Step
282 fExiting = false;
283 fExitedMother = false; // Full Step
284 }
285 }
286 //
287 // Search from top of history up through geo
288 // containing volume found:
289 // If on
290 // o OUTSIDE - Back up level, not/no longer
291 // o SURFACE and EXITING - Back up level, se
292 // else
293 // o containing volume found
294 //
295
296 while (notKnownContained) // Loop checking,
297 {
298 EVolume topVolumeType = fHistory.GetTopVol
299 if (topVolumeType!=kReplica && topVolumeTy
300 {
301 targetSolid = fHistory.GetTopVolume()->G
302 localPoint = fHistory.GetTopTransform().
303 insideCode = targetSolid->Inside(localPo
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::LocateGlob
317 << " Invoked Inside() for so
318 << ". Solid replied: " << solid
319 << " For local point p: " <<
320 }
321 #endif
322 }
323 else
324 {
325 if( topVolumeType == kReplica )
326 {
327 insideCode = freplicaNav.BackLocate(
328
329 // !CARE! if notKnownContained retur
330 // the containing placement volume o
331 // will result in the history being
332 // local point returned is the point
333 }
334 else
335 {
336 targetSolid = fHistory.GetTopVolume(
337 localPoint = fHistory.GetTopTransfor
338 G4ThreeVector localDirection =
339 fHistory.GetTopTransform().Transf
340 insideCode = fpExternalNav->Inside(t
341 }
342 }
343
344 // Point is inside current volume, break o
345 if ( insideCode == kInside ) { break; }
346
347 // Point is outside current volume, move u
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.GetTop
361 fBlockedReplicaNo = fHistory.GetTopRepli
362 fHistory.BackLevel();
363 fExiting = false;
364
365 if( noLevelsExited > 1 )
366 {
367 // The first transformation was done b
368 //
369 if(const auto *mRot = fBlockedPhysical
370 {
371 fGrandMotherExitNormal *= (*mRot).in
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 thi
383 // by using the direction
384 //
385 G4bool directionExiting = false;
386 G4ThreeVector localDirection =
387 fHistory.GetTopTransform().TransformAx
388
389 // Make sure localPoint in correct refer
390 // ( Was it already correct ? How ?
391 //
392 localPoint= fHistory.GetTopTransform().T
393 if ( fHistory.GetTopVolumeType() != kRep
394 {
395 G4ThreeVector normal = targetSolid->Su
396 directionExiting = normal.dot(localDir
397 isExiting = isExiting || directionExit
398 }
399 }
400
401 // Point is on a surface, but no longer ex
402 if ( !isExiting ) { break; }
403
404 ++noLevelsExited;
405
406 // Point is on the outer surface, leaving
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
415 fValidExitNormal = false;
416 fBlockedPhysicalVolume = fHistory.GetTopVo
417 fBlockedReplicaNo = fHistory.GetTopReplica
418 fHistory.BackLevel();
419
420 if( noLevelsExited > 1 )
421 {
422 // The first transformation was done by
423 //
424 const G4RotationMatrix* mRot =
425 fBlockedPhysicalVolume->GetRotation();
426 if( mRot != nullptr )
427 {
428 fGrandMotherExitNormal *= (*mRot).inve
429 fChangedGrandMotherRefFrame = true;
430 }
431 }
432 } // END while (notKnownContained)
433 //
434 // Search downwards until deepest containing
435 // blocking fBlockedPhysicalVolume/BlockedRe
436 //
437 // 3 Cases:
438 //
439 // o Parameterised daughters
440 // =>Must be one G4PVParameterised daughte
441 // o Positioned daughters & voxels
442 // o Positioned daughters & no voxels
443
444 noResult = true; // noResult should be rena
445 // something like enteredL
446 do
447 {
448 // Determine `type' of current mother volu
449 //
450 targetPhysical = fHistory.GetTopVolume();
451 if (targetPhysical == nullptr) { break; }
452 targetLogical = targetPhysical->GetLogical
453 switch( CharacteriseDaughters(targetLogica
454 {
455 case kNormal:
456 if ( targetLogical->GetVoxelHeader() !
457 {
458 noResult = GetVoxelNavigator().Level
459 fBl
460 fBl
461 glo
462 pGl
463 con
464 loc
465 }
466 else // do not u
467 {
468 noResult = fnormalNav.LevelLocate(fH
469 fB
470 fB
471 gl
472 pG
473 co
474 lo
475 }
476 break;
477 case kReplica:
478 noResult = freplicaNav.LevelLocate(fHi
479 fBl
480 fBl
481 glo
482 pGl
483 con
484 loc
485 break;
486 case kParameterised:
487 if( GetDaughtersRegularStructureId(tar
488 {
489 noResult = fparamNav.LevelLocate(fHi
490 fBl
491 fBl
492 glo
493 pGl
494 con
495 loc
496 }
497 else // Regular structure
498 {
499 noResult = fregularNav.LevelLocate(f
500 f
501 f
502 g
503 p
504 c
505 l
506 }
507 break;
508 case kExternal:
509 noResult = fpExternalNav->LevelLocate(
510
511
512
513
514
515
516 break;
517 }
518
519 // LevelLocate returns true if it finds a
520 // in which globalPoint is inside (or on t
521
522 if ( noResult )
523 {
524 // Entering a daughter after ascending
525 //
526 // The blocked volume is no longer valid
527 //
528 fBlockedPhysicalVolume = nullptr;
529 fBlockedReplicaNo = -1;
530
531 // fEntering should be false -- else blo
532 // fEnteredDaughter is used for ExitNorm
533 //
534 fEntering = false;
535 fEnteredDaughter = true;
536
537 if( fExitedMother )
538 {
539 G4VPhysicalVolume* enteredPhysical = f
540 const G4RotationMatrix* mRot = entered
541 if( mRot != nullptr )
542 {
543 // Go deeper, i.e. move 'down' in th
544 // Apply direct rotation, not invers
545 //
546 fGrandMotherExitNormal *= (*mRot);
547 fChangedGrandMotherRefFrame= true;
548 }
549 }
550
551 #ifdef G4DEBUG_NAVIGATION
552 if( fVerbose > 2 )
553 {
554 G4VPhysicalVolume* enteredPhysical =
555 G4cout << "*** G4Navigator::LocateGlo
556 G4cout << " Entering volume: " <<
557 << G4endl;
558 }
559 #endif
560 }
561 } while (noResult); // Loop checking, 07.10
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) { curPhysV
571 G4cout << " Return value = new volume =
572 G4cout << " ----- Upon exiting:" << G4e
573 PrintState();
574 if( fVerbose >= 5 )
575 {
576 G4cout << "Upon exiting LocateGlobalPoin
577 G4cout << " History = " << 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
592 // is updated in order to start the next st
593 // -> no check is performed whether pGlobalpoi
594 // original volume (this must be the case).
595 //
596 // Note: a direction could be added to the arg
597 // optional checking (via the old code b
598 // [ This would be done only in verbose
599 // *******************************************
600 //
601 void
602 G4Navigator::LocateGlobalPointWithinVolume(con
603 {
604 #ifdef G4DEBUG_NAVIGATION
605 assert( !fWasLimitedByGeometry );
606 // Check: Either step was not limited by a
607 // else the full step is no longer
608 #endif
609
610 fLastLocatedPointLocal = ComputeLocalPoint(
611 fLastTriedStepComputation = false;
612 fChangedGrandMotherRefFrame = false; // F
613
614 // For the case of Voxel (or Parameterised)
615 // Navigator must be messaged to update its
616
617 // Update the state of the Sub Navigators
618 // - in particular any voxel information th
619 //
620 G4VPhysicalVolume* motherPhysical = fHisto
621 G4LogicalVolume* motherLogical = mother
622
623 switch( CharacteriseDaughters(motherLogical
624 {
625 case kNormal:
626 GetVoxelNavigator().RelocateWithinVol
627 break;
628 case kParameterised:
629 fparamNav.RelocateWithinVolume( mothe
630 break;
631 case kReplica:
632 // Nothing to do
633 break;
634 case kExternal:
635 fpExternalNav->RelocateWithinVolume(
636
637 break;
638 }
639
640 // Reset the state variables
641 // - which would have been affected
642 // by the 'equivalent' call to LocateGl
643 // - who's values have been invalidated b
644 //
645 fBlockedPhysicalVolume = nullptr;
646 fBlockedReplicaNo = -1;
647 fEntering = false;
648 fEnteredDaughter = false; // Boundary not
649 fExiting = false;
650 fExitedMother = false; // Boundary not
651 }
652
653 // *******************************************
654 // SetSavedState
655 //
656 // Save the state, in case this is a parasitic
657 // Save fValidExitNormal, fExitNormal, fExitin
658 // fBlockedPhysicalVolume, fBlockedReplic
659 // *******************************************
660 //
661 void G4Navigator::SetSavedState()
662 {
663 // Note: the state of dependent objects is n
664 // ( This means that the full state is cha
665 // SetSavedState() and RestoreSavedState
666
667 fSaveState.sExitNormal = fExitNormal;
668 fSaveState.sValidExitNormal = fValidExitNorm
669 fSaveState.sExiting = fExiting;
670 fSaveState.sEntering = fEntering;
671
672 fSaveState.spBlockedPhysicalVolume = fBlocke
673 fSaveState.sBlockedReplicaNo = fBlockedRepli
674
675 fSaveState.sLastStepWasZero = static_cast<G4
676
677 fSaveState.sLocatedOutsideWorld = fLocatedOu
678 fSaveState.sLastLocatedPointLocal = fLastLoc
679 fSaveState.sEnteredDaughter = fEnteredDaught
680 fSaveState.sExitedMother = fExitedMother;
681 fSaveState.sWasLimitedByGeometry = fWasLimit
682
683 // Even the safety sphere - if you want to c
684 //
685 fSaveState.sPreviousSftOrigin = fPreviousSft
686 fSaveState.sPreviousSafety = fPreviousSafety
687 }
688
689 // *******************************************
690 // RestoreSavedState
691 //
692 // Restore the state (in Compute Step), in cas
693 // *******************************************
694 //
695 void G4Navigator::RestoreSavedState()
696 {
697 fExitNormal = fSaveState.sExitNormal;
698 fValidExitNormal = fSaveState.sValidExitNorm
699 fExiting = fSaveState.sExiting;
700 fEntering = fSaveState.sEntering;
701
702 fBlockedPhysicalVolume = fSaveState.spBlocke
703 fBlockedReplicaNo = fSaveState.sBlockedRepli
704
705 fLastStepWasZero = (fSaveState.sLastStepWasZ
706
707 fLocatedOutsideWorld = fSaveState.sLocatedOu
708 fLastLocatedPointLocal = fSaveState.sLastLoc
709 fEnteredDaughter = fSaveState.sEnteredDaught
710 fExitedMother = fSaveState.sExitedMother;
711 fWasLimitedByGeometry = fSaveState.sWasLimit
712
713 // The 'expected' behaviour is to restore th
714 fPreviousSftOrigin = fSaveState.sPreviousSft
715 fPreviousSafety = fSaveState.sPreviousSafety
716 }
717
718 // *******************************************
719 // ComputeStep
720 //
721 // Computes the next geometric Step: intersect
722 // mother and `daughter' volumes.
723 //
724 // NOTE:
725 //
726 // Flags on entry:
727 // --------------
728 // fValidExitNormal - Normal of exited volume
729 // coincident boundary)
730 // fExitNormal - Surface normal of exite
731 // fExiting - True if have exited sol
732 //
733 // fBlockedPhysicalVolume - Ptr to exited volu
734 // fBlockedReplicaNo - Replication no of exite
735 // fLastStepWasZero - True if last Step size
736 //
737 // Flags on exit:
738 // -------------
739 // fValidExitNormal - True if surface normal
740 // fExitNormal - Surface normal of exite
741 // reference system
742 // fExiting - True if exiting mother
743 // fEntering - True if entering `daugh
744 // fBlockedPhysicalVolume - Ptr to candidate (
745 // fBlockedReplicaNo - Replication no of candi
746 // fLastStepWasZero - True if this Step size
747 // *******************************************
748 //
749 G4double G4Navigator::ComputeStep( const G4Thr
750 const G4Thr
751 const G4dou
752 G4dou
753 {
754 #ifdef G4DEBUG_NAVIGATION
755 static G4ThreadLocal G4int sNavCScalls = 0;
756 ++sNavCScalls;
757 #endif
758
759 G4ThreeVector localDirection = ComputeLocalA
760 G4double Step = kInfinity;
761 G4VPhysicalVolume *motherPhysical = fHistor
762 G4LogicalVolume *motherLogical = motherPhysi
763
764 // All state relating to exiting normals mus
765 //
766 fExitNormalGlobalFrame = G4ThreeVector( 0.,
767 // Reset value - to erase its memory
768 fChangedGrandMotherRefFrame = false;
769 // Reset - used for local exit normal
770 fGrandMotherExitNormal = G4ThreeVector( 0.,
771 fCalculatedExitNormal = false;
772 // Reset for new step
773
774 #ifdef G4VERBOSE
775 if( fVerbose > 0 )
776 {
777 G4cout << "*** G4Navigator::ComputeStep: *
778 G4cout << " Volume = " << motherPhysica
779 << " - Proposed step length = " <<
780 << G4endl;
781 #ifdef G4DEBUG_NAVIGATION
782 if( fVerbose >= 2 )
783 {
784 G4cout << " Called with the arguments:
785 << " Globalpoint = " << std::set
786 << " Direction = " << std::set
787 if( fVerbose >= 4 )
788 {
789 G4cout << " ---- Upon entering : Stat
790 PrintState();
791 }
792 }
793 #endif
794 }
795 #endif
796
797 G4ThreeVector newLocalPoint = ComputeLocalPo
798
799 if( newLocalPoint != fLastLocatedPointLocal
800 {
801 // Check whether the relocation is within
802 //
803 G4ThreeVector oldLocalPoint = fLastLocated
804 G4double moveLenSq = (newLocalPoint-oldLoc
805
806 if ( moveLenSq >= fSqTol )
807 {
808 #ifdef G4VERBOSE
809 ComputeStepLog(pGlobalpoint, moveLenSq);
810 #endif
811 // Relocate the point within the same vo
812 //
813 LocateGlobalPointWithinVolume( pGlobalpo
814 }
815 }
816 if ( fHistory.GetTopVolumeType()!=kReplica )
817 {
818 switch( CharacteriseDaughters(motherLogica
819 {
820 case kNormal:
821 if ( motherLogical->GetVoxelHeader() !
822 {
823 Step = GetVoxelNavigator().ComputeSt
824 localDi
825 pCurren
826 pNewSaf
827 fHistor
828 fValidE
829 fExitNo
830 fExitin
831 fEnteri
832 &fBlock
833 fBlocke
834
835 }
836 else
837 {
838 if( motherPhysical->GetRegularStruct
839 {
840 Step = fnormalNav.ComputeStep(fLas
841 loca
842 pCur
843 pNew
844 fHis
845 fVal
846 fExi
847 fExi
848 fEnt
849 &fBl
850 fBlo
851 }
852 else // Regular (non-voxelised) str
853 {
854 LocateGlobalPointAndSetup( pGlobal
855 //
856 // if physical process limits the
857 // one given by ComputeStepSkippin
858 // point will be wrongly calculate
859
860 // There is a problem: when msc li
861 // assigned wrongly to phantom in
862 // of the container volume). Then
863 // reset the history topvolume to
864 //
865 if(fHistory.GetTopVolume()->GetReg
866 {
867 G4Exception("G4Navigator::Comput
868 "GeomNav1001", JustW
869 "Point is relocated in voxels,
870 Step = fnormalNav.ComputeStep(fL
871 lo
872 pC
873 pN
874 fH
875 fV
876 fE
877 fE
878 fE
879 &f
880 fB
881 }
882 else
883 {
884 Step = fregularNav.
885 ComputeStepSkippingEqualMat
886
887
888
889
890
891
892
893
894
895
896
897 }
898 }
899 }
900 break;
901 case kParameterised:
902 if( GetDaughtersRegularStructureId(mot
903 {
904 Step = fparamNav.ComputeStep(fLastLo
905 localDi
906 pCurren
907 pNewSaf
908 fHistor
909 fValidE
910 fExitNo
911 fExitin
912 fEnteri
913 &fBlock
914 fBlocke
915 }
916 else // Regular structure
917 {
918 Step = fregularNav.ComputeStep(fLast
919 local
920 pCurr
921 pNewS
922 fHist
923 fVali
924 fExit
925 fExit
926 fEnte
927 &fBlo
928 fBloc
929 }
930 break;
931 case kReplica:
932 G4Exception("G4Navigator::ComputeStep(
933 FatalException, "Not appli
934 break;
935 case kExternal:
936 Step = fpExternalNav->ComputeStep(fLas
937 loca
938 pCur
939 pNew
940 fHis
941 fVal
942 fExi
943 fExi
944 fEnt
945 &fBl
946 fBlo
947 break;
948 }
949 }
950 else
951 {
952 // In the case of a replica, it must handl
953 // edge/corner problem by itself
954 //
955 fExiting = fExitedMother;
956 Step = freplicaNav.ComputeStep(pGlobalpoin
957 pDirection,
958 fLastLocate
959 localDirect
960 pCurrentPro
961 pNewSafety,
962 fHistory,
963 fValidExitN
964 fCalculated
965 fExitNormal
966 fExiting,
967 fEntering,
968 &fBlockedPh
969 fBlockedRep
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 c
978 // - one, two (or a few) ca
979 // volumes
980 // - more than two is likel
981 // description or the Nav
982
983 // Rule of thumb: likely at an Edge if two c
984 // because at least two candi
985 // checked
986 //
987 fLocatedOnEdge = fLastStepWasZero && (Step
988 fLastStepWasZero = (Step<fMinStep);
989 if (fPushed) { fPushed = fLastStepWasZero;
990
991 // Handle large number of consecutive zero s
992 //
993 if ( fLastStepWasZero )
994 {
995 ++fNumberZeroSteps;
996
997 G4bool act = fNumberZeroSteps >= fActi
998 G4bool actAndReport = false;
999 G4bool abandon = fNumberZeroSteps >= fAban
1000 G4bool inform = false;
1001 #ifdef G4VERBOSE
1002 actAndReport = act && (!fPushed) && fWarn
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. P
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 ge
1024 << G4endl;
1025 message << " Track stuck, not moving
1026 << fNumberZeroSteps << " step
1027 << " Current phys volume: '
1028 << "'" << G4endl
1029 << " - at position : " << p
1030 << " in direction: " << p
1031 << " (local position: " <<
1032 << " (local direction: " <
1033 << " Previous phys volume: '
1034 << ( fLastMotherPhys != nullp
1035 << "'" << G4endl << G4endl;
1036 if( actAndReport || abandon )
1037 {
1038 message << " Likely geometry over
1039 << G4endl;
1040 }
1041 if( abandon ) // i.e. fNumberZeroStep
1042 {
1043 // Must kill this stuck track
1044 #ifdef G4VERBOSE
1045 if ( fWarnPush ) { CheckOverlapsIte
1046 #endif
1047 message << " Track *abandoned* due
1048 << " Event aborted. " << G4
1049 G4Exception("G4Navigator::ComputeSt
1050 EventMustBeAborted, mes
1051 }
1052 else
1053 {
1054 #ifdef G4VERBOSE
1055 if ( actAndReport ) // (!fPushed =
1056 {
1057 message << " *** Trying to get
1058 << " - expanding step to
1059 << " Potential ove
1060 G4Exception("G4Navigator::Comput
1061 JustWarning, message
1062 }
1063 #endif
1064 #ifdef G4DEBUG_NAVIGATION
1065 else
1066 {
1067 if( fNumberZeroSteps > 1 )
1068 {
1069 message << ", nav-comp-step ca
1070 << ", Step= " << Step
1071 G4cout << message.str();
1072 }
1073 }
1074 #endif
1075 } // end of else if ( abandon )
1076 } // end of if( actAndReport || abandon
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
1086 fExitedMother = fExiting;
1087
1088 fStepEndPoint = pGlobalpoint
1089 + std::min(Step,pCurrentPropo
1090 fLastStepEndPointLocal = fLastLocatedPointL
1091
1092 if( fExiting )
1093 {
1094 #ifdef G4DEBUG_NAVIGATION
1095 if( fVerbose > 2 )
1096 {
1097 G4cout << " At G4Nav CompStep End - if(
1098 << " fValidExitNormal = " << fVa
1099 G4cout << " fExitNormal= " << fExitNorm
1100 }
1101 #endif
1102
1103 if ( fValidExitNormal || fCalculatedExitN
1104 {
1105 // Convention: fExitNormal is in the 'g
1106 fGrandMotherExitNormal = fExitNormal;
1107 }
1108 else
1109 {
1110 // We must calculate the normal anyway
1111 //
1112 G4ThreeVector finalLocalPoint = fLastLo
1113 + localDi
1114
1115 if ( fHistory.GetTopVolumeType() != kR
1116 {
1117 // Find normal in the 'mother' coordi
1118 //
1119 G4ThreeVector exitNormalMotherFrame=
1120 motherLogical->GetSolid()->Surface
1121
1122 // Transform it to the 'grand-mother'
1123 //
1124 const G4RotationMatrix* mRot = mother
1125 if( mRot != nullptr )
1126 {
1127 fChangedGrandMotherRefFrame = true;
1128 fGrandMotherExitNormal = (*mRot).in
1129 }
1130 else
1131 {
1132 fGrandMotherExitNormal = exitNormal
1133 }
1134
1135 // Do not set fValidExitNormal -- thi
1136 // that the solid is convex!
1137 }
1138 else
1139 {
1140 fCalculatedExitNormal = false;
1141 //
1142 // Nothing can be done at this stage
1143 // Replica Navigation must have calcu
1144 // already.
1145 // Cases: mother is not convex, and e
1146
1147 #ifdef G4DEBUG_NAVIGATION
1148 G4ExceptionDescription desc;
1149
1150 desc << "Problem in ComputeStep: Rep
1151 << " valid exit Normal. " << G4e
1152 desc << " Do not know how calculate i
1153 desc << " Location = " << finalLo
1154 desc << " Volume name = " << motherP
1155 << " copy/replica No = " << mot
1156 G4Exception("G4Navigator::ComputeStep
1157 JustWarning, desc, "Norma
1158 #endif
1159 }
1160 }
1161
1162 if ( fHistory.GetTopVolumeType() != kRepl
1163 {
1164 fCalculatedExitNormal = true;
1165 }
1166
1167 // Now transform it to the global referen
1168 //
1169 if( fValidExitNormal || fCalculatedExitNo
1170 {
1171 auto depth = (G4int)fHistory.GetDepth(
1172 if( depth > 0 )
1173 {
1174 fExitNormalGlobalFrame = fHistory.Get
1175 .InverseTrans
1176 }
1177 else
1178 {
1179 fExitNormalGlobalFrame = fGrandMother
1180 }
1181 }
1182 else
1183 {
1184 fExitNormalGlobalFrame = G4ThreeVector(
1185 }
1186 }
1187
1188 if( (Step == pCurrentProposedStepLength) &&
1189 {
1190 // This if Step is not really limited by
1191 // The Navigator is obliged to return "in
1192 //
1193 Step = kInfinity;
1194 }
1195
1196 #ifdef G4VERBOSE
1197 if( fVerbose > 1 )
1198 {
1199 if( fVerbose >= 4 )
1200 {
1201 G4cout << " ----- Upon exiting :" <<
1202 PrintState();
1203 }
1204 G4cout << " Returned step= " << Step;
1205 if( fVerbose > 5 ) { G4cout << G4endl; }
1206 if( Step == kInfinity )
1207 {
1208 G4cout << " Requested step= " << pCurr
1209 if( fVerbose > 5) { G4cout << G4endl;
1210 }
1211 G4cout << " Safety = " << pNewSafety <<
1212 }
1213 #endif
1214
1215 fLastTriedStepComputation = true;
1216
1217 return Step;
1218 }
1219
1220 // ******************************************
1221 // CheckNextStep
1222 //
1223 // Compute the step without altering the navi
1224 // ******************************************
1225 //
1226 G4double G4Navigator::CheckNextStep( const G4
1227 const G4
1228 const G4
1229 G4
1230 {
1231 G4double step;
1232
1233 // Save the state, for this parasitic call
1234 //
1235 SetSavedState();
1236
1237 step = ComputeStep ( pGlobalpoint,
1238 pDirection,
1239 pCurrentProposedStepLe
1240 pNewSafety );
1241
1242 // It is a parasitic call, so attempt to re
1243 //
1244 RestoreSavedState();
1245 // NOTE: the state of the current subnaviga
1246 // ***> TODO: restore subnavigator state
1247 // if( last_located) Need
1248 // if( last_computed step) Need
1249
1250 return step;
1251 }
1252
1253 // ******************************************
1254 // ResetState
1255 //
1256 // Resets stack and minimum of navigator stat
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,
1275 fGrandMotherExitNormal = G4ThreeVector(0,0,
1276 fExitNormalGlobalFrame = G4ThreeVector(0,0,
1277
1278 fPreviousSftOrigin = G4ThreeVector(0,0,
1279 fPreviousSafety = 0.0;
1280
1281 fNumberZeroSteps = 0;
1282
1283 fBlockedPhysicalVolume = nullptr;
1284 fBlockedReplicaNo = -1;
1285
1286 fLastLocatedPointLocal = G4ThreeVector( kIn
1287 fLocatedOutsideWorld = false;
1288
1289 fLastMotherPhys = nullptr;
1290 }
1291
1292 // ******************************************
1293 // SetupHierarchy
1294 //
1295 // Renavigates & resets hierarchy described b
1296 // o Reset volumes
1297 // o Recompute transforms and/or solids of re
1298 // ******************************************
1299 //
1300 void G4Navigator::SetupHierarchy()
1301 {
1302 const auto depth = (G4int)fHistory.GetDept
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(fHi
1312 break;
1313 case kParameterised:
1314 G4VPhysicalVolume* current = fHistory
1315 G4int replicaNo = fHistory.GetReplica
1316 G4VPVParameterisation* pParam = curre
1317 G4VSolid* pSolid = pParam->ComputeSol
1318
1319 // Set up dimensions & transform in s
1320 //
1321 pSolid->ComputeDimensions(pParam, rep
1322 pParam->ComputeTransformation(replica
1323
1324 G4TouchableHistory* pTouchable = null
1325 if( pParam->IsNested() )
1326 {
1327 pTouchable= new G4TouchableHistory(
1328 pTouchable->MoveUpHistory(); // Mov
1329 // Adequate only if Nested at the
1330 // To extend to other cases:
1331 // pTouchable->MoveUpHistory(cdep
1332 // Move to the parent level of *C
1333 // Could replace this line and co
1334 // c-tor for History(levels to dr
1335 }
1336 // Set up the correct solid and mater
1337 //
1338 G4LogicalVolume* pLogical = current->
1339 pLogical->SetSolid( pSolid );
1340 pLogical->UpdateMaterial( pParam ->
1341 ComputeMaterial(replicaNo, current,
1342 delete pTouchable;
1343 break;
1344 }
1345 }
1346 }
1347
1348 // ******************************************
1349 // GetLocalExitNormal
1350 //
1351 // Obtains the Normal vector to a surface (in
1352 // pointing out of previous volume and into c
1353 // ******************************************
1354 //
1355 G4ThreeVector G4Navigator::GetLocalExitNormal
1356 {
1357 G4ThreeVector ExitNormal(0.,0.,0.);
1358 G4VSolid* currentSolid = nullptr;
1359 G4LogicalVolume* candidateLogical;
1360
1361 if ( fLastTriedStepComputation )
1362 {
1363 // use fLastLocatedPointLocal and next ca
1364 //
1365 G4ThreeVector nextSolidExitNormal(0.,0.,0
1366
1367 if( fEntering && (fBlockedPhysicalVolume!
1368 {
1369 candidateLogical = fBlockedPhysicalVolu
1370 if( candidateLogical != nullptr )
1371 {
1372 // fLastStepEndPointLocal is in the c
1373 // we need it in the daughter's coord
1374
1375 // The following code should also wor
1376 {
1377 // First transform fLastLocatedPoin
1378 // coordinates
1379 //
1380 G4AffineTransform MotherToDaughterT
1381 GetMotherToDaughterTransform( fBl
1382 fBl
1383 Vol
1384 G4ThreeVector daughterPointOwnLocal
1385 MotherToDaughterTransform.Transfo
1386
1387 // OK if it is a parameterised volu
1388 //
1389 EInside inSideIt;
1390 G4bool onSurface;
1391 G4double safety = -1.0;
1392 currentSolid = candidateLogical->Ge
1393 inSideIt = currentSolid->Inside(dau
1394 onSurface = (inSideIt == kSurface);
1395 if( !onSurface )
1396 {
1397 if( inSideIt == kOutside )
1398 {
1399 safety = (currentSolid->Distanc
1400 onSurface = safety < 100.0 * kC
1401 }
1402 else if (inSideIt == kInside )
1403 {
1404 safety = (currentSolid->Distanc
1405 onSurface = safety < 100.0 * kC
1406 }
1407 }
1408
1409 if( onSurface )
1410 {
1411 nextSolidExitNormal =
1412 currentSolid->SurfaceNormal(dau
1413
1414 // Entering the solid ==> opposit
1415 //
1416 // First flip ( ExitNormal = -nex
1417 // and then rotate the the norma
1418 ExitNormal = MotherToDaughterTran
1419 .InverseTransformAxis
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 surfac
1429 << " Point =
1430 << daughterPointOwnLoca
1431 << " Physical volume =
1432 << fBlockedPhysicalVolu
1433 << " Logical volume =
1434 << candidateLogical->Ge
1435 << " Solid =
1436 << " Type =
1437 << currentSolid->GetEnt
1438 << *currentSolid << G4e
1439 if( inSideIt == kOutside )
1440 {
1441 message << "Point is Outside.
1442 << " Safety (from ou
1443 }
1444 else // if( inSideIt == kInside
1445 {
1446 message << "Point is Inside.
1447 << " Safety (from in
1448 }
1449 G4Exception("G4Navigator::GetLo
1450 JustWarning, messag
1451 }
1452 #endif
1453 }
1454 *valid = onSurface; // was =tru
1455 }
1456 }
1457 }
1458 else if ( fExiting )
1459 {
1460 ExitNormal = fGrandMotherExitNormal;
1461 *valid = true;
1462 fCalculatedExitNormal = true; // Shoul
1463 }
1464 else // i.e. ( fBlockedPhysicalVolume =
1465 {
1466 *valid = false;
1467 G4Exception("G4Navigator::GetLocalExitN
1468 "GeomNav0003", JustWarning,
1469 "Incorrect call to GetLocal
1470 }
1471 }
1472 else // ( ! fLastTriedStepComputation ) i.
1473 {
1474 if ( EnteredDaughterVolume() )
1475 {
1476 G4VSolid* daughterSolid = fHistory.GetT
1477
1478 ExitNormal = -(daughterSolid->SurfaceNo
1479 if( std::fabs(ExitNormal.mag2()-1.0 ) >
1480 {
1481 G4ExceptionDescription desc;
1482 desc << " Parameters of solid: " << *
1483 << " Point for surface = " << fL
1484 G4Exception("G4Navigator::GetLocalExi
1485 "GeomNav0003", FatalExcep
1486 "Surface Normal returned
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. Exit
1500 {
1501 *valid = false;
1502 fCalculatedExitNormal = false;
1503 G4ExceptionDescription message;
1504 message << "Function called when *NOT
1505 message << "Exit Normal not calculate
1506 G4Exception("G4Navigator::GetLocalExi
1507 "GeomNav0003", JustWarnin
1508 }
1509 }
1510 }
1511 return ExitNormal;
1512 }
1513
1514 // ******************************************
1515 // GetMotherToDaughterTransform
1516 //
1517 // Obtains the mother to daughter affine tran
1518 // ******************************************
1519 //
1520 G4AffineTransform
1521 G4Navigator::GetMotherToDaughterTransform( G4
1522 G4
1523 EV
1524 {
1525 switch (enteringVolumeType)
1526 {
1527 case kNormal: // Nothing is needed to pr
1528 break; // It is stored already in
1529 case kReplica: // Sets the transform in t
1530 G4Exception("G4Navigator::GetMotherToDa
1531 "GeomNav0001", FatalExcepti
1532 "Method NOT Implemented yet
1533 break;
1534 case kParameterised:
1535 if( pEnteringPhysVol->GetRegularStructu
1536 {
1537 G4VPVParameterisation *pParam =
1538 pEnteringPhysVol->GetParameterisati
1539 G4VSolid* pSolid =
1540 pParam->ComputeSolid(enteringReplic
1541 pSolid->ComputeDimensions(pParam, ent
1542
1543 // Sets the transform in the Paramete
1544 //
1545 pParam->ComputeTransformation(enterin
1546
1547 // Set the correct solid and material
1548 //
1549 G4LogicalVolume* pLogical = pEntering
1550 pLogical->SetSolid( pSolid );
1551 }
1552 break;
1553 case kExternal:
1554 // Expect that nothing is needed to pre
1555 // It is stored already in the physical
1556 break;
1557 }
1558 return G4AffineTransform(pEnteringPhysVol->
1559 pEnteringPhysVol->
1560 }
1561
1562
1563 // ******************************************
1564 // GetLocalExitNormalAndCheck
1565 //
1566 // Obtains the Normal vector to a surface (in
1567 // pointing out of previous volume and into c
1568 // checks the current point against expected
1569 // ******************************************
1570 //
1571 G4ThreeVector
1572 G4Navigator::GetLocalExitNormalAndCheck(
1573 #ifdef G4DEBUG_NAVIGATION
1574 const G4ThreeVecto
1575 #else
1576 const G4ThreeVecto
1577 #endif
1578 G4bool* pVal
1579 {
1580 #ifdef G4DEBUG_NAVIGATION
1581 // Check Current point against expected 'lo
1582 //
1583 if ( fLastTriedStepComputation )
1584 {
1585 G4ThreeVector ExpectedBoundaryPointLocal;
1586
1587 const G4AffineTransform& GlobalToLocal =
1588 ExpectedBoundaryPointLocal =
1589 GlobalToLocal.TransformPoint( ExpectedB
1590
1591 // Add here: Comparison against expected
1592 // i.e. the endpoint of Comput
1593 }
1594 #endif
1595
1596 return GetLocalExitNormal( pValid );
1597 }
1598
1599 // ******************************************
1600 // GetGlobalExitNormal
1601 //
1602 // Obtains the Normal vector to a surface (in
1603 // pointing out of previous volume and into c
1604 // ******************************************
1605 //
1606 G4ThreeVector
1607 G4Navigator::GetGlobalExitNormal(const G4Thre
1608 G4bool
1609 {
1610 G4bool validNormal;
1611 G4ThreeVector localNormal, globalNormal;
1612
1613 G4bool usingStored = fCalculatedExitNormal
1614 ( fLastTriedStepComputation && fExitin
1615 ||
1616 ( !fLastTriedStepComputation
1617 && (IntersectPointGlobal-fStepEndPo
1618 // Calculated it 'just' before & t
1619 // but it did not move position
1620
1621 if( usingStored )
1622 {
1623 // This was computed in last call to Comp
1624 // and only if it arrived at boundary
1625 //
1626 globalNormal = fExitNormalGlobalFrame;
1627 G4double normMag2 = globalNormal.mag2();
1628 if( std::fabs ( normMag2 - 1.0 ) < perTho
1629 {
1630 *pNormalCalculated = true; // ComputeS
1631 // (fExitin
1632 }
1633 else
1634 {
1635 G4ExceptionDescription message;
1636 message.precision(10);
1637 message << " WARNING> Expected normal-
1638 << " i.e. a unit vector!" << G
1639 << " - but |normal| = " <<
1640 << " - and |normal|^2 = " <<
1641 << " which differs from 1.0 by
1642 << " n = " << fExitNormalGlo
1643 << " Global point: " << Inters
1644 << " Volume: " << fHistory.Get
1645 #ifdef G4VERBOSE
1646 G4LogicalVolume* candLog = fHistory.Ge
1647 if ( candLog != nullptr )
1648 {
1649 message << " Solid: " << candLog->Ge
1650 << ", Type: " << candLog->Ge
1651 << *candLog->GetSolid() << G
1652 }
1653 #endif
1654 message << "==========================
1655 << G4endl;
1656 G4int oldVerbose = fVerbose;
1657 fVerbose = 4;
1658 message << " State of Navigator: " <
1659 message << *this << G4endl;
1660 fVerbose = oldVerbose;
1661 message << "==========================
1662 << G4endl;
1663
1664 G4Exception("G4Navigator::GetGlobalExi
1665 "GeomNav0003",JustWarning,
1666 "Value obtained from stored glo
1667
1668 // (Re)Compute it now -- as either it
1669 //
1670 localNormal = GetLocalExitNormalAndChe
1671
1672 *pNormalCalculated = fCalculatedExitNo
1673 globalNormal = fHistory.GetTopTransfor
1674 .InverseTransformAxis(lo
1675 }
1676 }
1677 else
1678 {
1679 localNormal = GetLocalExitNormalAndCheck(
1680 *pNormalCalculated = fCalculatedExitNorma
1681
1682 #ifdef G4DEBUG_NAVIGATION
1683 usingStored = false;
1684
1685 if( (!validNormal) && !fCalculatedExitNor
1686 {
1687 G4ExceptionDescription edN;
1688 edN << " Calculated = " << fCalculated
1689 edN << " Entering= " << fEntering <<
1690 G4int oldVerbose = this->GetVerboseLeve
1691 this->SetVerboseLevel(4);
1692 edN << " State of Navigator: " << G4e
1693 edN << *this << G4endl;
1694 this->SetVerboseLevel( oldVerbose );
1695
1696 G4Exception("G4Navigator::GetGlobalExit
1697 "GeomNav0003", JustWarning,
1698 "LocalExitNormalAndCheck()
1699 }
1700 #endif
1701
1702 G4double localMag2 = localNormal.mag2();
1703 if( validNormal && (std::fabs(localMag2-
1704 {
1705 G4ExceptionDescription edN;
1706 edN.precision(10);
1707 edN << "G4Navigator::GetGlobalExitNorm
1708 << " Using Local Normal - from ca
1709 << G4endl
1710 << " Local Exit Normal : " << "
1711 << " vec = " << localNormal << G4e
1712 << " Global Exit Normal : " << "
1713 << " vec = " << globalNormal << G4
1714 << " Global point: " << Intersect
1715 edN << " Calculated It = " << fC
1716 << " Volume: " << fHistory.GetTop
1717 #ifdef G4VERBOSE
1718 G4LogicalVolume* candLog = fHistory.Ge
1719 if ( candLog != nullptr )
1720 {
1721 edN << " Solid: " << candLog->GetSo
1722 << ", Type: " << candLog->GetSol
1723 << *candLog->GetSolid();
1724 }
1725 #endif
1726 G4Exception("G4Navigator::GetGlobalExi
1727 "GeomNav0003",JustWarning,
1728 "Value obtained from new l
1729 localNormal = localNormal.unit(); // S
1730 }
1731 globalNormal = fHistory.GetTopTransform(
1732 .InverseTransformAxis(loca
1733 }
1734
1735 #ifdef G4DEBUG_NAVIGATION
1736 if( usingStored )
1737 {
1738 G4ThreeVector globalNormAgn;
1739
1740 localNormal = GetLocalExitNormalAndCheck(
1741
1742 globalNormAgn = fHistory.GetTopTransform(
1743 .InverseTransformAxis(loca
1744
1745 // Check the value computed against fExit
1746 G4ThreeVector diffNorm = globalNormAgn -
1747 if( diffNorm.mag2() > kToleranceNormalChe
1748 {
1749 G4ExceptionDescription edDfn;
1750 edDfn << "Found difference in normals i
1751 << "- when Get is called after Co
1752 edDfn << " Magnitude of diff = "
1753 edDfn << " Normal stored (Global)
1754 << G4endl;
1755 edDfn << " Global Computed from Local
1756 G4Exception("G4Navigator::GetGlobalExit
1757 JustWarning, edDfn);
1758 }
1759 }
1760 #endif
1761
1762 // Synchronise stored global exit normal as
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
1774 // ComputeStep.
1775 // ii) after (or at the end of) ComputeStep O
1776 // ******************************************
1777 //
1778 G4double G4Navigator::ComputeSafety( const G4
1779 const G4
1780 const G4
1781 {
1782 G4VPhysicalVolume *motherPhysical = fHisto
1783 G4double safety = 0.0;
1784
1785 G4double distEndpointSq = (pGlobalpoint-fSt
1786 G4bool stayedOnEndpoint = distEndpointSq <
1787 G4bool endpointOnSurface = fEnteredDaughter
1788
1789 G4bool onSurface = endpointOnSurface && sta
1790 if( ! onSurface )
1791 {
1792 safety= fpSafetyCalculator->SafetyInCurre
1793 // offload to G4SafetyCalculator - avoids
1794
1795 // Remember last safety origin & value
1796 //
1797 fPreviousSftOrigin = pGlobalpoint;
1798 fPreviousSafety = safety;
1799 // We overwrite the Safety 'sphere' - kee
1800 }
1801
1802 return safety;
1803 }
1804
1805 // ******************************************
1806 // CreateTouchableHistoryHandle
1807 // ******************************************
1808 //
1809 G4TouchableHandle G4Navigator::CreateTouchabl
1810 {
1811 return G4TouchableHandle( CreateTouchableHi
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 G4Navigat
1824 G4cout << " ValidExitNormal= " << fValid
1825 << " ExitNormal = " << fExitN
1826 << " Exiting = " << fExiti
1827 << " Entering = " << fEnter
1828 << " BlockedPhysicalVolume= " ;
1829 if (fBlockedPhysicalVolume==nullptr)
1830 {
1831 G4cout << "None";
1832 }
1833 else
1834 {
1835 G4cout << fBlockedPhysicalVolume->GetNa
1836 }
1837 G4cout << G4endl
1838 << " BlockedReplicaNo = " <<
1839 << " LastStepWasZero = " <<
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:Volu
1850 << std::setw( 9) << " ReplicaNo"
1851 << std::setw( 8) << " LastStepZero
1852 << G4endl;
1853 G4cout << "( " << std::setw(7) << fExitNo
1854 << ", " << std::setw(7) << fExitNo
1855 << ", " << std::setw(7) << fExitNo
1856 << std::setw( 5) << fValidExitNor
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)<< fBlockedPhysi
1863 G4cout << std::setw( 9) << fBlockedRepli
1864 << std::setw( 8) << fLastStepWasZ
1865 << G4endl;
1866 }
1867 if( fVerbose > 2 )
1868 {
1869 G4cout.precision(8);
1870 G4cout << " Current Localpoint = " << fLa
1871 G4cout << " PreviousSftOrigin = " << fPr
1872 G4cout << " PreviousSafety = " << fPr
1873 }
1874 G4cout.precision(oldcoutPrec);
1875 }
1876
1877 // ******************************************
1878 // ComputeStepLog
1879 // ******************************************
1880 //
1881 void G4Navigator::ComputeStepLog(const G4Thre
1882 G4doub
1883 {
1884 // The following checks only make sense if
1885 // than the tolerance.
1886
1887 const G4double fAccuracyForWarning = kCar
1888 fAccuracyForException = 1000
1889
1890 G4ThreeVector OriginalGlobalpoint = fHistor
1891 InverseTransfo
1892
1893 G4double shiftOriginSafSq = (fPreviousSftOr
1894
1895 // Check that the starting point of this st
1896 // within the isotropic safety sphere of th
1897 // to a accuracy/precision given by fAccur
1898 // If so give warning.
1899 // If it fails by more than fAccuracyForE
1900 //
1901 if( shiftOriginSafSq >= sqr(fPreviousSafety
1902 {
1903 G4double shiftOrigin = std::sqrt(shiftOri
1904 G4double diffShiftSaf = shiftOrigin - fPr
1905
1906 if( diffShiftSaf > fAccuracyForWarning )
1907 {
1908 G4long oldcoutPrec = G4cout.precision(8
1909 G4long oldcerrPrec = G4cerr.precision(1
1910 std::ostringstream message, suggestion;
1911 message << "Accuracy error or slightly
1912 << G4endl
1913 << " The Step's starting po
1914 << std::sqrt(moveLenSq)/mm << "
1915 << " since the last call to
1916 << " This has resulted in m
1917 << shiftOrigin/mm << " mm "
1918 << " from the last point at whi
1919 << " was calculated " << G4
1920 << " which is more than the
1921 << fPreviousSafety/mm << " mm
1922 << " This difference is "
1923 << diffShiftSaf/mm << " mm." <<
1924 << " The tolerated accuracy
1925 << fAccuracyForException/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
1933 << " - a process that has p
1934 << " larger than the current s
1935 << " - inaccuracy in the co
1936 suggestion << "We suggest that you "
1937 << " - find i) what part
1938 << " ii) through what part
1939 << " for example by r
1940 << G4endl
1941 << " /tracking/ver
1942 << " - check which proc
1943 << " this particle (and lo
1944 << G4endl
1945 << " - in case, create a
1946 << " of this event using:"
1947 << " /tracking/ver
1948 }
1949 G4Exception("G4Navigator::ComputeStep()
1950 "GeomNav1002", JustWarning,
1951 message, G4String(suggestio
1952 G4cout.precision(oldcoutPrec);
1953 G4cerr.precision(oldcerrPrec);
1954 }
1955 #ifdef G4DEBUG_NAVIGATION
1956 else
1957 {
1958 G4cerr << "WARNING - G4Navigator::Compu
1959 << " The Step's startin
1960 << std::sqrt(moveLenSq) << "," <
1961 << " which has taken it
1962 << " the current safety. " << G4
1963 }
1964 #endif
1965 }
1966 G4double safetyPlus = fPreviousSafety + fAc
1967 if ( shiftOriginSafSq > sqr(safetyPlus) )
1968 {
1969 std::ostringstream message;
1970 message << "May lead to a crash or unreli
1971 << " Position has shifted
1972 << " notifying the navigator !" <
1973 << " Tolerated safety: " <
1974 << " Computed shift : " <
1975 G4Exception("G4Navigator::ComputeStep()",
1976 JustWarning, message);
1977 }
1978 }
1979
1980 // ******************************************
1981 // CheckOverlapsIterative
1982 // ******************************************
1983 //
1984 G4bool G4Navigator::CheckOverlapsIterative(G4
1985 {
1986 // Check and report overlaps
1987 //
1988 G4bool foundOverlap = false;
1989 G4int nPoints = 300000, ntrials = 9, numO
1990 G4double trialLength = 1.0 * CLHEP::centim
1991 while ( ntrials-- > 0 && !foundOverlap )
1992 {
1993 if ( fVerbose > 1 )
1994 {
1995 G4cout << " ** Running overlap checks
1996 << vol->GetName()
1997 << " with length = " << trialLe
1998 }
1999 foundOverlap = vol->CheckOverlaps(nPoints
2000 fVerbos
2001 trialLength *= 0.1;
2002 if ( trialLength <= 1.0e-5 ) { numOverlap
2003 }
2004 return foundOverlap;
2005 }
2006
2007 // ******************************************
2008 // Operator <<
2009 // ******************************************
2010 //
2011 std::ostream& operator << (std::ostream &os,c
2012 {
2013 // Old version did only the following:
2014 // os << "Current History: " << G4endl << n
2015 // Old behaviour is recovered for fVerbose
2016
2017 // Adapted from G4Navigator::PrintState() c
2018
2019 G4long oldcoutPrec = os.precision(4);
2020 if( n.fVerbose >= 4 )
2021 {
2022 os << "The current state of G4Navigator i
2023 os << " ValidExitNormal= " << n.fValidEx
2024 << " ExitNormal = " << n.fExitNormal
2025 << " Exiting = " << n.fExiting
2026 << " Entering = " << n.fEntering
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.fBlo
2038 << " LastStepWasZero = " << n.fLa
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.fExitNorm
2053 << ", " << std::setw(7) << n.fExitNormal.
2054 << ", " << std::setw(7) << n.fExitNormal.
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.fBlockedPhysi
2062 os << std::setw( 9) << n.fBlockedReplica
2063 << std::setw( 8) << n.fLastStepWasZero
2064 << G4endl;
2065 }
2066 if( n.fVerbose > 2 )
2067 {
2068 os.precision(8);
2069 os << " Current Localpoint = " << n.fLast
2070 os << " PreviousSftOrigin = " << n.fPrev
2071 os << " PreviousSafety = " << n.fPrev
2072 }
2073 if( n.fVerbose > 3 || n.fVerbose == 0 )
2074 {
2075 os << "Current History: " << G4endl << n.
2076 }
2077
2078 os.precision(oldcoutPrec);
2079 return os;
2080 }
2081
2082 // ******************************************
2083 // SetVoxelNavigation: alternative navigator
2084 // ******************************************
2085 //
2086 void G4Navigator::SetVoxelNavigation(G4VoxelN
2087 {
2088 delete fpvoxelNav;
2089 fpvoxelNav = voxelNav;
2090 }
2091
2092 // ******************************************
2093 // InformLastStep: derived navigators can inf
2094 // used to update fLastStepWa
2095 // ******************************************
2096 void G4Navigator::InformLastStep(G4double la
2097 G4bool exit
2098 {
2099 G4bool zeroStep = ( lastStep == 0.0 );
2100 fLocatedOnEdge = fLastStepWasZero && zero
2101 fLastStepWasZero = zeroStep;
2102
2103 fExiting = exitsMotherVol;
2104 fEntering = entersDaughtVol;
2105 }
2106
2107 // ******************************************
2108 // SetExternalNavigation
2109 // ******************************************
2110 //
2111 void G4Navigator::SetExternalNavigation(G4VEx
2112 {
2113 fpExternalNav = externalNav;
2114 fpSafetyCalculator->SetExternalNavigation(e
2115 }
2116