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