Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/dna/management/src/G4ITNavigator1.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  1 //
  2 // ********************************************************************
  3 // * License and Disclaimer                                           *
  4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.                             *
 10 // *                                                                  *
 11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                                                  *
 18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // ********************************************************************
 25 //
 26 // 
 27 // class G4ITNavigator1 Implementation
 28 //
 29 // Original author: Paul Kent, July 95/96
 30 //
 31 // G4ITNavigator1 is a duplicate version of G4Navigator starting from Geant4.9.5
 32 // initially written by Paul Kent and colleagues.
 33 // The only difference resides in the way the information is saved and managed
 34 //
 35 // History:
 36 // - Created.                                  Paul Kent,     Jul 95/96
 37 // - Zero step protections                     J.A. / G.C.,   Nov  2004
 38 // - Added check mode                          G. Cosmo,      Mar  2004
 39 // - Made Navigator Abstract                   G. Cosmo,      Nov  2003
 40 // - G4ITNavigator1 created                     M.K.,          Nov  2012
 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::GetInstance()->GetSurfaceTolerance();
 75   fregularNav.SetNormalNavigation( &fnormalNav );
 76 
 77   fStepEndPoint = G4ThreeVector( kInfinity, kInfinity, kInfinity ); 
 78   fLastStepEndPointLocal = G4ThreeVector( kInfinity, kInfinity, kInfinity ); 
 79 
 80   fpVoxelSafety= new G4VoxelSafety();
 81   fpSaveState = nullptr;
 82 
 83   // this->SetVerboseLevel(3);
 84   // this->CheckMode(true);
 85 }
 86 
 87 // !>
 88 
 89 G4ITNavigator1::G4SaveNavigatorState::G4SaveNavigatorState()
 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,0);
102 
103     sPreviousSftOrigin     = G4ThreeVector(0,0,0);
104     sPreviousSafety        = 0.0;
105 
106     sNumberZeroSteps       = 0;
107 
108     spBlockedPhysicalVolume = nullptr;
109     sBlockedReplicaNo      = -1;
110 
111     sLastLocatedPointLocal = G4ThreeVector( kInfinity, -kInfinity, 0.0 );
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 G4ThreeVector &p,
130                                      const G4ThreeVector &direction,
131                                      const G4TouchableHistory &h)
132 {
133   ResetState();
134   fHistory = *h.GetHistory();
135   SetupHierarchy();
136   fLastTriedStepComputation= false;  // Redundant, but best
137   return LocateGlobalPointAndSetup(p, &direction, true, false);
138 }
139 
140 // ********************************************************************
141 // LocateGlobalPointAndSetup
142 //
143 // Locate the point in the hierarchy return 0 if outside
144 // The direction is required 
145 //    - if on an edge shared by more than two surfaces 
146 //      (to resolve likely looping in tracking)
147 //    - at initial location of a particle
148 //      (to resolve potential ambiguity at boundary)
149 // 
150 // Flags on exit: (comments to be completed)
151 // fEntering         - True if entering `daughter' volume (or replica)
152 //                     whether daughter of last mother directly 
153 //                     or daughter of that volume's ancestor.
154 // ********************************************************************
155 //
156 G4VPhysicalVolume* 
157 G4ITNavigator1::LocateGlobalPointAndSetup( const G4ThreeVector& globalPoint,
158                                         const G4ThreeVector* pGlobalDirection,
159                                         const G4bool relativeSearch,
160                                         const G4bool ignoreDirection )
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) || fLocatedOnEdge;
170 
171   fLastTriedStepComputation=   false;   
172   fChangedGrandMotherRefFrame= false;  // For local exit normal
173    
174   if( considerDirection && pGlobalDirection != nullptr )
175   {
176     globalDirection=*pGlobalDirection;
177   }
178 
179 
180 #ifdef G4VERBOSE
181   if( fVerbose > 2 )
182   {
183     G4long oldcoutPrec = G4cout.precision(8);
184     G4cout << "*** G4ITNavigator1::LocateGlobalPointAndSetup: ***" << G4endl;
185     G4cout << "    Called with arguments: " << G4endl
186            << "    Globalpoint = " << globalPoint << G4endl
187            << "    RelativeSearch = " << relativeSearch  << G4endl;
188     if( fVerbose == 4 )
189     {
190       G4cout << "    ----- Upon entering:" << G4endl;
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;   // Remember
207       fExitedMother = fExiting;       // Remember
208       if ( fExiting )
209       {
210         if ( fHistory.GetDepth() != 0u )
211         {
212           fBlockedPhysicalVolume = fHistory.GetTopVolume();
213           fBlockedReplicaNo = fHistory.GetTopReplicaNo();
214           fHistory.BackLevel();
215         }
216         else
217         {
218           fLastLocatedPointLocal = localPoint;
219           fLocatedOutsideWorld = true;
220           return nullptr;           // Have exited world volume
221         }
222         // A fix for the case where a volume is "entered" at an edge
223         // and a coincident surface exists outside it.
224         //  - This stops it from exiting further volumes and cycling
225         //  - However ReplicaNavigator treats this case itself
226         //
227         if ( fLocatedOnEdge && (VolumeType(fBlockedPhysicalVolume)!=kReplica ))
228         { 
229           fExiting= false;
230         }
231       }
232       else
233         if ( fEntering )
234         {
235           switch (VolumeType(fBlockedPhysicalVolume))
236           {
237             case kNormal:
238               fHistory.NewLevel(fBlockedPhysicalVolume, kNormal,
239                                 fBlockedPhysicalVolume->GetCopyNo());
240               break;
241             case kReplica:
242               freplicaNav.ComputeTransformation(fBlockedReplicaNo,
243                                                 fBlockedPhysicalVolume);
244               fHistory.NewLevel(fBlockedPhysicalVolume, kReplica,
245                                 fBlockedReplicaNo);
246               fBlockedPhysicalVolume->SetCopyNo(fBlockedReplicaNo);
247               break;
248             case kParameterised:
249               if( fBlockedPhysicalVolume->GetRegularStructureId() == 0 )
250               {
251                 G4VSolid *pSolid;
252                 G4VPVParameterisation *pParam;
253                 G4TouchableHistory parentTouchable( fHistory );
254                 pParam = fBlockedPhysicalVolume->GetParameterisation();
255                 pSolid = pParam->ComputeSolid(fBlockedReplicaNo,
256                                               fBlockedPhysicalVolume);
257                 pSolid->ComputeDimensions(pParam, fBlockedReplicaNo,
258                                           fBlockedPhysicalVolume);
259                 pParam->ComputeTransformation(fBlockedReplicaNo,
260                                               fBlockedPhysicalVolume);
261                 fHistory.NewLevel(fBlockedPhysicalVolume, kParameterised,
262                                   fBlockedReplicaNo);
263                 fBlockedPhysicalVolume->SetCopyNo(fBlockedReplicaNo);
264                 //
265                 // Set the correct solid and material in Logical Volume
266                 //
267                 G4LogicalVolume *pLogical;
268                 pLogical = fBlockedPhysicalVolume->GetLogicalVolume();
269                 pLogical->SetSolid( pSolid );
270                 pLogical->UpdateMaterial(pParam ->
271                   ComputeMaterial(fBlockedReplicaNo,
272                                   fBlockedPhysicalVolume, 
273                                   &parentTouchable));
274               }
275               break;
276            case kExternal:
277              G4Exception("G4ITNavigator1::LocateGlobalPointAndSetup()",
278                          "GeomNav0001", FatalException,
279                          "Not applicable for external volumes.");
280              break;
281           }
282           fEntering = false;
283           fBlockedPhysicalVolume = nullptr;
284           localPoint = fHistory.GetTopTransform().TransformPoint(globalPoint);
285           notKnownContained = false;
286         }
287     }
288     else
289     {
290       fBlockedPhysicalVolume = nullptr;
291       fEntering = false;
292       fEnteredDaughter = false;  // Full Step was not taken, did not enter
293       fExiting = false;
294       fExitedMother = false;     // Full Step was not taken, did not exit
295     }
296   }
297   //
298   // Search from top of history up through geometry until
299   // containing volume found:
300   // If on 
301   // o OUTSIDE - Back up level, not/no longer exiting volumes
302   // o SURFACE and EXITING - Back up level, setting new blocking no.s
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()->GetLogicalVolume()->GetSolid();
313       localPoint = fHistory.GetTopTransform().TransformPoint(globalPoint);
314       insideCode = targetSolid->Inside(localPoint);
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::LocateGlobalPointAndSetup(): ***" << G4endl
324                 << "    Invoked Inside() for solid: " << targetSolid->GetName()
325                 << ". Solid replied: " << solidResponse << G4endl
326                 << "    For local point p: " << localPoint << G4endl;
327       }
328 #endif
329     }
330     else
331     {
332       insideCode = freplicaNav.BackLocate(fHistory, globalPoint, localPoint,
333                                           fExiting, notKnownContained);
334       // !CARE! if notKnownContained returns false then the point is within
335       // the containing placement volume of the replica(s). If insidecode
336       // will result in the history being backed up one level, then the
337       // local point returned is the point in the system of this new level
338     }
339 
340 
341     if ( insideCode==kOutside )
342     {
343       noLevelsExited++; 
344       if ( fHistory.GetDepth() != 0u )
345       {
346         fBlockedPhysicalVolume = fHistory.GetTopVolume();
347         fBlockedReplicaNo = fHistory.GetTopReplicaNo();
348         fHistory.BackLevel();
349         fExiting = false;
350 
351         if( noLevelsExited > 1 )
352         {
353           // The first transformation was done by the sub-navigator
354           //
355           const G4RotationMatrix* mRot = fBlockedPhysicalVolume->GetRotation();
356           if( mRot != nullptr )
357           { 
358             fGrandMotherExitNormal *= (*mRot).inverse();
359             fChangedGrandMotherRefFrame= true;
360           }
361         }
362       }
363       else
364       {
365         fLastLocatedPointLocal = localPoint;
366         fLocatedOutsideWorld = true;
367           // No extra transformation for ExitNormal - is in frame of Top Volume
368         return nullptr;         // Have exited world volume
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 this level's volume
378           // by using the direction
379           //
380           G4bool directionExiting = false;
381           G4ThreeVector localDirection =
382               fHistory.GetTopTransform().TransformAxis(globalDirection);
383 
384           // Make sure localPoint in correct reference frame
385           //     ( Was it already correct ? How ? )
386           //
387           localPoint= fHistory.GetTopTransform().TransformPoint(globalPoint);
388           if ( fHistory.GetTopVolumeType()!=kReplica )
389           {
390             G4ThreeVector normal = targetSolid->SurfaceNormal(localPoint);
391             directionExiting = normal.dot(localDirection) > 0.0;
392             isExiting = isExiting || directionExiting;
393           }
394         }
395         if( isExiting )
396         {
397           noLevelsExited++; 
398           if ( fHistory.GetDepth() != 0u )
399           {
400             fBlockedPhysicalVolume = fHistory.GetTopVolume();
401             fBlockedReplicaNo = fHistory.GetTopReplicaNo();
402             fHistory.BackLevel();
403             //
404             // Still on surface but exited volume not necessarily convex
405             //
406             fValidExitNormal = false;
407 
408             if( noLevelsExited > 1 )
409             {
410               // The first transformation was done by the sub-navigator
411               //
412               const G4RotationMatrix* mRot=fBlockedPhysicalVolume->GetRotation();
413               if( mRot != nullptr )
414               { 
415                 fGrandMotherExitNormal *= (*mRot).inverse();
416                 fChangedGrandMotherRefFrame= true;
417               }
418             }
419           } 
420           else
421           {
422             fLastLocatedPointLocal = localPoint;
423             fLocatedOutsideWorld = true;
424               // No extra transformation for ExitNormal, is in frame of Top Vol
425             return nullptr;          // Have exited world volume
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 volume found,
440   // blocking fBlockedPhysicalVolume/BlockedReplicaNum
441   //
442   // 3 Cases:
443   //
444   // o Parameterised daughters
445   //   =>Must be one G4PVParameterised daughter & voxels
446   // o Positioned daughters & voxels
447   // o Positioned daughters & no voxels
448 
449   noResult = true;  // noResult should be renamed to 
450                     // something like enteredLevel, as that is its meaning.
451   do
452   {
453     // Determine `type' of current mother volume
454     //
455     targetPhysical = fHistory.GetTopVolume();
456     if (targetPhysical == nullptr) { break; }
457     targetLogical = targetPhysical->GetLogicalVolume();
458     switch( CharacteriseDaughters(targetLogical) )
459     {
460       case kNormal:
461         if ( targetLogical->GetVoxelHeader() != nullptr )  // use optimised navigation
462         {
463           noResult = fvoxelNav.LevelLocate(fHistory,
464                                            fBlockedPhysicalVolume,
465                                            fBlockedReplicaNo,
466                                            globalPoint,
467                                            pGlobalDirection,
468                                            considerDirection,
469                                            localPoint);
470         }
471         else                       // do not use optimised navigation
472         {
473           noResult = fnormalNav.LevelLocate(fHistory,
474                                             fBlockedPhysicalVolume,
475                                             fBlockedReplicaNo,
476                                             globalPoint,
477                                             pGlobalDirection,
478                                             considerDirection,
479                                             localPoint);
480         }
481         break;
482       case kReplica:
483         noResult = freplicaNav.LevelLocate(fHistory,
484                                            fBlockedPhysicalVolume,
485                                            fBlockedReplicaNo,
486                                            globalPoint,
487                                            pGlobalDirection,
488                                            considerDirection,
489                                            localPoint);
490         break;
491       case kParameterised:
492         if( GetDaughtersRegularStructureId(targetLogical) != 1 )
493         {
494           noResult = fparamNav.LevelLocate(fHistory,
495                                            fBlockedPhysicalVolume,
496                                            fBlockedReplicaNo,
497                                            globalPoint,
498                                            pGlobalDirection,
499                                            considerDirection,
500                                            localPoint);
501         }
502         else  // Regular structure
503         {
504           noResult = fregularNav.LevelLocate(fHistory,
505                                              fBlockedPhysicalVolume,
506                                              fBlockedReplicaNo,
507                                              globalPoint,
508                                              pGlobalDirection,
509                                              considerDirection,
510                                              localPoint);
511         }
512         break;
513       case kExternal:
514         G4Exception("G4ITNavigator1::LocateGlobalPointAndSetup()",
515                     "GeomNav0001", FatalException,
516                     "Not applicable for external volumes.");
517         break;
518     }
519 
520     // LevelLocate returns true if it finds a daughter volume 
521     // in which globalPoint is inside (or on the surface).
522 
523     if ( noResult )
524     {
525       // Entering a daughter after ascending
526       //
527       // The blocked volume is no longer valid - it was for another level
528       //
529       fBlockedPhysicalVolume = nullptr;
530       fBlockedReplicaNo = -1;
531 
532       // fEntering should be false -- else blockedVolume is assumed good.
533       // fEnteredDaughter is used for ExitNormal
534       //
535       fEntering = false;
536       fEnteredDaughter = true;
537 
538       if( fExitedMother )
539       {
540         G4VPhysicalVolume* enteredPhysical = fHistory.GetTopVolume();
541         const G4RotationMatrix* mRot = enteredPhysical->GetRotation();
542         if( mRot != nullptr )
543         { 
544           fGrandMotherExitNormal *= (*mRot).inverse();
545         }
546       }
547 
548 #ifdef G4DEBUG_NAVIGATION
549       if( fVerbose > 2 )
550       { 
551          G4VPhysicalVolume* enteredPhysical = fHistory.GetTopVolume();
552          G4cout << "*** G4ITNavigator1::LocateGlobalPointAndSetup() ***" << G4endl;
553          G4cout << "    Entering volume: " << enteredPhysical->GetName()
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)  { curPhysVol_Name = targetPhysical->GetName(); }
568     G4cout << "    Return value = new volume = " << curPhysVol_Name << G4endl;
569     G4cout << "    ----- Upon exiting:" << G4endl;
570     PrintState();
571     if( fVerbose == 5 )
572     {
573       G4cout << "Upon exiting LocateGlobalPointAndSetup():" << G4endl;
574       G4cout << "    History = " << G4endl << fHistory << G4endl << 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 and its subNavigators
589 //    is updated in order to start the next step at pGlobalpoint
590 // -> no check is performed whether pGlobalpoint is inside the 
591 //    original volume (this must be the case).
592 //
593 // Note: a direction could be added to the arguments, to aid in future
594 //       optional checking (via the old code below, flagged by OLD_LOCATE). 
595 //       [ This would be done only in verbose mode ]
596 // ********************************************************************
597 //
598 void
599 G4ITNavigator1::LocateGlobalPointWithinVolume(const G4ThreeVector& pGlobalpoint)
600 {  
601    fLastLocatedPointLocal = ComputeLocalPoint(pGlobalpoint);
602    fLastTriedStepComputation= false;
603    fChangedGrandMotherRefFrame= false;  //  Frame for Exit Normal
604 
605 #ifdef G4DEBUG_NAVIGATION
606    if( fVerbose > 2 )
607    { 
608      G4cout << "Entering LocateGlobalWithinVolume(): History = " << G4endl;
609      G4cout << fHistory << G4endl;
610    }
611 #endif
612 
613    // For the case of Voxel (or Parameterised) volume the respective 
614    // Navigator must be messaged to update its voxel information etc
615 
616    // Update the state of the Sub Navigators 
617    // - in particular any voxel information they store/cache
618    //
619    G4VPhysicalVolume*  motherPhysical = fHistory.GetTopVolume();
620    G4LogicalVolume*    motherLogical  = motherPhysical->GetLogicalVolume();
621    G4SmartVoxelHeader* pVoxelHeader   = motherLogical->GetVoxelHeader();
622 
623    if ( fHistory.GetTopVolumeType()!=kReplica )
624    {
625      switch( CharacteriseDaughters(motherLogical) )
626      {
627        case kNormal:
628          if ( pVoxelHeader != nullptr )
629          {
630            fvoxelNav.VoxelLocate( pVoxelHeader, fLastLocatedPointLocal );
631          }
632          break;
633        case kParameterised:
634          if( GetDaughtersRegularStructureId(motherLogical) != 1 )
635          {
636            // Resets state & returns voxel node
637            //
638            fparamNav.ParamVoxelLocate( pVoxelHeader, fLastLocatedPointLocal );
639          }
640          break;
641        case kReplica:
642          G4Exception("G4ITNavigator1::LocateGlobalPointWithinVolume()",
643                      "GeomNav0001", FatalException,
644                      "Not applicable for replicated volumes.");
645          break;
646        case kExternal:
647          G4Exception("G4ITNavigator1::LocateGlobalPointWithinVolume()",
648                      "GeomNav0001", FatalException,
649                      "Not applicable for external volumes.");
650          break;
651      }
652    }
653 
654    // Reset the state variables 
655    //   - which would have been affected
656    //     by the 'equivalent' call to LocateGlobalPointAndSetup
657    //   - who's values have been invalidated by the 'move'.
658    //
659    fBlockedPhysicalVolume = nullptr; 
660    fBlockedReplicaNo = -1;
661    fEntering = false;
662    fEnteredDaughter = false;  // Boundary not encountered, did not enter
663    fExiting = false;
664    fExitedMother = false;     // Boundary not encountered, did not exit
665 }
666 
667 // !>
668 G4ITNavigatorState_Lock1* G4ITNavigator1::GetNavigatorState()
669 {
670     SetSavedState();
671     return fpSaveState;
672 }
673 
674 void G4ITNavigator1::SetNavigatorState(G4ITNavigatorState_Lock1* navState)
675 {
676     fpSaveState = (G4SaveNavigatorState*) navState;
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 call
691 // Save fValidExitNormal, fExitNormal, fExiting, fEntering, 
692 //      fBlockedPhysicalVolume, fBlockedReplicaNo, fLastStepWasZero; 
693 // ********************************************************************
694 //
695 void G4ITNavigator1::SetSavedState()
696 {
697     // !>
698     // This check can be avoid if instead, at every first step of a track,
699     // the IT tracking uses NewNavigatorSate
700     // The normal tracking would just call once NewNavigatorState() before tracking
701 
702 //    if(fpSaveState == 0)
703 //        fpSaveState = new G4SaveNavigatorState;
704     // <!
705 
706   // fSaveExitNormal = fExitNormal;
707   fpSaveState->sExitNormal = fExitNormal;
708   fpSaveState->sValidExitNormal = fValidExitNormal;
709   fpSaveState->sExiting = fExiting;
710   fpSaveState->sEntering = fEntering;
711 
712   fpSaveState->spBlockedPhysicalVolume = fBlockedPhysicalVolume;
713   fpSaveState->sBlockedReplicaNo = fBlockedReplicaNo,
714 
715   fpSaveState->sLastStepWasZero = static_cast<G4int>(fLastStepWasZero);
716 
717   // !>
718   fpSaveState->sPreviousSftOrigin = fPreviousSftOrigin;
719   fpSaveState->sPreviousSafety = fPreviousSafety;
720   fpSaveState->sNumberZeroSteps = fNumberZeroSteps;
721   fpSaveState->sLocatedOnEdge = fLocatedOnEdge;
722   fpSaveState->sWasLimitedByGeometry= fWasLimitedByGeometry;
723   fpSaveState->sPushed=fPushed;
724   fpSaveState->sNumberZeroSteps=fNumberZeroSteps;
725   fpSaveState->sEnteredDaughter = fEnteredDaughter;
726   fpSaveState->sExitedMother = fExitedMother;
727 
728   fpSaveState->sLastLocatedPointLocal = fLastLocatedPointLocal;
729   fpSaveState->sLocatedOutsideWorld = fLocatedOutsideWorld;
730   // <!
731 }
732 
733 // ********************************************************************
734 // RestoreSavedState
735 //
736 // Restore the state (in Compute Step), in case this is a parasitic call
737 // ********************************************************************
738 //
739 void G4ITNavigator1::RestoreSavedState()
740 {
741   fExitNormal = fpSaveState->sExitNormal;
742   fValidExitNormal = fpSaveState->sValidExitNormal;
743   fExiting = fpSaveState->sExiting;
744   fEntering = fpSaveState->sEntering;
745 
746   fBlockedPhysicalVolume = fpSaveState->spBlockedPhysicalVolume;
747   fBlockedReplicaNo = fpSaveState->sBlockedReplicaNo,
748 
749   fLastStepWasZero = (fpSaveState->sLastStepWasZero != 0);
750 
751   // !>
752   fPreviousSftOrigin = fpSaveState->sPreviousSftOrigin ;
753   fPreviousSafety = fpSaveState->sPreviousSafety ;
754   fNumberZeroSteps = fpSaveState->sNumberZeroSteps ;
755   fLocatedOnEdge = fpSaveState->sLocatedOnEdge ;
756   fWasLimitedByGeometry = fpSaveState->sWasLimitedByGeometry;
757   fPushed = fpSaveState->sPushed;
758   fNumberZeroSteps = fpSaveState->sNumberZeroSteps;
759   fEnteredDaughter= fpSaveState->sEnteredDaughter ;
760   fExitedMother = fpSaveState->sExitedMother ;
761 
762   fLastLocatedPointLocal = fpSaveState->sLastLocatedPointLocal ;
763   fLocatedOutsideWorld = fpSaveState->sLocatedOutsideWorld;
764   // <!
765 }
766 // <!
767 
768 // ********************************************************************
769 // ComputeStep
770 //
771 // Computes the next geometric Step: intersections with current
772 // mother and `daughter' volumes.
773 //
774 // NOTE:
775 //
776 // Flags on entry:
777 // --------------
778 // fValidExitNormal  - Normal of exited volume is valid (convex, not a 
779 //                     coincident boundary)
780 // fExitNormal       - Surface normal of exited volume
781 // fExiting          - True if have exited solid
782 //
783 // fBlockedPhysicalVolume - Ptr to exited volume (or 0)
784 // fBlockedReplicaNo - Replication no of exited volume
785 // fLastStepWasZero  - True if last Step size was zero.
786 //
787 // Flags on exit:
788 // -------------
789 // fValidExitNormal  - True if surface normal of exited volume is valid
790 // fExitNormal       - Surface normal of exited volume rotated to mothers
791 //                    reference system
792 // fExiting          - True if exiting mother
793 // fEntering         - True if entering `daughter' volume (or replica)
794 // fBlockedPhysicalVolume - Ptr to candidate (entered) volume
795 // fBlockedReplicaNo - Replication no of candidate (entered) volume
796 // fLastStepWasZero  - True if this Step size was zero.
797 // ********************************************************************
798 //
799 G4double G4ITNavigator1::ComputeStep( const G4ThreeVector &pGlobalpoint,
800                                    const G4ThreeVector &pDirection,
801                                    const G4double pCurrentProposedStepLength,
802                                          G4double &pNewSafety)
803 {
804   G4ThreeVector localDirection = ComputeLocalAxis(pDirection);
805   G4double Step = kInfinity;
806   G4VPhysicalVolume  *motherPhysical = fHistory.GetTopVolume();
807   G4LogicalVolume *motherLogical = motherPhysical->GetLogicalVolume();
808 
809   // All state relating to exiting normals must be reset
810   //
811   fExitNormalGlobalFrame= G4ThreeVector( 0., 0., 0.);
812     // Reset value - to erase its memory
813   fChangedGrandMotherRefFrame= false;
814     // Reset - used for local exit normal
815   fGrandMotherExitNormal= G4ThreeVector( 0., 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: ***" << G4endl;
828     G4cout << "    Volume = " << motherPhysical->GetName() 
829            << " - Proposed step length = " << pCurrentProposedStepLength
830            << G4endl; 
831 #ifdef G4DEBUG_NAVIGATION
832     if( fVerbose >= 2 )
833     {
834       G4cout << "  Called with the arguments: " << G4endl
835              << "  Globalpoint = " << std::setw(25) << pGlobalpoint << G4endl
836              << "  Direction   = " << std::setw(25) << pDirection << G4endl;
837       if( fVerbose >= 4 )
838       {
839         G4cout << "  ---- Upon entering : State" << G4endl;
840         PrintState();
841       }
842     }
843 #endif
844   }
845 #endif
846 
847   G4ThreeVector newLocalPoint = ComputeLocalPoint(pGlobalpoint);
848   if( newLocalPoint != fLastLocatedPointLocal )
849   {
850     // Check whether the relocation is within safety
851     //
852     G4ThreeVector oldLocalPoint = fLastLocatedPointLocal;
853     G4double moveLenSq = (newLocalPoint-oldLocalPoint).mag2();
854 
855     if ( moveLenSq >= kCarTolerance*kCarTolerance )
856     {
857 #ifdef G4VERBOSE
858       ComputeStepLog(pGlobalpoint, moveLenSq);
859 #endif
860       // Relocate the point within the same volume
861       //
862       LocateGlobalPointWithinVolume( pGlobalpoint );
863       fLastTriedStepComputation= true;     // Ensure that this is set again !!
864     }
865   }
866   if ( fHistory.GetTopVolumeType()!=kReplica )
867   {
868     switch( CharacteriseDaughters(motherLogical) )
869     {
870       case kNormal:
871         if ( motherLogical->GetVoxelHeader() != nullptr )
872         {
873           Step = fvoxelNav.ComputeStep(fLastLocatedPointLocal,
874                                        localDirection,
875                                        pCurrentProposedStepLength,
876                                        pNewSafety,
877                                        fHistory,
878                                        fValidExitNormal,
879                                        fExitNormal,
880                                        fExiting,
881                                        fEntering,
882                                        &fBlockedPhysicalVolume,
883                                        fBlockedReplicaNo);
884       
885         }
886         else
887         {
888           if( motherPhysical->GetRegularStructureId() == 0 )
889           {
890             Step = fnormalNav.ComputeStep(fLastLocatedPointLocal,
891                                           localDirection,
892                                           pCurrentProposedStepLength,
893                                           pNewSafety,
894                                           fHistory,
895                                           fValidExitNormal,
896                                           fExitNormal,
897                                           fExiting,
898                                           fEntering,
899                                           &fBlockedPhysicalVolume,
900                                           fBlockedReplicaNo);
901           }
902           else  // Regular (non-voxelised) structure
903           {
904             LocateGlobalPointAndSetup( pGlobalpoint, &pDirection, true, true );
905             fLastTriedStepComputation= true; // Ensure that this is set again !!
906             //
907             // if physical process limits the step, the voxel will not be the
908             // one given by ComputeStepSkippingEqualMaterials() and the local
909             // point will be wrongly calculated.
910 
911             // There is a problem: when msc limits the step and the point is
912             // assigned wrongly to phantom in previous step (while it is out
913             // of the container volume). Then LocateGlobalPointAndSetup() has
914             // reset the history topvolume to world.
915             //
916             if(fHistory.GetTopVolume()->GetRegularStructureId() == 0 )
917             { 
918               G4Exception("G4ITNavigator1::ComputeStep()",
919                           "GeomNav1001", JustWarning,
920                 "Point is relocated in voxels, while it should be outside!");
921               Step = fnormalNav.ComputeStep(fLastLocatedPointLocal,
922                                             localDirection,
923                                             pCurrentProposedStepLength,
924                                             pNewSafety,
925                                             fHistory,
926                                             fValidExitNormal,
927                                             fExitNormal,
928                                             fExiting,
929                                             fEntering,
930                                             &fBlockedPhysicalVolume,
931                                             fBlockedReplicaNo);
932             }
933             else
934             {
935               Step = fregularNav.
936                    ComputeStepSkippingEqualMaterials(fLastLocatedPointLocal,
937                                                      localDirection,
938                                                      pCurrentProposedStepLength,
939                                                      pNewSafety,
940                                                      fHistory,
941                                                      fValidExitNormal,
942                                                      fExitNormal,
943                                                      fExiting,
944                                                      fEntering,
945                                                      &fBlockedPhysicalVolume,
946                                                      fBlockedReplicaNo,
947                                                      motherPhysical);
948             }
949           }
950         }
951         break;
952       case kParameterised:
953         if( GetDaughtersRegularStructureId(motherLogical) != 1 )
954         {
955           Step = fparamNav.ComputeStep(fLastLocatedPointLocal,
956                                        localDirection,
957                                        pCurrentProposedStepLength,
958                                        pNewSafety,
959                                        fHistory,
960                                        fValidExitNormal,
961                                        fExitNormal,
962                                        fExiting,
963                                        fEntering,
964                                        &fBlockedPhysicalVolume,
965                                        fBlockedReplicaNo);
966         }
967         else  // Regular structure
968         {
969           Step = fregularNav.ComputeStep(fLastLocatedPointLocal,
970                                          localDirection,
971                                          pCurrentProposedStepLength,
972                                          pNewSafety,
973                                          fHistory,
974                                          fValidExitNormal,
975                                          fExitNormal,
976                                          fExiting,
977                                          fEntering,
978                                          &fBlockedPhysicalVolume,
979                                          fBlockedReplicaNo);
980         }
981         break;
982       case kReplica:
983         G4Exception("G4ITNavigator1::ComputeStep()", "GeomNav0001",
984                     FatalException, "Not applicable for replicated volumes.");
985         break;
986       case kExternal:
987         G4Exception("G4ITNavigator1::ComputeStep()", "GeomNav0001",
988                     FatalException, "Not applicable for external volumes.");
989         break;
990     }
991   }
992   else
993   {
994     // In the case of a replica, it must handle the exiting
995     // edge/corner problem by itself
996     //
997     G4bool exitingReplica = fExitedMother;
998     G4bool calculatedExitNormal= false;
999     
1000     Step = freplicaNav.ComputeStep(pGlobalpoint,
1001                                    pDirection,
1002                                    fLastLocatedPointLocal,
1003                                    localDirection,
1004                                    pCurrentProposedStepLength,
1005                                    pNewSafety,
1006                                    fHistory,
1007                                    fValidExitNormal,
1008                                    calculatedExitNormal,
1009                                    fExitNormal,
1010                                    exitingReplica,
1011                                    fEntering,
1012                                    &fBlockedPhysicalVolume,
1013                                    fBlockedReplicaNo);
1014     fExiting= exitingReplica;                          
1015     fCalculatedExitNormal= calculatedExitNormal;
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 changing momentum at a boundary
1024   //                  - one, two (or a few) can occur at common edges between
1025   //                    volumes
1026   //                  - more than two is likely a problem in the geometry
1027   //                    description or the Navigation 
1028 
1029   // Rule of thumb: likely at an Edge if two consecutive steps are zero,
1030   //                because at least two candidate volumes must have been
1031   //                checked
1032   //
1033   fLocatedOnEdge   = fLastStepWasZero && (Step==0.0);
1034   fLastStepWasZero = (Step==0.0);
1035   if (fPushed)  { fPushed = fLastStepWasZero; }
1036 
1037   // Handle large number of consecutive zero steps
1038   //
1039   if ( fLastStepWasZero )
1040   {
1041     fNumberZeroSteps++;
1042 #ifdef G4DEBUG_NAVIGATION
1043     if( fNumberZeroSteps > 1 )
1044     {
1045        G4cout << "G4ITNavigator1::ComputeStep(): another zero step, # "
1046               << fNumberZeroSteps
1047               << " at " << pGlobalpoint
1048               << " in volume " << motherPhysical->GetName()
1049               << " nav-comp-step calls # " << sNavCScalls
1050               << G4endl;
1051     }
1052 #endif
1053     if( fNumberZeroSteps > fActionThreshold_NoZeroSteps-1 )
1054     {
1055        // Act to recover this stuck track. Pushing it along direction
1056        //
1057        Step += 100*kCarTolerance;
1058 #ifdef G4VERBOSE
1059        if ((!fPushed) && (fWarnPush))
1060        {
1061          std::ostringstream message;
1062          message << "Track stuck or not moving." << G4endl
1063                  << "          Track stuck, not moving for " 
1064                  << fNumberZeroSteps << " steps" << G4endl
1065                  << "          in volume -" << motherPhysical->GetName()
1066                  << "- at point " << pGlobalpoint << G4endl
1067                  << "          direction: " << pDirection << "." << G4endl
1068                  << "          Potential geometry or navigation problem !"
1069                  << G4endl
1070                  << "          Trying pushing it of " << Step << " mm ...";
1071          G4Exception("G4ITNavigator1::ComputeStep()", "GeomNav1002",
1072                      JustWarning, message, "Potential overlap in geometry!");
1073        }
1074 #endif
1075        fPushed = true;
1076     }
1077     if( fNumberZeroSteps > fAbandonThreshold_NoZeroSteps-1 )
1078     {
1079       // Must kill this stuck track
1080       //
1081       std::ostringstream message;
1082       message << "Stuck Track: potential geometry or navigation problem."
1083               << G4endl
1084               << "        Track stuck, not moving for " 
1085               << fNumberZeroSteps << " steps" << G4endl
1086               << "        in volume -" << motherPhysical->GetName()
1087               << "- at point " << pGlobalpoint << G4endl
1088               << "        direction: " << pDirection << ".";
1089       motherPhysical->CheckOverlaps(5000, 0.0);
1090       G4Exception("G4ITNavigator1::ComputeStep()", "GeomNav0003",
1091                   EventMustBeAborted, message);
1092     }
1093   }
1094   else
1095   {
1096     if (!fPushed)  fNumberZeroSteps = 0;
1097   }
1098 
1099   fEnteredDaughter = fEntering;   // I expect to enter a volume in this Step
1100   fExitedMother = fExiting;
1101 
1102   fStepEndPoint = pGlobalpoint + Step * pDirection; 
1103   fLastStepEndPointLocal = fLastLocatedPointLocal + Step * localDirection; 
1104 
1105   if( fExiting )
1106   {
1107 #ifdef G4DEBUG_NAVIGATION
1108     if( fVerbose > 2 )
1109     { 
1110       G4cout << " At G4Nav CompStep End - if(exiting) - fExiting= " << fExiting 
1111              << " fValidExitNormal = " << fValidExitNormal  << G4endl;
1112       G4cout << " fExitNormal= " << fExitNormal << G4endl;
1113     }
1114 #endif
1115 
1116     if(fValidExitNormal || fCalculatedExitNormal)
1117     {
1118       if (  fHistory.GetTopVolumeType()!=kReplica )
1119       {
1120         // Convention: fExitNormal is in the 'grand-mother' coordinate system
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 (in order to have it if requested)
1133       //
1134       G4ThreeVector finalLocalPoint =
1135             fLastLocatedPointLocal + localDirection*Step;
1136 
1137       if (  fHistory.GetTopVolumeType()!=kReplica )
1138       {
1139         // Find normal in the 'mother' coordinate system
1140         //
1141         G4ThreeVector exitNormalMotherFrame=
1142            motherLogical->GetSolid()->SurfaceNormal(finalLocalPoint);
1143         
1144         // Transform it to the 'grand-mother' coordinate system
1145         //
1146         const G4RotationMatrix* mRot = motherPhysical->GetRotation();
1147         if( mRot != nullptr )
1148         {
1149           fChangedGrandMotherRefFrame= true;           
1150           fGrandMotherExitNormal = (*mRot).inverse() * exitNormalMotherFrame;
1151         }
1152         else
1153         {
1154           fGrandMotherExitNormal = exitNormalMotherFrame;
1155         }
1156 
1157         // Do not set fValidExitNormal -- this signifies
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 currently - to solve this
1167         // Replica Navigation must have calculated the normal for this case
1168         // already.
1169         // Cases: mother is not convex, and exit is at previous replica level
1170 
1171 #ifdef G4DEBUG_NAVIGATION
1172         G4ExceptionDescription desc;
1173 
1174         desc << "Problem in ComputeStep:  Replica Navigation did not provide"
1175              << " valid exit Normal. " << G4endl;
1176         desc << " Do not know how calculate it in this case." << G4endl;
1177         desc << "  Location    = " << finalLocalPoint << G4endl;
1178         desc << "  Volume name = " << motherPhysical->GetName()
1179              << "  copy/replica No = " << motherPhysical->GetCopyNo() << G4endl;
1180         G4Exception("G4ITNavigator1::ComputeStep()", "GeomNav0003",
1181                     JustWarning, desc, "Normal not available for exiting.");
1182 #endif
1183       }
1184     }
1185 
1186     // Now transform it to the global reference frame !!
1187     //
1188     if( fValidExitNormal || fCalculatedExitNormal )
1189     {
1190       auto  depth = (G4int)fHistory.GetDepth();
1191       if( depth > 0 )
1192       {
1193         G4AffineTransform GrandMotherToGlobalTransf =
1194           fHistory.GetTransform(depth-1).Inverse();
1195         fExitNormalGlobalFrame =
1196           GrandMotherToGlobalTransf.TransformAxis( fGrandMotherExitNormal );
1197       }
1198       else
1199       {
1200         fExitNormalGlobalFrame= fGrandMotherExitNormal;
1201       }
1202     }
1203     else
1204     {
1205       fExitNormalGlobalFrame= G4ThreeVector( 0., 0., 0.);
1206     }
1207   }
1208   fStepEndPoint= pGlobalpoint+Step*pDirection; 
1209 
1210   if( (Step == pCurrentProposedStepLength) && (!fExiting) && (!fEntering) )
1211   {
1212     // This if Step is not really limited by the geometry.
1213     // The Navigator is obliged to return "infinity"
1214     //
1215     Step = kInfinity;
1216   }
1217 
1218 #ifdef G4VERBOSE
1219   if( fVerbose > 1 )
1220   {
1221     if( fVerbose >= 4 )
1222     {
1223       G4cout << "    ----- Upon exiting :" << G4endl;
1224       PrintState();
1225     }
1226     G4cout << "  Returned step= " << Step;
1227     if( fVerbose > 5 )   G4cout << G4endl;
1228     if( Step == kInfinity )
1229     {
1230        G4cout << " Requested step= " << pCurrentProposedStepLength ;
1231        if( fVerbose > 5) G4cout << G4endl;
1232     }
1233     G4cout << "  Safety = " << pNewSafety << G4endl;
1234   }
1235 #endif
1236 
1237   return Step;
1238 }
1239 
1240 // ********************************************************************
1241 // CheckNextStep
1242 //
1243 // Compute the step without altering the navigator state
1244 // ********************************************************************
1245 //
1246 G4double G4ITNavigator1::CheckNextStep( const G4ThreeVector& pGlobalpoint,
1247                                      const G4ThreeVector& pDirection,
1248                                      const G4double pCurrentProposedStepLength,
1249                                            G4double& pNewSafety)
1250 {
1251   G4double step;
1252 
1253   // Save the state, for this parasitic call
1254   //
1255   SetSavedState();
1256 
1257   step = ComputeStep ( pGlobalpoint, 
1258                        pDirection,
1259                        pCurrentProposedStepLength, 
1260                        pNewSafety ); 
1261 
1262   // If a parasitic call, then attempt to restore the key parts of the state
1263   //
1264   RestoreSavedState(); 
1265 
1266   return step; 
1267 }
1268 
1269 // ********************************************************************
1270 // ResetState
1271 //
1272 // Resets stack and minimum of navigator state `machine'
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,0);
1291   fGrandMotherExitNormal = G4ThreeVector(0,0,0);
1292   fExitNormalGlobalFrame = G4ThreeVector(0,0,0);
1293 
1294   fPreviousSftOrigin     = G4ThreeVector(0,0,0);
1295   fPreviousSafety        = 0.0; 
1296 
1297   fNumberZeroSteps       = 0;
1298     
1299   fBlockedPhysicalVolume = nullptr;
1300   fBlockedReplicaNo      = -1;
1301 
1302   fLastLocatedPointLocal = G4ThreeVector( kInfinity, -kInfinity, 0.0 ); 
1303   fLocatedOutsideWorld   = false;
1304 }
1305 
1306 // ********************************************************************
1307 // SetupHierarchy
1308 //
1309 // Renavigates & resets hierarchy described by current history
1310 // o Reset volumes
1311 // o Recompute transforms and/or solids of replicated/parameterised volumes
1312 // ********************************************************************
1313 //
1314 void G4ITNavigator1::SetupHierarchy()
1315 {
1316   G4int i;
1317   const auto  cdepth = (G4int)fHistory.GetDepth();
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(fHistory.GetReplicaNo(i), current);
1331         break;
1332       case kParameterised: {
1333           G4int replicaNo;
1334           pParam = current->GetParameterisation();
1335           replicaNo = fHistory.GetReplicaNo(i);
1336           pSolid = pParam->ComputeSolid(replicaNo, current);
1337 
1338           // Set up dimensions & transform in solid/physical volume
1339           //
1340           pSolid->ComputeDimensions(pParam, replicaNo, current);
1341           pParam->ComputeTransformation(replicaNo, current);
1342 
1343           G4TouchableHistory touchable( fHistory );
1344           touchable.MoveUpHistory();  // move up to the parent level
1345         
1346           // Set up the correct solid and material in Logical Volume
1347           //
1348           G4LogicalVolume *pLogical = current->GetLogicalVolume();
1349           pLogical->SetSolid( pSolid );
1350           pLogical->UpdateMaterial( pParam ->
1351             ComputeMaterial(replicaNo, current, &touchable) );
1352           break;
1353       }
1354       case kExternal:
1355         G4Exception("G4ITNavigator1::SetupHierarchy()",
1356                     "GeomNav0001", FatalException,
1357                     "Not applicable for external volumes.");
1358         break;
1359     }
1360   }
1361 }
1362 
1363 // ********************************************************************
1364 // GetLocalExitNormal
1365 //
1366 // Obtains the Normal vector to a surface (in local coordinates)
1367 // pointing out of previous volume and into current volume
1368 // ********************************************************************
1369 //
1370 G4ThreeVector G4ITNavigator1::GetLocalExitNormal( G4bool* valid )
1371 {
1372   G4ThreeVector    ExitNormal(0.,0.,0.);
1373   G4VSolid        *currentSolid=nullptr;
1374   G4LogicalVolume *candidateLogical;
1375   if ( fLastTriedStepComputation ) 
1376   {
1377     // use fLastLocatedPointLocal and next candidate volume
1378     //
1379     G4ThreeVector nextSolidExitNormal(0.,0.,0.);
1380 
1381     if( fEntering && (fBlockedPhysicalVolume!=nullptr) ) 
1382     { 
1383       candidateLogical= fBlockedPhysicalVolume->GetLogicalVolume();
1384       if( candidateLogical != nullptr ) 
1385       {
1386         // fLastStepEndPointLocal is in the coordinates of the mother
1387         // we need it in the daughter's coordinate system.
1388 
1389         // The following code should also work in case of Replica
1390         {
1391           // First transform fLastLocatedPointLocal to the new daughter
1392           // coordinates
1393           //
1394           G4AffineTransform MotherToDaughterTransform=
1395             GetMotherToDaughterTransform( fBlockedPhysicalVolume, 
1396                                           fBlockedReplicaNo,
1397                                           VolumeType(fBlockedPhysicalVolume) ); 
1398           G4ThreeVector daughterPointOwnLocal= 
1399             MotherToDaughterTransform.TransformPoint( fLastStepEndPointLocal ); 
1400 
1401           // OK if it is a parameterised volume
1402           //
1403           EInside  inSideIt; 
1404           G4bool   onSurface;
1405           G4double safety= -1.0; 
1406           currentSolid= candidateLogical->GetSolid(); 
1407           inSideIt  =  currentSolid->Inside(daughterPointOwnLocal); 
1408           onSurface =  (inSideIt == kSurface); 
1409           if( ! onSurface ) 
1410           {
1411             if( inSideIt == kOutside )
1412             {
1413               safety = (currentSolid->DistanceToIn(daughterPointOwnLocal)); 
1414               onSurface = safety < 100.0 * kCarTolerance; 
1415             }
1416             else if (inSideIt == kInside )
1417             {
1418               safety = (currentSolid->DistanceToOut(daughterPointOwnLocal)); 
1419               onSurface = safety < 100.0 * kCarTolerance; 
1420             }
1421           }
1422 
1423           if( onSurface ) 
1424           {
1425             nextSolidExitNormal =
1426               currentSolid->SurfaceNormal(daughterPointOwnLocal); 
1427  
1428             // Entering the solid ==> opposite
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 surface ! " << G4endl
1440                       << "  Point           = "
1441                       << daughterPointOwnLocal << G4endl 
1442                       << "  Physical volume = "
1443                       << fBlockedPhysicalVolume->GetName() << G4endl
1444                       << "  Logical volume  = "
1445                       << candidateLogical->GetName() << G4endl
1446                       << "  Solid           = " << currentSolid->GetName() 
1447                       << "  Type            = "
1448                       << currentSolid->GetEntityType() << G4endl
1449                       << *currentSolid << G4endl;
1450               if( inSideIt == kOutside )
1451               { 
1452                 message << "Point is Outside. " << G4endl
1453                         << "  Safety (from outside) = " << safety << G4endl;
1454               }
1455               else // if( inSideIt == kInside ) 
1456               {
1457                 message << "Point is Inside. " << G4endl
1458                         << "  Safety (from inside) = " << safety << G4endl;              
1459               }
1460               G4Exception("G4ITNavigator1::GetLocalExitNormal()", "GeomNav1001",
1461                           JustWarning, message);
1462             }
1463 #endif
1464           }
1465           *valid = onSurface;   //   was =true;
1466         }
1467       }
1468     }
1469     else if ( fExiting ) 
1470     {
1471       ExitNormal = fGrandMotherExitNormal;
1472       *valid = true;
1473       fCalculatedExitNormal= true;  // Should be true already
1474     }
1475     else  // i.e.  ( fBlockedPhysicalVolume == 0 )
1476     {
1477       *valid = false;
1478       G4Exception("G4ITNavigator1::GetLocalExitNormal()",
1479                   "GeomNav0003", JustWarning, 
1480                   "Incorrect call to GetLocalSurfaceNormal." );
1481     }
1482   }
1483   else //  ( ! fLastTriedStepComputation ) ie. last call was to Locate
1484   {
1485     if ( EnteredDaughterVolume() )
1486     {
1487       G4VSolid* daughterSolid =fHistory.GetTopVolume()->GetLogicalVolume()
1488                                                       ->GetSolid();
1489       ExitNormal= -(daughterSolid->SurfaceNormal(fLastLocatedPointLocal));
1490       if( std::fabs(ExitNormal.mag2()-1.0 ) > CLHEP::perMillion )
1491       {
1492         G4ExceptionDescription desc;
1493         desc << " Parameters of solid: " << *daughterSolid
1494              << " Point for surface = " << fLastLocatedPointLocal << std::endl;
1495         G4Exception("G4ITNavigator1::GetLocalExitNormal()",
1496                     "GeomNav0003", FatalException, desc,
1497                     "Surface Normal returned by Solid is not a Unit Vector." );
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. ExitNormal remains (0,0,0)
1511       { 
1512         *valid = false;
1513         fCalculatedExitNormal= false; 
1514         G4ExceptionDescription message; 
1515         message << "Function called when *NOT* at a Boundary." << G4endl;
1516         G4Exception("G4ITNavigator1::GetLocalExitNormal()",
1517                     "GeomNav0003", JustWarning, message); 
1518       }
1519     }
1520   }
1521   return ExitNormal;
1522 }
1523 
1524 // ********************************************************************
1525 // GetMotherToDaughterTransform
1526 //
1527 // Obtains the mother to daughter affine transformation
1528 // ********************************************************************
1529 //
1530 G4AffineTransform
1531 G4ITNavigator1::GetMotherToDaughterTransform( G4VPhysicalVolume *pEnteringPhysVol,   // not Const
1532                                            G4int   enteringReplicaNo,
1533                                            EVolume enteringVolumeType ) 
1534 {
1535   switch (enteringVolumeType)
1536   {
1537     case kNormal:  // Nothing is needed to prepare the transformation
1538       break;       // It is stored already in the physical volume (placement)
1539     case kReplica: // Sets the transform in the Replica - tbc
1540       G4Exception("G4ITNavigator1::GetMotherToDaughterTransform()",
1541                   "GeomNav0001", FatalException,
1542                   "Method NOT Implemented yet for replica volumes.");
1543       break;
1544     case kParameterised:
1545       if( pEnteringPhysVol->GetRegularStructureId() == 0 )
1546       {
1547         G4VPVParameterisation *pParam =
1548           pEnteringPhysVol->GetParameterisation();
1549         G4VSolid* pSolid =
1550           pParam->ComputeSolid(enteringReplicaNo, pEnteringPhysVol);
1551         pSolid->ComputeDimensions(pParam, enteringReplicaNo, pEnteringPhysVol);
1552 
1553         // Sets the transform in the Parameterisation
1554         //
1555         pParam->ComputeTransformation(enteringReplicaNo, pEnteringPhysVol);
1556 
1557         // Set the correct solid and material in Logical Volume
1558         //
1559         G4LogicalVolume* pLogical = pEnteringPhysVol->GetLogicalVolume();
1560         pLogical->SetSolid( pSolid );
1561       }
1562       break;
1563       case kExternal:
1564         G4Exception("G4ITNavigator1::GetMotherToDaughterTransform()",
1565                     "GeomNav0001", FatalException,
1566                     "Not applicable for external volumes.");
1567         break;
1568   }
1569   return G4AffineTransform(pEnteringPhysVol->GetRotation(), 
1570                            pEnteringPhysVol->GetTranslation()).Invert(); 
1571 }
1572 
1573 // ********************************************************************
1574 // GetLocalExitNormalAndCheck
1575 //
1576 // Obtains the Normal vector to a surface (in local coordinates)
1577 // pointing out of previous volume and into current volume, and
1578 // checks the current point against expected 'local' value.
1579 // ********************************************************************
1580 //
1581 G4ThreeVector G4ITNavigator1::
1582 GetLocalExitNormalAndCheck( 
1583 #ifdef G4DEBUG_NAVIGATION
1584                            const G4ThreeVector& ExpectedBoundaryPointGlobal,
1585 #else
1586                            const G4ThreeVector&,
1587 #endif
1588                                  G4bool*        pValid)
1589 {
1590 #ifdef G4DEBUG_NAVIGATION
1591   // Check Current point against expected 'local' value
1592   //
1593   if ( fLastTriedStepComputation ) 
1594   {
1595     G4ThreeVector ExpectedBoundaryPointLocal;
1596 
1597     const G4AffineTransform& GlobalToLocal= GetGlobalToLocalTransform(); 
1598     ExpectedBoundaryPointLocal =
1599       GlobalToLocal.TransformPoint( ExpectedBoundaryPointGlobal ); 
1600 
1601     // Add here:  Comparison against expected position,
1602     //            i.e. the endpoint of ComputeStep
1603   }
1604 #endif
1605   
1606   return GetLocalExitNormal( pValid); 
1607 }
1608 
1609 // ********************************************************************
1610 // GetGlobalExitNormal
1611 //
1612 // Obtains the Normal vector to a surface (in global coordinates)
1613 // pointing out of previous volume and into current volume
1614 // ********************************************************************
1615 //
1616 G4ThreeVector 
1617 G4ITNavigator1::GetGlobalExitNormal(const G4ThreeVector& IntersectPointGlobal,
1618                                        G4bool*        pNormalCalculated)
1619 {
1620   G4bool         validNormal;
1621   G4ThreeVector  localNormal, globalNormal;
1622 
1623   if( fLastTriedStepComputation && fExiting )  
1624   {
1625     // This was computed in ComputeStep -- and only on arrival at boundary
1626     //
1627     globalNormal = fExitNormalGlobalFrame; 
1628     *pNormalCalculated = true; // ComputeStep always computes it if Exiting
1629                                // (fExiting==true)
1630   }
1631   else
1632   {
1633     localNormal = GetLocalExitNormalAndCheck(IntersectPointGlobal,&validNormal);
1634     *pNormalCalculated = fCalculatedExitNormal;
1635 
1636 #ifdef G4DEBUG_NAVIGATION
1637     if( (!validNormal) && !fCalculatedExitNormal)
1638     {
1639       G4ExceptionDescription edN;
1640       edN << "  Calculated = " << fCalculatedExitNormal << G4endl;
1641       edN << "   Entering= "  << fEntering << G4endl;
1642       G4int oldVerbose= this->GetVerboseLevel();
1643       this->SetVerboseLevel(4);
1644       edN << "   State of Navigator: " << G4endl;
1645       edN << *this << G4endl;
1646       this->SetVerboseLevel( oldVerbose );
1647        
1648       G4Exception("G4ITNavigator1::GetGlobalExitNormal()",
1649                   "GeomNav0003", JustWarning, edN,
1650                   "LocalExitNormalAndCheck() did not calculate Normal.");
1651      }
1652 #endif
1653      
1654      G4double localMag2= localNormal.mag2();
1655      if( validNormal && (std::fabs(localMag2-1.0)) > CLHEP::perMillion )
1656      {
1657        G4ExceptionDescription edN;
1658 
1659        edN << "G4ITNavigator1::GetGlobalExitNormal: "
1660            << "  Using Local Normal - from call to GetLocalExitNormalAndCheck. "
1661            << G4endl
1662            << "  Local  Exit Normal = " << localNormal  << " || = "
1663            << std::sqrt(localMag2) << G4endl
1664            << "  Global Exit Normal = " << globalNormal << " || = "
1665            << globalNormal.mag() << G4endl;
1666        edN << "  Calculated It      = " << fCalculatedExitNormal << G4endl;
1667 
1668        G4Exception("G4ITNavigator1::GetGlobalExitNormal()",
1669                    "GeomNav0003",JustWarning, edN,
1670                    "Value obtained from new local *solid* is incorrect.");
1671        localNormal = localNormal.unit(); // Should we correct it ??
1672      }
1673      G4AffineTransform localToGlobal = GetLocalToGlobalTransform();
1674      globalNormal = localToGlobal.TransformAxis( localNormal );
1675   }
1676 
1677 #ifdef G4DEBUG_NAVIGATION
1678   // Temporary extra checks
1679   if( fLastTriedStepComputation && fExiting)
1680   {
1681     localNormal = GetLocalExitNormalAndCheck( IntersectPointGlobal, &validNormal);
1682     *pNormalCalculated = fCalculatedExitNormal;
1683 
1684     G4AffineTransform localToGlobal = GetLocalToGlobalTransform();
1685     globalNormal = localToGlobal.TransformAxis( localNormal );
1686     
1687     // Check the value computed against fExitNormalGlobalFrame
1688     G4ThreeVector diffNorm = globalNormal - fExitNormalGlobalFrame;
1689     if( diffNorm.mag2() > perMillion*CLHEP::perMillion)
1690     {
1691       G4ExceptionDescription edDfn;
1692       edDfn << "Found difference in normals in case of exiting mother "
1693             << "- when Get is called after ComputingStep " << G4endl;
1694       edDfn << "  Magnitude of diff =      " << diffNorm.mag() << G4endl;
1695       edDfn << "  Normal stored (Global)     = " << fExitNormalGlobalFrame
1696             << G4endl;
1697       edDfn << "  Global Computed from Local = " << globalNormal << G4endl;
1698       G4Exception("G4ITNavigator1::GetGlobalExitNormal()", "GeomNav0003",
1699                   JustWarning, edDfn);
1700     }
1701   }
1702 #endif
1703    
1704   return globalNormal;
1705 }
1706 
1707 // To make the new Voxel Safety the default, uncomment the next line
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 as the EndPoint of the
1715 //     ComputeStep.
1716 // ii) after (or at the end of) ComputeStep OR after the relocation.
1717 // ********************************************************************
1718 //
1719 G4double G4ITNavigator1::ComputeSafety( const G4ThreeVector &pGlobalpoint,
1720                                      const G4double pMaxLength,
1721                                      const G4bool keepState)
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::ComputeSafety: ***" << G4endl
1730            << "    Called at point: " << pGlobalpoint << G4endl;
1731 
1732     G4VPhysicalVolume  *motherPhysical = fHistory.GetTopVolume();
1733     G4cout << "    Volume = " << motherPhysical->GetName() 
1734            << " - Maximum length = " << pMaxLength << G4endl; 
1735     if( fVerbose >= 4 )
1736     {
1737        G4cout << "    ----- Upon entering Compute Safety:" << G4endl;
1738        PrintState();
1739     }
1740   }
1741 #endif
1742 
1743   if (keepState)  { SetSavedState(); }
1744 
1745   G4double distEndpointSq = (pGlobalpoint-fStepEndPoint).mag2(); 
1746   G4bool   stayedOnEndpoint  = distEndpointSq < kCarTolerance*kCarTolerance; 
1747   G4bool   endpointOnSurface = fEnteredDaughter || fExitedMother;
1748 
1749   if( !(endpointOnSurface && stayedOnEndpoint) )
1750   {
1751     // Pseudo-relocate to this point (updates voxel information only)
1752     //
1753     LocateGlobalPointWithinVolume( pGlobalpoint ); 
1754       // --->> DANGER: Side effects on sub-navigator voxel information <<---
1755       //       Could be replaced again by 'granular' calls to sub-navigator
1756       //       locates (similar side-effects, but faster.  
1757       //       Solutions:
1758       //        1) Re-locate (to where?)
1759       //        2) Insure that the methods using (G4ComputeStep?)
1760       //           does a relocation (if information is disturbed only ?)
1761 
1762 #ifdef G4DEBUG_NAVIGATION
1763     if( fVerbose >= 2 )
1764     {
1765       G4cout << "  G4ITNavigator1::ComputeSafety() relocates-in-volume to point: "
1766              << pGlobalpoint << G4endl;
1767     }
1768 #endif 
1769     G4VPhysicalVolume *motherPhysical = fHistory.GetTopVolume();
1770     G4LogicalVolume *motherLogical = motherPhysical->GetLogicalVolume();
1771     G4SmartVoxelHeader* pVoxelHeader = motherLogical->GetVoxelHeader();
1772     G4ThreeVector localPoint = ComputeLocalPoint(pGlobalpoint);
1773 
1774     if ( fHistory.GetTopVolumeType()!=kReplica )
1775     {
1776       switch(CharacteriseDaughters(motherLogical))
1777       {
1778         case kNormal:
1779           if ( pVoxelHeader != nullptr )
1780           {
1781 #ifdef G4NEW_SAFETY
1782             G4double safetyTwo = fpVoxelSafety->ComputeSafety(localPoint,
1783                                            *motherPhysical, pMaxLength);
1784             newSafety= safetyTwo;   // Faster and best available
1785 #else
1786             G4double safetyOldVoxel;
1787             safetyOldVoxel =
1788               fvoxelNav.ComputeSafety(localPoint,fHistory,pMaxLength);
1789             newSafety= safetyOldVoxel;
1790 #endif
1791           }
1792           else
1793           {
1794             newSafety=fnormalNav.ComputeSafety(localPoint,fHistory,pMaxLength);
1795           }
1796           break;
1797         case kParameterised:
1798           if( GetDaughtersRegularStructureId(motherLogical) != 1 )
1799           {
1800             newSafety=fparamNav.ComputeSafety(localPoint,fHistory,pMaxLength);
1801           }
1802           else  // Regular structure
1803           {
1804             newSafety=fregularNav.ComputeSafety(localPoint,fHistory,pMaxLength);
1805           }
1806           break;
1807         case kReplica:
1808           G4Exception("G4ITNavigator1::ComputeSafety()", "GeomNav0001",
1809                       FatalException, "Not applicable for replicated volumes.");
1810           break;
1811         case kExternal:
1812           G4Exception("G4ITNavigator1::ComputeSafety()",
1813                       "GeomNav0001", FatalException,
1814                       "Not applicable for external volumes.");
1815          break;
1816       }
1817     }
1818     else
1819     {
1820       newSafety = freplicaNav.ComputeSafety(pGlobalpoint, localPoint,
1821                                             fHistory, pMaxLength);
1822     }
1823   }
1824   else // if( endpointOnSurface && stayedOnEndpoint )
1825   {
1826 #ifdef G4DEBUG_NAVIGATION
1827     if( fVerbose >= 2 )
1828     {
1829       G4cout << "    G4ITNavigator1::ComputeSafety() finds that point - "
1830              << pGlobalpoint << " - is on surface " << G4endl; 
1831       if( fEnteredDaughter ) { G4cout << "   entered new daughter volume"; }
1832       if( fExitedMother )    { G4cout << "   and exited previous volume."; }
1833       G4cout << G4endl;
1834       G4cout << " EndPoint was = " << fStepEndPoint << G4endl;
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  " << G4endl;
1851     if( fVerbose > 2 )  { PrintState(); }
1852     G4cout << "    Returned value of Safety = " << newSafety << G4endl;
1853   }
1854   G4cout.precision(oldcoutPrec);
1855 #endif
1856 
1857   return newSafety;
1858 }
1859 
1860 // ********************************************************************
1861 // CreateTouchableHistoryHandle
1862 // ********************************************************************
1863 //
1864 G4TouchableHandle G4ITNavigator1::CreateTouchableHistoryHandle() const
1865 {
1866   return G4TouchableHandle( CreateTouchableHistory() );
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 G4ITNavigator1 is: " << G4endl;
1879     G4cout << "  ValidExitNormal= " << fValidExitNormal << G4endl
1880            << "  ExitNormal     = " << fExitNormal      << G4endl
1881            << "  Exiting        = " << fExiting         << G4endl
1882            << "  Entering       = " << fEntering        << G4endl
1883            << "  BlockedPhysicalVolume= " ;
1884     if (fBlockedPhysicalVolume==nullptr)
1885       G4cout << "None";
1886     else
1887       G4cout << fBlockedPhysicalVolume->GetName();
1888     G4cout << G4endl
1889            << "  BlockedReplicaNo     = " <<  fBlockedReplicaNo       << G4endl
1890            << "  LastStepWasZero      = " <<   fLastStepWasZero       << G4endl
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:Volume "  << " "   
1901            << std::setw( 9) << " ReplicaNo"        << " "  
1902            << std::setw( 8) << " LastStepZero  "   << " "   
1903            << G4endl;   
1904     G4cout << "( " << std::setw(7) << fExitNormal.x() 
1905            << ", " << std::setw(7) << fExitNormal.y()
1906            << ", " << std::setw(7) << fExitNormal.z() << " ) "
1907            << std::setw( 5)  << fValidExitNormal  << " "   
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)<< fBlockedPhysicalVolume->GetName();
1917     }
1918     G4cout << std::setw( 9)  << fBlockedReplicaNo  << " "
1919            << std::setw( 8)  << fLastStepWasZero   << " "
1920            << G4endl;   
1921   }
1922   if( fVerbose > 2 ) 
1923   {
1924     G4cout.precision(8);
1925     G4cout << " Current Localpoint = " << fLastLocatedPointLocal << G4endl;
1926     G4cout << " PreviousSftOrigin  = " << fPreviousSftOrigin << G4endl;
1927     G4cout << " PreviousSafety     = " << fPreviousSafety << G4endl; 
1928   }
1929   G4cout.precision(oldcoutPrec);
1930 }
1931 
1932 // ********************************************************************
1933 // ComputeStepLog
1934 // ********************************************************************
1935 //
1936 void G4ITNavigator1::ComputeStepLog(const G4ThreeVector& pGlobalpoint,
1937                                        G4double moveLenSq) const
1938 {
1939   //  The following checks only make sense if the move is larger
1940   //  than the tolerance.
1941 
1942   static const G4double fAccuracyForWarning   = kCarTolerance,
1943                         fAccuracyForException = 1000*kCarTolerance;
1944 
1945   G4ThreeVector OriginalGlobalpoint = fHistory.GetTopTransform().Inverse().
1946                                       TransformPoint(fLastLocatedPointLocal); 
1947 
1948   G4double shiftOriginSafSq = (fPreviousSftOrigin-pGlobalpoint).mag2();
1949 
1950   // Check that the starting point of this step is 
1951   // within the isotropic safety sphere of the last point
1952   // to a accuracy/precision  given by fAccuracyForWarning.
1953   //   If so give warning.
1954   //   If it fails by more than fAccuracyForException exit with error.
1955   //
1956   if( shiftOriginSafSq >= sqr(fPreviousSafety) )
1957   {
1958     G4double shiftOrigin = std::sqrt(shiftOriginSafSq);
1959     G4double diffShiftSaf = shiftOrigin - fPreviousSafety;
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 inaccurate position shift."
1967               << G4endl
1968               << "     The Step's starting point has moved " 
1969               << std::sqrt(moveLenSq)/mm << " mm " << G4endl
1970               << "     since the last call to a Locate method." << G4endl
1971               << "     This has resulted in moving " 
1972               << shiftOrigin/mm << " mm " 
1973               << " from the last point at which the safety " 
1974               << "     was calculated " << G4endl
1975               << "     which is more than the computed safety= " 
1976               << fPreviousSafety/mm << " mm  at that point." << G4endl
1977               << "     This difference is " 
1978               << diffShiftSaf/mm << " mm." << G4endl
1979               << "     The tolerated accuracy is "
1980               << fAccuracyForException/mm << " 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 to either " << G4endl
1988                << "    - a process that has proposed a displacement"
1989                << " larger than the current safety , or" << G4endl
1990                << "    - inaccuracy in the computation of the safety";
1991         suggestion << "We suggest that you " << G4endl
1992                    << "   - find i) what particle is being tracked, and "
1993                    << " ii) through what part of your geometry " << G4endl
1994                    << "      for example by re-running this event with "
1995                    << G4endl
1996                    << "         /tracking/verbose 1 "  << G4endl
1997                    << "    - check which processes you declare for"
1998                    << " this particle (and look at non-standard ones)"
1999                    << G4endl
2000                    << "   - in case, create a detailed logfile"
2001                    << " of this event using:" << G4endl
2002                    << "         /tracking/verbose 6 ";
2003       }
2004       G4Exception("G4ITNavigator1::ComputeStep()",
2005                   "GeomNav1002", JustWarning,
2006                   message, G4String(suggestion.str()));
2007       G4cout.precision(oldcoutPrec);
2008       G4cerr.precision(oldcerrPrec);
2009     }
2010 #ifdef G4DEBUG_NAVIGATION
2011     else
2012     {
2013       G4cerr << "WARNING - G4ITNavigator1::ComputeStep()" << G4endl
2014              << "          The Step's starting point has moved "
2015              << std::sqrt(moveLenSq) << "," << G4endl
2016              << "          which has taken it to the limit of"
2017              << " the current safety. " << G4endl;
2018     }
2019 #endif
2020   }
2021   G4double safetyPlus = fPreviousSafety + fAccuracyForException;
2022   if ( shiftOriginSafSq > sqr(safetyPlus) )
2023   {
2024     std::ostringstream message;
2025     message << "May lead to a crash or unreliable results." << G4endl
2026             << "        Position has shifted considerably without"
2027             << " notifying the navigator !" << G4endl
2028             << "        Tolerated safety: " << safetyPlus << G4endl
2029             << "        Computed shift  : " << shiftOriginSafSq;
2030     G4Exception("G4ITNavigator1::ComputeStep()", "GeomNav1002",
2031                 JustWarning, message);
2032   }
2033 }
2034 
2035 // ********************************************************************
2036 // Operator <<
2037 // ********************************************************************
2038 //
2039 std::ostream& operator << (std::ostream &os,const G4ITNavigator1 &n)
2040 {
2041   //  Old version did only the following:
2042   // os << "Current History: " << G4endl << n.fHistory;
2043   //  Old behaviour is recovered for fVerbose = 0
2044   
2045   // Adapted from G4ITNavigator1::PrintState() const
2046 
2047   G4long oldcoutPrec = os.precision(4);
2048   if( n.fVerbose >= 4 )
2049   {
2050     os << "The current state of G4ITNavigator1 is: " << G4endl;
2051     os << "  ValidExitNormal= " << n.fValidExitNormal << G4endl
2052     << "  ExitNormal     = " << n.fExitNormal      << G4endl
2053     << "  Exiting        = " << n.fExiting         << G4endl
2054     << "  Entering       = " << n.fEntering        << G4endl
2055     << "  BlockedPhysicalVolume= " ;
2056     if (n.fBlockedPhysicalVolume==nullptr)
2057       os << "None";
2058     else
2059       os << n.fBlockedPhysicalVolume->GetName();
2060     os << G4endl
2061     << "  BlockedReplicaNo     = " <<  n.fBlockedReplicaNo       << G4endl
2062     << "  LastStepWasZero      = " <<   n.fLastStepWasZero       << G4endl
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.fExitNormal.x()
2077     << ", " << std::setw(7) << n.fExitNormal.y()
2078     << ", " << std::setw(7) << n.fExitNormal.z() << " ) "
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.fBlockedPhysicalVolume->GetName(); }
2086     os << std::setw( 9)  << n.fBlockedReplicaNo  << " "
2087     << std::setw( 8)  << n.fLastStepWasZero   << " "
2088     << G4endl;
2089   }
2090   if( n.fVerbose > 2 )
2091   {
2092     os.precision(8);
2093     os << " Current Localpoint = " << n.fLastLocatedPointLocal << G4endl;
2094     os << " PreviousSftOrigin  = " << n.fPreviousSftOrigin << G4endl;
2095     os << " PreviousSafety     = " << n.fPreviousSafety << G4endl;
2096   }
2097   if( n.fVerbose > 3 || n.fVerbose == 0 )
2098   {
2099     os << "Current History: " << G4endl << n.fHistory;
2100   }
2101     
2102   os.precision(oldcoutPrec);
2103   return os;
2104 }
2105