| 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 //
27 // class G4ITNavigator2 Implementation
28 //
29 // Original author: Paul Kent, July 95/96
30 //
31 // G4ITNavigator2 is a duplicate version of G4
32 // initially written by Paul Kent and colleagu
33 // The only difference resides in the way the
34 //
35 // History:
36 // - Created.
37 // - Zero step protections
38 // - Added check mode
39 // - Made Navigator Abstract
40 // - G4ITNavigator2 created
41 // -------------------------------------------
42
43 #include <iomanip>
44
45 //#include "G4ITNavigator2.hh"
46 #include "G4ITNavigator.hh"
47 #include "G4ios.hh"
48 #include "G4SystemOfUnits.hh"
49 #include "G4GeometryTolerance.hh"
50 #include "G4VPhysicalVolume.hh"
51
52 //#define G4DEBUG_NAVIGATION 1
53 #include "G4VoxelSafety.hh"
54 #include "Randomize.hh"
55 #include "G4VoxelLimits.hh"
56
57 #define fHistory fpNavigatorState->fHistory
58 #define fLastTriedStepComputation fpNavigatorS
59 #define fEnteredDaughter fpNavigatorState->fEn
60 #define fExitedMother fpNavigatorState->fExite
61 #define fWasLimitedByGeometry fpNavigatorState
62 #define fStepEndPoint fpNavigatorState->fStepE
63 #define fLastStepEndPointLocal fpNavigatorStat
64 #define fPushed fpNavigatorState->fPushed
65 #define fLastTriedStepComputation fpNavigatorS
66 #define fEntering fpNavigatorState->fEntering
67 #define fExiting fpNavigatorState->fExiting
68 #define fBlockedPhysicalVolume fpNavigatorStat
69 #define fBlockedReplicaNo fpNavigatorState->fB
70 #define fLastLocatedPointLocal fpNavigatorStat
71 #define fLocatedOutsideWorld fpNavigatorState-
72 #define fValidExitNormal fpNavigatorState->fVa
73 #define fExitNormal fpNavigatorState->fExitNor
74 #define fGrandMotherExitNormal fpNavigatorStat
75 #define fChangedGrandMotherRefFrame fpNavigato
76 #define fExitNormalGlobalFrame fpNavigatorStat
77 #define fCalculatedExitNormal fpNavigatorState
78 #define fLastStepWasZero fpNavigatorState->fLa
79 #define fLocatedOnEdge fpNavigatorState->fLoca
80 #define fNumberZeroSteps fpNavigatorState->fNu
81 #define fPreviousSftOrigin fpNavigatorState->f
82 #define fPreviousSafety fpNavigatorState->fPre
83
84 //--------------------------------------------
85
86 namespace BoundingBox
87 {
88 enum Boundary
89 {
90 kMin,
91 kMax
92 };
93 }
94
95 // *******************************************
96 // Constructor
97 // *******************************************
98 //
99 G4ITNavigator2::G4ITNavigator2()
100 {
101 fActive= false;
102 fActionThreshold_NoZeroSteps = 1000;
103 fAbandonThreshold_NoZeroSteps = 2500;
104 kCarTolerance = G4GeometryTolerance::GetInst
105 fregularNav.SetNormalNavigation( &fnormalNav
106
107 fpNavigatorState = nullptr;
108 // fSaveState = 0;
109 fpVoxelSafety= new G4VoxelSafety();
110
111 // this->SetVerboseLevel(3);
112 // this->CheckMode(true);
113 }
114
115 // *******************************************
116 // Destructor
117 // *******************************************
118 //
119 G4ITNavigator2::~G4ITNavigator2()
120 { delete fpVoxelSafety; }
121
122 // *******************************************
123 // ResetHierarchyAndLocate
124 // *******************************************
125 //
126 G4VPhysicalVolume*
127 G4ITNavigator2::ResetHierarchyAndLocate(const
128 const G4T
129 const G4T
130 {
131 // ResetState();
132 fHistory = *h.GetHistory();
133 SetupHierarchy();
134 fLastTriedStepComputation= false; // Redund
135 return LocateGlobalPointAndSetup(p, &directi
136 }
137
138 // *******************************************
139 // LocateGlobalPointAndSetup
140 //
141 // Locate the point in the hierarchy return 0
142 // The direction is required
143 // - if on an edge shared by more than two
144 // (to resolve likely looping in tracking
145 // - at initial location of a particle
146 // (to resolve potential ambiguity at bou
147 //
148 // Flags on exit: (comments to be completed)
149 // fEntering - True if entering `daugh
150 // whether daughter of las
151 // or daughter of that vol
152 // *******************************************
153 //
154 G4VPhysicalVolume*
155 G4ITNavigator2::LocateGlobalPointAndSetup( con
156 const
157 const
158 const
159 {
160 CheckNavigatorStateIsValid();
161 G4bool notKnownContained=true, noResult;
162 G4VPhysicalVolume *targetPhysical;
163 G4LogicalVolume *targetLogical;
164 G4VSolid *targetSolid=nullptr;
165 G4ThreeVector localPoint, globalDirection;
166 EInside insideCode;
167
168 G4bool considerDirection = (!ignoreDirection
169
170 fLastTriedStepComputation= false;
171 fChangedGrandMotherRefFrame= false; // For
172
173 if( considerDirection && pGlobalDirection !=
174 {
175 globalDirection=*pGlobalDirection;
176 }
177
178
179 #ifdef G4VERBOSE
180 if( fVerbose > 2 )
181 {
182 G4long oldcoutPrec = G4cout.precision(8);
183 G4cout << "*** G4ITNavigator2::LocateGloba
184 G4cout << " Called with arguments: " <<
185 << " Globalpoint = " << globalPo
186 << " RelativeSearch = " << relat
187 if( fVerbose == 4 )
188 {
189 G4cout << " ----- Upon entering:" <<
190 PrintState();
191 }
192 G4cout.precision(oldcoutPrec);
193 }
194 #endif
195
196 G4int noLevelsExited=0 ;
197
198 if ( !relativeSearch )
199 {
200 ResetStackAndState();
201 }
202 else
203 {
204 if ( fWasLimitedByGeometry )
205 {
206 fWasLimitedByGeometry = false;
207 fEnteredDaughter = fEntering; // Remem
208 fExitedMother = fExiting; // Remem
209 if ( fExiting )
210 {
211 noLevelsExited++; // count this first
212
213 if ( fHistory.GetDepth() )
214 {
215 fBlockedPhysicalVolume = fHistory.Ge
216 fBlockedReplicaNo = fHistory.GetTopR
217 fHistory.BackLevel();
218 }
219 else
220 {
221 fLastLocatedPointLocal = localPoint;
222 fLocatedOutsideWorld = true;
223 fBlockedPhysicalVolume = nullptr;
224 fBlockedReplicaNo = -1;
225 fEntering = false; // No
226 fEnteredDaughter = false;
227 fExitedMother = true; // ??
228
229 return nullptr; // Have ex
230 }
231 // A fix for the case where a volume i
232 // and a coincident surface exists out
233 // - This stops it from exiting furth
234 // - However ReplicaNavigator treats
235 //
236 // assert( fBlockedPhysicalVolume!=0 )
237
238 // Expect to be on edge => on surface
239 //
240 if ( fLocatedOnEdge && (VolumeType(fBl
241 {
242 fExiting= false;
243 // Consider effect on Exit Normal !?
244 }
245 }
246 else
247 if ( fEntering )
248 {
249 // assert( fBlockedPhysicalVolume!=0
250
251 switch (VolumeType(fBlockedPhysicalV
252 {
253 case kNormal:
254 fHistory.NewLevel(fBlockedPhysic
255 fBlockedPhysic
256 break;
257 case kReplica:
258 freplicaNav.ComputeTransformatio
259
260 fHistory.NewLevel(fBlockedPhysic
261 fBlockedReplic
262 fBlockedPhysicalVolume->SetCopyN
263 break;
264 case kParameterised:
265 if( fBlockedPhysicalVolume->GetR
266 {
267 G4VSolid *pSolid;
268 G4VPVParameterisation *pParam;
269 G4TouchableHistory parentTouch
270 pParam = fBlockedPhysicalVolum
271 pSolid = pParam->ComputeSolid(
272
273 pSolid->ComputeDimensions(pPar
274 fBlo
275 pParam->ComputeTransformation(
276
277 fHistory.NewLevel(fBlockedPhys
278 fBlockedRepl
279 fBlockedPhysicalVolume->SetCop
280 //
281 // Set the correct solid and m
282 //
283 G4LogicalVolume *pLogical;
284 pLogical = fBlockedPhysicalVol
285 pLogical->SetSolid( pSolid );
286 pLogical->UpdateMaterial(pPara
287 ComputeMaterial(fBlockedRepl
288 fBlockedPhys
289 &parentTouch
290 }
291 break;
292 case kExternal:
293 G4Exception("G4ITNavigator2::Loca
294 "GeomNav0001", FatalE
295 "Not applicable for e
296 break;
297 }
298 fEntering = false;
299 fBlockedPhysicalVolume = nullptr;
300 localPoint = fHistory.GetTopTransfor
301 notKnownContained = false;
302 }
303 }
304 else
305 {
306 fBlockedPhysicalVolume = nullptr;
307 fEntering = false;
308 fEnteredDaughter = false; // Full Step
309 fExiting = false;
310 fExitedMother = false; // Full Step
311 }
312 }
313 //
314 // Search from top of history up through geo
315 // containing volume found:
316 // If on
317 // o OUTSIDE - Back up level, not/no longer
318 // o SURFACE and EXITING - Back up level, se
319 // else
320 // o containing volume found
321 //
322
323 while (notKnownContained)
324 {
325 if ( fHistory.GetTopVolumeType()!=kReplica
326 {
327 targetSolid = fHistory.GetTopVolume()->G
328 localPoint = fHistory.GetTopTransform().
329 insideCode = targetSolid->Inside(localPo
330 #ifdef G4VERBOSE
331 if(( fVerbose == 1 ) && ( fCheck ))
332 {
333 G4String solidResponse = "-kInside-";
334 if (insideCode == kOutside)
335 solidResponse = "-kOutside-";
336 else if (insideCode == kSurface)
337 solidResponse = "-kSurface-";
338 G4cout << "*** G4ITNavigator2::Locate
339 << " Invoked Inside() for s
340 << ". Solid replied: " << soli
341 << " For local point p: " <
342 }
343 #endif
344 }
345 else
346 {
347 insideCode = freplicaNav.BackLocate(fHis
348 fExi
349 // !CARE! if notKnownContained returns f
350 // the containing placement volume of th
351 // will result in the history being back
352 // local point returned is the point in
353 }
354
355
356 if ( insideCode==kOutside )
357 {
358 noLevelsExited++;
359 if ( fHistory.GetDepth() )
360 {
361 fBlockedPhysicalVolume = fHistory.GetT
362 fBlockedReplicaNo = fHistory.GetTopRep
363 fHistory.BackLevel();
364 fExiting = false;
365
366 if( noLevelsExited > 1 )
367 {
368 // The first transformation was done
369 //
370 const G4RotationMatrix* mRot = fBloc
371 if( mRot != nullptr )
372 {
373 fGrandMotherExitNormal *= (*mRot).
374 fChangedGrandMotherRefFrame= true;
375 }
376 }
377 }
378 else
379 {
380 fLastLocatedPointLocal = localPoint;
381 fLocatedOutsideWorld = true;
382 // No extra transformation for ExitN
383 return nullptr; // Have exited
384 }
385 }
386 else
387 if ( insideCode==kSurface )
388 {
389 G4bool isExiting = fExiting;
390 if( (!fExiting)&&considerDirection )
391 {
392 // Figure out whether we are exiting
393 // by using the direction
394 //
395 G4bool directionExiting = false;
396 G4ThreeVector localDirection =
397 fHistory.GetTopTransform().Trans
398
399 // Make sure localPoint in correct r
400 // ( Was it already correct ? Ho
401 //
402 localPoint= fHistory.GetTopTransform
403 if ( fHistory.GetTopVolumeType()!=kR
404 {
405 G4ThreeVector normal = targetSolid
406 directionExiting = normal.dot(loca
407 isExiting = isExiting || direction
408 }
409 }
410 if( isExiting )
411 {
412 noLevelsExited++;
413 if ( fHistory.GetDepth() )
414 {
415 fBlockedPhysicalVolume = fHistory.
416 fBlockedReplicaNo = fHistory.GetTo
417 fHistory.BackLevel();
418 //
419 // Still on surface but exited vol
420 //
421 fValidExitNormal = false;
422
423 if( noLevelsExited > 1 )
424 {
425 // The first transformation was
426 //
427 const G4RotationMatrix* mRot =
428 fBlockedPhysicalVolume->Ge
429 if( mRot != nullptr )
430 {
431 fGrandMotherExitNormal *= (*mR
432 fChangedGrandMotherRefFrame= t
433 }
434 }
435 }
436 else
437 {
438 fLastLocatedPointLocal = localPoin
439 fLocatedOutsideWorld = true;
440 // No extra transformation for E
441 return nullptr; // Have e
442 }
443 }
444 else
445 {
446 notKnownContained=false;
447 }
448 }
449 else
450 {
451 notKnownContained=false;
452 }
453 } // END while (notKnownContained)
454 //
455 // Search downwards until deepest containing
456 // blocking fBlockedPhysicalVolume/BlockedRe
457 //
458 // 3 Cases:
459 //
460 // o Parameterised daughters
461 // =>Must be one G4PVParameterised daughte
462 // o Positioned daughters & voxels
463 // o Positioned daughters & no voxels
464
465 noResult = true; // noResult should be rena
466 // something like enteredL
467 do
468 {
469 // Determine `type' of current mother volu
470 //
471 targetPhysical = fHistory.GetTopVolume();
472 if (targetPhysical == nullptr) { break; }
473 targetLogical = targetPhysical->GetLogical
474 switch( CharacteriseDaughters(targetLogica
475 {
476 case kNormal:
477 if ( targetLogical->GetVoxelHeader() !
478 {
479 noResult = fvoxelNav.LevelLocate(fHi
480 fBl
481 fBl
482 glo
483 pGl
484 con
485 loc
486 }
487 else // do not u
488 {
489 noResult = fnormalNav.LevelLocate(fH
490 fB
491 fB
492 gl
493 pG
494 co
495 lo
496 }
497 break;
498 case kReplica:
499 noResult = freplicaNav.LevelLocate(fHi
500 fBl
501 fBl
502 glo
503 pGl
504 con
505 loc
506 break;
507 case kParameterised:
508 if( GetDaughtersRegularStructureId(tar
509 {
510 noResult = fparamNav.LevelLocate(fHi
511 fBl
512 fBl
513 glo
514 pGl
515 con
516 loc
517 }
518 else // Regular structure
519 {
520 noResult = fregularNav.LevelLocate(f
521 f
522 f
523 g
524 p
525 c
526 l
527 }
528 break;
529 case kExternal:
530 G4Exception("G4ITNavigator2::LocateGlo
531 "GeomNav0001", FatalExcept
532 "Not applicable for extern
533 break;
534 }
535
536 // LevelLocate returns true if it finds a
537 // in which globalPoint is inside (or on t
538
539 if ( noResult )
540 {
541 // Entering a daughter after ascending
542 //
543 // The blocked volume is no longer valid
544 //
545 fBlockedPhysicalVolume = nullptr;
546 fBlockedReplicaNo = -1;
547
548 // fEntering should be false -- else blo
549 // fEnteredDaughter is used for ExitNorm
550 //
551 fEntering = false;
552 fEnteredDaughter = true;
553
554 if( fExitedMother )
555 {
556 G4VPhysicalVolume* enteredPhysical = f
557 const G4RotationMatrix* mRot = entered
558 if( mRot != nullptr )
559 {
560 // Go deeper, i.e. move 'down' in th
561 // Apply direct rotation, not invers
562 //
563 fGrandMotherExitNormal *= (*mRot);
564 fChangedGrandMotherRefFrame= true;
565 }
566 }
567
568 #ifdef G4DEBUG_NAVIGATION
569 if( fVerbose > 2 )
570 {
571 G4VPhysicalVolume* enteredPhysical =
572 G4cout << "*** G4ITNavigator2::Locate
573 G4cout << " Entering volume: " <<
574 << G4endl;
575 }
576 #endif
577 }
578 } while (noResult);
579
580 fLastLocatedPointLocal = localPoint;
581
582 #ifdef G4VERBOSE
583 if( fVerbose >= 4 )
584 {
585 G4long oldcoutPrec = G4cout.precision(8);
586 G4String curPhysVol_Name("None");
587 if (targetPhysical != nullptr) { curPhysV
588 G4cout << " Return value = new volume =
589 G4cout << " ----- Upon exiting:" << G4e
590 PrintState();
591 if( fVerbose >= 5 )
592 {
593 G4cout << "Upon exiting LocateGlobalPoin
594 G4cout << " History = " << G4endl <<
595 }
596 G4cout.precision(oldcoutPrec);
597 }
598 #endif
599
600 fLocatedOutsideWorld= false;
601
602 return targetPhysical;
603 }
604
605 // *******************************************
606 // LocateGlobalPointWithinVolume
607 //
608 // -> the state information of this Navigator
609 // is updated in order to start the next st
610 // -> no check is performed whether pGlobalpoi
611 // original volume (this must be the case).
612 //
613 // Note: a direction could be added to the arg
614 // optional checking (via the old code b
615 // [ This would be done only in verbose
616 // *******************************************
617 //
618 void
619 G4ITNavigator2::LocateGlobalPointWithinVolume(
620 {
621 CheckNavigatorStateIsValid();
622
623 #ifdef G4DEBUG_NAVIGATION
624 // Check: Either step was not limited by a
625 // or else the full step is no long
626 assert( !fWasLimitedByGeometry );
627 #endif
628
629 fLastLocatedPointLocal = ComputeLocalPoint(
630 fLastTriedStepComputation= false;
631 fChangedGrandMotherRefFrame= false; // Fr
632
633 #ifdef G4DEBUG_NAVIGATION
634 if( fVerbose > 2 )
635 {
636 G4cout << "Entering LocateGlobalWithinVol
637 G4cout << fHistory << G4endl;
638 }
639 #endif
640
641 // For the case of Voxel (or Parameterised)
642 // Navigator must be messaged to update its
643
644 // Update the state of the Sub Navigators
645 // - in particular any voxel information th
646 //
647 G4VPhysicalVolume* motherPhysical = fHisto
648 G4LogicalVolume* motherLogical = mother
649 G4SmartVoxelHeader* pVoxelHeader = mother
650
651 if ( fHistory.GetTopVolumeType()!=kReplica
652 {
653 switch( CharacteriseDaughters(motherLogic
654 {
655 case kNormal:
656 if ( pVoxelHeader != nullptr )
657 {
658 fvoxelNav.VoxelLocate( pVoxelHeader
659 }
660 break;
661 case kParameterised:
662 if( GetDaughtersRegularStructureId(mo
663 {
664 // Resets state & returns voxel nod
665 //
666 fparamNav.ParamVoxelLocate( pVoxelH
667 }
668 break;
669 case kReplica:
670 G4Exception("G4ITNavigator2::LocateGl
671 "GeomNav0001", FatalExcep
672 "Not applicable for repli
673 break;
674 case kExternal:
675 G4Exception("G4ITNavigator2::LocateGl
676 "GeomNav0001", FatalExcep
677 "Not applicable for exter
678 break;
679 }
680 }
681
682 // Reset the state variables
683 // - which would have been affected
684 // by the 'equivalent' call to LocateGl
685 // - who's values have been invalidated b
686 //
687 fBlockedPhysicalVolume = nullptr;
688 fBlockedReplicaNo = -1;
689 fEntering = false;
690 fEnteredDaughter = false; // Boundary not
691 fExiting = false;
692 fExitedMother = false; // Boundary not
693 }
694
695 // !>
696
697 void G4ITNavigator2::CheckNavigatorState() con
698 {
699 if(fpNavigatorState == nullptr)
700 {
701 G4ExceptionDescription exceptionDescri
702 exceptionDescription << "The navigator
703 exceptionDescription << "Either NewNav
704 exceptionDescription << "or the provid
705
706 G4Exception("G4ITNavigator::CheckNavig
707 "NavigatorStateNotValid",F
708 return;
709 }
710 }
711
712 G4ITNavigatorState_Lock2* G4ITNavigator2::GetN
713 {
714 return fpNavigatorState;
715 }
716
717 void G4ITNavigator2::SetNavigatorState(G4ITNav
718 {
719 fpNavigatorState = (G4NavigatorState*) nav
720 if(fpNavigatorState != nullptr) SetupHiera
721 // fpSaveState = (G4SaveNavigatorState*) na
722 // if(navState) RestoreSavedState();
723 }
724
725 void G4ITNavigator2::ResetNavigatorState()
726 {
727 fpNavigatorState = nullptr;
728 }
729
730 void G4ITNavigator2::NewNavigatorState()
731 {
732 fpNavigatorState = new G4NavigatorState();
733 if(fTopPhysical == nullptr)
734 {
735 G4ExceptionDescription exceptionDescri
736 exceptionDescription << "No World Volu
737
738 G4Exception("G4ITNavigator::NewNavigat
739 "NoWorldVolume",FatalExcep
740 return;
741 }
742
743 fHistory.SetFirstEntry(fTopPhysical );
744 SetupHierarchy();
745 }
746
747 void G4ITNavigator2::NewNavigatorState(const G
748 {
749 fpNavigatorState = new G4NavigatorState();
750 if(fTopPhysical == nullptr)
751 {
752 G4ExceptionDescription exceptionDescri
753 exceptionDescription << "No World Volu
754
755 G4Exception("G4ITNavigator::NewNavigat
756 "NoWorldVolume",FatalExcep
757 return;
758 }
759
760 fHistory = *h.GetHistory();
761 fLastTriedStepComputation= false; // Redu
762 SetupHierarchy();
763 }
764
765 G4VPhysicalVolume* G4ITNavigator2::NewNavigato
766
767 {
768 fpNavigatorState = new G4NavigatorState();
769
770 if(fTopPhysical == nullptr)
771 {
772 G4ExceptionDescription exceptionDescri
773 exceptionDescription << "No World Volu
774
775 G4Exception("G4ITNavigator::NewNavigat
776 "NoWorldVolume",FatalExcep
777 return nullptr;
778 }
779
780 fHistory.SetFirstEntry(fTopPhysical );
781 SetupHierarchy();
782 return LocateGlobalPointAndSetup(p, &direc
783 }
784
785 // *******************************************
786 // SetSavedState
787 //
788 // Save the state, in case this is a parasitic
789 // Save fValidExitNormal, fExitNormal, fExitin
790 // fBlockedPhysicalVolume, fBlockedReplic
791 // *******************************************
792 //
793 /*
794 void G4ITNavigator2::SetSavedState()
795 {
796 fSaveState = *fpNavigatorState;
797 }
798 */
799
800 /*
801 void G4ITNavigator2::SetSavedState()
802 {
803 // !>
804 // This check can be avoid if instead, at
805 // the IT tracking uses NewNavigatorSate
806 // The normal tracking would just call onc
807
808 // if(fpSaveState == 0)
809 // fpSaveState = new G4SaveNavigatorSta
810 // <!
811
812 // fSaveExitNormal = fExitNormal;
813 fpSaveState->sExitNormal = fExitNormal;
814 fpSaveState->sValidExitNormal = fValidExitNo
815 fpSaveState->sExiting = fExiting;
816 fpSaveState->sEntering = fEntering;
817
818 fpSaveState->spBlockedPhysicalVolume = fBloc
819 fpSaveState->sBlockedReplicaNo = fBlockedRep
820
821 fpSaveState->sLastStepWasZero = fLastStepWas
822
823 // !>
824 fpSaveState->sPreviousSftOrigin = fPreviousS
825 fpSaveState->sPreviousSafety = fPreviousSafe
826 fpSaveState->sNumberZeroSteps = fNumberZeroS
827 fpSaveState->sLocatedOnEdge = fLocatedOnEdge
828 fpSaveState->sWasLimitedByGeometry= fWasLimi
829 fpSaveState->sPushed=fPushed;
830 fpSaveState->sNumberZeroSteps=fNumberZeroSte
831 fpSaveState->sEnteredDaughter = fEnteredDaug
832 fpSaveState->sExitedMother = fExitedMother;
833
834 fpSaveState->sLastLocatedPointLocal = fLastL
835 fpSaveState->sLocatedOutsideWorld = fLocated
836 // <!
837 }
838 */
839 // *******************************************
840 // RestoreSavedState
841 //
842 // Restore the state (in Compute Step), in cas
843 // *******************************************
844 //
845 /*
846 void G4ITNavigator2::RestoreSavedState()
847 {
848 *fpNavigatorState = fSaveState;
849 }
850 */
851 /*
852 void G4ITNavigator2::RestoreSavedState()
853 {
854 fExitNormal = fpSaveState->sExitNormal;
855 fValidExitNormal = fpSaveState->sValidExitNo
856 fExiting = fpSaveState->sExiting;
857 fEntering = fpSaveState->sEntering;
858
859 fBlockedPhysicalVolume = fpSaveState->spBloc
860 fBlockedReplicaNo = fpSaveState->sBlockedRep
861
862 fLastStepWasZero = fpSaveState->sLastStepWas
863
864 // !>
865 fPreviousSftOrigin = fpSaveState->sPreviousS
866 fPreviousSafety = fpSaveState->sPreviousSafe
867 fNumberZeroSteps = fpSaveState->sNumberZeroS
868 fLocatedOnEdge = fpSaveState->sLocatedOnEdge
869 fWasLimitedByGeometry = fpSaveState->sWasLim
870 fPushed = fpSaveState->sPushed;
871 fNumberZeroSteps = fpSaveState->sNumberZeroS
872 fEnteredDaughter= fpSaveState->sEnteredDaugh
873 fExitedMother = fpSaveState->sExitedMother ;
874
875 fLastLocatedPointLocal = fpSaveState->sLastL
876 fLocatedOutsideWorld = fpSaveState->sLocated
877 // <!
878 }
879 */
880 // <!
881
882 // *******************************************
883 // ComputeStep
884 //
885 // Computes the next geometric Step: intersect
886 // mother and `daughter' volumes.
887 //
888 // NOTE:
889 //
890 // Flags on entry:
891 // --------------
892 // fValidExitNormal - Normal of exited volume
893 // coincident boundary)
894 // fExitNormal - Surface normal of exite
895 // fExiting - True if have exited sol
896 //
897 // fBlockedPhysicalVolume - Ptr to exited volu
898 // fBlockedReplicaNo - Replication no of exite
899 // fLastStepWasZero - True if last Step size
900 //
901 // Flags on exit:
902 // -------------
903 // fValidExitNormal - True if surface normal
904 // fExitNormal - Surface normal of exite
905 // reference system
906 // fExiting - True if exiting mother
907 // fEntering - True if entering `daugh
908 // fBlockedPhysicalVolume - Ptr to candidate (
909 // fBlockedReplicaNo - Replication no of candi
910 // fLastStepWasZero - True if this Step size
911 // *******************************************
912 //
913 G4double G4ITNavigator2::ComputeStep( const G4
914 const G4Thr
915 const G4dou
916 G4dou
917 {
918 CheckNavigatorStateIsValid();
919 G4ThreeVector localDirection = ComputeLocalA
920 G4double Step = kInfinity;
921 G4VPhysicalVolume *motherPhysical = fHistor
922 G4LogicalVolume *motherLogical = motherPhysi
923
924 // All state relating to exiting normals mus
925 //
926 fExitNormalGlobalFrame= G4ThreeVector( 0., 0
927 // Reset value - to erase its memory
928 fChangedGrandMotherRefFrame= false;
929 // Reset - used for local exit normal
930 fGrandMotherExitNormal= G4ThreeVector( 0., 0
931 fCalculatedExitNormal = false;
932 // Reset for new step
933
934 fLastTriedStepComputation= true;
935
936 #ifdef G4VERBOSE
937 if( fVerbose > 0 )
938 {
939 G4cout << "*** G4ITNavigator2::ComputeStep
940 G4cout << " Volume = " << motherPhysica
941 << " - Proposed step length = " <<
942 << G4endl;
943 #ifdef G4DEBUG_NAVIGATION
944 if( fVerbose >= 2 )
945 {
946 G4cout << " Called with the arguments:
947 << " Globalpoint = " << std::set
948 << " Direction = " << std::set
949 if( fVerbose >= 4 )
950 {
951 G4cout << " ---- Upon entering : Stat
952 PrintState();
953 }
954 }
955 #endif
956 }
957 #endif
958
959 G4ThreeVector newLocalPoint = ComputeLocalPo
960 if( newLocalPoint != fLastLocatedPointLocal
961 {
962 // Check whether the relocation is within
963 //
964 G4ThreeVector oldLocalPoint = fLastLocated
965 G4double moveLenSq = (newLocalPoint-oldLoc
966
967 if ( moveLenSq >= kCarTolerance*kCarTolera
968 {
969 #ifdef G4VERBOSE
970 ComputeStepLog(pGlobalpoint, moveLenSq);
971 #endif
972 // Relocate the point within the same vo
973 //
974 LocateGlobalPointWithinVolume( pGlobalpo
975 fLastTriedStepComputation= true; //
976 }
977 }
978 if ( fHistory.GetTopVolumeType()!=kReplica )
979 {
980 switch( CharacteriseDaughters(motherLogica
981 {
982 case kNormal:
983 if ( motherLogical->GetVoxelHeader() !
984 {
985 LocateGlobalPointWithinVolume(pGloba
986 Step = fvoxelNav.ComputeStep(fLastLo
987 localDi
988 pCurren
989 pNewSaf
990 fHistor
991 fValidE
992 fExitNo
993 fExitin
994 fEnteri
995 &fBlock
996 fBlocke
997
998 }
999 else
1000 {
1001 if( motherPhysical->GetRegularStruc
1002 {
1003 Step = fnormalNav.ComputeStep(fLa
1004 loc
1005 pCu
1006 pNe
1007 fHi
1008 fVa
1009 fEx
1010 fEx
1011 fEn
1012 &fB
1013 fBl
1014 }
1015 else // Regular (non-voxelised) st
1016 {
1017 LocateGlobalPointAndSetup( pGloba
1018 fLastTriedStepComputation= true;
1019 //
1020 // if physical process limits the
1021 // one given by ComputeStepSkippi
1022 // point will be wrongly calculat
1023
1024 // There is a problem: when msc l
1025 // assigned wrongly to phantom in
1026 // of the container volume). Then
1027 // reset the history topvolume to
1028 //
1029 if(fHistory.GetTopVolume()->GetRe
1030 {
1031 G4Exception("G4ITNavigator2::Co
1032 "GeomNav1001", Just
1033 "Point is relocated in voxels
1034 Step = fnormalNav.ComputeStep(f
1035 l
1036 p
1037 p
1038 f
1039 f
1040 f
1041 f
1042 f
1043 &
1044 f
1045 }
1046 else
1047 {
1048 Step = fregularNav.
1049 ComputeStepSkippingEqualMa
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061 }
1062 }
1063 }
1064 break;
1065 case kParameterised:
1066 if( GetDaughtersRegularStructureId(mo
1067 {
1068 Step = fparamNav.ComputeStep(fLastL
1069 localD
1070 pCurre
1071 pNewSa
1072 fHisto
1073 fValid
1074 fExitN
1075 fExiti
1076 fEnter
1077 &fBloc
1078 fBlock
1079 }
1080 else // Regular structure
1081 {
1082 Step = fregularNav.ComputeStep(fLas
1083 loca
1084 pCur
1085 pNew
1086 fHis
1087 fVal
1088 fExi
1089 fExi
1090 fEnt
1091 &fBl
1092 fBlo
1093 }
1094 break;
1095 case kReplica:
1096 G4Exception("G4ITNavigator2::ComputeS
1097 FatalException, "Not appl
1098 break;
1099 case kExternal:
1100 G4Exception("G4ITNavigator2::ComputeS
1101 FatalException, "Not appl
1102 break;
1103 }
1104 }
1105 else
1106 {
1107 // In the case of a replica, it must hand
1108 // edge/corner problem by itself
1109 //
1110 G4bool exitingReplica = fExitedMother;
1111 G4bool calculatedExitNormal;
1112 Step = freplicaNav.ComputeStep(pGlobalpoi
1113 pDirection
1114 fLastLocat
1115 localDirec
1116 pCurrentPr
1117 pNewSafety
1118 fHistory,
1119 fValidExit
1120 calculated
1121 fExitNorma
1122 exitingRep
1123 fEntering,
1124 &fBlockedP
1125 fBlockedRe
1126 fExiting= exitingReplica;
1127 fCalculatedExitNormal= calculatedExitNorm
1128 }
1129
1130
1131 // G4cout << " !!!! Step = " << Step << G4en
1132
1133 // Remember last safety origin & value.
1134 //
1135 fPreviousSftOrigin = pGlobalpoint;
1136 fPreviousSafety = pNewSafety;
1137
1138 // Count zero steps - one can occur due to
1139 // - one, two (or a few) c
1140 // volumes
1141 // - more than two is like
1142 // description or the Na
1143
1144 // Rule of thumb: likely at an Edge if two
1145 // because at least two cand
1146 // checked
1147 //
1148 fLocatedOnEdge = fLastStepWasZero && (Ste
1149 fLastStepWasZero = (Step==0.0);
1150 if (fPushed) { fPushed = fLastStepWasZero;
1151
1152 // Handle large number of consecutive zero
1153 //
1154 if ( fLastStepWasZero )
1155 {
1156 fNumberZeroSteps++;
1157 #ifdef G4DEBUG_NAVIGATION
1158 if( fNumberZeroSteps > 1 )
1159 {
1160 G4cout << "G4ITNavigator2::ComputeStep
1161 << fNumberZeroSteps
1162 << " at " << pGlobalpoint
1163 << " in volume " << motherPhysi
1164 << G4endl;
1165 }
1166 #endif
1167 if( fNumberZeroSteps > fActionThreshold_N
1168 {
1169 // Act to recover this stuck track. Pu
1170 //
1171 Step += 100*kCarTolerance;
1172 #ifdef G4VERBOSE
1173 if ((!fPushed) && (fWarnPush))
1174 {
1175 std::ostringstream message;
1176 message << "Track stuck or not movin
1177 << " Track stuck, n
1178 << fNumberZeroSteps << " ste
1179 << " in volume -" <
1180 << "- at point " << pGlobalp
1181 << " direction: " <
1182 << " Potential geom
1183 << G4endl
1184 << " Trying pushing
1185 G4Exception("G4ITNavigator2::Compute
1186 JustWarning, message, "P
1187 }
1188 #endif
1189 fPushed = true;
1190 }
1191 if( fNumberZeroSteps > fAbandonThreshold_
1192 {
1193 // Must kill this stuck track
1194 //
1195 std::ostringstream message;
1196 message << "Stuck Track: potential geom
1197 << G4endl
1198 << " Track stuck, not mo
1199 << fNumberZeroSteps << " steps"
1200 << " in volume -" << mot
1201 << "- at point " << pGlobalpoin
1202 << " direction: " << pDi
1203 motherPhysical->CheckOverlaps(5000, 0.0
1204 G4Exception("G4ITNavigator2::ComputeSte
1205 EventMustBeAborted, message
1206 }
1207 }
1208 else
1209 {
1210 if (!fPushed) fNumberZeroSteps = 0;
1211 }
1212
1213 fEnteredDaughter = fEntering; // I expect
1214 fExitedMother = fExiting;
1215
1216 fStepEndPoint = pGlobalpoint
1217 + std::min(Step,pCurrentPropo
1218 fLastStepEndPointLocal = fLastLocatedPointL
1219
1220 if( fExiting )
1221 {
1222 #ifdef G4DEBUG_NAVIGATION
1223 if( fVerbose > 2 )
1224 {
1225 G4cout << " At G4Nav CompStep End - if(
1226 << " fValidExitNormal = " << fVa
1227 G4cout << " fExitNormal= " << fExitNorm
1228 }
1229 #endif
1230
1231 if(fValidExitNormal || fCalculatedExitNor
1232 {
1233 if ( fHistory.GetTopVolumeType()!=kRep
1234 {
1235 // Convention: fExitNormal is in the
1236 //
1237 fGrandMotherExitNormal= fExitNormal;
1238 fCalculatedExitNormal= true;
1239 }
1240 else
1241 {
1242 fGrandMotherExitNormal = fExitNormal;
1243 }
1244 }
1245 else
1246 {
1247 // We must calculate the normal anyway
1248 //
1249 G4ThreeVector finalLocalPoint =
1250 fLastLocatedPointLocal + localDir
1251
1252 if ( fHistory.GetTopVolumeType()!=kRep
1253 {
1254 // Find normal in the 'mother' coordi
1255 //
1256 G4ThreeVector exitNormalMotherFrame=
1257 motherLogical->GetSolid()->Surface
1258
1259 // Transform it to the 'grand-mother'
1260 //
1261 const G4RotationMatrix* mRot = mother
1262 if( mRot != nullptr )
1263 {
1264 fChangedGrandMotherRefFrame= true;
1265 fGrandMotherExitNormal = (*mRot).in
1266 }
1267 else
1268 {
1269 fGrandMotherExitNormal = exitNormal
1270 }
1271
1272 // Do not set fValidExitNormal -- thi
1273 // that the solid is convex!
1274 //
1275 fCalculatedExitNormal= true;
1276 }
1277 else
1278 {
1279 fCalculatedExitNormal = false;
1280 //
1281 // Nothing can be done at this stage
1282 // Replica Navigation must have calcu
1283 // already.
1284 // Cases: mother is not convex, and e
1285
1286 #ifdef G4DEBUG_NAVIGATION
1287 G4ExceptionDescription desc;
1288
1289 desc << "Problem in ComputeStep: Rep
1290 << " valid exit Normal. " << G4e
1291 desc << " Do not know how calculate i
1292 desc << " Location = " << finalLo
1293 desc << " Volume name = " << motherP
1294 << " copy/replica No = " << mot
1295 G4Exception("G4ITNavigator2::ComputeS
1296 JustWarning, desc, "Norma
1297 #endif
1298 }
1299 }
1300
1301 // Now transform it to the global referen
1302 //
1303 if( fValidExitNormal || fCalculatedExitNo
1304 {
1305 auto depth= (G4int)fHistory.GetDepth()
1306 if( depth > 0 )
1307 {
1308 G4AffineTransform GrandMotherToGlobal
1309 fHistory.GetTransform(depth-1).Inve
1310 fExitNormalGlobalFrame =
1311 GrandMotherToGlobalTransf.Transform
1312 }
1313 else
1314 {
1315 fExitNormalGlobalFrame= fGrandMotherE
1316 }
1317 }
1318 else
1319 {
1320 fExitNormalGlobalFrame= G4ThreeVector(
1321 }
1322 }
1323
1324 if( (Step == pCurrentProposedStepLength) &&
1325 {
1326 // This if Step is not really limited by
1327 // The Navigator is obliged to return "in
1328 //
1329 Step = kInfinity;
1330 }
1331
1332 #ifdef G4VERBOSE
1333 if( fVerbose > 1 )
1334 {
1335 if( fVerbose >= 4 )
1336 {
1337 G4cout << " ----- Upon exiting :" <<
1338 PrintState();
1339 }
1340 G4cout << " Returned step= " << Step;
1341 if( fVerbose > 5 ) G4cout << G4endl;
1342 if( Step == kInfinity )
1343 {
1344 G4cout << " Requested step= " << pCurr
1345 if( fVerbose > 5) G4cout << G4endl;
1346 }
1347 G4cout << " Safety = " << pNewSafety <<
1348 }
1349 #endif
1350
1351 return Step;
1352 }
1353
1354 // ******************************************
1355 // CheckNextStep
1356 //
1357 // Compute the step without altering the navi
1358 // ******************************************
1359 //
1360 G4double G4ITNavigator2::CheckNextStep( const
1361 const
1362 const
1363
1364 {
1365 CheckNavigatorStateIsValid();
1366 G4double step;
1367
1368 // Save the state, for this parasitic call
1369 //
1370 //SetSavedState();
1371 // G4SaveNavigatorState savedState(fpNavigat
1372 G4NavigatorState savedState(*fpNavigatorSta
1373
1374 step = ComputeStep ( pGlobalpoint,
1375 pDirection,
1376 pCurrentProposedStepLe
1377 pNewSafety );
1378
1379 // If a parasitic call, then attempt to res
1380 //
1381 *fpNavigatorState = savedState;
1382 //RestoreSavedState();
1383 // NOTE: the state of the current subnaviga
1384 // ***> TODO: restore subnavigator state
1385 // if( last_located) Need
1386 // if( last_computed step) Need
1387
1388
1389 return step;
1390 }
1391
1392 // ******************************************
1393 // ResetState
1394 //
1395 // Resets stack and minimum of navigator stat
1396 // ******************************************
1397 //
1398 void G4ITNavigator2::ResetState()
1399 {
1400 fWasLimitedByGeometry = false;
1401 fEntering = false;
1402 fExiting = false;
1403 fLocatedOnEdge = false;
1404 fLastStepWasZero = false;
1405 fEnteredDaughter = false;
1406 fExitedMother = false;
1407 fPushed = false;
1408
1409 fValidExitNormal = false;
1410 fChangedGrandMotherRefFrame= false;
1411 fCalculatedExitNormal = false;
1412
1413 fExitNormal = G4ThreeVector(0,0,
1414 fGrandMotherExitNormal = G4ThreeVector(0,0,
1415 fExitNormalGlobalFrame = G4ThreeVector(0,0,
1416
1417 fPreviousSftOrigin = G4ThreeVector(0,0,
1418 fPreviousSafety = 0.0;
1419
1420 fNumberZeroSteps = 0;
1421
1422 fBlockedPhysicalVolume = nullptr;
1423 fBlockedReplicaNo = -1;
1424
1425 fLastLocatedPointLocal = G4ThreeVector( kIn
1426 fLocatedOutsideWorld = false;
1427 }
1428
1429 // ******************************************
1430 // SetupHierarchy
1431 //
1432 // Renavigates & resets hierarchy described b
1433 // o Reset volumes
1434 // o Recompute transforms and/or solids of re
1435 // ******************************************
1436 //
1437 void G4ITNavigator2::SetupHierarchy()
1438 {
1439 G4int i;
1440 const auto cdepth = (G4int)fHistory.GetDep
1441 G4VPhysicalVolume *current;
1442 G4VSolid *pSolid;
1443 G4VPVParameterisation *pParam;
1444
1445 for ( i=1; i<=cdepth; i++ )
1446 {
1447 current = fHistory.GetVolume(i);
1448 switch ( fHistory.GetVolumeType(i) )
1449 {
1450 case kNormal:
1451 break;
1452 case kReplica:
1453 freplicaNav.ComputeTransformation(fHi
1454 break;
1455 case kParameterised:
1456 {
1457 G4int replicaNo;
1458 pParam = current->GetParameterisation
1459 replicaNo = fHistory.GetReplicaNo(i);
1460 pSolid = pParam->ComputeSolid(replica
1461
1462 // Set up dimensions & transform in s
1463 //
1464 pSolid->ComputeDimensions(pParam, rep
1465 pParam->ComputeTransformation(replica
1466
1467 G4TouchableHistory *pTouchable= nullp
1468 if( pParam->IsNested() )
1469 {
1470 pTouchable= new G4TouchableHistory
1471 pTouchable->MoveUpHistory(); // Mo
1472 // Adequate only if Nested at th
1473 // To extend to other cases:
1474 // pTouchable->MoveUpHistory(cdept
1475 // Move to the parent level of *
1476 // Could replace this line and c
1477 // c-tor for History(levels to d
1478 }
1479 // Set up the correct solid and mater
1480 //
1481 G4LogicalVolume *pLogical = current->
1482 pLogical->SetSolid( pSolid );
1483 pLogical->UpdateMaterial( pParam ->
1484 ComputeMaterial(replicaNo, current,
1485 delete pTouchable;
1486 }
1487 break;
1488
1489 case kExternal:
1490 G4Exception("G4ITNavigator2::SetupHie
1491 "GeomNav0001", FatalExcep
1492 "Not applicable for exter
1493 break;
1494
1495 }
1496 }
1497 }
1498
1499 // ******************************************
1500 // GetLocalExitNormal
1501 //
1502 // Obtains the Normal vector to a surface (in
1503 // pointing out of previous volume and into c
1504 // ******************************************
1505 //
1506 G4ThreeVector G4ITNavigator2::GetLocalExitNor
1507 {
1508 CheckNavigatorStateIsValid();
1509 G4ThreeVector ExitNormal(0.,0.,0.);
1510 G4VSolid *currentSolid=nullptr;
1511 G4LogicalVolume *candidateLogical;
1512 if ( fLastTriedStepComputation )
1513 {
1514 // use fLastLocatedPointLocal and next ca
1515 //
1516 G4ThreeVector nextSolidExitNormal(0.,0.,0
1517
1518 if( fEntering && (fBlockedPhysicalVolume!
1519 {
1520 candidateLogical= fBlockedPhysicalVolum
1521 if( candidateLogical != nullptr )
1522 {
1523 // fLastStepEndPointLocal is in the c
1524 // we need it in the daughter's coord
1525
1526 // The following code should also wor
1527 {
1528 // First transform fLastLocatedPoin
1529 // coordinates
1530 //
1531 G4AffineTransform MotherToDaughterT
1532 GetMotherToDaughterTransform( fBl
1533 fBl
1534 Vol
1535 G4ThreeVector daughterPointOwnLocal
1536 MotherToDaughterTransform.Transfo
1537
1538 // OK if it is a parameterised volu
1539 //
1540 EInside inSideIt;
1541 G4bool onSurface;
1542 G4double safety= -1.0;
1543 currentSolid= candidateLogical->Get
1544 inSideIt = currentSolid->Inside(d
1545 onSurface = (inSideIt == kSurface)
1546 if( ! onSurface )
1547 {
1548 if( inSideIt == kOutside )
1549 {
1550 safety = (currentSolid->Distanc
1551 onSurface = safety < 100.0 * kC
1552 }
1553 else if (inSideIt == kInside )
1554 {
1555 safety = (currentSolid->Distanc
1556 onSurface = safety < 100.0 * kC
1557 }
1558 }
1559
1560 if( onSurface )
1561 {
1562 nextSolidExitNormal =
1563 currentSolid->SurfaceNormal(dau
1564
1565 // Entering the solid ==> opposit
1566 //
1567 ExitNormal = -nextSolidExitNormal
1568 fCalculatedExitNormal= true;
1569 }
1570 else
1571 {
1572 #ifdef G4VERBOSE
1573 if(( fVerbose == 1 ) && ( fCheck
1574 {
1575 std::ostringstream message;
1576 message << "Point not on surfac
1577 << " Point =
1578 << daughterPointOwnLoca
1579 << " Physical volume =
1580 << fBlockedPhysicalVolu
1581 << " Logical volume =
1582 << candidateLogical->Ge
1583 << " Solid =
1584 << " Type =
1585 << currentSolid->GetEnt
1586 << *currentSolid << G4e
1587 if( inSideIt == kOutside )
1588 {
1589 message << "Point is Outside.
1590 << " Safety (from ou
1591 }
1592 else // if( inSideIt == kInside
1593 {
1594 message << "Point is Inside.
1595 << " Safety (from in
1596 }
1597 G4Exception("G4ITNavigator2::Ge
1598 JustWarning, messag
1599 }
1600 #endif
1601 }
1602 *valid = onSurface; // was =tru
1603 }
1604 }
1605 }
1606 else if ( fExiting )
1607 {
1608 ExitNormal = fGrandMotherExitNormal;
1609 *valid = true;
1610 fCalculatedExitNormal= true; // Should
1611 }
1612 else // i.e. ( fBlockedPhysicalVolume =
1613 {
1614 *valid = false;
1615 G4Exception("G4ITNavigator2::GetLocalEx
1616 "GeomNav0003", JustWarning,
1617 "Incorrect call to GetLocal
1618 }
1619 }
1620 else // ( ! fLastTriedStepComputation ) ie
1621 {
1622 if ( EnteredDaughterVolume() )
1623 {
1624 G4VSolid* daughterSolid =fHistory.GetTo
1625
1626 ExitNormal= -(daughterSolid->SurfaceNor
1627 if( std::fabs(ExitNormal.mag2()-1.0 ) >
1628 {
1629 G4ExceptionDescription desc;
1630 desc << " Parameters of solid: " << *
1631 << " Point for surface = " << fL
1632 G4Exception("G4ITNavigator2::GetLocal
1633 "GeomNav0003", FatalExcep
1634 "Surface Normal returned
1635 }
1636 fCalculatedExitNormal= true;
1637 *valid = true;
1638 }
1639 else
1640 {
1641 if( fExitedMother )
1642 {
1643 ExitNormal = fGrandMotherExitNormal;
1644 *valid = true;
1645 fCalculatedExitNormal= true;
1646 }
1647 else // We are not at a boundary. Exit
1648 {
1649 *valid = false;
1650 fCalculatedExitNormal= false;
1651 G4ExceptionDescription message;
1652 message << "Function called when *NOT
1653 G4Exception("G4ITNavigator2::GetLocal
1654 "GeomNav0003", JustWarnin
1655 }
1656 }
1657 }
1658 return ExitNormal;
1659 }
1660
1661 // ******************************************
1662 // GetMotherToDaughterTransform
1663 //
1664 // Obtains the mother to daughter affine tran
1665 // ******************************************
1666 //
1667 G4AffineTransform
1668 G4ITNavigator2::GetMotherToDaughterTransform(
1669 G4
1670 EV
1671 {
1672 CheckNavigatorStateIsValid();
1673 switch (enteringVolumeType)
1674 {
1675 case kNormal: // Nothing is needed to pr
1676 break; // It is stored already in
1677 case kReplica: // Sets the transform in t
1678 G4Exception("G4ITNavigator2::GetMotherT
1679 "GeomNav0001", FatalExcepti
1680 "Method NOT Implemented yet
1681 break;
1682 case kParameterised:
1683 if( pEnteringPhysVol->GetRegularStructu
1684 {
1685 G4VPVParameterisation *pParam =
1686 pEnteringPhysVol->GetParameterisati
1687 G4VSolid* pSolid =
1688 pParam->ComputeSolid(enteringReplic
1689 pSolid->ComputeDimensions(pParam, ent
1690
1691 // Sets the transform in the Paramete
1692 //
1693 pParam->ComputeTransformation(enterin
1694
1695 // Set the correct solid and material
1696 //
1697 G4LogicalVolume* pLogical = pEntering
1698 pLogical->SetSolid( pSolid );
1699 }
1700 break;
1701 case kExternal:
1702 G4Exception("G4ITNavigator2::GetMothe
1703 "GeomNav0001", FatalExcep
1704 "Not applicable for exter
1705 break;
1706 }
1707 return G4AffineTransform(pEnteringPhysVol->
1708 pEnteringPhysVol->
1709 }
1710
1711 // ******************************************
1712 // GetLocalExitNormalAndCheck
1713 //
1714 // Obtains the Normal vector to a surface (in
1715 // pointing out of previous volume and into c
1716 // checks the current point against expected
1717 // ******************************************
1718 //
1719 G4ThreeVector G4ITNavigator2::
1720 GetLocalExitNormalAndCheck(
1721 #ifdef G4DEBUG_NAVIGATION
1722 const G4ThreeVecto
1723 #else
1724 const G4ThreeVecto
1725 #endif
1726 G4bool*
1727 {
1728 CheckNavigatorStateIsValid();
1729 #ifdef G4DEBUG_NAVIGATION
1730 // Check Current point against expected 'lo
1731 //
1732 if ( fLastTriedStepComputation )
1733 {
1734 G4ThreeVector ExpectedBoundaryPointLocal;
1735
1736 const G4AffineTransform& GlobalToLocal= G
1737 ExpectedBoundaryPointLocal =
1738 GlobalToLocal.TransformPoint( ExpectedB
1739
1740 // Add here: Comparison against expected
1741 // i.e. the endpoint of Comput
1742 }
1743 #endif
1744
1745 return GetLocalExitNormal( pValid);
1746 }
1747
1748 // ******************************************
1749 // GetGlobalExitNormal
1750 //
1751 // Obtains the Normal vector to a surface (in
1752 // pointing out of previous volume and into c
1753 // ******************************************
1754 //
1755 G4ThreeVector
1756 G4ITNavigator2::GetGlobalExitNormal(const G4T
1757 G4bool
1758 {
1759 CheckNavigatorStateIsValid();
1760 G4bool validNormal;
1761 G4ThreeVector localNormal, globalNormal;
1762 G4bool usingStored;
1763
1764 usingStored=
1765 fCalculatedExitNormal &&
1766 ( ( fLastTriedStepComputation && fExiti
1767 ||
1768 (!fLastTriedStepComputation
1769 && (IntersectPointGlobal-fStepEndPo
1770 < (10.0*kCarTolerance*kCarTolera
1771 ) // Calculated it 'just' before & t
1772 // but it did not move position
1773 );
1774
1775 if( usingStored )
1776 {
1777 // This was computed in ComputeStep -- an
1778 //
1779 globalNormal = fExitNormalGlobalFrame;
1780 G4double normMag2 = globalNormal.mag2();
1781 if( std::fabs ( normMag2 - 1.0 ) < perMil
1782 {
1783 *pNormalCalculated = true; // ComputeS
1784 // (fExitin
1785 }
1786 else
1787 {
1788 G4ExceptionDescription message;
1789
1790 message << " ERROR> Expected normal-gl
1791 << " - but |normal| = " << s
1792 << " - and |normal|^ = " <<
1793 << " which differs from 1.0 by
1794 << " n = " << fExitNormalGlo
1795 message << "==========================
1796 << G4endl;
1797 G4int oldVerbose = fVerbose;
1798 fVerbose=4;
1799 message << " State of Navigator: " <
1800 message << *this << G4endl;
1801 fVerbose = oldVerbose;
1802 message << "==========================
1803 << G4endl;
1804
1805 G4Exception("G4ITNavigator2::GetGlobal
1806 "GeomNav0003",JustWarning,
1807 "Value obtained from stored glo
1808
1809 // (Re)Compute it now -- as either it
1810 //
1811 localNormal = GetLocalExitNormalAndChe
1812
1813 *pNormalCalculated = fCalculatedExitNo
1814
1815 G4AffineTransform localToGlobal = GetL
1816 globalNormal = localToGlobal.Transform
1817 }
1818 }
1819 else
1820 {
1821 localNormal = GetLocalExitNormalAndCheck(
1822 *pNormalCalculated = fCalculatedExitNorma
1823
1824 #ifdef G4DEBUG_NAVIGATION
1825 usingStored= false;
1826
1827 if( (!validNormal) && !fCalculatedExitNor
1828 {
1829 G4ExceptionDescription edN;
1830 edN << " Calculated = " << fCalculated
1831 edN << " Entering= " << fEntering <<
1832 G4int oldVerbose= this->GetVerboseLevel
1833 this->SetVerboseLevel(4);
1834 edN << " State of Navigator: " << G4e
1835 edN << *this << G4endl;
1836 this->SetVerboseLevel( oldVerbose );
1837
1838 G4Exception("G4ITNavigator2::GetGlobalE
1839 "GeomNav0003", JustWarning,
1840 "LocalExitNormalAndCheck()
1841 }
1842 #endif
1843
1844 G4double localMag2= localNormal.mag2();
1845 if( validNormal && (std::fabs(localMag2-
1846 {
1847 G4ExceptionDescription edN;
1848
1849 edN << "G4ITNavigator2::GetGlobalExitN
1850 << " Using Local Normal - from ca
1851 << G4endl
1852 << " Local Exit Normal : " << "
1853 << " vec = " << localNormal << G4e
1854 << " Global Exit Normal : " << "
1855 << " vec = " << globalNormal << G4
1856 edN << " Calculated It = " << fC
1857
1858 G4Exception("G4ITNavigator2::GetGlobal
1859 "GeomNav0003",JustWarning,
1860 "Value obtained from new l
1861 localNormal = localNormal.unit(); // S
1862 }
1863 G4AffineTransform localToGlobal = GetLoc
1864 globalNormal = localToGlobal.TransformAx
1865 }
1866
1867 #ifdef G4DEBUG_NAVIGATION
1868 if( usingStored )
1869 {
1870 G4ThreeVector globalNormAgn;
1871
1872 localNormal= GetLocalExitNormalAndCheck(I
1873
1874 G4AffineTransform localToGlobal = GetLoca
1875 globalNormAgn = localToGlobal.TransformAx
1876
1877 // Check the value computed against fExit
1878 G4ThreeVector diffNorm = globalNormAgn -
1879 if( diffNorm.mag2() > perMillion*CLHEP::p
1880 {
1881 G4ExceptionDescription edDfn;
1882 edDfn << "Found difference in normals i
1883 << "- when Get is called after Co
1884 edDfn << " Magnitude of diff = "
1885 edDfn << " Normal stored (Global)
1886 << G4endl;
1887 edDfn << " Global Computed from Local
1888 G4Exception("G4ITNavigator::GetGlobalEx
1889 JustWarning, edDfn);
1890 }
1891 }
1892 #endif
1893
1894 return globalNormal;
1895 }
1896
1897 // To make the new Voxel Safety the default,
1898 #define G4NEW_SAFETY 1
1899
1900 // ******************************************
1901 // ComputeSafety
1902 //
1903 // It assumes that it will be
1904 // i) called at the Point in the same volume
1905 // ComputeStep.
1906 // ii) after (or at the end of) ComputeStep O
1907 // ******************************************
1908 //
1909 G4double G4ITNavigator2::ComputeSafety( const
1910 const G4
1911 const G4
1912 {
1913 CheckNavigatorStateIsValid();
1914 G4double newSafety = 0.0;
1915
1916 #ifdef G4DEBUG_NAVIGATION
1917 G4int oldcoutPrec = G4cout.precision(8);
1918 if( fVerbose > 0 )
1919 {
1920 G4cout << "*** G4ITNavigator2::ComputeSaf
1921 << " Called at point: " << pGlo
1922
1923 G4VPhysicalVolume *motherPhysical = fHis
1924 G4cout << " Volume = " << motherPhysic
1925 << " - Maximum length = " << pMaxL
1926 if( fVerbose >= 4 )
1927 {
1928 G4cout << " ----- Upon entering Com
1929 PrintState();
1930 }
1931 }
1932 #endif
1933
1934 G4SaveNavigatorState* savedState(nullptr);
1935
1936 G4double distEndpointSq = (pGlobalpoint-fSt
1937 G4bool stayedOnEndpoint = distEndpointSq
1938 G4bool endpointOnSurface = fEnteredDaught
1939
1940 if( endpointOnSurface && stayedOnEndpoint )
1941 {
1942 #ifdef G4DEBUG_NAVIGATION
1943 if( fVerbose >= 2 )
1944 {
1945 G4cout << " G4Navigator::ComputeSa
1946 << pGlobalpoint << " - is on surface
1947 if( fEnteredDaughter ) { G4cout << "
1948 if( fExitedMother ) { G4cout << "
1949 G4cout << G4endl;
1950 G4cout << " EndPoint was = " << fStep
1951 }
1952 #endif
1953 newSafety = 0.0;
1954 // return newSafety;
1955 }
1956 else // if( !(endpointOnSurface && stayedOn
1957 {
1958
1959 if (keepState)
1960 {
1961 // SetSavedState();
1962 savedState = new G4SaveNavigatorState(fpN
1963 }
1964
1965 // Pseudo-relocate to this point (updates
1966 //
1967 LocateGlobalPointWithinVolume( pGlobalpoi
1968 // --->> DANGER: Side effects on sub-na
1969 // Could be replaced again by 'gr
1970 // locates (similar side-effects,
1971 // Solutions:
1972 // 1) Re-locate (to where?)
1973 // 2) Insure that the methods us
1974 // does a relocation (if info
1975
1976 #ifdef G4DEBUG_NAVIGATION
1977 if( fVerbose >= 2 )
1978 {
1979 G4cout << " G4ITNavigator2::ComputeSaf
1980 << pGlobalpoint << G4endl;
1981 }
1982 #endif
1983 G4VPhysicalVolume *motherPhysical = fHist
1984 G4LogicalVolume *motherLogical = motherPh
1985 G4SmartVoxelHeader* pVoxelHeader = mother
1986 G4ThreeVector localPoint = ComputeLocalPo
1987
1988 if ( fHistory.GetTopVolumeType()!=kReplic
1989 {
1990 switch(CharacteriseDaughters(motherLogi
1991 {
1992 case kNormal:
1993 if ( pVoxelHeader != nullptr )
1994 {
1995 #ifdef G4NEW_SAFETY
1996 G4double safetyTwo = fpVoxelSafet
1997 *m
1998 newSafety= safetyTwo; // Faster
1999 #else
2000 G4double safetyOldVoxel;
2001 LocateGlobalPointWithinVolume(pGl
2002 safetyOldVoxel =
2003 fvoxelNav.ComputeSafety(localPo
2004 newSafety= safetyOldVoxel;
2005 #endif
2006 }
2007 else
2008 {
2009 newSafety=fnormalNav.ComputeSafet
2010 }
2011 break;
2012 case kParameterised:
2013 if( GetDaughtersRegularStructureId(
2014 {
2015 newSafety=fparamNav.ComputeSafety
2016 }
2017 else // Regular structure
2018 {
2019 newSafety=fregularNav.ComputeSafe
2020 }
2021 break;
2022 case kReplica:
2023 G4Exception("G4ITNavigator2::Comput
2024 FatalException, "Not ap
2025 break;
2026 case kExternal:
2027 G4Exception("G4ITNavigator2::Comput
2028 "GeomNav0001", FatalExc
2029 "Not applicable for ext
2030 break;
2031 }
2032 }
2033 else
2034 {
2035 newSafety = freplicaNav.ComputeSafety(p
2036 f
2037 }
2038
2039 if (keepState)
2040 {
2041 *fpNavigatorState = *savedState;
2042 delete savedState;
2043 // RestoreSavedState();
2044 // This now overwrites the values of the
2045 }
2046
2047 // Remember last safety origin & value
2048 //
2049 // We overwrite the Safety 'sphere' - kee
2050 fPreviousSftOrigin = pGlobalpoint;
2051 fPreviousSafety = newSafety;
2052 }
2053
2054 #ifdef G4DEBUG_NAVIGATION
2055 if( fVerbose > 1 )
2056 {
2057 G4cout << " ---- Exiting ComputeSafety
2058 if( fVerbose > 2 ) { PrintState(); }
2059 G4cout << " Returned value of Safety =
2060 }
2061 G4cout.precision(oldcoutPrec);
2062 #endif
2063
2064 return newSafety;
2065 }
2066
2067
2068 // ******************************************
2069 // RecheckDistanceToCurrentBoundary
2070 //
2071 // Trial method for checking potential displa
2072 // Check position aDisplacedGlobalpoint, to s
2073 // current volume (mother).
2074 // If in mother, check distance to boundary a
2075 // If in entering daughter, check distance ba
2076 // NOTE:
2077 // Can be called only after ComputeStep() is
2078 // Deals only with current volume (and potent
2079 // Caveats
2080 // First VERSION: Does not consider daughter
2081 // neither for checking whether it has
2082 // nor for determining distance to exit
2083 // ******************************************
2084 //
2085 G4bool G4ITNavigator2::RecheckDistanceToCurre
2086 const G4ThreeVector &aDi
2087 const G4ThreeVector &aNe
2088 const G4double ProposedM
2089 G4double *prDistance,
2090 G4double *prNewSafety) c
2091 {
2092 G4ThreeVector localPosition = ComputeLocal
2093 G4ThreeVector localDirection = ComputeLocal
2094 // G4double Step = kInfinity;
2095
2096 G4bool validExitNormal;
2097 G4ThreeVector exitNormal;
2098 // Check against mother solid
2099 G4VPhysicalVolume *motherPhysical = fHisto
2100 G4LogicalVolume *motherLogical = motherPhys
2101
2102 #ifdef CHECK_ORDER_OF_METHODS
2103 if( ! fLastTriedStepComputation )
2104 {
2105 G4Exception("G4Navigator::RecheckDistanc
2106 "GeomNav0001", FatalExceptio
2107 "Method must be called after ComputeStep
2108 }
2109 #endif
2110
2111 EInside locatedDaug; // = kUndefined;
2112 G4double daughterStep= DBL_MAX;
2113 G4double daughterSafety= DBL_MAX;
2114
2115 if( fEnteredDaughter )
2116 {
2117 if( motherLogical->CharacteriseDaughters
2118
2119 // Track arrived at boundary of a daught
2120 // the last call of ComputeStep().
2121 // In case the proposed displaced point
2122 // it must backtrack at least to the e
2123 // NOTE: No check is made against other
2124 // assumed that the proposed displacem
2125 // this is not needed.
2126
2127 // Must check boundary of current daught
2128 //
2129 G4VPhysicalVolume *candPhysical= fBlocke
2130 G4LogicalVolume *candLogical= candPhysic
2131 G4VSolid *candSolid= candLogica
2132
2133 G4AffineTransform nextLevelTrf(candPhysi
2134 candPhysi
2135
2136 G4ThreeVector dgPosition= nextLevelTrf.
2137 G4ThreeVector dgDirection= nextLevelTrf.
2138 locatedDaug = candSolid->Inside(dgPositi
2139
2140 if( locatedDaug == kInside )
2141 {
2142 // Reverse direction - and find first
2143 // Must backtrack
2144 G4double distanceBackOut =
2145 candSolid->DistanceToOut(dgPosition
2146 - dgDirect
2147 true, &val
2148 daughterStep= - distanceBackOut;
2149 // No check is made whether the par
2150 // at this location without encount
2151 // a different psurface of the curr
2152 if( prNewSafety != nullptr )
2153 {
2154 daughterSafety= candSolid->Distanc
2155 }
2156 }
2157 else
2158 {
2159 if( locatedDaug == kOutside )
2160 {
2161 // See whether it still intersect
2162 //
2163 daughterStep= candSolid->Distanc
2164
2165 if( prNewSafety != nullptr )
2166 {
2167 daughterSafety= candSolid->Dist
2168 }
2169 }
2170 else
2171 {
2172 // The point remains on the surfac
2173 //
2174 daughterStep= 0.0;
2175 daughterSafety= 0.0;
2176 }
2177 }
2178
2179 // If trial point is in daughter (or on
2180 // answer, the rest is not relevant
2181 //
2182 if( locatedDaug != kOutside )
2183 {
2184 *prDistance= daughterStep;
2185 if( prNewSafety != nullptr ) { *prNe
2186 return true;
2187 }
2188 // If ever extended, so that some type o
2189 // this would change
2190 }
2191
2192 G4VSolid *motherSolid= motherLogical->GetSo
2193
2194 G4double motherStep= DBL_MAX, motherSafety=
2195
2196 // Check distance to boundary of mother
2197 //
2198 if ( fHistory.GetTopVolumeType()!=kReplica
2199 {
2200 EInside locatedMoth = motherSolid->Insid
2201
2202 if( locatedMoth == kInside )
2203 {
2204 motherSafety= motherSolid->DistanceTo
2205 if( ProposedMove >= motherSafety )
2206 {
2207 motherStep= motherSolid->DistanceT
2208
2209 true, &validExit
2210 }
2211 else
2212 {
2213 motherStep= ProposedMove;
2214 }
2215 }
2216 else if( locatedMoth == kOutside)
2217 {
2218 motherSafety= motherSolid->DistanceTo
2219 if( ProposedMove >= motherSafety )
2220 {
2221 motherStep= - motherSolid->Distan
2222
2223 }
2224 }
2225 else
2226 {
2227 motherSafety= 0.0;
2228 *prDistance= 0.0; // On surface - n
2229 if( prNewSafety != nullptr ) { *prNe
2230 return false;
2231 }
2232 }
2233 else
2234 {
2235 return false;
2236 }
2237
2238 *prDistance= std::min( motherStep, daughte
2239 if( prNewSafety != nullptr )
2240 {
2241 *prNewSafety= std::min( motherSafety, da
2242 }
2243 return true;
2244 }
2245
2246
2247 // ******************************************
2248 // CreateTouchableHistoryHandle
2249 // ******************************************
2250 //
2251 G4TouchableHandle G4ITNavigator2::CreateTouch
2252 {
2253 CheckNavigatorStateIsValid();
2254 return G4TouchableHandle( CreateTouchableHi
2255 }
2256
2257 // ******************************************
2258 // PrintState
2259 // ******************************************
2260 //
2261 void G4ITNavigator2::PrintState() const
2262 {
2263 CheckNavigatorStateIsValid();
2264 G4long oldcoutPrec = G4cout.precision(4);
2265 if( fVerbose >= 4 )
2266 {
2267 G4cout << "The current state of G4Navigat
2268 G4cout << " ValidExitNormal= " << fValid
2269 << " ExitNormal = " << fExitN
2270 << " Exiting = " << fExiti
2271 << " Entering = " << fEnter
2272 << " BlockedPhysicalVolume= " ;
2273 if (fBlockedPhysicalVolume==nullptr)
2274 G4cout << "None";
2275 else
2276 G4cout << fBlockedPhysicalVolume->GetNa
2277 G4cout << G4endl
2278 << " BlockedReplicaNo = " <<
2279 << " LastStepWasZero = " <<
2280 << G4endl;
2281 }
2282 if( ( 1 < fVerbose) && (fVerbose < 4) )
2283 {
2284 G4cout << G4endl; // Make sure to line up
2285 G4cout << std::setw(30) << " ExitNormal "
2286 << std::setw( 5) << " Valid "
2287 << std::setw( 9) << " Exiting "
2288 << std::setw( 9) << " Entering"
2289 << std::setw(15) << " Blocked:Volu
2290 << std::setw( 9) << " ReplicaNo"
2291 << std::setw( 8) << " LastStepZero
2292 << G4endl;
2293 G4cout << "( " << std::setw(7) << fExitNo
2294 << ", " << std::setw(7) << fExitNo
2295 << ", " << std::setw(7) << fExitNo
2296 << std::setw( 5) << fValidExitNor
2297 << std::setw( 9) << fExiting
2298 << std::setw( 9) << fEntering
2299 if ( fBlockedPhysicalVolume==nullptr )
2300 {
2301 G4cout << std::setw(15) << "None";
2302 }
2303 else
2304 {
2305 G4cout << std::setw(15)<< fBlockedPhysi
2306 }
2307 G4cout << std::setw( 9) << fBlockedRepli
2308 << std::setw( 8) << fLastStepWasZ
2309 << G4endl;
2310 }
2311 if( fVerbose > 2 )
2312 {
2313 G4cout.precision(8);
2314 G4cout << " Current Localpoint = " << fLa
2315 G4cout << " PreviousSftOrigin = " << fPr
2316 G4cout << " PreviousSafety = " << fPr
2317 }
2318 G4cout.precision(oldcoutPrec);
2319 }
2320
2321 // ******************************************
2322 // ComputeStepLog
2323 // ******************************************
2324 //
2325 void G4ITNavigator2::ComputeStepLog(const G4T
2326 G4doub
2327 {
2328 CheckNavigatorStateIsValid();
2329 // The following checks only make sense if
2330 // than the tolerance.
2331
2332 static const G4double fAccuracyForWarning
2333 fAccuracyForException
2334
2335 G4ThreeVector OriginalGlobalpoint = fHistor
2336 Transfo
2337
2338 G4double shiftOriginSafSq = (fPreviousSftOr
2339
2340 // Check that the starting point of this st
2341 // within the isotropic safety sphere of th
2342 // to a accuracy/precision given by fAccur
2343 // If so give warning.
2344 // If it fails by more than fAccuracyForE
2345 //
2346 if( shiftOriginSafSq >= sqr(fPreviousSafety
2347 {
2348 G4double shiftOrigin = std::sqrt(shiftOri
2349 G4double diffShiftSaf = shiftOrigin - fPr
2350
2351 if( diffShiftSaf > fAccuracyForWarning )
2352 {
2353 G4long oldcoutPrec= G4cout.precision(8)
2354 G4long oldcerrPrec= G4cerr.precision(10
2355 std::ostringstream message, suggestion;
2356 message << "Accuracy error or slightly
2357 << G4endl
2358 << " The Step's starting po
2359 << std::sqrt(moveLenSq)/mm << "
2360 << " since the last call to
2361 << " This has resulted in m
2362 << shiftOrigin/mm << " mm "
2363 << " from the last point at whi
2364 << " was calculated " << G4
2365 << " which is more than the
2366 << fPreviousSafety/mm << " mm
2367 << " This difference is "
2368 << diffShiftSaf/mm << " mm." <<
2369 << " The tolerated accuracy
2370 << fAccuracyForException/mm <<
2371
2372 suggestion << " ";
2373 static G4ThreadLocal G4int warnNow = 0;
2374 if( ((++warnNow % 100) == 1) )
2375 {
2376 message << G4endl
2377 << " This problem can be due
2378 << " - a process that has p
2379 << " larger than the current s
2380 << " - inaccuracy in the co
2381 suggestion << "We suggest that you "
2382 << " - find i) what part
2383 << " ii) through what part
2384 << " for example by r
2385 << G4endl
2386 << " /tracking/ver
2387 << " - check which proc
2388 << " this particle (and lo
2389 << G4endl
2390 << " - in case, create a
2391 << " of this event using:"
2392 << " /tracking/ver
2393 }
2394 G4Exception("G4ITNavigator2::ComputeSte
2395 "GeomNav1002", JustWarning,
2396 message, G4String(suggestio
2397 G4cout.precision(oldcoutPrec);
2398 G4cerr.precision(oldcerrPrec);
2399 }
2400 #ifdef G4DEBUG_NAVIGATION
2401 else
2402 {
2403 G4cerr << "WARNING - G4ITNavigator2::Co
2404 << " The Step's startin
2405 << std::sqrt(moveLenSq) << "," <
2406 << " which has taken it
2407 << " the current safety. " << G4
2408 }
2409 #endif
2410 }
2411 G4double safetyPlus = fPreviousSafety + fAc
2412 if ( shiftOriginSafSq > sqr(safetyPlus) )
2413 {
2414 std::ostringstream message;
2415 message << "May lead to a crash or unreli
2416 << " Position has shifted
2417 << " notifying the navigator !" <
2418 << " Tolerated safety: " <
2419 << " Computed shift : " <
2420 G4Exception("G4ITNavigator2::ComputeStep(
2421 JustWarning, message);
2422 }
2423 }
2424
2425 // ******************************************
2426 // Operator <<
2427 // ******************************************
2428 //
2429 std::ostream& operator << (std::ostream &os,c
2430 {
2431 // Old version did only the following:
2432 // os << "Current History: " << G4endl << n
2433 // Old behaviour is recovered for fVerbose
2434
2435 // Adapted from G4ITNavigator2::PrintState(
2436
2437 G4long oldcoutPrec = os.precision(4);
2438 if( n.fVerbose >= 4 )
2439 {
2440 os << "The current state of G4ITNavigator
2441 os << " ValidExitNormal= " << n.fValidEx
2442 << " ExitNormal = " << n.fExitNormal
2443 << " Exiting = " << n.fExiting
2444 << " Entering = " << n.fEntering
2445 << " BlockedPhysicalVolume= " ;
2446
2447 if (n.fBlockedPhysicalVolume==nullptr)
2448 {
2449 os << "None";
2450 }
2451 else
2452 {
2453 os << n.fBlockedPhysicalVolume->GetName
2454 }
2455
2456 os << G4endl
2457 << " BlockedReplicaNo = " << n.fBlo
2458 << " LastStepWasZero = " << n.fLa
2459 << G4endl;
2460 }
2461 if( ( 1 < n.fVerbose) && (n.fVerbose < 4) )
2462 {
2463 os << G4endl; // Make sure to line up
2464 os << std::setw(30) << " ExitNormal " <<
2465 << std::setw( 5) << " Valid " << "
2466 << std::setw( 9) << " Exiting " << "
2467 << std::setw( 9) << " Entering" << "
2468 << std::setw(15) << " Blocked:Volume " <
2469 << std::setw( 9) << " ReplicaNo" <
2470 << std::setw( 8) << " LastStepZero " <
2471 << G4endl;
2472 os << "( " << std::setw(7) << n.fExitNorm
2473 << ", " << std::setw(7) << n.fExitNormal.
2474 << ", " << std::setw(7) << n.fExitNormal.
2475 << std::setw( 5) << n.fValidExitNormal
2476 << std::setw( 9) << n.fExiting
2477 << std::setw( 9) << n.fEntering
2478
2479 if ( n.fBlockedPhysicalVolume==nullptr )
2480 { os << std::setw(15) << "None"; }
2481 else
2482 { os << std::setw(15)<< n.fBlockedPhysica
2483
2484 os << std::setw( 9) << n.fBlockedReplica
2485 << std::setw( 8) << n.fLastStepWasZero
2486 << G4endl;
2487 }
2488 if( n.fVerbose > 2 )
2489 {
2490 os.precision(8);
2491 os << " Current Localpoint = " << n.fLast
2492 os << " PreviousSftOrigin = " << n.fPrev
2493 os << " PreviousSafety = " << n.fPrev
2494 }
2495 if( n.fVerbose > 3 || n.fVerbose == 0 )
2496 {
2497 os << "Current History: " << G4endl << n.
2498 }
2499
2500 os.precision(oldcoutPrec);
2501 return os;
2502 }
2503
2504 //-------------------------------------------
2505
2506 EInside
2507 G4ITNavigator2::InsideCurrentVolume(const G4T
2508 {
2509 const G4AffineTransform& transform = GetGlo
2510 G4ThreeVector localPoint(transform.Transfor
2511
2512 G4VSolid* solid = fHistory.GetTopVolume()->
2513
2514 return solid->Inside(localPoint);
2515 }
2516
2517 //-------------------------------------------
2518
2519 void G4ITNavigator2::GetRandomInCurrentVolume
2520 {
2521 const G4AffineTransform& local2Global = Get
2522 G4VSolid* solid = fHistory.GetTopVolume()->
2523
2524 G4VoxelLimits voxelLimits; // Defaults to
2525 G4double vmin, vmax;
2526 G4AffineTransform dummy;
2527 std::vector<std::vector<G4double> > fExtend
2528
2529 solid->CalculateExtent(kXAxis,voxelLimits,d
2530 fExtend[kXAxis][BoundingBox::kMin] = vmin;
2531 fExtend[kXAxis][BoundingBox::kMax] = vmax;
2532
2533 solid->CalculateExtent(kYAxis,voxelLimits,d
2534 fExtend[kYAxis][BoundingBox::kMin] = vmin;
2535 fExtend[kYAxis][BoundingBox::kMax] = vmax;
2536
2537 solid->CalculateExtent(kZAxis,voxelLimits,d
2538 fExtend[kZAxis][BoundingBox::kMin] = vmin;
2539 fExtend[kZAxis][BoundingBox::kMax] = vmax;
2540
2541 G4ThreeVector rndmPos;
2542
2543 while(true)
2544 {
2545 for(G4int i = 0 ; i < 3 ; ++i)
2546 {
2547 G4double min = fExtend[i][BoundingBox::
2548 G4double max = fExtend[i][BoundingBox::
2549 rndmPos[i] = G4UniformRand()*(max-min)+
2550 }
2551
2552 if(solid->Inside(rndmPos) == kInside) bre
2553 }
2554
2555 _rndmPoint = local2Global.TransformPoint(rn
2556 }
2557
2558