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