Geant4  10.01.p02
G4Navigator.cc
Go to the documentation of this file.
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 // $Id: G4Navigator.cc 2012-11-01 japost Exp $
28 // GEANT4 tag $ Name: $
29 //
30 // class G4Navigator Implementation
31 //
32 // Original author: Paul Kent, July 95/96
33 // Responsible 1996-present: John Apostolakis
34 // Co-maintainer: Gabriele Cosmo
35 // Additional revisions by: Pedro Arce, Vladimir Grichine
36 // --------------------------------------------------------------------
37 
38 #include <iomanip>
39 
40 #include "G4Navigator.hh"
41 #include "G4ios.hh"
42 #include "G4SystemOfUnits.hh"
43 #include "G4GeometryTolerance.hh"
44 #include "G4VPhysicalVolume.hh"
45 
46 #include "G4VoxelSafety.hh"
47 
48 // ********************************************************************
49 // Constructor
50 // ********************************************************************
51 //
53  : fWasLimitedByGeometry(false), fVerbose(0),
54  fTopPhysical(0), fCheck(false), fPushed(false), fWarnPush(true)
55 {
56  fActive= false;
58 
60  // Initialises also all
61  // - exit / entry flags
62  // - flags & variables for exit normals
63  // - zero step counters
64  // - blocked volume
65 
68 
71 
74 
76 }
77 
78 // ********************************************************************
79 // Destructor
80 // ********************************************************************
81 //
83 {
84  delete fpVoxelSafety;
85 }
86 
87 // ********************************************************************
88 // ResetHierarchyAndLocate
89 // ********************************************************************
90 //
93  const G4ThreeVector &direction,
94  const G4TouchableHistory &h)
95 {
96  ResetState();
97  fHistory = *h.GetHistory();
99  fLastTriedStepComputation= false; // Redundant, but best
100  return LocateGlobalPointAndSetup(p, &direction, true, false);
101 }
102 
103 // ********************************************************************
104 // LocateGlobalPointAndSetup
105 //
106 // Locate the point in the hierarchy return 0 if outside
107 // The direction is required
108 // - if on an edge shared by more than two surfaces
109 // (to resolve likely looping in tracking)
110 // - at initial location of a particle
111 // (to resolve potential ambiguity at boundary)
112 //
113 // Flags on exit: (comments to be completed)
114 // fEntering - True if entering `daughter' volume (or replica)
115 // whether daughter of last mother directly
116 // or daughter of that volume's ancestor.
117 // fExiting - True if exited 'mother' volume
118 // (always ? - how about if going back down ? - tbc)
119 // ********************************************************************
120 //
123  const G4ThreeVector* pGlobalDirection,
124  const G4bool relativeSearch,
125  const G4bool ignoreDirection )
126 {
127  G4bool notKnownContained=true, noResult;
128  G4VPhysicalVolume *targetPhysical;
129  G4LogicalVolume *targetLogical;
130  G4VSolid *targetSolid=0;
131  G4ThreeVector localPoint, globalDirection;
132  EInside insideCode;
133 
134  G4bool considerDirection = (!ignoreDirection) || fLocatedOnEdge;
135 
137  fChangedGrandMotherRefFrame= false; // For local exit normal
138 
139  if( considerDirection && pGlobalDirection != 0 )
140  {
141  globalDirection=*pGlobalDirection;
142  }
143 
144 #ifdef G4VERBOSE
145  if( fVerbose > 2 )
146  {
147  G4int oldcoutPrec = G4cout.precision(8);
148  G4cout << "*** G4Navigator::LocateGlobalPointAndSetup: ***" << G4endl;
149  G4cout << " Called with arguments: " << G4endl
150  << " Globalpoint = " << globalPoint << G4endl
151  << " RelativeSearch = " << relativeSearch << G4endl;
152  if( fVerbose >= 4 )
153  {
154  G4cout << " ----- Upon entering:" << G4endl;
155  PrintState();
156  }
157  G4cout.precision(oldcoutPrec);
158  }
159 #endif
160 
161  G4int noLevelsExited=0 ;
162  G4int noLevelsEntered= 0;
163 
164  if ( !relativeSearch )
165  {
167  }
168  else
169  {
170  if ( fWasLimitedByGeometry )
171  {
172  fWasLimitedByGeometry = false;
173  fEnteredDaughter = fEntering; // Remember
174  fExitedMother = fExiting; // Remember
175  if ( fExiting )
176  {
177  noLevelsExited++; // count this first level entered too
178 
179  if ( fHistory.GetDepth() )
180  {
184  }
185  else
186  {
187  fLastLocatedPointLocal = localPoint;
188  fLocatedOutsideWorld = true;
189  fBlockedPhysicalVolume = 0; // to be sure
190  fBlockedReplicaNo = -1;
191  fEntering = false; // No longer
192  fEnteredDaughter = false;
193  fExitedMother = true; // ??
194 
195  return 0; // Have exited world volume
196  }
197  // A fix for the case where a volume is "entered" at an edge
198  // and a coincident surface exists outside it.
199  // - This stops it from exiting further volumes and cycling
200  // - However ReplicaNavigator treats this case itself
201  //
202  // assert( fBlockedPhysicalVolume!=0 );
203 
204  // Expect to be on edge => on surface
205  //
207  {
208  fExiting= false;
209  // Consider effect on Exit Normal !?
210  }
211  }
212  else
213  if ( fEntering )
214  {
215  // assert( fBlockedPhysicalVolume!=0 );
216 
217  noLevelsEntered++; // count the first level entered too
218 
220  {
221  case kNormal:
224  break;
225  case kReplica:
231  break;
232  case kParameterised:
234  {
235  G4VSolid *pSolid;
236  G4VPVParameterisation *pParam;
237  G4TouchableHistory parentTouchable( fHistory );
239  pSolid = pParam->ComputeSolid(fBlockedReplicaNo,
241  pSolid->ComputeDimensions(pParam, fBlockedReplicaNo,
248  //
249  // Set the correct solid and material in Logical Volume
250  //
251  G4LogicalVolume *pLogical;
253  pLogical->SetSolid( pSolid );
254  pLogical->UpdateMaterial(pParam ->
255  ComputeMaterial(fBlockedReplicaNo,
257  &parentTouchable));
258  }
259  break;
260  }
261  fEntering = false;
263  localPoint = fHistory.GetTopTransform().TransformPoint(globalPoint);
264  notKnownContained = false;
265  }
266  }
267  else
268  {
270  fEntering = false;
271  fEnteredDaughter = false; // Full Step was not taken, did not enter
272  fExiting = false;
273  fExitedMother = false; // Full Step was not taken, did not exit
274  }
275  }
276  //
277  // Search from top of history up through geometry until
278  // containing volume found:
279  // If on
280  // o OUTSIDE - Back up level, not/no longer exiting volumes
281  // o SURFACE and EXITING - Back up level, setting new blocking no.s
282  // else
283  // o containing volume found
284  //
285 
286  while (notKnownContained)
287  {
289  {
290  targetSolid = fHistory.GetTopVolume()->GetLogicalVolume()->GetSolid();
291  localPoint = fHistory.GetTopTransform().TransformPoint(globalPoint);
292  insideCode = targetSolid->Inside(localPoint);
293 #ifdef G4VERBOSE
294  if(( fVerbose == 1 ) && ( fCheck ))
295  {
296  G4String solidResponse = "-kInside-";
297  if (insideCode == kOutside)
298  solidResponse = "-kOutside-";
299  else if (insideCode == kSurface)
300  solidResponse = "-kSurface-";
301  G4cout << "*** G4Navigator::LocateGlobalPointAndSetup(): ***" << G4endl
302  << " Invoked Inside() for solid: " << targetSolid->GetName()
303  << ". Solid replied: " << solidResponse << G4endl
304  << " For local point p: " << localPoint << G4endl;
305  }
306 #endif
307  }
308  else
309  {
310  insideCode = freplicaNav.BackLocate(fHistory, globalPoint, localPoint,
311  fExiting, notKnownContained);
312  // !CARE! if notKnownContained returns false then the point is within
313  // the containing placement volume of the replica(s). If insidecode
314  // will result in the history being backed up one level, then the
315  // local point returned is the point in the system of this new level
316  }
317 
318 
319  if ( insideCode==kOutside )
320  {
321  noLevelsExited++;
322  if ( fHistory.GetDepth() )
323  {
327  fExiting = false;
328 
329  if( noLevelsExited > 1 )
330  {
331  // The first transformation was done by the sub-navigator
332  //
334  if( mRot )
335  {
336  fGrandMotherExitNormal *= (*mRot).inverse();
338  }
339  }
340  }
341  else
342  {
343  fLastLocatedPointLocal = localPoint;
344  fLocatedOutsideWorld = true;
345  // No extra transformation for ExitNormal - is in frame of Top Volume
346  return 0; // Have exited world volume
347  }
348  }
349  else
350  if ( insideCode==kSurface )
351  {
352  G4bool isExiting = fExiting;
353  if( (!fExiting)&&considerDirection )
354  {
355  // Figure out whether we are exiting this level's volume
356  // by using the direction
357  //
358  G4bool directionExiting = false;
359  G4ThreeVector localDirection =
360  fHistory.GetTopTransform().TransformAxis(globalDirection);
361 
362  // Make sure localPoint in correct reference frame
363  // ( Was it already correct ? How ? )
364  //
365  localPoint= fHistory.GetTopTransform().TransformPoint(globalPoint);
367  {
368  G4ThreeVector normal = targetSolid->SurfaceNormal(localPoint);
369  directionExiting = normal.dot(localDirection) > 0.0;
370  isExiting = isExiting || directionExiting;
371  }
372  }
373  if( isExiting )
374  {
375  noLevelsExited++;
376  if ( fHistory.GetDepth() )
377  {
381  //
382  // Still on surface but exited volume not necessarily convex
383  //
384  fValidExitNormal = false;
385 
386  if( noLevelsExited > 1 )
387  {
388  // The first transformation was done by the sub-navigator
389  //
390  const G4RotationMatrix* mRot =
392  if( mRot )
393  {
394  fGrandMotherExitNormal *= (*mRot).inverse();
396  }
397  }
398  }
399  else
400  {
401  fLastLocatedPointLocal = localPoint;
402  fLocatedOutsideWorld = true;
403  // No extra transformation for ExitNormal, is in frame of Top Vol
404  return 0; // Have exited world volume
405  }
406  }
407  else
408  {
409  notKnownContained=false;
410  }
411  }
412  else
413  {
414  notKnownContained=false;
415  }
416  } // END while (notKnownContained)
417  //
418  // Search downwards until deepest containing volume found,
419  // blocking fBlockedPhysicalVolume/BlockedReplicaNum
420  //
421  // 3 Cases:
422  //
423  // o Parameterised daughters
424  // =>Must be one G4PVParameterised daughter & voxels
425  // o Positioned daughters & voxels
426  // o Positioned daughters & no voxels
427 
428  noResult = true; // noResult should be renamed to
429  // something like enteredLevel, as that is its meaning.
430  do
431  {
432  // Determine `type' of current mother volume
433  //
434  targetPhysical = fHistory.GetTopVolume();
435  if (!targetPhysical) { break; }
436  targetLogical = targetPhysical->GetLogicalVolume();
437  switch( CharacteriseDaughters(targetLogical) )
438  {
439  case kNormal:
440  if ( targetLogical->GetVoxelHeader() ) // use optimised navigation
441  {
442  noResult = fvoxelNav.LevelLocate(fHistory,
445  globalPoint,
446  pGlobalDirection,
447  considerDirection,
448  localPoint);
449  }
450  else // do not use optimised navigation
451  {
452  noResult = fnormalNav.LevelLocate(fHistory,
455  globalPoint,
456  pGlobalDirection,
457  considerDirection,
458  localPoint);
459  }
460  break;
461  case kReplica:
462  noResult = freplicaNav.LevelLocate(fHistory,
465  globalPoint,
466  pGlobalDirection,
467  considerDirection,
468  localPoint);
469  break;
470  case kParameterised:
471  if( GetDaughtersRegularStructureId(targetLogical) != 1 )
472  {
473  noResult = fparamNav.LevelLocate(fHistory,
476  globalPoint,
477  pGlobalDirection,
478  considerDirection,
479  localPoint);
480  }
481  else // Regular structure
482  {
483  noResult = fregularNav.LevelLocate(fHistory,
486  globalPoint,
487  pGlobalDirection,
488  considerDirection,
489  localPoint);
490  }
491  break;
492  }
493 
494  // LevelLocate returns true if it finds a daughter volume
495  // in which globalPoint is inside (or on the surface).
496 
497  if ( noResult )
498  {
499  noLevelsEntered++;
500 
501  // Entering a daughter after ascending
502  //
503  // The blocked volume is no longer valid - it was for another level
504  //
506  fBlockedReplicaNo = -1;
507 
508  // fEntering should be false -- else blockedVolume is assumed good.
509  // fEnteredDaughter is used for ExitNormal
510  //
511  fEntering = false;
512  fEnteredDaughter = true;
513 
514  if( fExitedMother )
515  {
516  G4VPhysicalVolume* enteredPhysical = fHistory.GetTopVolume();
517  const G4RotationMatrix* mRot = enteredPhysical->GetRotation();
518  if( mRot )
519  {
520  // Go deeper, i.e. move 'down' in the hierarchy
521  // Apply direct rotation, not inverse
522  //
523  fGrandMotherExitNormal *= (*mRot);
525  }
526  }
527 
528 #ifdef G4DEBUG_NAVIGATION
529  if( fVerbose > 2 )
530  {
531  G4VPhysicalVolume* enteredPhysical = fHistory.GetTopVolume();
532  G4cout << "*** G4Navigator::LocateGlobalPointAndSetup() ***" << G4endl;
533  G4cout << " Entering volume: " << enteredPhysical->GetName()
534  << G4endl;
535  }
536 #endif
537  }
538  } while (noResult);
539 
540  fLastLocatedPointLocal = localPoint;
541 
542 #ifdef G4VERBOSE
543  if( fVerbose >= 4 )
544  {
545  G4int oldcoutPrec = G4cout.precision(8);
546  G4String curPhysVol_Name("None");
547  if (targetPhysical) { curPhysVol_Name = targetPhysical->GetName(); }
548  G4cout << " Return value = new volume = " << curPhysVol_Name << G4endl;
549  G4cout << " ----- Upon exiting:" << G4endl;
550  PrintState();
551  if( fVerbose >= 5 )
552  {
553  G4cout << "Upon exiting LocateGlobalPointAndSetup():" << G4endl;
554  G4cout << " History = " << G4endl << fHistory << G4endl << G4endl;
555  }
556  G4cout.precision(oldcoutPrec);
557  }
558 #endif
559 
560  fLocatedOutsideWorld= false;
561 
562  return targetPhysical;
563 }
564 
565 // ********************************************************************
566 // LocateGlobalPointWithinVolume
567 //
568 // -> the state information of this Navigator and its subNavigators
569 // is updated in order to start the next step at pGlobalpoint
570 // -> no check is performed whether pGlobalpoint is inside the
571 // original volume (this must be the case).
572 //
573 // Note: a direction could be added to the arguments, to aid in future
574 // optional checking (via the old code below, flagged by OLD_LOCATE).
575 // [ This would be done only in verbose mode ]
576 // ********************************************************************
577 //
578 void
580 {
581 #ifdef G4DEBUG_NAVIGATION
582  // Check: Either step was not limited by a boundary
583  // or else the full step is no longer being taken
584  assert( !fWasLimitedByGeometry );
585 #endif
586 
589  fChangedGrandMotherRefFrame= false; // Frame for Exit Normal
590 
591 #ifdef G4DEBUG_NAVIGATION
592  if( fVerbose > 2 )
593  {
594  G4cout << "Entering LocateGlobalWithinVolume(): History = " << G4endl;
595  G4cout << fHistory << G4endl;
596  }
597 #endif
598 
599  // For the case of Voxel (or Parameterised) volume the respective
600  // Navigator must be messaged to update its voxel information etc
601 
602  // Update the state of the Sub Navigators
603  // - in particular any voxel information they store/cache
604  //
605  G4VPhysicalVolume* motherPhysical = fHistory.GetTopVolume();
606  G4LogicalVolume* motherLogical = motherPhysical->GetLogicalVolume();
607  G4SmartVoxelHeader* pVoxelHeader = motherLogical->GetVoxelHeader();
608 
610  {
611  switch( CharacteriseDaughters(motherLogical) )
612  {
613  case kNormal:
614  if ( pVoxelHeader )
615  {
617  }
618  break;
619  case kParameterised:
620  if( GetDaughtersRegularStructureId(motherLogical) != 1 )
621  {
622  // Resets state & returns voxel node
623  //
625  }
626  break;
627  case kReplica:
628  G4Exception("G4Navigator::LocateGlobalPointWithinVolume()",
629  "GeomNav0001", FatalException,
630  "Not applicable for replicated volumes.");
631  break;
632  }
633  }
634 
635  // Reset the state variables
636  // - which would have been affected
637  // by the 'equivalent' call to LocateGlobalPointAndSetup
638  // - who's values have been invalidated by the 'move'.
639  //
641  fBlockedReplicaNo = -1;
642  fEntering = false;
643  fEnteredDaughter = false; // Boundary not encountered, did not enter
644  fExiting = false;
645  fExitedMother = false; // Boundary not encountered, did not exit
646 }
647 
648 // ********************************************************************
649 // SetSavedState
650 //
651 // Save the state, in case this is a parasitic call
652 // Save fValidExitNormal, fExitNormal, fExiting, fEntering,
653 // fBlockedPhysicalVolume, fBlockedReplicaNo, fLastStepWasZero;
654 // ********************************************************************
655 //
657 {
658  // Note: the state of dependent objects is not currently saved.
659  // ( This means that the full state is changed by calls between
660  // SetSavedState() and RestoreSavedState();
661 
666 
669 
671 
677 
678  // Even the safety sphere - if you want to change it do it explicitly!
679  //
682 }
683 
684 // ********************************************************************
685 // RestoreSavedState
686 //
687 // Restore the state (in Compute Step), in case this is a parasitic call
688 // ********************************************************************
689 //
691 {
696 
699 
701 
707 
708  // The 'expected' behaviour is to restore these too (fix 2014.05.26)
711 }
712 
713 // ********************************************************************
714 // ComputeStep
715 //
716 // Computes the next geometric Step: intersections with current
717 // mother and `daughter' volumes.
718 //
719 // NOTE:
720 //
721 // Flags on entry:
722 // --------------
723 // fValidExitNormal - Normal of exited volume is valid (convex, not a
724 // coincident boundary)
725 // fExitNormal - Surface normal of exited volume
726 // fExiting - True if have exited solid
727 //
728 // fBlockedPhysicalVolume - Ptr to exited volume (or 0)
729 // fBlockedReplicaNo - Replication no of exited volume
730 // fLastStepWasZero - True if last Step size was zero.
731 //
732 // Flags on exit:
733 // -------------
734 // fValidExitNormal - True if surface normal of exited volume is valid
735 // fExitNormal - Surface normal of exited volume rotated to mothers
736 // reference system
737 // fExiting - True if exiting mother
738 // fEntering - True if entering `daughter' volume (or replica)
739 // fBlockedPhysicalVolume - Ptr to candidate (entered) volume
740 // fBlockedReplicaNo - Replication no of candidate (entered) volume
741 // fLastStepWasZero - True if this Step size was zero.
742 // ********************************************************************
743 //
745  const G4ThreeVector &pDirection,
746  const G4double pCurrentProposedStepLength,
747  G4double &pNewSafety)
748 {
749  G4ThreeVector localDirection = ComputeLocalAxis(pDirection);
751  G4VPhysicalVolume *motherPhysical = fHistory.GetTopVolume();
752  G4LogicalVolume *motherLogical = motherPhysical->GetLogicalVolume();
753 
754  // All state relating to exiting normals must be reset
755  //
757  // Reset value - to erase its memory
759  // Reset - used for local exit normal
760  fGrandMotherExitNormal= G4ThreeVector( 0., 0., 0.);
761  fCalculatedExitNormal = false;
762  // Reset for new step
763 
764  static G4ThreadLocal G4int sNavCScalls=0;
765  sNavCScalls++;
766 
768 
769 #ifdef G4VERBOSE
770  if( fVerbose > 0 )
771  {
772  G4cout << "*** G4Navigator::ComputeStep: ***" << G4endl;
773  G4cout << " Volume = " << motherPhysical->GetName()
774  << " - Proposed step length = " << pCurrentProposedStepLength
775  << G4endl;
776 #ifdef G4DEBUG_NAVIGATION
777  if( fVerbose >= 2 )
778  {
779  G4cout << " Called with the arguments: " << G4endl
780  << " Globalpoint = " << std::setw(25) << pGlobalpoint << G4endl
781  << " Direction = " << std::setw(25) << pDirection << G4endl;
782  if( fVerbose >= 4 )
783  {
784  G4cout << " ---- Upon entering : State" << G4endl;
785  PrintState();
786  }
787  }
788 #endif
789  }
790 #endif
791 
792  G4ThreeVector newLocalPoint = ComputeLocalPoint(pGlobalpoint);
793  if( newLocalPoint != fLastLocatedPointLocal )
794  {
795  // Check whether the relocation is within safety
796  //
797  G4ThreeVector oldLocalPoint = fLastLocatedPointLocal;
798  G4double moveLenSq = (newLocalPoint-oldLocalPoint).mag2();
799 
800  if ( moveLenSq >= kCarTolerance*kCarTolerance )
801  {
802 #ifdef G4VERBOSE
803  ComputeStepLog(pGlobalpoint, moveLenSq);
804 #endif
805  // Relocate the point within the same volume
806  //
807  LocateGlobalPointWithinVolume( pGlobalpoint );
808  fLastTriedStepComputation= true; // Ensure that this is set again !!
809  }
810  }
812  {
813  switch( CharacteriseDaughters(motherLogical) )
814  {
815  case kNormal:
816  if ( motherLogical->GetVoxelHeader() )
817  {
819  localDirection,
820  pCurrentProposedStepLength,
821  pNewSafety,
822  fHistory,
824  fExitNormal,
825  fExiting,
826  fEntering,
829 
830  }
831  else
832  {
833  if( motherPhysical->GetRegularStructureId() == 0 )
834  {
836  localDirection,
837  pCurrentProposedStepLength,
838  pNewSafety,
839  fHistory,
841  fExitNormal,
842  fExiting,
843  fEntering,
846  }
847  else // Regular (non-voxelised) structure
848  {
849  LocateGlobalPointAndSetup( pGlobalpoint, &pDirection, true, true );
850  fLastTriedStepComputation= true; // Ensure that this is set again !!
851  //
852  // if physical process limits the step, the voxel will not be the
853  // one given by ComputeStepSkippingEqualMaterials() and the local
854  // point will be wrongly calculated.
855 
856  // There is a problem: when msc limits the step and the point is
857  // assigned wrongly to phantom in previous step (while it is out
858  // of the container volume). Then LocateGlobalPointAndSetup() has
859  // reset the history topvolume to world.
860  //
862  {
863  G4Exception("G4Navigator::ComputeStep()",
864  "GeomNav1001", JustWarning,
865  "Point is relocated in voxels, while it should be outside!");
867  localDirection,
868  pCurrentProposedStepLength,
869  pNewSafety,
870  fHistory,
872  fExitNormal,
873  fExiting,
874  fEntering,
877  }
878  else
879  {
880  Step = fregularNav.
881  ComputeStepSkippingEqualMaterials(fLastLocatedPointLocal,
882  localDirection,
883  pCurrentProposedStepLength,
884  pNewSafety,
885  fHistory,
887  fExitNormal,
888  fExiting,
889  fEntering,
892  motherPhysical);
893  }
894  }
895  }
896  break;
897  case kParameterised:
898  if( GetDaughtersRegularStructureId(motherLogical) != 1 )
899  {
901  localDirection,
902  pCurrentProposedStepLength,
903  pNewSafety,
904  fHistory,
906  fExitNormal,
907  fExiting,
908  fEntering,
911  }
912  else // Regular structure
913  {
915  localDirection,
916  pCurrentProposedStepLength,
917  pNewSafety,
918  fHistory,
920  fExitNormal,
921  fExiting,
922  fEntering,
925  }
926  break;
927  case kReplica:
928  G4Exception("G4Navigator::ComputeStep()", "GeomNav0001",
929  FatalException, "Not applicable for replicated volumes.");
930  break;
931  }
932  }
933  else
934  {
935  // In the case of a replica, it must handle the exiting
936  // edge/corner problem by itself
937  //
938  G4bool exitingReplica = fExitedMother;
939  G4bool calculatedExitNormal;
940  Step = freplicaNav.ComputeStep(pGlobalpoint,
941  pDirection,
943  localDirection,
944  pCurrentProposedStepLength,
945  pNewSafety,
946  fHistory,
948  calculatedExitNormal,
949  fExitNormal,
950  exitingReplica,
951  fEntering,
954  fExiting= exitingReplica;
955  fCalculatedExitNormal= calculatedExitNormal;
956  }
957 
958  // Remember last safety origin & value.
959  //
960  fPreviousSftOrigin = pGlobalpoint;
961  fPreviousSafety = pNewSafety;
962 
963  // Count zero steps - one can occur due to changing momentum at a boundary
964  // - one, two (or a few) can occur at common edges between
965  // volumes
966  // - more than two is likely a problem in the geometry
967  // description or the Navigation
968 
969  // Rule of thumb: likely at an Edge if two consecutive steps are zero,
970  // because at least two candidate volumes must have been
971  // checked
972  //
973  fLocatedOnEdge = fLastStepWasZero && (Step==0.0);
974  fLastStepWasZero = (Step==0.0);
975  if (fPushed) { fPushed = fLastStepWasZero; }
976 
977  // Handle large number of consecutive zero steps
978  //
979  if ( fLastStepWasZero )
980  {
982 #ifdef G4DEBUG_NAVIGATION
983  if( fNumberZeroSteps > 1 )
984  {
985  G4cout << "G4Navigator::ComputeStep(): another zero step, # "
987  << " at " << pGlobalpoint
988  << " in volume " << motherPhysical->GetName()
989  << " nav-comp-step calls # " << sNavCScalls
990  << G4endl;
991  }
992 #endif
994  {
995  // Act to recover this stuck track. Pushing it along direction
996  //
997  Step += 100*kCarTolerance;
998 #ifdef G4VERBOSE
999  if ((!fPushed) && (fWarnPush))
1000  {
1001  std::ostringstream message;
1002  message << "Track stuck or not moving." << G4endl
1003  << " Track stuck, not moving for "
1004  << fNumberZeroSteps << " steps" << G4endl
1005  << " in volume -" << motherPhysical->GetName()
1006  << "- at point " << pGlobalpoint << G4endl
1007  << " direction: " << pDirection << "." << G4endl
1008  << " Potential geometry or navigation problem !"
1009  << G4endl
1010  << " Trying pushing it of " << Step << " mm ...";
1011  G4Exception("G4Navigator::ComputeStep()", "GeomNav1002",
1012  JustWarning, message, "Potential overlap in geometry!");
1013  }
1014 #endif
1015  fPushed = true;
1016  }
1018  {
1019  // Must kill this stuck track
1020  //
1021  std::ostringstream message;
1022  message << "Stuck Track: potential geometry or navigation problem."
1023  << G4endl
1024  << " Track stuck, not moving for "
1025  << fNumberZeroSteps << " steps" << G4endl
1026  << " in volume -" << motherPhysical->GetName()
1027  << "- at point " << pGlobalpoint << G4endl
1028  << " direction: " << pDirection << ".";
1029  motherPhysical->CheckOverlaps(5000, false);
1030  G4Exception("G4Navigator::ComputeStep()", "GeomNav0003",
1031  EventMustBeAborted, message);
1032  }
1033  }
1034  else
1035  {
1036  if (!fPushed) fNumberZeroSteps = 0;
1037  }
1038 
1039  fEnteredDaughter = fEntering; // I expect to enter a volume in this Step
1041 
1042  fStepEndPoint = pGlobalpoint
1043  + std::min(Step,pCurrentProposedStepLength) * pDirection;
1044  fLastStepEndPointLocal = fLastLocatedPointLocal + Step * localDirection;
1045 
1046  if( fExiting )
1047  {
1048 #ifdef G4DEBUG_NAVIGATION
1049  if( fVerbose > 2 )
1050  {
1051  G4cout << " At G4Nav CompStep End - if(exiting) - fExiting= " << fExiting
1052  << " fValidExitNormal = " << fValidExitNormal << G4endl;
1053  G4cout << " fExitNormal= " << fExitNormal << G4endl;
1054  }
1055 #endif
1056 
1058  {
1060  {
1061  // Convention: fExitNormal is in the 'grand-mother' coordinate system
1062  //
1064  fCalculatedExitNormal= true;
1065  }
1066  else
1067  {
1069  }
1070  }
1071  else
1072  {
1073  // We must calculate the normal anyway (in order to have it if requested)
1074  //
1075  G4ThreeVector finalLocalPoint =
1076  fLastLocatedPointLocal + localDirection*Step;
1077 
1079  {
1080  // Find normal in the 'mother' coordinate system
1081  //
1082  G4ThreeVector exitNormalMotherFrame=
1083  motherLogical->GetSolid()->SurfaceNormal(finalLocalPoint);
1084 
1085  // Transform it to the 'grand-mother' coordinate system
1086  //
1087  const G4RotationMatrix* mRot = motherPhysical->GetRotation();
1088  if( mRot )
1089  {
1091  fGrandMotherExitNormal = (*mRot).inverse() * exitNormalMotherFrame;
1092  }
1093  else
1094  {
1095  fGrandMotherExitNormal = exitNormalMotherFrame;
1096  }
1097 
1098  // Do not set fValidExitNormal -- this signifies
1099  // that the solid is convex!
1100  //
1101  fCalculatedExitNormal= true;
1102  }
1103  else
1104  {
1105  fCalculatedExitNormal = false;
1106  //
1107  // Nothing can be done at this stage currently - to solve this
1108  // Replica Navigation must have calculated the normal for this case
1109  // already.
1110  // Cases: mother is not convex, and exit is at previous replica level
1111 
1112 #ifdef G4DEBUG_NAVIGATION
1114 
1115  desc << "Problem in ComputeStep: Replica Navigation did not provide"
1116  << " valid exit Normal. " << G4endl;
1117  desc << " Do not know how calculate it in this case." << G4endl;
1118  desc << " Location = " << finalLocalPoint << G4endl;
1119  desc << " Volume name = " << motherPhysical->GetName()
1120  << " copy/replica No = " << motherPhysical->GetCopyNo() << G4endl;
1121  G4Exception("G4Navigator::ComputeStep()", "GeomNav0003",
1122  JustWarning, desc, "Normal not available for exiting.");
1123 #endif
1124  }
1125  }
1126 
1127  // Now transform it to the global reference frame !!
1128  //
1130  {
1131  G4int depth= fHistory.GetDepth();
1132  if( depth > 0 )
1133  {
1134  G4AffineTransform GrandMotherToGlobalTransf =
1135  fHistory.GetTransform(depth-1).Inverse();
1137  GrandMotherToGlobalTransf.TransformAxis( fGrandMotherExitNormal );
1138  }
1139  else
1140  {
1142  }
1143  }
1144  else
1145  {
1146  fExitNormalGlobalFrame= G4ThreeVector( 0., 0., 0.);
1147  }
1148  }
1149 
1150  if( (Step == pCurrentProposedStepLength) && (!fExiting) && (!fEntering) )
1151  {
1152  // This if Step is not really limited by the geometry.
1153  // The Navigator is obliged to return "infinity"
1154  //
1155  Step = kInfinity;
1156  }
1157 
1158 #ifdef G4VERBOSE
1159  if( fVerbose > 1 )
1160  {
1161  if( fVerbose >= 4 )
1162  {
1163  G4cout << " ----- Upon exiting :" << G4endl;
1164  PrintState();
1165  }
1166  G4cout << " Returned step= " << Step;
1167  if( fVerbose > 5 ) G4cout << G4endl;
1168  if( Step == kInfinity )
1169  {
1170  G4cout << " Requested step= " << pCurrentProposedStepLength ;
1171  if( fVerbose > 5) G4cout << G4endl;
1172  }
1173  G4cout << " Safety = " << pNewSafety << G4endl;
1174  }
1175 #endif
1176 
1177  return Step;
1178 }
1179 
1180 // ********************************************************************
1181 // CheckNextStep
1182 //
1183 // Compute the step without altering the navigator state
1184 // ********************************************************************
1185 //
1187  const G4ThreeVector& pDirection,
1188  const G4double pCurrentProposedStepLength,
1189  G4double& pNewSafety)
1190 {
1191  G4double step;
1192 
1193  // Save the state, for this parasitic call
1194  //
1195  SetSavedState();
1196 
1197  step = ComputeStep ( pGlobalpoint,
1198  pDirection,
1199  pCurrentProposedStepLength,
1200  pNewSafety );
1201 
1202  // It is a parasitic call, so attempt to restore the key parts of the state
1203  //
1204  RestoreSavedState();
1205  // NOTE: the state of the current subnavigator is NOT restored.
1206  // ***> TODO: restore subnavigator state
1207  // if( last_located) Need Position of last location
1208  // if( last_computed step) Need Endposition of last step
1209 
1210  return step;
1211 }
1212 
1213 // ********************************************************************
1214 // ResetState
1215 //
1216 // Resets stack and minimum of navigator state `machine'
1217 // ********************************************************************
1218 //
1220 {
1221  fWasLimitedByGeometry = false;
1222  fEntering = false;
1223  fExiting = false;
1224  fLocatedOnEdge = false;
1225  fLastStepWasZero = false;
1226  fEnteredDaughter = false;
1227  fExitedMother = false;
1228  fPushed = false;
1229 
1230  fValidExitNormal = false;
1232  fCalculatedExitNormal = false;
1233 
1234  fExitNormal = G4ThreeVector(0,0,0);
1237 
1239  fPreviousSafety = 0.0;
1240 
1241  fNumberZeroSteps = 0;
1242 
1244  fBlockedReplicaNo = -1;
1245 
1247  fLocatedOutsideWorld = false;
1248 }
1249 
1250 // ********************************************************************
1251 // SetupHierarchy
1252 //
1253 // Renavigates & resets hierarchy described by current history
1254 // o Reset volumes
1255 // o Recompute transforms and/or solids of replicated/parameterised volumes
1256 // ********************************************************************
1257 //
1259 {
1260  G4int i;
1261  const G4int cdepth = fHistory.GetDepth();
1262  G4VPhysicalVolume *current;
1263  G4VSolid *pSolid;
1264  G4VPVParameterisation *pParam;
1265 
1266  for ( i=1; i<=cdepth; i++ )
1267  {
1268  current = fHistory.GetVolume(i);
1269  switch ( fHistory.GetVolumeType(i) )
1270  {
1271  case kNormal:
1272  break;
1273  case kReplica:
1275  break;
1276  case kParameterised:
1277  G4int replicaNo;
1278  pParam = current->GetParameterisation();
1279  replicaNo = fHistory.GetReplicaNo(i);
1280  pSolid = pParam->ComputeSolid(replicaNo, current);
1281 
1282  // Set up dimensions & transform in solid/physical volume
1283  //
1284  pSolid->ComputeDimensions(pParam, replicaNo, current);
1285  pParam->ComputeTransformation(replicaNo, current);
1286 
1287  G4TouchableHistory *pTouchable= 0;
1288  if( pParam->IsNested() )
1289  {
1290  pTouchable= new G4TouchableHistory( fHistory );
1291  pTouchable->MoveUpHistory(); // Move up to the parent level
1292  // Adequate only if Nested at the Branch level (last)
1293  // To extend to other cases:
1294  // pTouchable->MoveUpHistory(cdepth-i-1);
1295  // Move to the parent level of *Current* level
1296  // Could replace this line and constructor with a revised
1297  // c-tor for History(levels to drop)
1298  }
1299  // Set up the correct solid and material in Logical Volume
1300  //
1301  G4LogicalVolume *pLogical = current->GetLogicalVolume();
1302  pLogical->SetSolid( pSolid );
1303  pLogical->UpdateMaterial( pParam ->
1304  ComputeMaterial(replicaNo, current, pTouchable) );
1305  delete pTouchable;
1306  break;
1307  }
1308  }
1309 }
1310 
1311 // ********************************************************************
1312 // GetLocalExitNormal
1313 //
1314 // Obtains the Normal vector to a surface (in local coordinates)
1315 // pointing out of previous volume and into current volume
1316 // ********************************************************************
1317 //
1319 {
1320  G4ThreeVector ExitNormal(0.,0.,0.);
1321  G4VSolid *currentSolid=0;
1322  G4LogicalVolume *candidateLogical;
1324  {
1325  // use fLastLocatedPointLocal and next candidate volume
1326  //
1327  G4ThreeVector nextSolidExitNormal(0.,0.,0.);
1328 
1329  if( fEntering && (fBlockedPhysicalVolume!=0) )
1330  {
1331  candidateLogical= fBlockedPhysicalVolume->GetLogicalVolume();
1332  if( candidateLogical )
1333  {
1334  // fLastStepEndPointLocal is in the coordinates of the mother
1335  // we need it in the daughter's coordinate system.
1336 
1337  // The following code should also work in case of Replica
1338  {
1339  // First transform fLastLocatedPointLocal to the new daughter
1340  // coordinates
1341  //
1342  G4AffineTransform MotherToDaughterTransform=
1346  G4ThreeVector daughterPointOwnLocal=
1347  MotherToDaughterTransform.TransformPoint( fLastStepEndPointLocal );
1348 
1349  // OK if it is a parameterised volume
1350  //
1351  EInside inSideIt;
1352  G4bool onSurface;
1353  G4double safety= -1.0;
1354  currentSolid= candidateLogical->GetSolid();
1355  inSideIt = currentSolid->Inside(daughterPointOwnLocal);
1356  onSurface = (inSideIt == kSurface);
1357  if( ! onSurface )
1358  {
1359  if( inSideIt == kOutside )
1360  {
1361  safety = (currentSolid->DistanceToIn(daughterPointOwnLocal));
1362  onSurface = safety < 100.0 * kCarTolerance;
1363  }
1364  else if (inSideIt == kInside )
1365  {
1366  safety = (currentSolid->DistanceToOut(daughterPointOwnLocal));
1367  onSurface = safety < 100.0 * kCarTolerance;
1368  }
1369  }
1370 
1371  if( onSurface )
1372  {
1373  nextSolidExitNormal =
1374  currentSolid->SurfaceNormal(daughterPointOwnLocal);
1375 
1376  // Entering the solid ==> opposite
1377  //
1378  ExitNormal = -nextSolidExitNormal;
1379  fCalculatedExitNormal= true;
1380  }
1381  else
1382  {
1383 #ifdef G4VERBOSE
1384  if(( fVerbose == 1 ) && ( fCheck ))
1385  {
1386  std::ostringstream message;
1387  message << "Point not on surface ! " << G4endl
1388  << " Point = "
1389  << daughterPointOwnLocal << G4endl
1390  << " Physical volume = "
1392  << " Logical volume = "
1393  << candidateLogical->GetName() << G4endl
1394  << " Solid = " << currentSolid->GetName()
1395  << " Type = "
1396  << currentSolid->GetEntityType() << G4endl
1397  << *currentSolid << G4endl;
1398  if( inSideIt == kOutside )
1399  {
1400  message << "Point is Outside. " << G4endl
1401  << " Safety (from outside) = " << safety << G4endl;
1402  }
1403  else // if( inSideIt == kInside )
1404  {
1405  message << "Point is Inside. " << G4endl
1406  << " Safety (from inside) = " << safety << G4endl;
1407  }
1408  G4Exception("G4Navigator::GetLocalExitNormal()", "GeomNav1001",
1409  JustWarning, message);
1410  }
1411 #endif
1412  }
1413  *valid = onSurface; // was =true;
1414  }
1415  }
1416  }
1417  else if ( fExiting )
1418  {
1419  ExitNormal = fGrandMotherExitNormal;
1420  *valid = true;
1421  fCalculatedExitNormal= true; // Should be true already
1422  }
1423  else // i.e. ( fBlockedPhysicalVolume == 0 )
1424  {
1425  *valid = false;
1426  G4Exception("G4Navigator::GetLocalExitNormal()",
1427  "GeomNav0003", JustWarning,
1428  "Incorrect call to GetLocalSurfaceNormal." );
1429  }
1430  }
1431  else // ( ! fLastTriedStepComputation ) ie. last call was to Locate
1432  {
1433  if ( EnteredDaughterVolume() )
1434  {
1435  G4VSolid* daughterSolid =fHistory.GetTopVolume()->GetLogicalVolume()
1436  ->GetSolid();
1437  ExitNormal= -(daughterSolid->SurfaceNormal(fLastLocatedPointLocal));
1438  if( std::fabs(ExitNormal.mag2()-1.0 ) > CLHEP::perMillion )
1439  {
1441  desc << " Parameters of solid: " << *daughterSolid
1442  << " Point for surface = " << fLastLocatedPointLocal << std::endl;
1443  G4Exception("G4Navigator::GetLocalExitNormal()",
1444  "GeomNav0003", FatalException, desc,
1445  "Surface Normal returned by Solid is not a Unit Vector." );
1446  }
1447  fCalculatedExitNormal= true;
1448  *valid = true;
1449  }
1450  else
1451  {
1452  if( fExitedMother )
1453  {
1454  ExitNormal = fGrandMotherExitNormal;
1455  *valid = true;
1456  fCalculatedExitNormal= true;
1457  }
1458  else // We are not at a boundary. ExitNormal remains (0,0,0)
1459  {
1460  *valid = false;
1461  fCalculatedExitNormal= false;
1462  G4ExceptionDescription message;
1463  message << "Function called when *NOT* at a Boundary." << G4endl;
1464  message << "Exit Normal not calculated." << G4endl;
1465  G4Exception("G4Navigator::GetLocalExitNormal()",
1466  "GeomNav0003", JustWarning, message);
1467  }
1468  }
1469  }
1470  return ExitNormal;
1471 }
1472 
1473 // ********************************************************************
1474 // GetMotherToDaughterTransform
1475 //
1476 // Obtains the mother to daughter affine transformation
1477 // ********************************************************************
1478 //
1481  G4int enteringReplicaNo,
1482  EVolume enteringVolumeType )
1483 {
1484  switch (enteringVolumeType)
1485  {
1486  case kNormal: // Nothing is needed to prepare the transformation
1487  break; // It is stored already in the physical volume (placement)
1488  case kReplica: // Sets the transform in the Replica - tbc
1489  G4Exception("G4Navigator::GetMotherToDaughterTransform()",
1490  "GeomNav0001", FatalException,
1491  "Method NOT Implemented yet for replica volumes.");
1492  break;
1493  case kParameterised:
1494  if( pEnteringPhysVol->GetRegularStructureId() == 0 )
1495  {
1496  G4VPVParameterisation *pParam =
1497  pEnteringPhysVol->GetParameterisation();
1498  G4VSolid* pSolid =
1499  pParam->ComputeSolid(enteringReplicaNo, pEnteringPhysVol);
1500  pSolid->ComputeDimensions(pParam, enteringReplicaNo, pEnteringPhysVol);
1501 
1502  // Sets the transform in the Parameterisation
1503  //
1504  pParam->ComputeTransformation(enteringReplicaNo, pEnteringPhysVol);
1505 
1506  // Set the correct solid and material in Logical Volume
1507  //
1508  G4LogicalVolume* pLogical = pEnteringPhysVol->GetLogicalVolume();
1509  pLogical->SetSolid( pSolid );
1510  }
1511  break;
1512  }
1513  return G4AffineTransform(pEnteringPhysVol->GetRotation(),
1514  pEnteringPhysVol->GetTranslation()).Invert();
1515 }
1516 
1517 
1518 // ********************************************************************
1519 // GetLocalExitNormalAndCheck
1520 //
1521 // Obtains the Normal vector to a surface (in local coordinates)
1522 // pointing out of previous volume and into current volume, and
1523 // checks the current point against expected 'local' value.
1524 // ********************************************************************
1525 //
1528 #ifdef G4DEBUG_NAVIGATION
1529  const G4ThreeVector& ExpectedBoundaryPointGlobal,
1530 #else
1531  const G4ThreeVector&,
1532 #endif
1533  G4bool* pValid)
1534 {
1535 #ifdef G4DEBUG_NAVIGATION
1536  // Check Current point against expected 'local' value
1537  //
1539  {
1540  G4ThreeVector ExpectedBoundaryPointLocal;
1541 
1542  const G4AffineTransform& GlobalToLocal= GetGlobalToLocalTransform();
1543  ExpectedBoundaryPointLocal =
1544  GlobalToLocal.TransformPoint( ExpectedBoundaryPointGlobal );
1545 
1546  // Add here: Comparison against expected position,
1547  // i.e. the endpoint of ComputeStep
1548  }
1549 #endif
1550 
1551  return GetLocalExitNormal( pValid);
1552 }
1553 
1554 
1555 // ********************************************************************
1556 // GetGlobalExitNormal
1557 //
1558 // Obtains the Normal vector to a surface (in global coordinates)
1559 // pointing out of previous volume and into current volume
1560 // ********************************************************************
1561 //
1564  G4bool* pNormalCalculated)
1565 {
1566  G4bool validNormal;
1567  G4ThreeVector localNormal, globalNormal;
1568  G4bool usingStored;
1569 
1570  usingStored=
1572  ( ( fLastTriedStepComputation && fExiting) // Just calculated it
1573  || // No locate in between
1575  && (IntersectPointGlobal-fStepEndPoint).mag2()
1576  < (10.0*kCarTolerance*kCarTolerance)
1577  ) // Calculated it 'just' before & then called locate
1578  // but it did not move position
1579  );
1580 
1581  if( usingStored )
1582  {
1583  // This was computed in last call to ComputeStep
1584  // and only if it arrived at boundary
1585  //
1586  globalNormal = fExitNormalGlobalFrame;
1587  G4double normMag2 = globalNormal.mag2();
1588  if( std::fabs ( normMag2 - 1.0 ) < perMillion ) // Value is good
1589  {
1590  *pNormalCalculated = true; // ComputeStep always computes it if Exiting
1591  // (fExiting==true)
1592  }
1593  else
1594  {
1595  G4ExceptionDescription message;
1596  message.precision(10);
1597  message << " WARNING> Expected normal-global-frame to be valid, "
1598  << " i.e. a unit vector!" << G4endl
1599  << " - but |normal| = " << std::sqrt(normMag2)
1600  << " - and |normal|^2 = " << normMag2 << G4endl
1601  << " which differs from 1.0 by " << normMag2 - 1.0 << G4endl
1602  << " n = " << fExitNormalGlobalFrame << G4endl;
1603  message << "============================================================"
1604  << G4endl;
1605  G4int oldVerbose = fVerbose;
1606  fVerbose=4;
1607  message << " State of Navigator: " << G4endl;
1608  message << *this << G4endl;
1609  fVerbose = oldVerbose;
1610  message << "============================================================"
1611  << G4endl;
1612 
1613  G4Exception("G4Navigator::GetGlobalExitNormal()",
1614  "GeomNav0003",JustWarning, message,
1615  "Value obtained from stored global-normal is not a unit vector.");
1616 
1617  // (Re)Compute it now -- as either it was not computed, or it is wrong.
1618  //
1619  localNormal = GetLocalExitNormalAndCheck(IntersectPointGlobal,
1620  &validNormal);
1621  *pNormalCalculated = fCalculatedExitNormal;
1622 
1623  G4AffineTransform localToGlobal = GetLocalToGlobalTransform();
1624  globalNormal = localToGlobal.TransformAxis( localNormal );
1625  }
1626  }
1627  else
1628  {
1629  localNormal = GetLocalExitNormalAndCheck(IntersectPointGlobal,&validNormal);
1630  *pNormalCalculated = fCalculatedExitNormal;
1631 
1632 #ifdef G4DEBUG_NAVIGATION
1633  usingStored= false;
1634 
1635  if( (!validNormal) && !fCalculatedExitNormal)
1636  {
1638  edN << " Calculated = " << fCalculatedExitNormal << G4endl;
1639  edN << " Entering= " << fEntering << G4endl;
1640  G4int oldVerbose= this->GetVerboseLevel();
1641  this->SetVerboseLevel(4);
1642  edN << " State of Navigator: " << G4endl;
1643  edN << *this << G4endl;
1644  this->SetVerboseLevel( oldVerbose );
1645 
1646  G4Exception("G4Navigator::GetGlobalExitNormal()",
1647  "GeomNav0003", JustWarning, edN,
1648  "LocalExitNormalAndCheck() did not calculate Normal.");
1649  }
1650 #endif
1651 
1652  G4double localMag2= localNormal.mag2();
1653  if( validNormal && (std::fabs(localMag2-1.0)) > CLHEP::perMillion )
1654  {
1656 
1657  edN << "G4Navigator::GetGlobalExitNormal: "
1658  << " Using Local Normal - from call to GetLocalExitNormalAndCheck. "
1659  << G4endl
1660  << " Local Exit Normal : " << " || = " << std::sqrt(localMag2)
1661  << " vec = " << localNormal << G4endl
1662  << " Global Exit Normal : " << " || = " << globalNormal.mag()
1663  << " vec = " << globalNormal << G4endl;
1664  edN << " Calculated It = " << fCalculatedExitNormal << G4endl;
1665 
1666  G4Exception("G4Navigator::GetGlobalExitNormal()",
1667  "GeomNav0003",JustWarning, edN,
1668  "Value obtained from new local *solid* is incorrect.");
1669  localNormal = localNormal.unit(); // Should we correct it ??
1670  }
1671  G4AffineTransform localToGlobal = GetLocalToGlobalTransform();
1672  globalNormal = localToGlobal.TransformAxis( localNormal );
1673  }
1674 
1675 #ifdef G4DEBUG_NAVIGATION
1676  if( usingStored )
1677  {
1678  G4ThreeVector globalNormAgn;
1679 
1680  localNormal= GetLocalExitNormalAndCheck(IntersectPointGlobal, &validNormal);
1681 
1682  G4AffineTransform localToGlobal = GetLocalToGlobalTransform();
1683  globalNormAgn = localToGlobal.TransformAxis( localNormal );
1684 
1685  // Check the value computed against fExitNormalGlobalFrame
1686  G4ThreeVector diffNorm = globalNormAgn - fExitNormalGlobalFrame;
1687  if( diffNorm.mag2() > perMillion*CLHEP::perMillion)
1688  {
1689  G4ExceptionDescription edDfn;
1690  edDfn << "Found difference in normals in case of exiting mother "
1691  << "- when Get is called after ComputingStep " << G4endl;
1692  edDfn << " Magnitude of diff = " << diffNorm.mag() << G4endl;
1693  edDfn << " Normal stored (Global) = " << fExitNormalGlobalFrame
1694  << G4endl;
1695  edDfn << " Global Computed from Local = " << globalNormAgn << G4endl;
1696  G4Exception("G4Navigator::GetGlobalExitNormal()", "GeomNav0003",
1697  JustWarning, edDfn);
1698  }
1699  }
1700 #endif
1701 
1702  return globalNormal;
1703 }
1704 
1705 // To make the new Voxel Safety the default, uncomment the next line
1706 #define G4NEW_SAFETY 1
1707 
1708 // ********************************************************************
1709 // ComputeSafety
1710 //
1711 // It assumes that it will be
1712 // i) called at the Point in the same volume as the EndPoint of the
1713 // ComputeStep.
1714 // ii) after (or at the end of) ComputeStep OR after the relocation.
1715 // ********************************************************************
1716 //
1718  const G4double pMaxLength,
1719  const G4bool keepState)
1720 {
1721  G4double newSafety = 0.0;
1722 
1723 #ifdef G4DEBUG_NAVIGATION
1724  G4int oldcoutPrec = G4cout.precision(8);
1725  if( fVerbose > 0 )
1726  {
1727  G4cout << "*** G4Navigator::ComputeSafety: ***" << G4endl
1728  << " Called at point: " << pGlobalpoint << G4endl;
1729 
1730  G4VPhysicalVolume *motherPhysical = fHistory.GetTopVolume();
1731  G4cout << " Volume = " << motherPhysical->GetName()
1732  << " - Maximum length = " << pMaxLength << G4endl;
1733  if( fVerbose >= 4 )
1734  {
1735  G4cout << " ----- Upon entering Compute Safety:" << G4endl;
1736  PrintState();
1737  }
1738  }
1739 #endif
1740 
1741  G4double distEndpointSq = (pGlobalpoint-fStepEndPoint).mag2();
1742  G4bool stayedOnEndpoint = distEndpointSq < kCarTolerance*kCarTolerance;
1743  G4bool endpointOnSurface = fEnteredDaughter || fExitedMother;
1744 
1745  if( endpointOnSurface && stayedOnEndpoint )
1746  {
1747 #ifdef G4DEBUG_NAVIGATION
1748  if( fVerbose >= 2 )
1749  {
1750  G4cout << " G4Navigator::ComputeSafety() finds that point - "
1751  << pGlobalpoint << " - is on surface " << G4endl;
1752  if( fEnteredDaughter ) { G4cout << " entered new daughter volume"; }
1753  if( fExitedMother ) { G4cout << " and exited previous volume."; }
1754  G4cout << G4endl;
1755  G4cout << " EndPoint was = " << fStepEndPoint << G4endl;
1756  }
1757 #endif
1758  newSafety = 0.0;
1759  // return newSafety;
1760  }
1761  else // if( !(endpointOnSurface && stayedOnEndpoint) )
1762  {
1763  if (keepState) { SetSavedState(); }
1764 
1765  // Pseudo-relocate to this point (updates voxel information only)
1766  //
1767  LocateGlobalPointWithinVolume( pGlobalpoint );
1768  // --->> DANGER: Side effects on sub-navigator voxel information <<---
1769  // Could be replaced again by 'granular' calls to sub-navigator
1770  // locates (similar side-effects, but faster.
1771  // Solutions:
1772  // 1) Re-locate (to where?)
1773  // 2) Insure that the methods using (G4ComputeStep?)
1774  // does a relocation (if information is disturbed only ?)
1775 
1776 #ifdef G4DEBUG_NAVIGATION
1777  if( fVerbose >= 2 )
1778  {
1779  G4cout << " G4Navigator::ComputeSafety() relocates-in-volume to point: "
1780  << pGlobalpoint << G4endl;
1781  }
1782 #endif
1783  G4VPhysicalVolume *motherPhysical = fHistory.GetTopVolume();
1784  G4LogicalVolume *motherLogical = motherPhysical->GetLogicalVolume();
1785  G4SmartVoxelHeader* pVoxelHeader = motherLogical->GetVoxelHeader();
1786  G4ThreeVector localPoint = ComputeLocalPoint(pGlobalpoint);
1787 
1789  {
1790  switch(CharacteriseDaughters(motherLogical))
1791  {
1792  case kNormal:
1793  if ( pVoxelHeader )
1794  {
1795 #ifdef G4NEW_SAFETY
1796  G4double safetyTwo = fpVoxelSafety->ComputeSafety(localPoint,
1797  *motherPhysical, pMaxLength);
1798  newSafety= safetyTwo; // Faster and best available
1799 #else
1800  G4double safetyOldVoxel;
1801  safetyOldVoxel =
1802  fvoxelNav.ComputeSafety(localPoint,fHistory,pMaxLength);
1803  newSafety= safetyOldVoxel;
1804 #endif
1805  }
1806  else
1807  {
1808  newSafety=fnormalNav.ComputeSafety(localPoint,fHistory,pMaxLength);
1809  }
1810  break;
1811  case kParameterised:
1812  if( GetDaughtersRegularStructureId(motherLogical) != 1 )
1813  {
1814  newSafety=fparamNav.ComputeSafety(localPoint,fHistory,pMaxLength);
1815  }
1816  else // Regular structure
1817  {
1818  newSafety=fregularNav.ComputeSafety(localPoint,fHistory,pMaxLength);
1819  }
1820  break;
1821  case kReplica:
1822  G4Exception("G4Navigator::ComputeSafety()", "GeomNav0001",
1823  FatalException, "Not applicable for replicated volumes.");
1824  break;
1825  }
1826  }
1827  else
1828  {
1829  newSafety = freplicaNav.ComputeSafety(pGlobalpoint, localPoint,
1830  fHistory, pMaxLength);
1831  }
1832 
1833  if (keepState)
1834  {
1836  // This now overwrites the values of the Safety 'sphere' (correction)
1837  }
1838 
1839  // Remember last safety origin & value
1840  //
1841  // We overwrite the Safety 'sphere' - keeping old behaviour
1842  fPreviousSftOrigin = pGlobalpoint;
1843  fPreviousSafety = newSafety;
1844  }
1845 
1846 #ifdef G4DEBUG_NAVIGATION
1847  if( fVerbose > 1 )
1848  {
1849  G4cout << " ---- Exiting ComputeSafety " << G4endl;
1850  if( fVerbose > 2 ) { PrintState(); }
1851  G4cout << " Returned value of Safety = " << newSafety << G4endl;
1852  }
1853  G4cout.precision(oldcoutPrec);
1854 #endif
1855 
1856  return newSafety;
1857 }
1858 
1859 
1860 // ********************************************************************
1861 // RecheckDistanceToCurrentBoundary
1862 //
1863 // Trial method for checking potential displacement for MS
1864 // Check position aDisplacedGlobalpoint, to see whether it is in the
1865 // current volume (mother).
1866 // If in mother, check distance to boundary along aNewDirection.
1867 // If in entering daughter, check distance back to boundary.
1868 // NOTE:
1869 // Can be called only after ComputeStep() is called - before ReLocation
1870 // Deals only with current volume (and potentially entered)
1871 // Caveats
1872 // First VERSION: Does not consider daughter volumes if it remained in mother
1873 // neither for checking whether it has exited current (mother) volume,
1874 // nor for determining distance to exit this (mother) volume.
1875 // ********************************************************************
1876 //
1878  const G4ThreeVector &aDisplacedGlobalPoint,
1879  const G4ThreeVector &aNewDirection,
1880  const G4double ProposedMove,
1881  G4double *prDistance,
1882  G4double *prNewSafety) const
1883 {
1884  G4ThreeVector localPosition = ComputeLocalPoint(aDisplacedGlobalPoint);
1885  G4ThreeVector localDirection = ComputeLocalAxis(aNewDirection);
1886  // G4double Step = kInfinity;
1887 
1888  G4bool validExitNormal;
1889  G4ThreeVector exitNormal;
1890  // Check against mother solid
1891  G4VPhysicalVolume *motherPhysical = fHistory.GetTopVolume();
1892  G4LogicalVolume *motherLogical = motherPhysical->GetLogicalVolume();
1893 
1894 #ifdef CHECK_ORDER_OF_METHODS
1896  {
1897  G4Exception("G4Navigator::RecheckDistanceToCurrentBoundary()",
1898  "GeomNav0001", FatalException,
1899  "Method must be called after ComputeStep(), before call to LocateMethod.");
1900  }
1901 #endif
1902 
1903  EInside locatedDaug; // = kUndefined;
1904  G4double daughterStep= DBL_MAX;
1905  G4double daughterSafety= DBL_MAX;
1906 
1907  if( fEnteredDaughter )
1908  {
1909  if( motherLogical->CharacteriseDaughters() ==kReplica ) { return false; }
1910 
1911  // Track arrived at boundary of a daughter volume at
1912  // the last call of ComputeStep().
1913  // In case the proposed displaced point is inside this daughter,
1914  // it must backtrack at least to the entry point.
1915  // NOTE: No check is made against other daughter volumes. It is
1916  // assumed that the proposed displacement is small enough that
1917  // this is not needed.
1918 
1919  // Must check boundary of current daughter
1920  //
1921  G4VPhysicalVolume *candPhysical= fBlockedPhysicalVolume;
1922  G4LogicalVolume *candLogical= candPhysical->GetLogicalVolume();
1923  G4VSolid *candSolid= candLogical->GetSolid();
1924 
1925  G4AffineTransform nextLevelTrf(candPhysical->GetRotation(),
1926  candPhysical->GetTranslation());
1927 
1928  G4ThreeVector dgPosition= nextLevelTrf.TransformPoint(localPosition);
1929  G4ThreeVector dgDirection= nextLevelTrf.TransformAxis(localDirection);
1930  locatedDaug = candSolid->Inside(dgPosition);
1931 
1932  if( locatedDaug == kInside )
1933  {
1934  // Reverse direction - and find first exit. ( Is it valid?)
1935  // Must backtrack
1936  G4double distanceBackOut =
1937  candSolid->DistanceToOut(dgPosition,
1938  - dgDirection, // Reverse direction
1939  true, &validExitNormal, &exitNormal);
1940  daughterStep= - distanceBackOut;
1941  // No check is made whether the particle could have arrived at
1942  // at this location without encountering another volume or
1943  // a different psurface of the current volume
1944  if( prNewSafety )
1945  {
1946  daughterSafety= candSolid->DistanceToOut(dgPosition);
1947  }
1948  }
1949  else
1950  {
1951  if( locatedDaug == kOutside )
1952  {
1953  // See whether it still intersects it
1954  //
1955  daughterStep= candSolid->DistanceToIn(dgPosition,
1956  dgDirection);
1957  if( prNewSafety )
1958  {
1959  daughterSafety= candSolid->DistanceToIn(dgPosition);
1960  }
1961  }
1962  else
1963  {
1964  // The point remains on the surface of candidate solid
1965  //
1966  daughterStep= 0.0;
1967  daughterSafety= 0.0;
1968  }
1969  }
1970 
1971  // If trial point is in daughter (or on its surface) we have the
1972  // answer, the rest is not relevant
1973  //
1974  if( locatedDaug != kOutside )
1975  {
1976  *prDistance= daughterStep;
1977  if( prNewSafety ) { *prNewSafety= daughterSafety; }
1978  return true;
1979  }
1980  // If ever extended, so that some type of mother cut daughter,
1981  // this would change
1982  }
1983 
1984  G4VSolid *motherSolid= motherLogical->GetSolid();
1985 
1986  G4double motherStep= DBL_MAX, motherSafety= DBL_MAX;
1987 
1988  // Check distance to boundary of mother
1989  //
1991  {
1992  EInside locatedMoth = motherSolid->Inside(localPosition);
1993 
1994  if( locatedMoth == kInside )
1995  {
1996  motherSafety= motherSolid->DistanceToOut(localPosition);
1997  if( ProposedMove >= motherSafety )
1998  {
1999  motherStep= motherSolid->DistanceToOut(localPosition,
2000  localDirection,
2001  true, &validExitNormal, &exitNormal);
2002  }
2003  else
2004  {
2005  motherStep= ProposedMove;
2006  }
2007  }
2008  else if( locatedMoth == kOutside)
2009  {
2010  motherSafety= motherSolid->DistanceToIn(localPosition);
2011  if( ProposedMove >= motherSafety )
2012  {
2013  motherStep= - motherSolid->DistanceToIn(localPosition,
2014  -localDirection);
2015  }
2016  }
2017  else
2018  {
2019  motherSafety= 0.0;
2020  *prDistance= 0.0; // On surface - no move // motherStep;
2021  if( prNewSafety ) { *prNewSafety= motherSafety; }
2022  return false;
2023  }
2024  }
2025  else
2026  {
2027  return false;
2028  }
2029 
2030  *prDistance= std::min( motherStep, daughterStep );
2031  if( prNewSafety )
2032  {
2033  *prNewSafety= std::min( motherSafety, daughterSafety );
2034  }
2035  return true;
2036 }
2037 
2038 
2039 // ********************************************************************
2040 // CreateTouchableHistoryHandle
2041 // ********************************************************************
2042 //
2044 {
2046 }
2047 
2048 // ********************************************************************
2049 // PrintState
2050 // ********************************************************************
2051 //
2053 {
2054  G4int oldcoutPrec = G4cout.precision(4);
2055  if( fVerbose >= 4 )
2056  {
2057  G4cout << "The current state of G4Navigator is: " << G4endl;
2058  G4cout << " ValidExitNormal= " << fValidExitNormal // << G4endl
2059  << " ExitNormal = " << fExitNormal // << G4endl
2060  << " Exiting = " << fExiting // << G4endl
2061  << " Entering = " << fEntering // << G4endl
2062  << " BlockedPhysicalVolume= " ;
2063  if (fBlockedPhysicalVolume==0)
2064  G4cout << "None";
2065  else
2067  G4cout << G4endl
2068  << " BlockedReplicaNo = " << fBlockedReplicaNo // << G4endl
2069  << " LastStepWasZero = " << fLastStepWasZero // << G4endl
2070  << G4endl;
2071  }
2072  if( ( 1 < fVerbose) && (fVerbose < 4) )
2073  {
2074  G4cout << G4endl; // Make sure to line up
2075  G4cout << std::setw(30) << " ExitNormal " << " "
2076  << std::setw( 5) << " Valid " << " "
2077  << std::setw( 9) << " Exiting " << " "
2078  << std::setw( 9) << " Entering" << " "
2079  << std::setw(15) << " Blocked:Volume " << " "
2080  << std::setw( 9) << " ReplicaNo" << " "
2081  << std::setw( 8) << " LastStepZero " << " "
2082  << G4endl;
2083  G4cout << "( " << std::setw(7) << fExitNormal.x()
2084  << ", " << std::setw(7) << fExitNormal.y()
2085  << ", " << std::setw(7) << fExitNormal.z() << " ) "
2086  << std::setw( 5) << fValidExitNormal << " "
2087  << std::setw( 9) << fExiting << " "
2088  << std::setw( 9) << fEntering << " ";
2089  if ( fBlockedPhysicalVolume==0 )
2090  G4cout << std::setw(15) << "None";
2091  else
2092  G4cout << std::setw(15)<< fBlockedPhysicalVolume->GetName();
2093  G4cout << std::setw( 9) << fBlockedReplicaNo << " "
2094  << std::setw( 8) << fLastStepWasZero << " "
2095  << G4endl;
2096  }
2097  if( fVerbose > 2 )
2098  {
2099  G4cout.precision(8);
2100  G4cout << " Current Localpoint = " << fLastLocatedPointLocal << G4endl;
2101  G4cout << " PreviousSftOrigin = " << fPreviousSftOrigin << G4endl;
2102  G4cout << " PreviousSafety = " << fPreviousSafety << G4endl;
2103  }
2104  G4cout.precision(oldcoutPrec);
2105 }
2106 
2107 // ********************************************************************
2108 // ComputeStepLog
2109 // ********************************************************************
2110 //
2112  G4double moveLenSq) const
2113 {
2114  // The following checks only make sense if the move is larger
2115  // than the tolerance.
2116 
2117  static const G4double fAccuracyForWarning = kCarTolerance,
2118  fAccuracyForException = 1000*kCarTolerance;
2119 
2120  G4ThreeVector OriginalGlobalpoint = fHistory.GetTopTransform().Inverse().
2121  TransformPoint(fLastLocatedPointLocal);
2122 
2123  G4double shiftOriginSafSq = (fPreviousSftOrigin-pGlobalpoint).mag2();
2124 
2125  // Check that the starting point of this step is
2126  // within the isotropic safety sphere of the last point
2127  // to a accuracy/precision given by fAccuracyForWarning.
2128  // If so give warning.
2129  // If it fails by more than fAccuracyForException exit with error.
2130  //
2131  if( shiftOriginSafSq >= sqr(fPreviousSafety) )
2132  {
2133  G4double shiftOrigin = std::sqrt(shiftOriginSafSq);
2134  G4double diffShiftSaf = shiftOrigin - fPreviousSafety;
2135 
2136  if( diffShiftSaf > fAccuracyForWarning )
2137  {
2138  G4int oldcoutPrec= G4cout.precision(8);
2139  G4int oldcerrPrec= G4cerr.precision(10);
2140  std::ostringstream message, suggestion;
2141  message << "Accuracy error or slightly inaccurate position shift."
2142  << G4endl
2143  << " The Step's starting point has moved "
2144  << std::sqrt(moveLenSq)/mm << " mm " << G4endl
2145  << " since the last call to a Locate method." << G4endl
2146  << " This has resulted in moving "
2147  << shiftOrigin/mm << " mm "
2148  << " from the last point at which the safety "
2149  << " was calculated " << G4endl
2150  << " which is more than the computed safety= "
2151  << fPreviousSafety/mm << " mm at that point." << G4endl
2152  << " This difference is "
2153  << diffShiftSaf/mm << " mm." << G4endl
2154  << " The tolerated accuracy is "
2155  << fAccuracyForException/mm << " mm.";
2156 
2157  suggestion << " ";
2158  static G4ThreadLocal G4int warnNow = 0;
2159  if( ((++warnNow % 100) == 1) )
2160  {
2161  message << G4endl
2162  << " This problem can be due to either " << G4endl
2163  << " - a process that has proposed a displacement"
2164  << " larger than the current safety , or" << G4endl
2165  << " - inaccuracy in the computation of the safety";
2166  suggestion << "We suggest that you " << G4endl
2167  << " - find i) what particle is being tracked, and "
2168  << " ii) through what part of your geometry " << G4endl
2169  << " for example by re-running this event with "
2170  << G4endl
2171  << " /tracking/verbose 1 " << G4endl
2172  << " - check which processes you declare for"
2173  << " this particle (and look at non-standard ones)"
2174  << G4endl
2175  << " - in case, create a detailed logfile"
2176  << " of this event using:" << G4endl
2177  << " /tracking/verbose 6 ";
2178  }
2179  G4Exception("G4Navigator::ComputeStep()",
2180  "GeomNav1002", JustWarning,
2181  message, G4String(suggestion.str()));
2182  G4cout.precision(oldcoutPrec);
2183  G4cerr.precision(oldcerrPrec);
2184  }
2185 #ifdef G4DEBUG_NAVIGATION
2186  else
2187  {
2188  G4cerr << "WARNING - G4Navigator::ComputeStep()" << G4endl
2189  << " The Step's starting point has moved "
2190  << std::sqrt(moveLenSq) << "," << G4endl
2191  << " which has taken it to the limit of"
2192  << " the current safety. " << G4endl;
2193  }
2194 #endif
2195  }
2196  G4double safetyPlus = fPreviousSafety + fAccuracyForException;
2197  if ( shiftOriginSafSq > sqr(safetyPlus) )
2198  {
2199  std::ostringstream message;
2200  message << "May lead to a crash or unreliable results." << G4endl
2201  << " Position has shifted considerably without"
2202  << " notifying the navigator !" << G4endl
2203  << " Tolerated safety: " << safetyPlus << G4endl
2204  << " Computed shift : " << shiftOriginSafSq;
2205  G4Exception("G4Navigator::ComputeStep()", "GeomNav1002",
2206  JustWarning, message);
2207  }
2208 }
2209 
2210 // ********************************************************************
2211 // Operator <<
2212 // ********************************************************************
2213 //
2214 std::ostream& operator << (std::ostream &os,const G4Navigator &n)
2215 {
2216  // Old version did only the following:
2217  // os << "Current History: " << G4endl << n.fHistory;
2218  // Old behaviour is recovered for fVerbose = 0
2219 
2220  // Adapted from G4Navigator::PrintState() const
2221 
2222  G4int oldcoutPrec = os.precision(4);
2223  if( n.fVerbose >= 4 )
2224  {
2225  os << "The current state of G4Navigator is: " << G4endl;
2226  os << " ValidExitNormal= " << n.fValidExitNormal << G4endl
2227  << " ExitNormal = " << n.fExitNormal << G4endl
2228  << " Exiting = " << n.fExiting << G4endl
2229  << " Entering = " << n.fEntering << G4endl
2230  << " BlockedPhysicalVolume= " ;
2231  if (n.fBlockedPhysicalVolume==0)
2232  os << "None";
2233  else
2234  os << n.fBlockedPhysicalVolume->GetName();
2235  os << G4endl
2236  << " BlockedReplicaNo = " << n.fBlockedReplicaNo << G4endl
2237  << " LastStepWasZero = " << n.fLastStepWasZero << G4endl
2238  << G4endl;
2239  }
2240  if( ( 1 < n.fVerbose) && (n.fVerbose < 4) )
2241  {
2242  os << G4endl; // Make sure to line up
2243  os << std::setw(30) << " ExitNormal " << " "
2244  << std::setw( 5) << " Valid " << " "
2245  << std::setw( 9) << " Exiting " << " "
2246  << std::setw( 9) << " Entering" << " "
2247  << std::setw(15) << " Blocked:Volume " << " "
2248  << std::setw( 9) << " ReplicaNo" << " "
2249  << std::setw( 8) << " LastStepZero " << " "
2250  << G4endl;
2251  os << "( " << std::setw(7) << n.fExitNormal.x()
2252  << ", " << std::setw(7) << n.fExitNormal.y()
2253  << ", " << std::setw(7) << n.fExitNormal.z() << " ) "
2254  << std::setw( 5) << n.fValidExitNormal << " "
2255  << std::setw( 9) << n.fExiting << " "
2256  << std::setw( 9) << n.fEntering << " ";
2257  if ( n.fBlockedPhysicalVolume==0 )
2258  { os << std::setw(15) << "None"; }
2259  else
2260  { os << std::setw(15)<< n.fBlockedPhysicalVolume->GetName(); }
2261  os << std::setw( 9) << n.fBlockedReplicaNo << " "
2262  << std::setw( 8) << n.fLastStepWasZero << " "
2263  << G4endl;
2264  }
2265  if( n.fVerbose > 2 )
2266  {
2267  os.precision(8);
2268  os << " Current Localpoint = " << n.fLastLocatedPointLocal << G4endl;
2269  os << " PreviousSftOrigin = " << n.fPreviousSftOrigin << G4endl;
2270  os << " PreviousSafety = " << n.fPreviousSafety << G4endl;
2271  }
2272  if( n.fVerbose > 3 || n.fVerbose == 0 )
2273  {
2274  os << "Current History: " << G4endl << n.fHistory;
2275  }
2276 
2277  os.precision(oldcoutPrec);
2278  return os;
2279 }
G4String GetName() const
G4ReplicaNavigation freplicaNav
Definition: G4Navigator.hh:503
G4bool fWarnPush
Definition: G4Navigator.hh:495
G4ParameterisedNavigation fparamNav
Definition: G4Navigator.hh:502
G4SmartVoxelHeader * GetVoxelHeader() const
G4bool fValidExitNormal
Definition: G4Navigator.hh:423
G4VPhysicalVolume * GetTopVolume() const
G4VoxelNavigation fvoxelNav
Definition: G4Navigator.hh:501
G4VPhysicalVolume * fBlockedPhysicalVolume
Definition: G4Navigator.hh:415
virtual void ResetState()
G4bool fLastTriedStepComputation
Definition: G4Navigator.hh:400
virtual G4bool LevelLocate(G4NavigationHistory &history, const G4VPhysicalVolume *blockedVol, const G4int blockedNum, const G4ThreeVector &globalPoint, const G4ThreeVector *globalDirection, const G4bool pLocatedOnEdge, G4ThreeVector &localPoint)
G4VoxelSafety * fpVoxelSafety
Definition: G4Navigator.hh:505
G4double ComputeStep(const G4ThreeVector &globalPoint, const G4ThreeVector &globalDirection, const G4double currentProposedStepLength, G4double &newSafety, G4NavigationHistory &history, G4bool &validExitNormal, G4ThreeVector &exitNormal, G4bool &exiting, G4bool &entering, G4VPhysicalVolume *(*pBlockedPhysical), G4int &blockedReplicaNo)
G4bool fCheck
Definition: G4Navigator.hh:493
const G4ThreeVector & GetTranslation() const
static const G4double kInfinity
Definition: geomdefs.hh:42
virtual G4bool IsNested() const
#define fWasLimitedByGeometry
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
void UpdateMaterial(G4Material *pMaterial)
G4String GetName() const
CLHEP::Hep3Vector G4ThreeVector
CLHEP::HepRotation G4RotationMatrix
G4ThreeVector fLastStepEndPointLocal
Definition: G4Navigator.hh:388
virtual void ComputeDimensions(G4VPVParameterisation *p, const G4int n, const G4VPhysicalVolume *pRep)
Definition: G4VSolid.cc:137
G4AffineTransform Inverse() const
G4int GetVerboseLevel() const
G4ThreeVector fGrandMotherExitNormal
Definition: G4Navigator.hh:427
G4ThreeVector ComputeLocalAxis(const G4ThreeVector &pVec) const
G4bool fPushed
Definition: G4Navigator.hh:495
virtual G4TouchableHistoryHandle CreateTouchableHistoryHandle() const
virtual G4VSolid * ComputeSolid(const G4int, G4VPhysicalVolume *)
void ResetStackAndState()
G4int GetDepth() const
G4int fNumberZeroSteps
Definition: G4Navigator.hh:447
G4double GetSurfaceTolerance() const
G4VPhysicalVolume * spBlockedPhysicalVolume
Definition: G4Navigator.hh:471
virtual G4double ComputeStep(const G4ThreeVector &pGlobalPoint, const G4ThreeVector &pDirection, const G4double pCurrentProposedStepLength, G4double &pNewSafety)
Definition: G4Navigator.cc:744
G4int fVerbose
Definition: G4Navigator.hh:392
G4double ComputeSafety(const G4ThreeVector &localPoint, const G4VPhysicalVolume &currentPhysical, G4double maxLength=DBL_MAX)
void PrintState() const
void SetSolid(G4VSolid *pSolid)
G4int fAbandonThreshold_NoZeroSteps
Definition: G4Navigator.hh:451
G4SmartVoxelNode * ParamVoxelLocate(G4SmartVoxelHeader *pHead, const G4ThreeVector &localPoint)
G4bool fEnteredDaughter
Definition: G4Navigator.hh:372
EVolume GetVolumeType(G4int n) const
G4double ComputeStep(const G4ThreeVector &globalPoint, const G4ThreeVector &globalDirection, const G4ThreeVector &localPoint, const G4ThreeVector &localDirection, const G4double currentProposedStepLength, G4double &newSafety, G4NavigationHistory &history, G4bool &validExitNormal, G4bool &calculatedExitNormal, G4ThreeVector &exitNormal, G4bool &exiting, G4bool &entering, G4VPhysicalVolume *(*pBlockedPhysical), G4int &blockedReplicaNo)
G4bool fLocatedOnEdge
Definition: G4Navigator.hh:445
const G4AffineTransform GetLocalToGlobalTransform() const
virtual G4GeometryType GetEntityType() const =0
#define G4ThreadLocal
Definition: tls.hh:89
G4ReferenceCountedHandle< G4TouchableHistory > G4TouchableHistoryHandle
G4bool LevelLocate(G4NavigationHistory &history, const G4VPhysicalVolume *blockedVol, const G4int blockedNum, const G4ThreeVector &globalPoint, const G4ThreeVector *globalDirection, const G4bool pLocatedOnEdge, G4ThreeVector &localPoint)
void ComputeTransformation(const G4int replicaNo, G4VPhysicalVolume *pVol, G4ThreeVector &point) const
int G4int
Definition: G4Types.hh:78
G4bool fLastStepWasZero
Definition: G4Navigator.hh:442
virtual G4VPhysicalVolume * ResetHierarchyAndLocate(const G4ThreeVector &point, const G4ThreeVector &direction, const G4TouchableHistory &h)
Definition: G4Navigator.cc:92
G4ThreeVector fExitNormalGlobalFrame
Definition: G4Navigator.hh:431
virtual G4ThreeVector GetLocalExitNormalAndCheck(const G4ThreeVector &point, G4bool *valid)
G4int GetTopReplicaNo() const
void SetVerboseLevel(G4int level)
#define fPushed
G4double ComputeSafety(const G4ThreeVector &globalpoint, const G4NavigationHistory &history, const G4double pMaxLength=DBL_MAX)
static double normal(HepRandomEngine *eptr)
Definition: RandPoisson.cc:77
G4bool fWasLimitedByGeometry
Definition: G4Navigator.hh:382
G4NormalNavigation fnormalNav
Definition: G4Navigator.hh:500
G4bool fEntering
Definition: G4Navigator.hh:406
void SetNormalNavigation(G4NormalNavigation *fnormnav)
G4double ComputeSafety(const G4ThreeVector &localPoint, const G4NavigationHistory &history, const G4double pProposedMaxLength=DBL_MAX)
G4GLOB_DLL std::ostream G4cout
const G4String & GetName() const
virtual EInside Inside(const G4ThreeVector &p) const =0
virtual G4ThreeVector GetLocalExitNormal(G4bool *valid)
G4double ComputeStep(const G4ThreeVector &globalPoint, const G4ThreeVector &globalDirection, const G4double currentProposedStepLength, G4double &newSafety, G4NavigationHistory &history, G4bool &validExitNormal, G4ThreeVector &exitNormal, G4bool &exiting, G4bool &entering, G4VPhysicalVolume *(*pBlockedPhysical), G4int &blockedReplicaNo)
G4ThreeVector fPreviousSftOrigin
Definition: G4Navigator.hh:454
G4double CheckNextStep(const G4ThreeVector &pGlobalPoint, const G4ThreeVector &pDirection, const G4double pCurrentProposedStepLength, G4double &pNewSafety)
bool G4bool
Definition: G4Types.hh:79
EVolume GetTopVolumeType() const
virtual G4int GetRegularStructureId() const =0
G4NavigationHistory fHistory
Definition: G4Navigator.hh:368
virtual G4double ComputeSafety(const G4ThreeVector &globalpoint, const G4NavigationHistory &history, const G4double pMaxLength=DBL_MAX)
virtual G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const =0
virtual G4VPVParameterisation * GetParameterisation() const =0
virtual void SetCopyNo(G4int CopyNo)=0
G4int GetReplicaNo(G4int n) const
G4bool LevelLocate(G4NavigationHistory &history, const G4VPhysicalVolume *blockedVol, const G4int blockedNum, const G4ThreeVector &globalPoint, const G4ThreeVector *globalDirection, const G4bool pLocatedOnEdge, G4ThreeVector &localPoint)
G4bool fExitedMother
Definition: G4Navigator.hh:378
EVolume CharacteriseDaughters() const
virtual G4bool RecheckDistanceToCurrentBoundary(const G4ThreeVector &pGlobalPoint, const G4ThreeVector &pDirection, const G4double CurrentProposedStepLength, G4double *prDistance, G4double *prNewSafety=0) const
G4bool EnteredDaughterVolume() const
virtual G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const =0
G4int fBlockedReplicaNo
Definition: G4Navigator.hh:416
#define G4DEBUG_NAVIGATION
const G4int n
EVolume CharacteriseDaughters(const G4LogicalVolume *pLog) const
G4ThreeVector fLastLocatedPointLocal
Definition: G4Navigator.hh:418
G4double fPreviousSafety
Definition: G4Navigator.hh:455
G4bool fLocatedOutsideWorld
Definition: G4Navigator.hh:420
void SetSavedState()
Definition: G4Navigator.cc:656
G4bool fActive
Definition: G4Navigator.hh:397
G4ThreeVector ComputeLocalPoint(const G4ThreeVector &rGlobPoint) const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4double kCarTolerance
Definition: G4Navigator.hh:361
static const double perMillion
Definition: G4SIunits.hh:298
void NewLevel(G4VPhysicalVolume *pNewMother, EVolume vType=kNormal, G4int nReplica=-1)
G4ThreeVector TransformPoint(const G4ThreeVector &vec) const
G4TouchableHistory * CreateTouchableHistory() const
G4double ComputeStep(const G4ThreeVector &localPoint, const G4ThreeVector &localDirection, const G4double currentProposedStepLength, G4double &newSafety, G4NavigationHistory &history, G4bool &validExitNormal, G4ThreeVector &exitNormal, G4bool &exiting, G4bool &entering, G4VPhysicalVolume *(*pBlockedPhysical), G4int &blockedReplicaNo)
G4AffineTransform GetMotherToDaughterTransform(G4VPhysicalVolume *dVolume, G4int dReplicaNo, EVolume dVolumeType)
G4LogicalVolume * GetLogicalVolume() const
const G4AffineTransform & GetTransform(G4int n) const
virtual G4ThreeVector GetGlobalExitNormal(const G4ThreeVector &point, G4bool *valid)
G4bool fExiting
Definition: G4Navigator.hh:406
EInside
Definition: geomdefs.hh:58
G4int GetDaughtersRegularStructureId(const G4LogicalVolume *pLog) const
virtual void ComputeTransformation(const G4int, G4VPhysicalVolume *) const =0
G4ThreeVector TransformAxis(const G4ThreeVector &axis) const
virtual G4int GetCopyNo() const =0
const G4AffineTransform & GetGlobalToLocalTransform() const
virtual G4VPhysicalVolume * LocateGlobalPointAndSetup(const G4ThreeVector &point, const G4ThreeVector *direction=0, const G4bool pRelativeSearch=true, const G4bool ignoreDirection=true)
Definition: G4Navigator.cc:122
struct G4Navigator::G4SaveNavigatorState fSaveState
const G4RotationMatrix * GetRotation() const
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
Definition: Step.hh:41
G4SmartVoxelNode * VoxelLocate(G4SmartVoxelHeader *pHead, const G4ThreeVector &localPoint)
G4double ComputeSafety(const G4ThreeVector &localPoint, const G4NavigationHistory &history, const G4double pProposedMaxLength=DBL_MAX)
#define G4endl
Definition: G4ios.hh:61
const G4AffineTransform & GetTopTransform() const
std::ostream & operator<<(std::ostream &os, const G4Navigator &n)
G4bool fCalculatedExitNormal
Definition: G4Navigator.hh:433
G4ThreeVector fStepEndPoint
Definition: G4Navigator.hh:385
const G4NavigationHistory * GetHistory() const
G4ThreeVector fExitNormal
Definition: G4Navigator.hh:424
G4int fActionThreshold_NoZeroSteps
Definition: G4Navigator.hh:449
T sqr(const T &x)
Definition: templates.hh:145
G4bool fChangedGrandMotherRefFrame
Definition: G4Navigator.hh:429
EVolume
Definition: geomdefs.hh:68
G4int MoveUpHistory(G4int num_levels=1)
G4double ComputeSafety(const G4ThreeVector &globalPoint, const G4ThreeVector &localPoint, G4NavigationHistory &history, const G4double pProposedMaxLength=DBL_MAX)
G4bool LevelLocate(G4NavigationHistory &history, const G4VPhysicalVolume *blockedVol, const G4int blockedNum, const G4ThreeVector &globalPoint, const G4ThreeVector *globalDirection, const G4bool pLocatedOnEdge, G4ThreeVector &localPoint)
double G4double
Definition: G4Types.hh:76
G4VPhysicalVolume * GetVolume(G4int n) const
G4bool LevelLocate(G4NavigationHistory &history, const G4VPhysicalVolume *blockedVol, const G4int blockedNum, const G4ThreeVector &globalPoint, const G4ThreeVector *globalDirection, const G4bool pLocatedOnEdge, G4ThreeVector &localPoint)
virtual G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=false, G4bool *validNorm=0, G4ThreeVector *n=0) const =0
virtual G4double ComputeSafety(const G4ThreeVector &globalpoint, const G4double pProposedMaxLength=DBL_MAX, const G4bool keepState=true)
virtual G4double ComputeStep(const G4ThreeVector &globalPoint, const G4ThreeVector &globalDirection, const G4double currentProposedStepLength, G4double &newSafety, G4NavigationHistory &history, G4bool &validExitNormal, G4ThreeVector &exitNormal, G4bool &exiting, G4bool &entering, G4VPhysicalVolume *(*pBlockedPhysical), G4int &blockedReplicaNo)
G4RegularNavigation fregularNav
Definition: G4Navigator.hh:504
EInside BackLocate(G4NavigationHistory &history, const G4ThreeVector &globalPoint, G4ThreeVector &localPoint, const G4bool &exiting, G4bool &notKnownInside) const
#define DBL_MAX
Definition: templates.hh:83
virtual void SetupHierarchy()
static const double mm
Definition: G4SIunits.hh:102
virtual G4bool CheckOverlaps(G4int res=1000, G4double tol=0., G4bool verbose=true, G4int errMax=1)
EVolume VolumeType(const G4VPhysicalVolume *pVol) const
static G4GeometryTolerance * GetInstance()
virtual void LocateGlobalPointWithinVolume(const G4ThreeVector &position)
Definition: G4Navigator.cc:579
G4VSolid * GetSolid() const
void ComputeStepLog(const G4ThreeVector &pGlobalpoint, G4double moveLenSq) const
void RestoreSavedState()
Definition: G4Navigator.cc:690
virtual ~G4Navigator()
Definition: G4Navigator.cc:82
G4GLOB_DLL std::ostream G4cerr