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