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