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