Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ITNavigator.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // $Id: G4ITNavigator.cc 64376 2012-10-31 16:39:40Z gcosmo $
27 //
28 // class G4ITNavigator Implementation
29 //
30 // Original author: Paul Kent, July 95/96
31 //
32 // G4ITNavigator is a duplicate version of G4Navigator starting from Geant4.9.5
33 // initially written by Paul Kent and colleagues.
34 // The only difference resides in the way the information is saved and managed
35 //
36 // History:
37 // - Created. Paul Kent, Jul 95/96
38 // - Zero step protections J.A. / G.C., Nov 2004
39 // - Added check mode G. Cosmo, Mar 2004
40 // - Made Navigator Abstract G. Cosmo, Nov 2003
41 // - G4ITNavigator created M.K., Nov 2012
42 // --------------------------------------------------------------------
43 
44 #include <iomanip>
45 
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 
54 // ********************************************************************
55 // Constructor
56 // ********************************************************************
57 //
59  : fWasLimitedByGeometry(false), fVerbose(0),
60  fTopPhysical(0), fCheck(false), fPushed(false), fWarnPush(true)
61 {
62  fActive= false;
63  fLastTriedStepComputation= false;
65 
66  fActionThreshold_NoZeroSteps = 10;
67  fAbandonThreshold_NoZeroSteps = 25;
68 
70  fregularNav.SetNormalNavigation( &fnormalNav );
71 
72  fStepEndPoint = G4ThreeVector( kInfinity, kInfinity, kInfinity );
73  fLastStepEndPointLocal = G4ThreeVector( kInfinity, kInfinity, kInfinity );
74 
75  fpSaveState = 0;
76 
77  // this->SetVerboseLevel(3);
78  // this->CheckMode(true);
79 }
80 
81 // !>
82 
83 G4ITNavigator::G4SaveNavigatorState::G4SaveNavigatorState() : G4ITNavigatorState_Lock()
84 {
85  sWasLimitedByGeometry = false;
86  sEntering = false;
87  sExiting = false;
88  sLocatedOnEdge = false;
89  sLastStepWasZero = false;
90  sEnteredDaughter = false;
91  sExitedMother = false;
92  sPushed = false;
93 
94  sValidExitNormal = false;
95  sExitNormal = G4ThreeVector(0,0,0);
96 
97  sPreviousSftOrigin = G4ThreeVector(0,0,0);
98  sPreviousSafety = 0.0;
99 
100  sNumberZeroSteps = 0;
101 
102  spBlockedPhysicalVolume = 0;
103  sBlockedReplicaNo = -1;
104 
105  sLastLocatedPointLocal = G4ThreeVector( kInfinity, -kInfinity, 0.0 );
106  sLocatedOutsideWorld = false;
107 }
108 
109 // <!
110 
111 // ********************************************************************
112 // Destructor
113 // ********************************************************************
114 //
116 {;}
117 
118 // ********************************************************************
119 // ResetHierarchyAndLocate
120 // ********************************************************************
121 //
124  const G4ThreeVector &direction,
125  const G4TouchableHistory &h)
126 {
127  ResetState();
128  fHistory = *h.GetHistory();
129  SetupHierarchy();
130  fLastTriedStepComputation= false; // Redundant, but best
131  return LocateGlobalPointAndSetup(p, &direction, true, false);
132 }
133 
134 // ********************************************************************
135 // LocateGlobalPointAndSetup
136 //
137 // Locate the point in the hierarchy return 0 if outside
138 // The direction is required
139 // - if on an edge shared by more than two surfaces
140 // (to resolve likely looping in tracking)
141 // - at initial location of a particle
142 // (to resolve potential ambiguity at boundary)
143 //
144 // Flags on exit: (comments to be completed)
145 // fEntering - True if entering `daughter' volume (or replica)
146 // whether daughter of last mother directly
147 // or daughter of that volume's ancestor.
148 // ********************************************************************
149 //
152  const G4ThreeVector* pGlobalDirection,
153  const G4bool relativeSearch,
154  const G4bool ignoreDirection )
155 {
156  G4bool notKnownContained=true, noResult;
157  G4VPhysicalVolume *targetPhysical;
158  G4LogicalVolume *targetLogical;
159  G4VSolid *targetSolid=0;
160  G4ThreeVector localPoint, globalDirection;
161  EInside insideCode;
162 
163  G4bool considerDirection = (!ignoreDirection) || fLocatedOnEdge;
164  fLastTriedStepComputation= false;
165 
166  if( considerDirection && pGlobalDirection != 0 )
167  {
168  globalDirection=*pGlobalDirection;
169  }
170 
171 #ifdef G4VERBOSE
172  if( fVerbose > 2 )
173  {
174  G4int oldcoutPrec = G4cout.precision(8);
175  G4cout << "*** G4ITNavigator::LocateGlobalPointAndSetup: ***" << G4endl;
176  G4cout << " Called with arguments: " << G4endl
177  << " Globalpoint = " << globalPoint << G4endl
178  << " RelativeSearch = " << relativeSearch << G4endl;
179  if( fVerbose == 4 )
180  {
181  G4cout << " ----- Upon entering:" << G4endl;
182  PrintState();
183  }
184  G4cout.precision(oldcoutPrec);
185  }
186 #endif
187 
188  if ( !relativeSearch )
189  {
191  }
192  else
193  {
194  if ( fWasLimitedByGeometry )
195  {
196  fWasLimitedByGeometry = false;
197  fEnteredDaughter = fEntering; // Remember
198  fExitedMother = fExiting; // Remember
199  if ( fExiting )
200  {
201  if ( fHistory.GetDepth() )
202  {
203  fBlockedPhysicalVolume = fHistory.GetTopVolume();
204  fBlockedReplicaNo = fHistory.GetTopReplicaNo();
206  }
207  else
208  {
209  fLastLocatedPointLocal = localPoint;
210  fLocatedOutsideWorld = true;
211  return 0; // Have exited world volume
212  }
213  // A fix for the case where a volume is "entered" at an edge
214  // and a coincident surface exists outside it.
215  // - This stops it from exiting further volumes and cycling
216  // - However ReplicaNavigator treats this case itself
217  //
218  if ( fLocatedOnEdge && (VolumeType(fBlockedPhysicalVolume)!=kReplica ))
219  {
220  fExiting= false;
221  }
222  }
223  else
224  if ( fEntering )
225  {
226  switch (VolumeType(fBlockedPhysicalVolume))
227  {
228  case kNormal:
229  fHistory.NewLevel(fBlockedPhysicalVolume, kNormal,
230  fBlockedPhysicalVolume->GetCopyNo());
231  break;
232  case kReplica:
233  freplicaNav.ComputeTransformation(fBlockedReplicaNo,
234  fBlockedPhysicalVolume);
235  fHistory.NewLevel(fBlockedPhysicalVolume, kReplica,
236  fBlockedReplicaNo);
237  fBlockedPhysicalVolume->SetCopyNo(fBlockedReplicaNo);
238  break;
239  case kParameterised:
240  if( fBlockedPhysicalVolume->GetRegularStructureId() == 0 )
241  {
242  G4VSolid *pSolid;
243  G4VPVParameterisation *pParam;
244  G4TouchableHistory parentTouchable( fHistory );
245  pParam = fBlockedPhysicalVolume->GetParameterisation();
246  pSolid = pParam->ComputeSolid(fBlockedReplicaNo,
247  fBlockedPhysicalVolume);
248  pSolid->ComputeDimensions(pParam, fBlockedReplicaNo,
249  fBlockedPhysicalVolume);
250  pParam->ComputeTransformation(fBlockedReplicaNo,
251  fBlockedPhysicalVolume);
252  fHistory.NewLevel(fBlockedPhysicalVolume, kParameterised,
253  fBlockedReplicaNo);
254  fBlockedPhysicalVolume->SetCopyNo(fBlockedReplicaNo);
255  //
256  // Set the correct solid and material in Logical Volume
257  //
258  G4LogicalVolume *pLogical;
259  pLogical = fBlockedPhysicalVolume->GetLogicalVolume();
260  pLogical->SetSolid( pSolid );
261  pLogical->UpdateMaterial(pParam ->
262  ComputeMaterial(fBlockedReplicaNo,
263  fBlockedPhysicalVolume,
264  &parentTouchable));
265  }
266  break;
267  }
268  fEntering = false;
269  fBlockedPhysicalVolume = 0;
270  localPoint = fHistory.GetTopTransform().TransformPoint(globalPoint);
271  notKnownContained = false;
272  }
273  }
274  else
275  {
276  fBlockedPhysicalVolume = 0;
277  fEntering = false;
278  fEnteredDaughter = false; // Full Step was not taken, did not enter
279  fExiting = false;
280  fExitedMother = false; // Full Step was not taken, did not exit
281  }
282  }
283  //
284  // Search from top of history up through geometry until
285  // containing volume found:
286  // If on
287  // o OUTSIDE - Back up level, not/no longer exiting volumes
288  // o SURFACE and EXITING - Back up level, setting new blocking no.s
289  // else
290  // o containing volume found
291  //
292  while (notKnownContained)
293  {
295  {
296  targetSolid = fHistory.GetTopVolume()->GetLogicalVolume()->GetSolid();
297  localPoint = fHistory.GetTopTransform().TransformPoint(globalPoint);
298  insideCode = targetSolid->Inside(localPoint);
299 #ifdef G4VERBOSE
300  if(( fVerbose == 1 ) && ( fCheck ))
301  {
302  G4String solidResponse = "-kInside-";
303  if (insideCode == kOutside)
304  solidResponse = "-kOutside-";
305  else if (insideCode == kSurface)
306  solidResponse = "-kSurface-";
307  G4cout << "*** G4ITNavigator::LocateGlobalPointAndSetup(): ***" << G4endl
308  << " Invoked Inside() for solid: " << targetSolid->GetName()
309  << ". Solid replied: " << solidResponse << G4endl
310  << " For local point p: " << localPoint << G4endl;
311  }
312 #endif
313  }
314  else
315  {
316  insideCode = freplicaNav.BackLocate(fHistory, globalPoint, localPoint,
317  fExiting, notKnownContained);
318  // !CARE! if notKnownContained returns false then the point is within
319  // the containing placement volume of the replica(s). If insidecode
320  // will result in the history being backed up one level, then the
321  // local point returned is the point in the system of this new level
322  }
323  if ( insideCode==kOutside )
324  {
325  if ( fHistory.GetDepth() )
326  {
327  fBlockedPhysicalVolume = fHistory.GetTopVolume();
328  fBlockedReplicaNo = fHistory.GetTopReplicaNo();
330  fExiting = false;
331  }
332  else
333  {
334  fLastLocatedPointLocal = localPoint;
335  fLocatedOutsideWorld = true;
336  return 0; // Have exited world volume
337  }
338  }
339  else
340  if ( insideCode==kSurface )
341  {
342  G4bool isExiting = fExiting;
343  if( (!fExiting)&&considerDirection )
344  {
345  // Figure out whether we are exiting this level's volume
346  // by using the direction
347  //
348  G4bool directionExiting = false;
349  G4ThreeVector localDirection =
350  fHistory.GetTopTransform().TransformAxis(globalDirection);
352  {
353  G4ThreeVector normal = targetSolid->SurfaceNormal(localPoint);
354  directionExiting = normal.dot(localDirection) > 0.0;
355  isExiting = isExiting || directionExiting;
356  }
357  }
358  if( isExiting )
359  {
360  if ( fHistory.GetDepth() )
361  {
362  fBlockedPhysicalVolume = fHistory.GetTopVolume();
363  fBlockedReplicaNo = fHistory.GetTopReplicaNo();
365  //
366  // Still on surface but exited volume not necessarily convex
367  //
368  fValidExitNormal = false;
369  }
370  else
371  {
372  fLastLocatedPointLocal = localPoint;
373  fLocatedOutsideWorld = true;
374  return 0; // Have exited world volume
375  }
376  }
377  else
378  {
379  notKnownContained=false;
380  }
381  }
382  else
383  {
384  notKnownContained=false;
385  }
386  } // END while (notKnownContained)
387  //
388  // Search downwards until deepest containing volume found,
389  // blocking fBlockedPhysicalVolume/BlockedReplicaNum
390  //
391  // 3 Cases:
392  //
393  // o Parameterised daughters
394  // =>Must be one G4PVParameterised daughter & voxels
395  // o Positioned daughters & voxels
396  // o Positioned daughters & no voxels
397 
398  noResult = true; // noResult should be renamed to
399  // something like enteredLevel, as that is its meaning.
400  do
401  {
402  // Determine `type' of current mother volume
403  //
404  targetPhysical = fHistory.GetTopVolume();
405  if (!targetPhysical) { break; }
406  targetLogical = targetPhysical->GetLogicalVolume();
407  switch( CharacteriseDaughters(targetLogical) )
408  {
409  case kNormal:
410  if ( targetLogical->GetVoxelHeader() ) // use optimised navigation
411  {
412  noResult = fvoxelNav.LevelLocate(fHistory,
413  fBlockedPhysicalVolume,
414  fBlockedReplicaNo,
415  globalPoint,
416  pGlobalDirection,
417  considerDirection,
418  localPoint);
419  }
420  else // do not use optimised navigation
421  {
422  noResult = fnormalNav.LevelLocate(fHistory,
423  fBlockedPhysicalVolume,
424  fBlockedReplicaNo,
425  globalPoint,
426  pGlobalDirection,
427  considerDirection,
428  localPoint);
429  }
430  break;
431  case kReplica:
432  noResult = freplicaNav.LevelLocate(fHistory,
433  fBlockedPhysicalVolume,
434  fBlockedReplicaNo,
435  globalPoint,
436  pGlobalDirection,
437  considerDirection,
438  localPoint);
439  break;
440  case kParameterised:
441  if( GetDaughtersRegularStructureId(targetLogical) != 1 )
442  {
443  noResult = fparamNav.LevelLocate(fHistory,
444  fBlockedPhysicalVolume,
445  fBlockedReplicaNo,
446  globalPoint,
447  pGlobalDirection,
448  considerDirection,
449  localPoint);
450  }
451  else // Regular structure
452  {
453  noResult = fregularNav.LevelLocate(fHistory,
454  fBlockedPhysicalVolume,
455  fBlockedReplicaNo,
456  globalPoint,
457  pGlobalDirection,
458  considerDirection,
459  localPoint);
460  }
461  break;
462  }
463 
464  // LevelLocate returns true if it finds a daughter volume
465  // in which globalPoint is inside (or on the surface).
466 
467  if ( noResult )
468  {
469  // Entering a daughter after ascending
470  //
471  // The blocked volume is no longer valid - it was for another level
472  //
473  fBlockedPhysicalVolume = 0;
474  fBlockedReplicaNo = -1;
475 
476  // fEntering should be false -- else blockedVolume is assumed good.
477  // fEnteredDaughter is used for ExitNormal
478  //
479  fEntering = false;
480  fEnteredDaughter = true;
481 #ifdef G4DEBUG_NAVIGATION
482  if( fVerbose > 2 )
483  {
484  G4VPhysicalVolume* enteredPhysical = fHistory.GetTopVolume();
485  G4cout << "*** G4ITNavigator::LocateGlobalPointAndSetup() ***" << G4endl;
486  G4cout << " Entering volume: " << enteredPhysical->GetName()
487  << G4endl;
488  }
489 #endif
490  }
491  } while (noResult);
492 
493  fLastLocatedPointLocal = localPoint;
494 
495 #ifdef G4VERBOSE
496  if( fVerbose == 4 )
497  {
498  G4int oldcoutPrec = G4cout.precision(8);
499  G4String curPhysVol_Name("None");
500  if (targetPhysical) { curPhysVol_Name = targetPhysical->GetName(); }
501  G4cout << " Return value = new volume = " << curPhysVol_Name << G4endl;
502  G4cout << " ----- Upon exiting:" << G4endl;
503  PrintState();
504 #ifdef G4DEBUG_NAVIGATION
505  G4cout << "Upon exiting LocateGlobalPointAndSetup():" << G4endl;
506  G4cout << " History = " << G4endl << fHistory << G4endl << G4endl;
507 #endif
508  G4cout.precision(oldcoutPrec);
509  }
510 #endif
511 
512  fLocatedOutsideWorld= false;
513 
514  return targetPhysical;
515 }
516 
517 // ********************************************************************
518 // LocateGlobalPointWithinVolume
519 //
520 // -> the state information of this Navigator and its subNavigators
521 // is updated in order to start the next step at pGlobalpoint
522 // -> no check is performed whether pGlobalpoint is inside the
523 // original volume (this must be the case).
524 //
525 // Note: a direction could be added to the arguments, to aid in future
526 // optional checking (via the old code below, flagged by OLD_LOCATE).
527 // [ This would be done only in verbose mode ]
528 // ********************************************************************
529 //
530 void
532 {
533  fLastLocatedPointLocal = ComputeLocalPoint(pGlobalpoint);
534  fLastTriedStepComputation= false;
535 
536 #ifdef G4DEBUG_NAVIGATION
537  if( fVerbose > 2 )
538  {
539  G4cout << "Entering LocateGlobalWithinVolume(): History = " << G4endl;
540  G4cout << fHistory << G4endl;
541  }
542 #endif
543 
544  // For the case of Voxel (or Parameterised) volume the respective
545  // Navigator must be messaged to update its voxel information etc
546 
547  // Update the state of the Sub Navigators
548  // - in particular any voxel information they store/cache
549  //
550  G4VPhysicalVolume* motherPhysical = fHistory.GetTopVolume();
551  G4LogicalVolume* motherLogical = motherPhysical->GetLogicalVolume();
552  G4SmartVoxelHeader* pVoxelHeader = motherLogical->GetVoxelHeader();
553 
555  {
556  switch( CharacteriseDaughters(motherLogical) )
557  {
558  case kNormal:
559  if ( pVoxelHeader )
560  {
561  fvoxelNav.VoxelLocate( pVoxelHeader, fLastLocatedPointLocal );
562  }
563  break;
564  case kParameterised:
565  if( GetDaughtersRegularStructureId(motherLogical) != 1 )
566  {
567  // Resets state & returns voxel node
568  //
569  fparamNav.ParamVoxelLocate( pVoxelHeader, fLastLocatedPointLocal );
570  }
571  break;
572  case kReplica:
573  G4Exception("G4ITNavigator::LocateGlobalPointWithinVolume()",
574  "GeomNav0001", FatalException,
575  "Not applicable for replicated volumes.");
576  break;
577  }
578  }
579 
580  // Reset the state variables
581  // - which would have been affected
582  // by the 'equivalent' call to LocateGlobalPointAndSetup
583  // - who's values have been invalidated by the 'move'.
584  //
585  fBlockedPhysicalVolume = 0;
586  fBlockedReplicaNo = -1;
587  fEntering = false;
588  fEnteredDaughter = false; // Boundary not encountered, did not enter
589  fExiting = false;
590  fExitedMother = false; // Boundary not encountered, did not exit
591 }
592 
593 // !>
595 {
596  SetSavedState();
597  return fpSaveState;
598 }
599 
601 {
602  fpSaveState = (G4SaveNavigatorState*) navState;
603  if(navState) RestoreSavedState();
604 }
605 
607 {
608  fpSaveState = new G4SaveNavigatorState();
609  ResetState();
610 }
611 
612 
613 // ********************************************************************
614 // SetSavedState
615 //
616 // Save the state, in case this is a parasitic call
617 // Save fValidExitNormal, fExitNormal, fExiting, fEntering,
618 // fBlockedPhysicalVolume, fBlockedReplicaNo, fLastStepWasZero;
619 // ********************************************************************
620 //
622 {
623  // !>
624  // This check can be avoid if instead, at every first step of a track,
625  // the IT tracking uses NewNavigatorSate
626  // The normal tracking would just call once NewNavigatorState() before tracking
627 
628 // if(fpSaveState == 0)
629 // fpSaveState = new G4SaveNavigatorState;
630  // <!
631 
632  // fSaveExitNormal = fExitNormal;
633  fpSaveState->sExitNormal = fExitNormal;
634  fpSaveState->sValidExitNormal = fValidExitNormal;
635  fpSaveState->sExiting = fExiting;
636  fpSaveState->sEntering = fEntering;
637 
638  fpSaveState->spBlockedPhysicalVolume = fBlockedPhysicalVolume;
639  fpSaveState->sBlockedReplicaNo = fBlockedReplicaNo,
640 
641  fpSaveState->sLastStepWasZero = fLastStepWasZero;
642 
643  // !>
644  fpSaveState->sPreviousSftOrigin = fPreviousSftOrigin;
645  fpSaveState->sPreviousSafety = fPreviousSafety;
646  fpSaveState->sNumberZeroSteps = fNumberZeroSteps;
647  fpSaveState->sLocatedOnEdge = fLocatedOnEdge;
648  fpSaveState->sWasLimitedByGeometry= fWasLimitedByGeometry;
649  fpSaveState->sPushed=fPushed;
650  fpSaveState->sNumberZeroSteps=fNumberZeroSteps;
651  fpSaveState->sEnteredDaughter = fEnteredDaughter;
652  fpSaveState->sExitedMother = fExitedMother;
653 
654  fpSaveState->sLastLocatedPointLocal = fLastLocatedPointLocal;
655  fpSaveState->sLocatedOutsideWorld = fLocatedOutsideWorld;
656  // <!
657 }
658 
659 // ********************************************************************
660 // RestoreSavedState
661 //
662 // Restore the state (in Compute Step), in case this is a parasitic call
663 // ********************************************************************
664 //
666 {
667  fExitNormal = fpSaveState->sExitNormal;
668  fValidExitNormal = fpSaveState->sValidExitNormal;
669  fExiting = fpSaveState->sExiting;
670  fEntering = fpSaveState->sEntering;
671 
672  fBlockedPhysicalVolume = fpSaveState->spBlockedPhysicalVolume;
673  fBlockedReplicaNo = fpSaveState->sBlockedReplicaNo,
674 
675  fLastStepWasZero = fpSaveState->sLastStepWasZero;
676 
677  // !>
678  fPreviousSftOrigin = fpSaveState->sPreviousSftOrigin ;
679  fPreviousSafety = fpSaveState->sPreviousSafety ;
680  fNumberZeroSteps = fpSaveState->sNumberZeroSteps ;
681  fLocatedOnEdge = fpSaveState->sLocatedOnEdge ;
682  fWasLimitedByGeometry = fpSaveState->sWasLimitedByGeometry;
683  fPushed = fpSaveState->sPushed;
684  fNumberZeroSteps = fpSaveState->sNumberZeroSteps;
685  fEnteredDaughter= fpSaveState->sEnteredDaughter ;
686  fExitedMother = fpSaveState->sExitedMother ;
687 
688  fLastLocatedPointLocal = fpSaveState->sLastLocatedPointLocal ;
689  fLocatedOutsideWorld = fpSaveState->sLocatedOutsideWorld;
690  // <!
691 }
692 // <!
693 
694 // ********************************************************************
695 // ComputeStep
696 //
697 // Computes the next geometric Step: intersections with current
698 // mother and `daughter' volumes.
699 //
700 // NOTE:
701 //
702 // Flags on entry:
703 // --------------
704 // fValidExitNormal - Normal of exited volume is valid (convex, not a
705 // coincident boundary)
706 // fExitNormal - Surface normal of exited volume
707 // fExiting - True if have exited solid
708 //
709 // fBlockedPhysicalVolume - Ptr to exited volume (or 0)
710 // fBlockedReplicaNo - Replication no of exited volume
711 // fLastStepWasZero - True if last Step size was zero.
712 //
713 // Flags on exit:
714 // -------------
715 // fValidExitNormal - True if surface normal of exited volume is valid
716 // fExitNormal - Surface normal of exited volume rotated to mothers
717 // reference system
718 // fExiting - True if exiting mother
719 // fEntering - True if entering `daughter' volume (or replica)
720 // fBlockedPhysicalVolume - Ptr to candidate (entered) volume
721 // fBlockedReplicaNo - Replication no of candidate (entered) volume
722 // fLastStepWasZero - True if this Step size was zero.
723 // ********************************************************************
724 //
726  const G4ThreeVector &pDirection,
727  const G4double pCurrentProposedStepLength,
728  G4double &pNewSafety)
729 {
730  G4ThreeVector localDirection = ComputeLocalAxis(pDirection);
731  G4double Step = kInfinity;
732  G4VPhysicalVolume *motherPhysical = fHistory.GetTopVolume();
733  G4LogicalVolume *motherLogical = motherPhysical->GetLogicalVolume();
734 
735  static G4int sNavCScalls=0;
736  sNavCScalls++;
737 
738  fLastTriedStepComputation= true;
739 
740 #ifdef G4VERBOSE
741  if( fVerbose > 0 )
742  {
743  G4cout << "*** G4ITNavigator::ComputeStep: ***" << G4endl;
744  G4cout << " Volume = " << motherPhysical->GetName()
745  << " - Proposed step length = " << pCurrentProposedStepLength
746  << G4endl;
747 #ifdef G4DEBUG_NAVIGATION
748  if( fVerbose >= 4 )
749  {
750  G4cout << " Called with the arguments: " << G4endl
751  << " Globalpoint = " << std::setw(25) << pGlobalpoint << G4endl
752  << " Direction = " << std::setw(25) << pDirection << G4endl;
753  G4cout << " ---- Upon entering :" << G4endl;
754  PrintState();
755  }
756 #endif
757  }
758 #endif
759 
760  G4ThreeVector newLocalPoint = ComputeLocalPoint(pGlobalpoint);
761  if( newLocalPoint != fLastLocatedPointLocal )
762  {
763  // Check whether the relocation is within safety
764  //
765  G4ThreeVector oldLocalPoint = fLastLocatedPointLocal;
766  G4double moveLenSq = (newLocalPoint-oldLocalPoint).mag2();
767 
768  if ( moveLenSq >= kCarTolerance*kCarTolerance )
769  {
770 #ifdef G4VERBOSE
771  ComputeStepLog(pGlobalpoint, moveLenSq);
772 #endif
773  // Relocate the point within the same volume
774  //
775  LocateGlobalPointWithinVolume( pGlobalpoint );
776  fLastTriedStepComputation= true; // Ensure that this is set again !!
777  }
778  }
780  {
781  switch( CharacteriseDaughters(motherLogical) )
782  {
783  case kNormal:
784  if ( motherLogical->GetVoxelHeader() )
785  {
786  Step = fvoxelNav.ComputeStep(fLastLocatedPointLocal,
787  localDirection,
788  pCurrentProposedStepLength,
789  pNewSafety,
790  fHistory,
791  fValidExitNormal,
792  fExitNormal,
793  fExiting,
794  fEntering,
795  &fBlockedPhysicalVolume,
796  fBlockedReplicaNo);
797 
798  }
799  else
800  {
801  if( motherPhysical->GetRegularStructureId() == 0 )
802  {
803  Step = fnormalNav.ComputeStep(fLastLocatedPointLocal,
804  localDirection,
805  pCurrentProposedStepLength,
806  pNewSafety,
807  fHistory,
808  fValidExitNormal,
809  fExitNormal,
810  fExiting,
811  fEntering,
812  &fBlockedPhysicalVolume,
813  fBlockedReplicaNo);
814  }
815  else // Regular (non-voxelised) structure
816  {
817  LocateGlobalPointAndSetup( pGlobalpoint, &pDirection, true, true );
818  fLastTriedStepComputation= true; // Ensure that this is set again !!
819  //
820  // if physical process limits the step, the voxel will not be the
821  // one given by ComputeStepSkippingEqualMaterials() and the local
822  // point will be wrongly calculated.
823 
824  // There is a problem: when msc limits the step and the point is
825  // assigned wrongly to phantom in previous step (while it is out
826  // of the container volume). Then LocateGlobalPointAndSetup() has
827  // reset the history topvolume to world.
828  //
830  {
831  G4Exception("G4ITNavigator::ComputeStep()",
832  "GeomNav1001", JustWarning,
833  "Point is relocated in voxels, while it should be outside!");
834  Step = fnormalNav.ComputeStep(fLastLocatedPointLocal,
835  localDirection,
836  pCurrentProposedStepLength,
837  pNewSafety,
838  fHistory,
839  fValidExitNormal,
840  fExitNormal,
841  fExiting,
842  fEntering,
843  &fBlockedPhysicalVolume,
844  fBlockedReplicaNo);
845  }
846  else
847  {
848  Step = fregularNav.
849  ComputeStepSkippingEqualMaterials(fLastLocatedPointLocal,
850  localDirection,
851  pCurrentProposedStepLength,
852  pNewSafety,
853  fHistory,
854  fValidExitNormal,
855  fExitNormal,
856  fExiting,
857  fEntering,
858  &fBlockedPhysicalVolume,
859  fBlockedReplicaNo,
860  motherPhysical);
861  }
862  }
863  }
864  break;
865  case kParameterised:
866  if( GetDaughtersRegularStructureId(motherLogical) != 1 )
867  {
868  Step = fparamNav.ComputeStep(fLastLocatedPointLocal,
869  localDirection,
870  pCurrentProposedStepLength,
871  pNewSafety,
872  fHistory,
873  fValidExitNormal,
874  fExitNormal,
875  fExiting,
876  fEntering,
877  &fBlockedPhysicalVolume,
878  fBlockedReplicaNo);
879  }
880  else // Regular structure
881  {
882  Step = fregularNav.ComputeStep(fLastLocatedPointLocal,
883  localDirection,
884  pCurrentProposedStepLength,
885  pNewSafety,
886  fHistory,
887  fValidExitNormal,
888  fExitNormal,
889  fExiting,
890  fEntering,
891  &fBlockedPhysicalVolume,
892  fBlockedReplicaNo);
893  }
894  break;
895  case kReplica:
896  G4Exception("G4ITNavigator::ComputeStep()", "GeomNav0001",
897  FatalException, "Not applicable for replicated volumes.");
898  break;
899  }
900  }
901  else
902  {
903  // In the case of a replica, it must handle the exiting
904  // edge/corner problem by itself
905  //
906  G4bool exitingReplica = fExitedMother;
907  G4bool calculatedExitNormal= false;
908 
909  Step = freplicaNav.ComputeStep(pGlobalpoint,
910  pDirection,
911  fLastLocatedPointLocal,
912  localDirection,
913  pCurrentProposedStepLength,
914  pNewSafety,
915  fHistory,
916  fValidExitNormal,
917  calculatedExitNormal,
918  fExitNormal,
919  exitingReplica,
920  fEntering,
921  &fBlockedPhysicalVolume,
922  fBlockedReplicaNo);
923  fExiting= exitingReplica; // still ok to set it ??
924  }
925 
926  // Remember last safety origin & value.
927  //
928  fPreviousSftOrigin = pGlobalpoint;
929  fPreviousSafety = pNewSafety;
930 
931  // Count zero steps - one can occur due to changing momentum at a boundary
932  // - one, two (or a few) can occur at common edges between
933  // volumes
934  // - more than two is likely a problem in the geometry
935  // description or the Navigation
936 
937  // Rule of thumb: likely at an Edge if two consecutive steps are zero,
938  // because at least two candidate volumes must have been
939  // checked
940  //
941  fLocatedOnEdge = fLastStepWasZero && (Step==0.0);
942  fLastStepWasZero = (Step==0.0);
943  if (fPushed) fPushed = fLastStepWasZero;
944 
945  // Handle large number of consecutive zero steps
946  //
947  if ( fLastStepWasZero )
948  {
949  fNumberZeroSteps++;
950 #ifdef G4DEBUG_NAVIGATION
951  if( fNumberZeroSteps > 1 )
952  {
953  G4cout << "G4ITNavigator::ComputeStep(): another zero step, # "
954  << fNumberZeroSteps
955  << " at " << pGlobalpoint
956  << " in volume " << motherPhysical->GetName()
957  << " nav-comp-step calls # " << sNavCScalls
958  << G4endl;
959  }
960 #endif
961  if( fNumberZeroSteps > fActionThreshold_NoZeroSteps-1 )
962  {
963  // Act to recover this stuck track. Pushing it along direction
964  //
965  Step += 100*kCarTolerance;
966 #ifdef G4VERBOSE
967  if ((!fPushed) && (fWarnPush))
968  {
969  std::ostringstream message;
970  message << "Track stuck or not moving." << G4endl
971  << " Track stuck, not moving for "
972  << fNumberZeroSteps << " steps" << G4endl
973  << " in volume -" << motherPhysical->GetName()
974  << "- at point " << pGlobalpoint << G4endl
975  << " direction: " << pDirection << "." << G4endl
976  << " Potential geometry or navigation problem !"
977  << G4endl
978  << " Trying pushing it of " << Step << " mm ...";
979  G4Exception("G4ITNavigator::ComputeStep()", "GeomNav1002",
980  JustWarning, message, "Potential overlap in geometry!");
981  }
982 #endif
983  fPushed = true;
984  }
985  if( fNumberZeroSteps > fAbandonThreshold_NoZeroSteps-1 )
986  {
987  // Must kill this stuck track
988  //
989  std::ostringstream message;
990  message << "Stuck Track: potential geometry or navigation problem."
991  << G4endl
992  << " Track stuck, not moving for "
993  << fNumberZeroSteps << " steps" << G4endl
994  << " in volume -" << motherPhysical->GetName()
995  << "- at point " << pGlobalpoint << G4endl
996  << " direction: " << pDirection << ".";
997  motherPhysical->CheckOverlaps(5000, false);
998  G4Exception("G4ITNavigator::ComputeStep()", "GeomNav0003",
999  EventMustBeAborted, message);
1000  }
1001  }
1002  else
1003  {
1004  if (!fPushed) fNumberZeroSteps = 0;
1005  }
1006 
1007  fEnteredDaughter = fEntering; // I expect to enter a volume in this Step
1008  fExitedMother = fExiting;
1009 
1010  fStepEndPoint = pGlobalpoint + Step * pDirection;
1011  fLastStepEndPointLocal = fLastLocatedPointLocal + Step * localDirection;
1012 
1013  if( fExiting )
1014  {
1015 #ifdef G4DEBUG_NAVIGATION
1016  if( fVerbose > 2 )
1017  {
1018  G4cout << " At G4Nav CompStep End - if(exiting) - fExiting= " << fExiting
1019  << " fValidExitNormal = " << fValidExitNormal << G4endl;
1020  G4cout << " fExitNormal= " << fExitNormal << G4endl;
1021  }
1022 #endif
1023 
1024  if(fValidExitNormal)
1025  {
1026  // Convention: fExitNormal is in the 'grand-mother' coordinate system
1027  //
1028  fGrandMotherExitNormal= fExitNormal;
1029  }
1030  else
1031  {
1032  // We must calculate the normal anyway (in order to have it if requested)
1033  //
1034  G4ThreeVector finalLocalPoint =
1035  fLastLocatedPointLocal + localDirection*Step;
1036 
1037  // Now fGrandMotherExitNormal is in the 'grand-mother' coordinate system
1038  //
1039  fGrandMotherExitNormal =
1040  motherLogical->GetSolid()->SurfaceNormal(finalLocalPoint);
1041 
1042  const G4RotationMatrix* mRot = motherPhysical->GetRotation();
1043  if( mRot )
1044  {
1045  fGrandMotherExitNormal *= (*mRot).inverse();
1046  }
1047  // Do not set fValidExitNormal -- this signifies that the solid is convex!
1048  }
1049  }
1050  fStepEndPoint= pGlobalpoint+Step*pDirection;
1051 
1052  if( (Step == pCurrentProposedStepLength) && (!fExiting) && (!fEntering) )
1053  {
1054  // This if Step is not really limited by the geometry.
1055  // The Navigator is obliged to return "infinity"
1056  //
1057  Step = kInfinity;
1058  }
1059 
1060 #ifdef G4VERBOSE
1061  if( fVerbose > 1 )
1062  {
1063  if( fVerbose >= 4 )
1064  {
1065  G4cout << " ----- Upon exiting :" << G4endl;
1066  PrintState();
1067  }
1068  G4cout <<" Returned step = " << Step << G4endl;
1069  if( Step == kInfinity )
1070  {
1071  G4cout << " Original proposed step = "
1072  << pCurrentProposedStepLength << G4endl;
1073  }
1074  G4cout << " Safety = " << pNewSafety << G4endl;
1075  }
1076 #endif
1077 
1078  return Step;
1079 }
1080 
1081 // ********************************************************************
1082 // CheckNextStep
1083 //
1084 // Compute the step without altering the navigator state
1085 // ********************************************************************
1086 //
1088  const G4ThreeVector& pDirection,
1089  const G4double pCurrentProposedStepLength,
1090  G4double& pNewSafety)
1091 {
1092  G4double step;
1093 
1094  // Save the state, for this parasitic call
1095  //
1096  SetSavedState();
1097 
1098  step = ComputeStep ( pGlobalpoint,
1099  pDirection,
1100  pCurrentProposedStepLength,
1101  pNewSafety );
1102 
1103  // If a parasitic call, then attempt to restore the key parts of the state
1104  //
1105  RestoreSavedState();
1106 
1107  return step;
1108 }
1109 
1110 // ********************************************************************
1111 // ResetState
1112 //
1113 // Resets stack and minimum of navigator state `machine'
1114 // ********************************************************************
1115 //
1117 {
1118  fWasLimitedByGeometry = false;
1119  fEntering = false;
1120  fExiting = false;
1121  fLocatedOnEdge = false;
1122  fLastStepWasZero = false;
1123  fEnteredDaughter = false;
1124  fExitedMother = false;
1125  fPushed = false;
1126 
1127  fValidExitNormal = false;
1128  fExitNormal = G4ThreeVector(0,0,0);
1129 
1130  fPreviousSftOrigin = G4ThreeVector(0,0,0);
1131  fPreviousSafety = 0.0;
1132 
1133  fNumberZeroSteps = 0;
1134 
1135  fBlockedPhysicalVolume = 0;
1136  fBlockedReplicaNo = -1;
1137 
1138  fLastLocatedPointLocal = G4ThreeVector( kInfinity, -kInfinity, 0.0 );
1139  fLocatedOutsideWorld = false;
1140 }
1141 
1142 // ********************************************************************
1143 // SetupHierarchy
1144 //
1145 // Renavigates & resets hierarchy described by current history
1146 // o Reset volumes
1147 // o Recompute transforms and/or solids of replicated/parameterised volumes
1148 // ********************************************************************
1149 //
1151 {
1152  G4int i;
1153  const G4int cdepth = fHistory.GetDepth();
1154  G4VPhysicalVolume *current;
1155  G4VSolid *pSolid;
1156  G4VPVParameterisation *pParam;
1157 
1158  for ( i=1; i<=cdepth; i++ )
1159  {
1160  current = fHistory.GetVolume(i);
1161  switch ( fHistory.GetVolumeType(i) )
1162  {
1163  case kNormal:
1164  break;
1165  case kReplica:
1166  freplicaNav.ComputeTransformation(fHistory.GetReplicaNo(i), current);
1167  break;
1168  case kParameterised:
1169  G4int replicaNo;
1170  pParam = current->GetParameterisation();
1171  replicaNo = fHistory.GetReplicaNo(i);
1172  pSolid = pParam->ComputeSolid(replicaNo, current);
1173 
1174  // Set up dimensions & transform in solid/physical volume
1175  //
1176  pSolid->ComputeDimensions(pParam, replicaNo, current);
1177  pParam->ComputeTransformation(replicaNo, current);
1178 
1179  G4TouchableHistory touchable( fHistory );
1180  touchable.MoveUpHistory(); // move up to the parent level
1181 
1182  // Set up the correct solid and material in Logical Volume
1183  //
1184  G4LogicalVolume *pLogical = current->GetLogicalVolume();
1185  pLogical->SetSolid( pSolid );
1186  pLogical->UpdateMaterial( pParam ->
1187  ComputeMaterial(replicaNo, current, &touchable) );
1188  break;
1189  }
1190  }
1191 }
1192 
1193 // ********************************************************************
1194 // GetLocalExitNormal
1195 //
1196 // Obtains the Normal vector to a surface (in local coordinates)
1197 // pointing out of previous volume and into current volume
1198 // ********************************************************************
1199 //
1201 {
1202  G4ThreeVector ExitNormal(0.,0.,0.);
1203  G4VSolid *currentSolid=0;
1204  G4LogicalVolume *candidateLogical;
1205  if ( fLastTriedStepComputation )
1206  {
1207  // use fLastLocatedPointLocal
1208  // and next candidate volume
1209  G4ThreeVector nextSolidExitNormal(0.,0.,0.);
1210 
1211  if( fEntering && (fBlockedPhysicalVolume!=0) )
1212  {
1213  candidateLogical= fBlockedPhysicalVolume->GetLogicalVolume();
1214  if( candidateLogical )
1215  {
1216  // fLastStepEndPointLocal is in the coordinates of the mother
1217  // we need it in the daughter's coordinate system.
1218 
1219  if( CharacteriseDaughters(candidateLogical) != kReplica )
1220  {
1221  // First transform fLastLocatedPointLocal to the new daughter
1222  // coordinates
1223  G4AffineTransform MotherToDaughterTransform=
1224  GetMotherToDaughterTransform( fBlockedPhysicalVolume,
1225  fBlockedReplicaNo,
1226  VolumeType(fBlockedPhysicalVolume) );
1227  G4ThreeVector daughterPointOwnLocal=
1228  MotherToDaughterTransform.TransformPoint( fLastStepEndPointLocal );
1229 
1230  // OK if it is a parameterised volume
1231  //
1232  EInside inSideIt;
1233  G4bool onSurface;
1234  G4double safety= -1.0;
1235  currentSolid= candidateLogical->GetSolid();
1236  inSideIt = currentSolid->Inside(daughterPointOwnLocal);
1237  onSurface = (inSideIt == kSurface);
1238  if( ! onSurface )
1239  {
1240  if( inSideIt == kOutside )
1241  {
1242  safety = (currentSolid->DistanceToIn(daughterPointOwnLocal));
1243  onSurface = safety < 100.0 * kCarTolerance;
1244  }
1245  else if (inSideIt == kInside )
1246  {
1247  safety = (currentSolid->DistanceToOut(daughterPointOwnLocal));
1248  onSurface = safety < 100.0 * kCarTolerance;
1249  }
1250  }
1251 
1252  if( onSurface )
1253  {
1254  nextSolidExitNormal =
1255  currentSolid->SurfaceNormal(daughterPointOwnLocal);
1256 
1257  // Entering the solid ==> opposite
1258  //
1259  ExitNormal = -nextSolidExitNormal;
1260  }
1261  else
1262  {
1263 #ifdef G4VERBOSE
1264  if(( fVerbose == 1 ) && ( fCheck ))
1265  {
1266  std::ostringstream message;
1267  message << "Point not on surface ! " << G4endl
1268  << " Point = "
1269  << daughterPointOwnLocal << G4endl
1270  << " Physical volume = "
1271  << fBlockedPhysicalVolume->GetName() << G4endl
1272  << " Logical volume = "
1273  << candidateLogical->GetName() << G4endl
1274  << " Solid = " << currentSolid->GetName()
1275  << " Type = "
1276  << currentSolid->GetEntityType() << G4endl
1277  << *currentSolid << G4endl;
1278  if( inSideIt == kOutside )
1279  {
1280  message << "Point is Outside. " << G4endl
1281  << " Safety (from outside) = " << safety << G4endl;
1282  }
1283  else // if( inSideIt == kInside )
1284  {
1285  message << "Point is Inside. " << G4endl
1286  << " Safety (from inside) = " << safety << G4endl;
1287  }
1288  G4Exception("G4ITNavigator::GetLocalExitNormal()", "GeomNav1001",
1289  JustWarning, message);
1290  }
1291 #endif
1292  }
1293  *valid = onSurface; // was =true;
1294  }
1295  else
1296  {
1297  *valid = false; // TODO: Need Separate code for replica!!!!
1298 #ifdef G4DEBUG_NAVIGATION
1299  G4Exception("G4ITNavigator::GetLocalExitNormal()", "GeomNav0001",
1300  FatalException,
1301  "Local normal not (yet) available for replica volumes.");
1302 #endif
1303  }
1304  }
1305  }
1306  else if ( fExiting )
1307  {
1308  ExitNormal = fGrandMotherExitNormal;
1309  *valid = true;
1310  }
1311  else // ie ( fBlockedPhysicalVolume == 0 )
1312  {
1313  *valid = false;
1314  }
1315  }
1316  else
1317  {
1318  if ( EnteredDaughterVolume() )
1319  {
1320  ExitNormal= -(fHistory.GetTopVolume()->GetLogicalVolume()->
1321  GetSolid()->SurfaceNormal(fLastLocatedPointLocal));
1322  *valid = true;
1323  }
1324  else
1325  {
1326  if( fExitedMother )
1327  {
1328  ExitNormal = fGrandMotherExitNormal;
1329  *valid = true;
1330  }
1331  else // We are not at a boundary. ExitNormal remains (0,0,0)
1332  {
1333  *valid = false;
1334  }
1335  }
1336  }
1337  return ExitNormal;
1338 }
1339 
1340 // ********************************************************************
1341 // GetMotherToDaughterTransform
1342 //
1343 // Obtains the mother to daughter affine transformation
1344 // ********************************************************************
1345 //
1348  G4int enteringReplicaNo,
1349  EVolume enteringVolumeType )
1350 {
1351  switch (enteringVolumeType)
1352  {
1353  case kNormal: // Nothing is needed to prepare the transformation
1354  break; // It is stored already in the physical volume (placement)
1355  case kReplica: // Sets the transform in the Replica - tbc
1356  G4Exception("G4ITNavigator::GetMotherToDaughterTransform()",
1357  "GeomNav0001", FatalException,
1358  "Method NOT Implemented yet for replica volumes.");
1359  break;
1360  case kParameterised:
1361  if( pEnteringPhysVol->GetRegularStructureId() == 0 )
1362  {
1363  G4VPVParameterisation *pParam =
1364  pEnteringPhysVol->GetParameterisation();
1365  G4VSolid* pSolid =
1366  pParam->ComputeSolid(enteringReplicaNo, pEnteringPhysVol);
1367  pSolid->ComputeDimensions(pParam, enteringReplicaNo, pEnteringPhysVol);
1368 
1369  // Sets the transform in the Parameterisation
1370  //
1371  pParam->ComputeTransformation(enteringReplicaNo, pEnteringPhysVol);
1372 
1373  // Set the correct solid and material in Logical Volume
1374  //
1375  G4LogicalVolume* pLogical = pEnteringPhysVol->GetLogicalVolume();
1376  pLogical->SetSolid( pSolid );
1377  }
1378  break;
1379  }
1380  return G4AffineTransform(pEnteringPhysVol->GetRotation(),
1381  pEnteringPhysVol->GetTranslation()).Invert();
1382 }
1383 
1384 // ********************************************************************
1385 // GetLocalExitNormalAndCheck
1386 //
1387 // Obtains the Normal vector to a surface (in local coordinates)
1388 // pointing out of previous volume and into current volume, and
1389 // checks the current point against expected 'local' value.
1390 // ********************************************************************
1391 //
1393 GetLocalExitNormalAndCheck(const G4ThreeVector& ExpectedBoundaryPointGlobal,
1394  G4bool* pValid)
1395 {
1396  G4ThreeVector ExpectedBoundaryPointLocal;
1397 
1398  // Check Current point against expected 'local' value
1399  //
1400  if ( fLastTriedStepComputation )
1401  {
1402  const G4AffineTransform& GlobalToLocal= GetGlobalToLocalTransform();
1403  ExpectedBoundaryPointLocal =
1404  GlobalToLocal.TransformPoint( ExpectedBoundaryPointGlobal );
1405  }
1406 
1407  return GetLocalExitNormal( pValid);
1408 }
1409 
1410 // ********************************************************************
1411 // GetGlobalExitNormal
1412 //
1413 // Obtains the Normal vector to a surface (in global coordinates)
1414 // pointing out of previous volume and into current volume
1415 // ********************************************************************
1416 //
1419  G4bool* pValidNormal)
1420 {
1421  G4bool validNormal;
1422  G4ThreeVector localNormal, globalNormal;
1423 
1424  localNormal = GetLocalExitNormalAndCheck( IntersectPointGlobal, &validNormal);
1425  *pValidNormal = validNormal;
1426  G4AffineTransform localToGlobal = GetLocalToGlobalTransform();
1427  globalNormal = localToGlobal.TransformAxis( localNormal );
1428 
1429  return globalNormal;
1430 }
1431 
1432 // ********************************************************************
1433 // ComputeSafety
1434 //
1435 // It assumes that it will be
1436 // i) called at the Point in the same volume as the EndPoint of the
1437 // ComputeStep.
1438 // ii) after (or at the end of) ComputeStep OR after the relocation.
1439 // ********************************************************************
1440 //
1442  const G4double pMaxLength,
1443  const G4bool keepState)
1444 {
1445  G4double newSafety = 0.0;
1446 
1447 #ifdef G4DEBUG_NAVIGATION
1448  G4int oldcoutPrec = G4cout.precision(8);
1449  if( fVerbose > 0 )
1450  {
1451  G4cout << "*** G4ITNavigator::ComputeSafety: ***" << G4endl
1452  << " Called at point: " << pGlobalpoint << G4endl;
1453 
1454  G4VPhysicalVolume *motherPhysical = fHistory.GetTopVolume();
1455  G4cout << " Volume = " << motherPhysical->GetName()
1456  << " - Maximum length = " << pMaxLength << G4endl;
1457  if( fVerbose >= 4 )
1458  {
1459  G4cout << " ----- Upon entering Compute Safety:" << G4endl;
1460  PrintState();
1461  }
1462  }
1463 #endif
1464 
1465  if (keepState) { SetSavedState(); }
1466  // fLastTriedStepComputation= true; -- this method is NOT computing the Step size
1467 
1468  G4double distEndpointSq = (pGlobalpoint-fStepEndPoint).mag2();
1469  G4bool stayedOnEndpoint = distEndpointSq < kCarTolerance*kCarTolerance;
1470  G4bool endpointOnSurface = fEnteredDaughter || fExitedMother;
1471 
1472  if( !(endpointOnSurface && stayedOnEndpoint) )
1473  {
1474  // Pseudo-relocate to this point (updates voxel information only)
1475  //
1476  LocateGlobalPointWithinVolume( pGlobalpoint );
1477  // --->> Danger: Side effects on sub-navigator voxel information <<---
1478  // Could be replaced again by 'granular' calls to sub-navigator
1479  // locates (similar side-effects, but faster.
1480  // Solutions:
1481  // 1) Re-locate (to where?)
1482  // 2) Insure that the methods using (G4ComputeStep?)
1483  // does a relocation (if information is disturbed only ?)
1484 
1485 #ifdef G4DEBUG_NAVIGATION
1486  if( fVerbose >= 2 )
1487  {
1488  G4cout << " G4ITNavigator::ComputeSafety() relocates-in-volume to point: "
1489  << pGlobalpoint << G4endl;
1490  }
1491 #endif
1492  G4VPhysicalVolume *motherPhysical = fHistory.GetTopVolume();
1493  G4LogicalVolume *motherLogical = motherPhysical->GetLogicalVolume();
1494  G4SmartVoxelHeader* pVoxelHeader = motherLogical->GetVoxelHeader();
1495  G4ThreeVector localPoint = ComputeLocalPoint(pGlobalpoint);
1496 
1498  {
1499  switch(CharacteriseDaughters(motherLogical))
1500  {
1501  case kNormal:
1502  if ( pVoxelHeader )
1503  {
1504  newSafety=fvoxelNav.ComputeSafety(localPoint,fHistory,pMaxLength);
1505  }
1506  else
1507  {
1508  newSafety=fnormalNav.ComputeSafety(localPoint,fHistory,pMaxLength);
1509  }
1510  break;
1511  case kParameterised:
1512  if( GetDaughtersRegularStructureId(motherLogical) != 1 )
1513  {
1514  newSafety = fparamNav.ComputeSafety(localPoint,fHistory,pMaxLength);
1515  }
1516  else // Regular structure
1517  {
1518  newSafety = fregularNav.ComputeSafety(localPoint,fHistory,pMaxLength);
1519  }
1520  break;
1521  case kReplica:
1522  G4Exception("G4ITNavigator::ComputeSafety()", "NotApplicable",
1523  FatalException, "Not applicable for replicated volumes.");
1524  break;
1525  }
1526  }
1527  else
1528  {
1529  newSafety = freplicaNav.ComputeSafety(pGlobalpoint, localPoint,
1530  fHistory, pMaxLength);
1531  }
1532  }
1533  else // if( endpointOnSurface && stayedOnEndpoint )
1534  {
1535 #ifdef G4DEBUG_NAVIGATION
1536  if( fVerbose >= 2 )
1537  {
1538  G4cout << " G4ITNavigator::ComputeSafety() finds that point - "
1539  << pGlobalpoint << " - is on surface " << G4endl;
1540  if( fEnteredDaughter ) { G4cout << " entered new daughter volume"; }
1541  if( fExitedMother ) { G4cout << " and exited previous volume."; }
1542  G4cout << G4endl;
1543  G4cout << " EndPoint was = " << fStepEndPoint << G4endl;
1544  }
1545 #endif
1546  newSafety = 0.0;
1547  }
1548 
1549  // Remember last safety origin & value
1550  //
1551  fPreviousSftOrigin = pGlobalpoint;
1552  fPreviousSafety = newSafety;
1553 
1554  if (keepState) { RestoreSavedState(); }
1555 
1556 #ifdef G4DEBUG_NAVIGATION
1557  if( fVerbose > 1 )
1558  {
1559  G4cout << " ---- Exiting ComputeSafety " << G4endl;
1560  if( fVerbose > 2 ) { PrintState(); }
1561  G4cout << " Returned value of Safety = " << newSafety << G4endl;
1562  }
1563  G4cout.precision(oldcoutPrec);
1564 #endif
1565 
1566  return newSafety;
1567 }
1568 
1569 // ********************************************************************
1570 // CreateTouchableHistoryHandle
1571 // ********************************************************************
1572 //
1574 {
1576 }
1577 
1578 // ********************************************************************
1579 // PrintState
1580 // ********************************************************************
1581 //
1583 {
1584  G4int oldcoutPrec = G4cout.precision(4);
1585  if( fVerbose == 4 )
1586  {
1587  G4cout << "The current state of G4ITNavigator is: " << G4endl;
1588  G4cout << " ValidExitNormal= " << fValidExitNormal << G4endl
1589  << " ExitNormal = " << fExitNormal << G4endl
1590  << " Exiting = " << fExiting << G4endl
1591  << " Entering = " << fEntering << G4endl
1592  << " BlockedPhysicalVolume= " ;
1593  if (fBlockedPhysicalVolume==0)
1594  G4cout << "None";
1595  else
1596  G4cout << fBlockedPhysicalVolume->GetName();
1597  G4cout << G4endl
1598  << " BlockedReplicaNo = " << fBlockedReplicaNo << G4endl
1599  << " LastStepWasZero = " << fLastStepWasZero << G4endl
1600  << G4endl;
1601  }
1602  if( ( 1 < fVerbose) && (fVerbose < 4) )
1603  {
1604  G4cout << std::setw(30) << " ExitNormal " << " "
1605  << std::setw( 5) << " Valid " << " "
1606  << std::setw( 9) << " Exiting " << " "
1607  << std::setw( 9) << " Entering" << " "
1608  << std::setw(15) << " Blocked:Volume " << " "
1609  << std::setw( 9) << " ReplicaNo" << " "
1610  << std::setw( 8) << " LastStepZero " << " "
1611  << G4endl;
1612  G4cout << "( " << std::setw(7) << fExitNormal.x()
1613  << ", " << std::setw(7) << fExitNormal.y()
1614  << ", " << std::setw(7) << fExitNormal.z() << " ) "
1615  << std::setw( 5) << fValidExitNormal << " "
1616  << std::setw( 9) << fExiting << " "
1617  << std::setw( 9) << fEntering << " ";
1618  if ( fBlockedPhysicalVolume==0 )
1619  G4cout << std::setw(15) << "None";
1620  else
1621  G4cout << std::setw(15)<< fBlockedPhysicalVolume->GetName();
1622  G4cout << std::setw( 9) << fBlockedReplicaNo << " "
1623  << std::setw( 8) << fLastStepWasZero << " "
1624  << G4endl;
1625  }
1626  if( fVerbose > 2 )
1627  {
1628  G4cout.precision(8);
1629  G4cout << " Current Localpoint = " << fLastLocatedPointLocal << G4endl;
1630  G4cout << " PreviousSftOrigin = " << fPreviousSftOrigin << G4endl;
1631  G4cout << " PreviousSafety = " << fPreviousSafety << G4endl;
1632  }
1633  G4cout.precision(oldcoutPrec);
1634 }
1635 
1636 // ********************************************************************
1637 // ComputeStepLog
1638 // ********************************************************************
1639 //
1640 void G4ITNavigator::ComputeStepLog(const G4ThreeVector& pGlobalpoint,
1641  G4double moveLenSq) const
1642 {
1643  // The following checks only make sense if the move is larger
1644  // than the tolerance.
1645 
1646  static const G4double fAccuracyForWarning = kCarTolerance,
1647  fAccuracyForException = 1000*kCarTolerance;
1648 
1649  G4ThreeVector OriginalGlobalpoint = fHistory.GetTopTransform().Inverse().
1650  TransformPoint(fLastLocatedPointLocal);
1651 
1652  G4double shiftOriginSafSq = (fPreviousSftOrigin-pGlobalpoint).mag2();
1653 
1654  // Check that the starting point of this step is
1655  // within the isotropic safety sphere of the last point
1656  // to a accuracy/precision given by fAccuracyForWarning.
1657  // If so give warning.
1658  // If it fails by more than fAccuracyForException exit with error.
1659  //
1660  if( shiftOriginSafSq >= sqr(fPreviousSafety) )
1661  {
1662  G4double shiftOrigin = std::sqrt(shiftOriginSafSq);
1663  G4double diffShiftSaf = shiftOrigin - fPreviousSafety;
1664 
1665  if( diffShiftSaf > fAccuracyForWarning )
1666  {
1667  G4int oldcoutPrec= G4cout.precision(8);
1668  G4int oldcerrPrec= G4cerr.precision(10);
1669  std::ostringstream message, suggestion;
1670  message << "Accuracy error or slightly inaccurate position shift."
1671  << G4endl
1672  << " The Step's starting point has moved "
1673  << std::sqrt(moveLenSq)/mm << " mm " << G4endl
1674  << " since the last call to a Locate method." << G4endl
1675  << " This has resulted in moving "
1676  << shiftOrigin/mm << " mm "
1677  << " from the last point at which the safety "
1678  << " was calculated " << G4endl
1679  << " which is more than the computed safety= "
1680  << fPreviousSafety/mm << " mm at that point." << G4endl
1681  << " This difference is "
1682  << diffShiftSaf/mm << " mm." << G4endl
1683  << " The tolerated accuracy is "
1684  << fAccuracyForException/mm << " mm.";
1685 
1686  suggestion << " ";
1687  static G4int warnNow = 0;
1688  if( ((++warnNow % 100) == 1) )
1689  {
1690  message << G4endl
1691  << " This problem can be due to either " << G4endl
1692  << " - a process that has proposed a displacement"
1693  << " larger than the current safety , or" << G4endl
1694  << " - inaccuracy in the computation of the safety";
1695  suggestion << "We suggest that you " << G4endl
1696  << " - find i) what particle is being tracked, and "
1697  << " ii) through what part of your geometry " << G4endl
1698  << " for example by re-running this event with "
1699  << G4endl
1700  << " /tracking/verbose 1 " << G4endl
1701  << " - check which processes you declare for"
1702  << " this particle (and look at non-standard ones)"
1703  << G4endl
1704  << " - in case, create a detailed logfile"
1705  << " of this event using:" << G4endl
1706  << " /tracking/verbose 6 ";
1707  }
1708  G4Exception("G4ITNavigator::ComputeStep()",
1709  "GeomNav1002", JustWarning,
1710  message, G4String(suggestion.str()));
1711  G4cout.precision(oldcoutPrec);
1712  G4cerr.precision(oldcerrPrec);
1713  }
1714 #ifdef G4DEBUG_NAVIGATION
1715  else
1716  {
1717  G4cerr << "WARNING - G4ITNavigator::ComputeStep()" << G4endl
1718  << " The Step's starting point has moved "
1719  << std::sqrt(moveLenSq) << "," << G4endl
1720  << " which has taken it to the limit of"
1721  << " the current safety. " << G4endl;
1722  }
1723 #endif
1724  }
1725  G4double safetyPlus = fPreviousSafety + fAccuracyForException;
1726  if ( shiftOriginSafSq > sqr(safetyPlus) )
1727  {
1728  std::ostringstream message;
1729  message << "May lead to a crash or unreliable results." << G4endl
1730  << " Position has shifted considerably without"
1731  << " notifying the navigator !" << G4endl
1732  << " Tolerated safety: " << safetyPlus << G4endl
1733  << " Computed shift : " << shiftOriginSafSq;
1734  G4Exception("G4ITNavigator::ComputeStep()", "GeomNav1002",
1735  JustWarning, message);
1736  }
1737 }
1738 
1739 // ********************************************************************
1740 // Operator <<
1741 // ********************************************************************
1742 //
1743 std::ostream& operator << (std::ostream &os,const G4ITNavigator &n)
1744 {
1745  os << "Current History: " << G4endl << n.fHistory;
1746  return os;
1747 }