source: trunk/source/geometry/navigation/src/G4Navigator.cc @ 1358

Last change on this file since 1358 was 1347, checked in by garnier, 14 years ago

geant4 tag 9.4
File size: 51.0 KB
Line 
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,v 1.46 2010/11/15 14:03:27 gcosmo Exp $
28// GEANT4 tag $ Name:  $
29//
30// class G4Navigator Implementation
31//
32// Original author: Paul Kent, July 95/96
33//
34// --------------------------------------------------------------------
35
36#include "G4Navigator.hh"
37#include "G4ios.hh"
38#include <iomanip>
39
40#include "G4GeometryTolerance.hh"
41#include "G4VPhysicalVolume.hh"
42
43// ********************************************************************
44// Constructor
45// ********************************************************************
46//
47G4Navigator::G4Navigator()
48  : fWasLimitedByGeometry(false), fVerbose(0),
49    fTopPhysical(0), fCheck(false), fPushed(false), fWarnPush(true)
50{
51  fActive= false; 
52  ResetStackAndState();
53
54  fActionThreshold_NoZeroSteps  = 10; 
55  fAbandonThreshold_NoZeroSteps = 25; 
56
57  kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
58  fregularNav.SetNormalNavigation( &fnormalNav );
59
60  fStepEndPoint = G4ThreeVector( kInfinity, kInfinity, kInfinity ); 
61}
62
63// ********************************************************************
64// Destructor
65// ********************************************************************
66//
67G4Navigator::~G4Navigator()
68{;}
69
70// ********************************************************************
71// ResetHierarchyAndLocate
72// ********************************************************************
73//
74G4VPhysicalVolume*
75G4Navigator::ResetHierarchyAndLocate(const G4ThreeVector &p,
76                                     const G4ThreeVector &direction,
77                                     const G4TouchableHistory &h)
78{
79  ResetState();
80  fHistory = *h.GetHistory();
81  SetupHierarchy();
82  return LocateGlobalPointAndSetup(p, &direction, true, false);
83}
84
85// ********************************************************************
86// LocateGlobalPointAndSetup
87//
88// Locate the point in the hierarchy return 0 if outside
89// The direction is required
90//    - if on an edge shared by more than two surfaces
91//      (to resolve likely looping in tracking)
92//    - at initial location of a particle
93//      (to resolve potential ambiguity at boundary)
94//
95// Flags on exit: (comments to be completed)
96// fEntering         - True if entering `daughter' volume (or replica)
97//                     whether daughter of last mother directly
98//                     or daughter of that volume's ancestor.
99// ********************************************************************
100//
101G4VPhysicalVolume* 
102G4Navigator::LocateGlobalPointAndSetup( const G4ThreeVector& globalPoint,
103                                        const G4ThreeVector* pGlobalDirection,
104                                        const G4bool relativeSearch,
105                                        const G4bool ignoreDirection )
106{
107  G4bool notKnownContained=true, noResult;
108  G4VPhysicalVolume *targetPhysical;
109  G4LogicalVolume *targetLogical;
110  G4VSolid *targetSolid=0;
111  G4ThreeVector localPoint, globalDirection;
112  EInside insideCode;
113 
114  G4bool considerDirection = (!ignoreDirection) || fLocatedOnEdge;
115 
116  if( considerDirection && pGlobalDirection != 0 )
117  {
118    globalDirection=*pGlobalDirection;
119  }
120
121#ifdef G4DEBUG_NAVIGATION
122  if( fVerbose > 2 )
123  {
124    G4cout << "Upon entering LocateGlobalPointAndSetup():" << G4endl;
125    G4cout << "    History = " << G4endl << fHistory << G4endl << G4endl;
126  }
127#endif
128
129#ifdef G4VERBOSE
130  if( fVerbose > 2 )
131  {
132    G4int oldcoutPrec = G4cout.precision(8);
133    G4cout << "*** G4Navigator::LocateGlobalPointAndSetup: ***" << G4endl; 
134    G4cout << "    Called with arguments: " << G4endl
135           << "    Globalpoint = " << globalPoint << G4endl
136           << "    RelativeSearch = " << relativeSearch  << G4endl;
137    if( fVerbose == 4 )
138    {
139      G4cout << "    ----- Upon entering:" << G4endl;
140      PrintState();
141    }
142    G4cout.precision(oldcoutPrec);
143  }
144#endif
145
146  if ( !relativeSearch )
147  {
148    ResetStackAndState();
149  }
150  else
151  {
152    if ( fWasLimitedByGeometry )
153    {
154      fWasLimitedByGeometry = false;
155      fEnteredDaughter = fEntering;   // Remember
156      fExitedMother = fExiting;       // Remember
157      if ( fExiting )
158      {
159        if ( fHistory.GetDepth() )
160        {
161          fBlockedPhysicalVolume = fHistory.GetTopVolume();
162          fBlockedReplicaNo = fHistory.GetTopReplicaNo();
163          fHistory.BackLevel();
164        }
165        else
166        {
167          fLastLocatedPointLocal = localPoint;
168          fLocatedOutsideWorld = true;
169          return 0;           // Have exited world volume
170        }
171        // A fix for the case where a volume is "entered" at an edge
172        // and a coincident surface exists outside it.
173        //  - This stops it from exiting further volumes and cycling
174        //  - However ReplicaNavigator treats this case itself
175        //
176        if ( fLocatedOnEdge && (VolumeType(fBlockedPhysicalVolume)!=kReplica ))
177        { 
178          fExiting= false;
179        }
180      }
181      else
182        if ( fEntering )
183        {
184          switch (VolumeType(fBlockedPhysicalVolume))
185          {
186            case kNormal:
187              fHistory.NewLevel(fBlockedPhysicalVolume, kNormal,
188                                fBlockedPhysicalVolume->GetCopyNo());
189              break;
190            case kReplica:
191              freplicaNav.ComputeTransformation(fBlockedReplicaNo,
192                                                fBlockedPhysicalVolume);
193              fHistory.NewLevel(fBlockedPhysicalVolume, kReplica,
194                                fBlockedReplicaNo);
195              fBlockedPhysicalVolume->SetCopyNo(fBlockedReplicaNo);
196              break;
197            case kParameterised:
198              if( fBlockedPhysicalVolume->GetRegularStructureId() == 0 )
199              {
200                G4VSolid *pSolid;
201                G4VPVParameterisation *pParam;
202                G4TouchableHistory parentTouchable( fHistory );
203                pParam = fBlockedPhysicalVolume->GetParameterisation();
204                pSolid = pParam->ComputeSolid(fBlockedReplicaNo,
205                                              fBlockedPhysicalVolume);
206                pSolid->ComputeDimensions(pParam, fBlockedReplicaNo,
207                                          fBlockedPhysicalVolume);
208                pParam->ComputeTransformation(fBlockedReplicaNo,
209                                              fBlockedPhysicalVolume);
210                fHistory.NewLevel(fBlockedPhysicalVolume, kParameterised,
211                                  fBlockedReplicaNo);
212                fBlockedPhysicalVolume->SetCopyNo(fBlockedReplicaNo);
213                //
214                // Set the correct solid and material in Logical Volume
215                //
216                G4LogicalVolume *pLogical;
217                pLogical = fBlockedPhysicalVolume->GetLogicalVolume();
218                pLogical->SetSolid( pSolid );
219                pLogical->UpdateMaterial(pParam ->
220                  ComputeMaterial(fBlockedReplicaNo,
221                                  fBlockedPhysicalVolume, 
222                                  &parentTouchable));
223              }
224              break;
225          }
226          fEntering = false;
227          fBlockedPhysicalVolume = 0;
228          localPoint = fHistory.GetTopTransform().TransformPoint(globalPoint);
229          notKnownContained = false;
230        }
231    }
232    else
233    {
234      fBlockedPhysicalVolume = 0;
235      fEntering = false;
236      fEnteredDaughter = false;  // Full Step was not taken, did not enter
237      fExiting = false;
238      fExitedMother = false;     // Full Step was not taken, did not exit
239    }
240  }
241  //
242  // Search from top of history up through geometry until
243  // containing volume found:
244  // If on
245  // o OUTSIDE - Back up level, not/no longer exiting volumes
246  // o SURFACE and EXITING - Back up level, setting new blocking no.s
247  // else
248  // o containing volume found
249  //
250  while (notKnownContained)
251  {
252    if ( fHistory.GetTopVolumeType()!=kReplica )
253    {
254      targetSolid = fHistory.GetTopVolume()->GetLogicalVolume()->GetSolid();
255      localPoint = fHistory.GetTopTransform().TransformPoint(globalPoint);
256      insideCode = targetSolid->Inside(localPoint);
257#ifdef G4VERBOSE
258      if(( fVerbose == 1 ) && ( fCheck ))
259      {
260         G4String solidResponse = "-kInside-";
261         if (insideCode == kOutside)
262           solidResponse = "-kOutside-";
263         else if (insideCode == kSurface)
264           solidResponse = "-kSurface-";
265         G4cout << "*** G4Navigator::LocateGlobalPointAndSetup(): ***" << G4endl
266                << "    Invoked Inside() for solid: " << targetSolid->GetName()
267                << ". Solid replied: " << solidResponse << G4endl
268                << "    For local point p: " << localPoint << G4endl;
269      }
270#endif
271    }
272    else
273    {
274      insideCode = freplicaNav.BackLocate(fHistory, globalPoint, localPoint,
275                                          fExiting, notKnownContained);
276      // !CARE! if notKnownContained returns false then the point is within
277      // the containing placement volume of the replica(s). If insidecode
278      // will result in the history being backed up one level, then the
279      // local point returned is the point in the system of this new level
280    }
281    if ( insideCode==kOutside )
282    {
283      if ( fHistory.GetDepth() )
284      {
285        fBlockedPhysicalVolume = fHistory.GetTopVolume();
286        fBlockedReplicaNo = fHistory.GetTopReplicaNo();
287        fHistory.BackLevel();
288        fExiting = false;
289      }
290      else
291      {
292        fLastLocatedPointLocal = localPoint;
293        fLocatedOutsideWorld = true;
294        return 0;         // Have exited world volume
295      }
296    }
297    else
298      if ( insideCode==kSurface )
299      {
300        G4bool isExiting = fExiting;
301        if( (!fExiting)&&considerDirection )
302        {
303          // Figure out whether we are exiting this level's volume
304          // by using the direction
305          //
306          G4bool directionExiting = false;
307          G4ThreeVector localDirection =
308              fHistory.GetTopTransform().TransformAxis(globalDirection);
309          if ( fHistory.GetTopVolumeType()!=kReplica )
310          {
311            G4ThreeVector normal = targetSolid->SurfaceNormal(localPoint);
312            directionExiting = normal.dot(localDirection) > 0.0;
313            isExiting = isExiting || directionExiting;
314          }
315        }
316        if( isExiting )
317        {
318          if ( fHistory.GetDepth() )
319          {
320            fBlockedPhysicalVolume = fHistory.GetTopVolume();
321            fBlockedReplicaNo = fHistory.GetTopReplicaNo();
322            fHistory.BackLevel();
323            //
324            // Still on surface but exited volume not necessarily convex
325            //
326            fValidExitNormal = false;
327          } 
328          else
329          {
330            fLastLocatedPointLocal = localPoint;
331            fLocatedOutsideWorld = true;
332            return 0;          // Have exited world volume
333          }
334        }
335        else
336        {
337          notKnownContained=false;
338        }
339      }
340      else
341      {
342        notKnownContained=false;
343      }
344  }  // END while (notKnownContained)
345  //
346  // Search downwards until deepest containing volume found,
347  // blocking fBlockedPhysicalVolume/BlockedReplicaNum
348  //
349  // 3 Cases:
350  //
351  // o Parameterised daughters
352  //   =>Must be one G4PVParameterised daughter & voxels
353  // o Positioned daughters & voxels
354  // o Positioned daughters & no voxels
355
356  noResult = true;  // noResult should be renamed to
357                    // something like enteredLevel, as that is its meaning.
358  do
359  {
360    // Determine `type' of current mother volume
361    //
362    targetPhysical = fHistory.GetTopVolume();
363    if (!targetPhysical) { break; }
364    targetLogical = targetPhysical->GetLogicalVolume();
365    switch( CharacteriseDaughters(targetLogical) )
366    {
367      case kNormal:
368        if ( targetLogical->GetVoxelHeader() )  // use optimised navigation
369        {
370          noResult = fvoxelNav.LevelLocate(fHistory,
371                                           fBlockedPhysicalVolume,
372                                           fBlockedReplicaNo,
373                                           globalPoint,
374                                           pGlobalDirection,
375                                           considerDirection,
376                                           localPoint);
377        }
378        else                       // do not use optimised navigation
379        {
380          noResult = fnormalNav.LevelLocate(fHistory,
381                                            fBlockedPhysicalVolume,
382                                            fBlockedReplicaNo,
383                                            globalPoint,
384                                            pGlobalDirection,
385                                            considerDirection,
386                                            localPoint);
387        }
388        break;
389      case kReplica:
390        noResult = freplicaNav.LevelLocate(fHistory,
391                                           fBlockedPhysicalVolume,
392                                           fBlockedReplicaNo,
393                                           globalPoint,
394                                           pGlobalDirection,
395                                           considerDirection,
396                                           localPoint);
397        break;
398      case kParameterised:
399        if( GetDaughtersRegularStructureId(targetLogical) != 1 )
400        {
401          noResult = fparamNav.LevelLocate(fHistory,
402                                           fBlockedPhysicalVolume,
403                                           fBlockedReplicaNo,
404                                           globalPoint,
405                                           pGlobalDirection,
406                                           considerDirection,
407                                           localPoint);
408        }
409        else  // Regular structure
410        {
411          noResult = fregularNav.LevelLocate(fHistory,
412                                             fBlockedPhysicalVolume,
413                                             fBlockedReplicaNo,
414                                             globalPoint,
415                                             pGlobalDirection,
416                                             considerDirection,
417                                             localPoint);
418        }
419        break;
420    }
421
422    // LevelLocate returns true if it finds a daughter volume
423    // in which globalPoint is inside (or on the surface).
424
425    if ( noResult )
426    {
427      // Entering a daughter after ascending
428      //
429      // The blocked volume is no longer valid - it was for another level
430      //
431      fBlockedPhysicalVolume = 0;
432      fBlockedReplicaNo = -1;
433
434      // fEntering should be false -- else blockedVolume is assumed good.
435      // fEnteredDaughter is used for ExitNormal
436      //
437      fEntering = false;
438      fEnteredDaughter = true;
439#ifdef G4DEBUG_NAVIGATION
440      if( fVerbose > 2 )
441      { 
442         G4VPhysicalVolume* enteredPhysical = fHistory.GetTopVolume();
443         G4cout << "*** G4Navigator::LocateGlobalPointAndSetup() ***" << G4endl; 
444         G4cout << "    Entering volume: " << enteredPhysical->GetName()
445                << G4endl;
446      }
447#endif
448    }
449  } while (noResult);
450
451  fLastLocatedPointLocal = localPoint;
452
453#ifdef G4VERBOSE
454  if( fVerbose == 4 )
455  {
456    G4int oldcoutPrec = G4cout.precision(8);
457    G4String curPhysVol_Name("None");
458    if (targetPhysical)  { curPhysVol_Name = targetPhysical->GetName(); }
459    G4cout << "    Return value = new volume = " << curPhysVol_Name << G4endl;
460    G4cout << "    ----- Upon exiting:" << G4endl;
461    PrintState();
462#ifdef G4DEBUG_NAVIGATION
463    G4cout << "Upon exiting LocateGlobalPointAndSetup():" << G4endl;
464    G4cout << "    History = " << G4endl << fHistory << G4endl << G4endl;
465#endif
466    G4cout.precision(oldcoutPrec);
467  }
468#endif
469
470  fLocatedOutsideWorld= false;
471
472  return targetPhysical;
473}
474
475// ********************************************************************
476// LocateGlobalPointWithinVolume
477//
478// -> the state information of this Navigator and its subNavigators
479//    is updated in order to start the next step at pGlobalpoint
480// -> no check is performed whether pGlobalpoint is inside the
481//    original volume (this must be the case).
482//
483// Note: a direction could be added to the arguments, to aid in future
484//       optional checking (via the old code below, flagged by OLD_LOCATE).
485//       [ This would be done only in verbose mode ]
486// ********************************************************************
487//
488void
489G4Navigator::LocateGlobalPointWithinVolume(const G4ThreeVector& pGlobalpoint)
490{ 
491   fLastLocatedPointLocal = ComputeLocalPoint(pGlobalpoint);
492
493#ifdef G4DEBUG_NAVIGATION
494   if( fVerbose > 2 )
495   { 
496     G4cout << "Entering LocateGlobalWithinVolume(): History = " << G4endl;
497     G4cout << fHistory << G4endl;
498   }
499#endif
500
501   // For the case of Voxel (or Parameterised) volume the respective
502   // Navigator must be messaged to update its voxel information etc
503
504   // Update the state of the Sub Navigators
505   // - in particular any voxel information they store/cache
506   //
507   G4VPhysicalVolume*  motherPhysical = fHistory.GetTopVolume();
508   G4LogicalVolume*    motherLogical  = motherPhysical->GetLogicalVolume();
509   G4SmartVoxelHeader* pVoxelHeader   = motherLogical->GetVoxelHeader();
510
511   if ( fHistory.GetTopVolumeType()!=kReplica )
512   {
513     switch( CharacteriseDaughters(motherLogical) )
514     {
515       case kNormal:
516         if ( pVoxelHeader )
517         {
518           fvoxelNav.VoxelLocate( pVoxelHeader, fLastLocatedPointLocal );
519         }
520         break;
521       case kParameterised:
522         if( GetDaughtersRegularStructureId(motherLogical) != 1 )
523         {
524           // Resets state & returns voxel node
525           //
526           fparamNav.ParamVoxelLocate( pVoxelHeader, fLastLocatedPointLocal );
527         }
528         break;
529       case kReplica:
530         G4Exception("G4Navigator::LocateGlobalPointWithinVolume()",
531                     "NotApplicable", FatalException,
532                     "Not applicable for replicated volumes.");
533         break;
534     }
535   }
536
537   // Reset the state variables
538   //   - which would have been affected
539   //     by the 'equivalent' call to LocateGlobalPointAndSetup
540   //   - who's values have been invalidated by the 'move'.
541   //
542   fBlockedPhysicalVolume = 0; 
543   fBlockedReplicaNo = -1;
544   fEntering = false;
545   fEnteredDaughter = false;  // Boundary not encountered, did not enter
546   fExiting = false;
547   fExitedMother = false;     // Boundary not encountered, did not exit
548}
549
550// ********************************************************************
551// SetSavedState
552//
553// Save the state, in case this is a parasitic call
554// Save fValidExitNormal, fExitNormal, fExiting, fEntering,
555//      fBlockedPhysicalVolume, fBlockedReplicaNo, fLastStepWasZero;
556// ********************************************************************
557//
558void G4Navigator::SetSavedState()
559{
560  // fSaveExitNormal = fExitNormal;
561  fSaveState.sExitNormal = fExitNormal;
562  fSaveState.sValidExitNormal = fValidExitNormal;
563  fSaveState.sExiting = fExiting;
564  fSaveState.sEntering = fEntering;
565
566  fSaveState.spBlockedPhysicalVolume = fBlockedPhysicalVolume;
567  fSaveState.sBlockedReplicaNo = fBlockedReplicaNo, 
568
569  fSaveState.sLastStepWasZero = fLastStepWasZero; 
570}
571
572// ********************************************************************
573// RestoreSavedState
574//
575// Restore the state (in Compute Step), in case this is a parasitic call
576// ********************************************************************
577//
578void G4Navigator::RestoreSavedState()
579{
580  fExitNormal = fSaveState.sExitNormal;
581  fValidExitNormal = fSaveState.sValidExitNormal;
582  fExiting = fSaveState.sExiting;
583  fEntering = fSaveState.sEntering;
584
585  fBlockedPhysicalVolume = fSaveState.spBlockedPhysicalVolume;
586  fBlockedReplicaNo = fSaveState.sBlockedReplicaNo, 
587
588  fLastStepWasZero = fSaveState.sLastStepWasZero; 
589}
590
591// ********************************************************************
592// ComputeStep
593//
594// Computes the next geometric Step: intersections with current
595// mother and `daughter' volumes.
596//
597// NOTE:
598//
599// Flags on entry:
600// --------------
601// fValidExitNormal  - Normal of exited volume is valid (convex, not a
602//                     coincident boundary)
603// fExitNormal       - Surface normal of exited volume
604// fExiting          - True if have exited solid
605//
606// fBlockedPhysicalVolume - Ptr to exited volume (or 0)
607// fBlockedReplicaNo - Replication no of exited volume
608// fLastStepWasZero  - True if last Step size was zero.
609//
610// Flags on exit:
611// -------------
612// fValidExitNormal  - True if surface normal of exited volume is valid
613// fExitNormal       - Surface normal of exited volume rotated to mothers
614//                    reference system
615// fExiting          - True if exiting mother
616// fEntering         - True if entering `daughter' volume (or replica)
617// fBlockedPhysicalVolume - Ptr to candidate (entered) volume
618// fBlockedReplicaNo - Replication no of candidate (entered) volume
619// fLastStepWasZero  - True if this Step size was zero.
620// ********************************************************************
621//
622G4double G4Navigator::ComputeStep( const G4ThreeVector &pGlobalpoint,
623                                   const G4ThreeVector &pDirection,
624                                   const G4double pCurrentProposedStepLength,
625                                         G4double &pNewSafety)
626{
627  G4ThreeVector localDirection = ComputeLocalAxis(pDirection);
628  G4double Step = kInfinity;
629  G4VPhysicalVolume  *motherPhysical = fHistory.GetTopVolume();
630  G4LogicalVolume *motherLogical = motherPhysical->GetLogicalVolume();
631
632  static G4int sNavCScalls=0;
633  sNavCScalls++;
634
635#ifdef G4VERBOSE
636  if( fVerbose > 0 )
637  {
638    G4cout << "*** G4Navigator::ComputeStep: ***" << G4endl; 
639    G4cout << "    Volume = " << motherPhysical->GetName() 
640           << " - Proposed step length = " << pCurrentProposedStepLength
641           << G4endl; 
642#ifdef G4DEBUG_NAVIGATION
643    if( fVerbose >= 4 ) 
644    {
645      G4cout << "  Called with the arguments: " << G4endl
646             << "  Globalpoint = " << std::setw(25) << pGlobalpoint << G4endl
647             << "  Direction   = " << std::setw(25) << pDirection << G4endl;
648      G4cout << "  ---- Upon entering :" << G4endl;
649      PrintState();
650    }
651#endif
652  }
653#endif
654
655  G4ThreeVector newLocalPoint = ComputeLocalPoint(pGlobalpoint);
656  if( newLocalPoint != fLastLocatedPointLocal )
657  {
658    // Check whether the relocation is within safety
659    //
660    G4ThreeVector oldLocalPoint = fLastLocatedPointLocal;
661    G4double moveLenSq = (newLocalPoint-oldLocalPoint).mag2();
662
663    if ( moveLenSq >= kCarTolerance*kCarTolerance )
664    {
665#ifdef G4VERBOSE
666      ComputeStepLog(pGlobalpoint, moveLenSq);
667#endif
668      // Relocate the point within the same volume
669      //
670      LocateGlobalPointWithinVolume( pGlobalpoint );
671    }
672  }
673  if ( fHistory.GetTopVolumeType()!=kReplica )
674  {
675    switch( CharacteriseDaughters(motherLogical) )
676    {
677      case kNormal:
678        if ( motherLogical->GetVoxelHeader() )
679        {
680          Step = fvoxelNav.ComputeStep(fLastLocatedPointLocal,
681                                       localDirection,
682                                       pCurrentProposedStepLength,
683                                       pNewSafety,
684                                       fHistory,
685                                       fValidExitNormal,
686                                       fExitNormal,
687                                       fExiting,
688                                       fEntering,
689                                       &fBlockedPhysicalVolume,
690                                       fBlockedReplicaNo);
691     
692        }
693        else
694        {
695          if( motherPhysical->GetRegularStructureId() == 0 )
696          {
697            Step = fnormalNav.ComputeStep(fLastLocatedPointLocal,
698                                          localDirection,
699                                          pCurrentProposedStepLength,
700                                          pNewSafety,
701                                          fHistory,
702                                          fValidExitNormal,
703                                          fExitNormal,
704                                          fExiting,
705                                          fEntering,
706                                          &fBlockedPhysicalVolume,
707                                          fBlockedReplicaNo);
708          }
709          else  // Regular (non-voxelised) structure
710          {
711            LocateGlobalPointAndSetup( pGlobalpoint, &pDirection, true, true );
712            //
713            // if physical process limits the step, the voxel will not be the
714            // one given by ComputeStepSkippingEqualMaterials() and the local
715            // point will be wrongly calculated.
716
717            // There is a problem: when msc limits the step and the point is
718            // assigned wrongly to phantom in previous step (while it is out
719            // of the container volume). Then LocateGlobalPointAndSetup() has
720            // reset the history topvolume to world.
721            //
722            if(fHistory.GetTopVolume()->GetRegularStructureId() == 0 )
723            { 
724              G4Exception("G4Navigator::ComputeStep()",
725                          "Bad-location-of-point", JustWarning,
726                "Point is relocated in voxels, while it should be outside!");
727              Step = fnormalNav.ComputeStep(fLastLocatedPointLocal,
728                                            localDirection,
729                                            pCurrentProposedStepLength,
730                                            pNewSafety,
731                                            fHistory,
732                                            fValidExitNormal,
733                                            fExitNormal,
734                                            fExiting,
735                                            fEntering,
736                                            &fBlockedPhysicalVolume,
737                                            fBlockedReplicaNo);
738            }
739            else
740            {
741              Step = fregularNav.
742                   ComputeStepSkippingEqualMaterials(fLastLocatedPointLocal,
743                                                     localDirection,
744                                                     pCurrentProposedStepLength,
745                                                     pNewSafety,
746                                                     fHistory,
747                                                     fValidExitNormal,
748                                                     fExitNormal,
749                                                     fExiting,
750                                                     fEntering,
751                                                     &fBlockedPhysicalVolume,
752                                                     fBlockedReplicaNo,
753                                                     motherPhysical);
754            }
755          }
756        }
757        break;
758      case kParameterised:
759        if( GetDaughtersRegularStructureId(motherLogical) != 1 )
760        {
761          Step = fparamNav.ComputeStep(fLastLocatedPointLocal,
762                                       localDirection,
763                                       pCurrentProposedStepLength,
764                                       pNewSafety,
765                                       fHistory,
766                                       fValidExitNormal,
767                                       fExitNormal,
768                                       fExiting,
769                                       fEntering,
770                                       &fBlockedPhysicalVolume,
771                                       fBlockedReplicaNo);
772        }
773        else  // Regular structure
774        {
775          Step = fregularNav.ComputeStep(fLastLocatedPointLocal,
776                                         localDirection,
777                                         pCurrentProposedStepLength,
778                                         pNewSafety,
779                                         fHistory,
780                                         fValidExitNormal,
781                                         fExitNormal,
782                                         fExiting,
783                                         fEntering,
784                                         &fBlockedPhysicalVolume,
785                                         fBlockedReplicaNo);
786        }
787        break;
788      case kReplica:
789        G4Exception("G4Navigator::ComputeStep()", "NotApplicable",
790                    FatalException, "Not applicable for replicated volumes.");
791        break;
792    }
793  }
794  else
795  {
796    // In the case of a replica, it must handle the exiting
797    // edge/corner problem by itself
798    //
799    G4bool exitingReplica = fExitedMother;
800    Step = freplicaNav.ComputeStep(pGlobalpoint,
801                                   pDirection,
802                                   fLastLocatedPointLocal,
803                                   localDirection,
804                                   pCurrentProposedStepLength,
805                                   pNewSafety,
806                                   fHistory,
807                                   fValidExitNormal,
808                                   fExitNormal,
809                                   exitingReplica,
810                                   fEntering,
811                                   &fBlockedPhysicalVolume,
812                                   fBlockedReplicaNo);
813    fExiting= exitingReplica;                          // still ok to set it ??
814  }
815
816  // Remember last safety origin & value.
817  //
818  fPreviousSftOrigin = pGlobalpoint;
819  fPreviousSafety = pNewSafety; 
820
821  // Count zero steps - one can occur due to changing momentum at a boundary
822  //                  - one, two (or a few) can occur at common edges between
823  //                    volumes
824  //                  - more than two is likely a problem in the geometry
825  //                    description or the Navigation
826
827  // Rule of thumb: likely at an Edge if two consecutive steps are zero,
828  //                because at least two candidate volumes must have been
829  //                checked
830  //
831  fLocatedOnEdge   = fLastStepWasZero && (Step==0.0);
832  fLastStepWasZero = (Step==0.0);
833  if (fPushed)  fPushed = fLastStepWasZero;
834
835  // Handle large number of consecutive zero steps
836  //
837  if ( fLastStepWasZero )
838  {
839    fNumberZeroSteps++;
840#ifdef G4DEBUG_NAVIGATION
841    if( fNumberZeroSteps > 1 )
842    {
843       G4cout << "G4Navigator::ComputeStep(): another zero step, # "
844              << fNumberZeroSteps
845              << " at " << pGlobalpoint
846              << " in volume " << motherPhysical->GetName()
847              << " nav-comp-step calls # " << sNavCScalls
848              << G4endl;
849    }
850#endif
851    if( fNumberZeroSteps > fActionThreshold_NoZeroSteps-1 )
852    {
853       // Act to recover this stuck track. Pushing it along direction
854       //
855       Step += 100*kCarTolerance;
856#ifdef G4VERBOSE
857       if ((!fPushed) && (fWarnPush))
858       {
859         G4cerr << "WARNING - G4Navigator::ComputeStep()" << G4endl
860                << "          Track stuck, not moving for " 
861                << fNumberZeroSteps << " steps" << G4endl
862                << "          in volume -" << motherPhysical->GetName()
863                << "- at point " << pGlobalpoint << G4endl
864                << "          direction: " << pDirection << "." << G4endl
865                << "          Potential geometry or navigation problem !"
866                << G4endl
867                << "          Trying pushing it of " << Step << " mm ..."
868                << G4endl;
869       }
870#endif
871       fPushed = true;
872    }
873    if( fNumberZeroSteps > fAbandonThreshold_NoZeroSteps-1 )
874    {
875      // Must kill this stuck track
876      //
877      G4cerr << "ERROR - G4Navigator::ComputeStep()" << G4endl
878             << "        Track stuck, not moving for " 
879             << fNumberZeroSteps << " steps" << G4endl
880             << "        in volume -" << motherPhysical->GetName()
881             << "- at point " << pGlobalpoint << G4endl
882             << "        direction: " << pDirection << "." << G4endl;
883      motherPhysical->CheckOverlaps(5000, false);
884      G4Exception("G4Navigator::ComputeStep()",
885                  "StuckTrack", EventMustBeAborted, 
886                  "Stuck Track: potential geometry or navigation problem.");
887    }
888  }
889  else
890  {
891    if (!fPushed)  fNumberZeroSteps = 0;
892  }
893
894  fEnteredDaughter = fEntering;   // I expect to enter a volume in this Step
895  fExitedMother = fExiting;
896
897  if( fExiting )
898  {
899#ifdef G4DEBUG_NAVIGATION
900    if( fVerbose > 2 )
901    { 
902      G4cout << " At G4Nav CompStep End - if(exiting) - fExiting= " << fExiting
903             << " fValidExitNormal = " << fValidExitNormal  << G4endl;
904      G4cout << " fExitNormal= " << fExitNormal << G4endl;
905    }
906#endif
907
908    if(fValidExitNormal)
909    {
910      // Convention: fExitNormal is in the 'grand-mother' coordinate system
911      //
912      fGrandMotherExitNormal= fExitNormal;
913    }
914    else
915    { 
916      // We must calculate the normal anyway (in order to have it if requested)
917      //
918      G4ThreeVector finalLocalPoint =
919        fLastLocatedPointLocal + localDirection*Step;
920
921      // Now fGrandMotherExitNormal is in the 'grand-mother' coordinate system
922      //
923      fGrandMotherExitNormal =
924        motherLogical->GetSolid()->SurfaceNormal(finalLocalPoint);
925
926      const G4RotationMatrix* mRot = motherPhysical->GetRotation();
927      if( mRot )
928      { 
929        fGrandMotherExitNormal *= (*mRot).inverse();
930      }
931    }
932  }
933  fStepEndPoint= pGlobalpoint+Step*pDirection; 
934
935  if( (Step == pCurrentProposedStepLength) && (!fExiting) && (!fEntering) )
936  {
937    // This if Step is not really limited by the geometry.
938    // The Navigator is obliged to return "infinity"
939    //
940    Step = kInfinity;
941  }
942
943#ifdef G4VERBOSE
944  if( fVerbose > 1 )
945  {
946    if( fVerbose >= 4 )
947    {
948      G4cout << "    ----- Upon exiting :" << G4endl;
949      PrintState();
950    }
951    G4cout <<"    Returned step = " << Step << G4endl;
952    if( Step == kInfinity )
953    {
954      G4cout << "    Original proposed step = "
955             << pCurrentProposedStepLength << G4endl;
956    }
957    G4cout << "    Safety = " << pNewSafety << G4endl;
958  }
959#endif
960
961  return Step;
962}
963
964// ********************************************************************
965// CheckNextStep
966//
967// Compute the step without altering the navigator state
968// ********************************************************************
969//
970G4double G4Navigator::CheckNextStep( const G4ThreeVector& pGlobalpoint,
971                                     const G4ThreeVector& pDirection,
972                                     const G4double pCurrentProposedStepLength,
973                                           G4double& pNewSafety)
974{
975  G4double step;
976
977  // Save the state, for this parasitic call
978  //
979  SetSavedState();
980
981  step = ComputeStep ( pGlobalpoint, 
982                       pDirection,
983                       pCurrentProposedStepLength, 
984                       pNewSafety ); 
985
986  // If a parasitic call, then attempt to restore the key parts of the state
987  //
988  RestoreSavedState(); 
989
990  return step; 
991}
992
993// ********************************************************************
994// ResetState
995//
996// Resets stack and minimum of navigator state `machine'
997// ********************************************************************
998//
999void G4Navigator::ResetState()
1000{
1001  fWasLimitedByGeometry  = false;
1002  fEntering              = false;
1003  fExiting               = false;
1004  fLocatedOnEdge         = false;
1005  fLastStepWasZero       = false;
1006  fEnteredDaughter       = false;
1007  fExitedMother          = false;
1008  fPushed                = false;
1009
1010  fValidExitNormal       = false;
1011  fExitNormal            = G4ThreeVector(0,0,0);
1012
1013  fPreviousSftOrigin     = G4ThreeVector(0,0,0);
1014  fPreviousSafety        = 0.0; 
1015
1016  fNumberZeroSteps       = 0;
1017   
1018  fBlockedPhysicalVolume = 0;
1019  fBlockedReplicaNo      = -1;
1020
1021  fLastLocatedPointLocal = G4ThreeVector( kInfinity, -kInfinity, 0.0 ); 
1022  fLocatedOutsideWorld   = false;
1023}
1024
1025// ********************************************************************
1026// SetupHierarchy
1027//
1028// Renavigates & resets hierarchy described by current history
1029// o Reset volumes
1030// o Recompute transforms and/or solids of replicated/parameterised volumes
1031// ********************************************************************
1032//
1033void G4Navigator::SetupHierarchy()
1034{
1035  G4int i;
1036  const G4int cdepth = fHistory.GetDepth();
1037  G4VPhysicalVolume *mother, *current;
1038  G4VSolid *pSolid;
1039  G4VPVParameterisation *pParam;
1040
1041  for ( i=1; i<=cdepth; i++ )
1042  {
1043    current = fHistory.GetVolume(i);
1044    switch ( fHistory.GetVolumeType(i) )
1045    {
1046      case kNormal:
1047        break;
1048      case kReplica:
1049        freplicaNav.ComputeTransformation(fHistory.GetReplicaNo(i), current);
1050        break;
1051      case kParameterised:
1052        G4int replicaNo;
1053        pParam = current->GetParameterisation();
1054        replicaNo = fHistory.GetReplicaNo(i);
1055        pSolid = pParam->ComputeSolid(replicaNo, current);
1056
1057        // Set up dimensions & transform in solid/physical volume
1058        //
1059        pSolid->ComputeDimensions(pParam, replicaNo, current);
1060        pParam->ComputeTransformation(replicaNo, current);
1061
1062        G4TouchableHistory touchable( fHistory );
1063        touchable.MoveUpHistory();  // move up to the parent level
1064     
1065        // Set up the correct solid and material in Logical Volume
1066        //
1067        G4LogicalVolume *pLogical = current->GetLogicalVolume();
1068        pLogical->SetSolid( pSolid );
1069        pLogical->UpdateMaterial( pParam ->
1070          ComputeMaterial(replicaNo, current, &touchable) );
1071        break;
1072    }
1073    mother = current;
1074  }
1075}
1076
1077// ********************************************************************
1078// GetLocalExitNormal
1079//
1080// Obtains the Normal vector to a surface (in local coordinates)
1081// pointing out of previous volume and into current volume
1082// ********************************************************************
1083//
1084G4ThreeVector G4Navigator::GetLocalExitNormal( G4bool* valid )
1085{
1086  G4ThreeVector ExitNormal(0.,0.,0.);
1087
1088  if ( EnteredDaughterVolume() )
1089  {
1090    ExitNormal= -(fHistory.GetTopVolume()->GetLogicalVolume()->
1091                  GetSolid()->SurfaceNormal(fLastLocatedPointLocal));
1092    *valid = true;
1093  }
1094  else
1095  {
1096    if( fExitedMother )
1097    {
1098      ExitNormal = fGrandMotherExitNormal;
1099      *valid = true;
1100    }
1101    else
1102    {
1103      // We are not at a boundary.
1104      // ExitNormal remains (0,0,0)
1105      //
1106      *valid = false;
1107    }
1108  }
1109  return ExitNormal;
1110}
1111
1112// ********************************************************************
1113// ComputeSafety
1114//
1115// It assumes that it will be
1116//  i) called at the Point in the same volume as the EndPoint of the
1117//     ComputeStep.
1118// ii) after (or at the end of) ComputeStep OR after the relocation.
1119// ********************************************************************
1120//
1121G4double G4Navigator::ComputeSafety( const G4ThreeVector &pGlobalpoint,
1122                                     const G4double pMaxLength,
1123                                     const G4bool keepState)
1124{
1125  G4double newSafety = 0.0;
1126
1127#ifdef G4DEBUG_NAVIGATION
1128  G4int oldcoutPrec = G4cout.precision(8);
1129  if( fVerbose > 0 )
1130  {
1131    G4cout << "*** G4Navigator::ComputeSafety: ***" << G4endl
1132           << "    Called at point: " << pGlobalpoint << G4endl;
1133
1134    G4VPhysicalVolume  *motherPhysical = fHistory.GetTopVolume();
1135    G4cout << "    Volume = " << motherPhysical->GetName() 
1136           << " - Maximum length = " << pMaxLength << G4endl; 
1137    if( fVerbose >= 4 )
1138    {
1139       G4cout << "    ----- Upon entering Compute Safety:" << G4endl;
1140       PrintState();
1141    }
1142  }
1143#endif
1144
1145  if (keepState)  { SetSavedState(); }
1146
1147  G4double distEndpointSq = (pGlobalpoint-fStepEndPoint).mag2(); 
1148  G4bool   stayedOnEndpoint  = distEndpointSq < kCarTolerance*kCarTolerance; 
1149  G4bool   endpointOnSurface = fEnteredDaughter || fExitedMother;
1150
1151  if( !(endpointOnSurface && stayedOnEndpoint) )
1152  {
1153    // Pseudo-relocate to this point (updates voxel information only)
1154    //
1155    LocateGlobalPointWithinVolume( pGlobalpoint );
1156      // --->> Danger: Side effects on sub-navigator voxel information <<---
1157      //       Could be replaced again by 'granular' calls to sub-navigator
1158      //       locates (similar side-effects, but faster. 
1159      //       Solutions:
1160      //        1) Re-locate (to where?)
1161      //        2) Insure that the methods using (G4ComputeStep?)
1162      //           does a relocation (if information is disturbed only ?)
1163
1164#ifdef G4DEBUG_NAVIGATION
1165    if( fVerbose >= 2 )
1166    {
1167      G4cout << "  G4Navigator::ComputeSafety() relocates-in-volume to point: "
1168             << pGlobalpoint << G4endl;
1169    }
1170#endif
1171    G4VPhysicalVolume *motherPhysical = fHistory.GetTopVolume();
1172    G4LogicalVolume *motherLogical = motherPhysical->GetLogicalVolume();
1173    G4SmartVoxelHeader* pVoxelHeader = motherLogical->GetVoxelHeader();
1174    G4ThreeVector localPoint = ComputeLocalPoint(pGlobalpoint);
1175
1176    if ( fHistory.GetTopVolumeType()!=kReplica )
1177    {
1178      switch(CharacteriseDaughters(motherLogical))
1179      {
1180        case kNormal:
1181          if ( pVoxelHeader )
1182          {
1183            newSafety=fvoxelNav.ComputeSafety(localPoint,fHistory,pMaxLength);
1184          }
1185          else
1186          {
1187            newSafety=fnormalNav.ComputeSafety(localPoint,fHistory,pMaxLength);
1188          }
1189          break;
1190        case kParameterised:
1191          if( GetDaughtersRegularStructureId(motherLogical) != 1 )
1192          {
1193            newSafety = fparamNav.ComputeSafety(localPoint,fHistory,pMaxLength);
1194          }
1195          else  // Regular structure
1196          {
1197            newSafety = fregularNav.ComputeSafety(localPoint,fHistory,pMaxLength);
1198          }
1199          break;
1200        case kReplica:
1201          G4Exception("G4Navigator::ComputeSafety()", "NotApplicable",
1202                      FatalException, "Not applicable for replicated volumes.");
1203          break;
1204      }
1205    }
1206    else
1207    {
1208      newSafety = freplicaNav.ComputeSafety(pGlobalpoint, localPoint,
1209                                            fHistory, pMaxLength);
1210    }
1211  }
1212  else // if( endpointOnSurface && stayedOnEndpoint )
1213  {
1214#ifdef G4DEBUG_NAVIGATION
1215    if( fVerbose >= 2 )
1216    {
1217      G4cout << "    G4Navigator::ComputeSafety() finds that point - " 
1218             << pGlobalpoint << " - is on surface " << G4endl; 
1219      if( fEnteredDaughter ) { G4cout << "   entered new daughter volume"; }
1220      if( fExitedMother )    { G4cout << "   and exited previous volume."; }
1221      G4cout << G4endl;
1222      G4cout << " EndPoint was = " << fStepEndPoint << G4endl;
1223    } 
1224#endif
1225    newSafety = 0.0; 
1226  }
1227
1228  // Remember last safety origin & value
1229  //
1230  fPreviousSftOrigin = pGlobalpoint;
1231  fPreviousSafety = newSafety; 
1232
1233  if (keepState)  { RestoreSavedState(); }
1234
1235#ifdef G4DEBUG_NAVIGATION
1236  if( fVerbose > 1 )
1237  {
1238    G4cout << "   ---- Exiting ComputeSafety  " << G4endl;
1239    if( fVerbose > 2 )  { PrintState(); }
1240    G4cout << "    Returned value of Safety = " << newSafety << G4endl;
1241  }
1242  G4cout.precision(oldcoutPrec);
1243#endif
1244
1245  return newSafety;
1246}
1247
1248// ********************************************************************
1249// CreateTouchableHistoryHandle
1250// ********************************************************************
1251//
1252G4TouchableHistoryHandle G4Navigator::CreateTouchableHistoryHandle() const
1253{
1254  return G4TouchableHistoryHandle( CreateTouchableHistory() );
1255}
1256
1257// ********************************************************************
1258// PrintState
1259// ********************************************************************
1260//
1261void  G4Navigator::PrintState() const
1262{
1263  G4int oldcoutPrec = G4cout.precision(4);
1264  if( fVerbose == 4 )
1265  {
1266    G4cout << "The current state of G4Navigator is: " << G4endl;
1267    G4cout << "  ValidExitNormal= " << fValidExitNormal << G4endl
1268           << "  ExitNormal     = " << fExitNormal      << G4endl
1269           << "  Exiting        = " << fExiting         << G4endl
1270           << "  Entering       = " << fEntering        << G4endl
1271           << "  BlockedPhysicalVolume= " ;
1272    if (fBlockedPhysicalVolume==0)
1273      G4cout << "None";
1274    else
1275      G4cout << fBlockedPhysicalVolume->GetName();
1276    G4cout << G4endl
1277           << "  BlockedReplicaNo     = " <<  fBlockedReplicaNo       << G4endl
1278           << "  LastStepWasZero      = " <<   fLastStepWasZero       << G4endl
1279           << G4endl;   
1280  }
1281  if( ( 1 < fVerbose) && (fVerbose < 4) )
1282  {
1283    G4cout << std::setw(30) << " ExitNormal "  << " "     
1284           << std::setw( 5) << " Valid "       << " "     
1285           << std::setw( 9) << " Exiting "     << " "     
1286           << std::setw( 9) << " Entering"     << " " 
1287           << std::setw(15) << " Blocked:Volume "  << " "   
1288           << std::setw( 9) << " ReplicaNo"        << " " 
1289           << std::setw( 8) << " LastStepZero  "   << " "   
1290           << G4endl;   
1291    G4cout << "( " << std::setw(7) << fExitNormal.x() 
1292           << ", " << std::setw(7) << fExitNormal.y()
1293           << ", " << std::setw(7) << fExitNormal.z() << " ) "
1294           << std::setw( 5)  << fValidExitNormal  << " "   
1295           << std::setw( 9)  << fExiting          << " "
1296           << std::setw( 9)  << fEntering         << " ";
1297    if ( fBlockedPhysicalVolume==0 )
1298      G4cout << std::setw(15) << "None";
1299    else
1300      G4cout << std::setw(15)<< fBlockedPhysicalVolume->GetName();
1301      G4cout << std::setw( 9)  << fBlockedReplicaNo  << " "
1302             << std::setw( 8)  << fLastStepWasZero   << " "
1303             << G4endl;   
1304  }
1305  if( fVerbose > 2 ) 
1306  {
1307    G4cout.precision(8);
1308    G4cout << " Current Localpoint = " << fLastLocatedPointLocal << G4endl;
1309    G4cout << " PreviousSftOrigin  = " << fPreviousSftOrigin << G4endl;
1310    G4cout << " PreviousSafety     = " << fPreviousSafety << G4endl; 
1311  }
1312  G4cout.precision(oldcoutPrec);
1313}
1314
1315// ********************************************************************
1316// ComputeStepLog
1317// ********************************************************************
1318//
1319void G4Navigator::ComputeStepLog(const G4ThreeVector& pGlobalpoint,
1320                                       G4double moveLenSq) const
1321{
1322  //  The following checks only make sense if the move is larger
1323  //  than the tolerance.
1324
1325  static const G4double fAccuracyForWarning   = kCarTolerance,
1326                        fAccuracyForException = 1000*kCarTolerance;
1327
1328  G4ThreeVector OriginalGlobalpoint = fHistory.GetTopTransform().Inverse().
1329                                      TransformPoint(fLastLocatedPointLocal); 
1330
1331  G4double shiftOriginSafSq = (fPreviousSftOrigin-pGlobalpoint).mag2();
1332
1333  // Check that the starting point of this step is
1334  // within the isotropic safety sphere of the last point
1335  // to a accuracy/precision  given by fAccuracyForWarning.
1336  //   If so give warning.
1337  //   If it fails by more than fAccuracyForException exit with error.
1338  //
1339  if( shiftOriginSafSq >= sqr(fPreviousSafety) )
1340  {
1341    G4double shiftOrigin = std::sqrt(shiftOriginSafSq);
1342    G4double diffShiftSaf = shiftOrigin - fPreviousSafety;
1343
1344    if( diffShiftSaf > fAccuracyForWarning )
1345    {
1346      G4int oldcoutPrec= G4cout.precision(8);
1347      G4int oldcerrPrec= G4cerr.precision(10);
1348      G4Exception("G4Navigator::ComputeStep()",
1349                  "UnexpectedPositionShift", JustWarning,
1350                  "Accuracy error or slightly inaccurate position shift.");
1351      G4cerr << "     The Step's starting point has moved " 
1352             << std::sqrt(moveLenSq)/mm << " mm " << G4endl
1353             << "     since the last call to a Locate method." << G4endl;
1354      G4cerr << "     This has resulted in moving " 
1355             << shiftOrigin/mm << " mm " 
1356             << " from the last point at which the safety " 
1357             << "     was calculated " << G4endl;
1358      G4cerr << "     which is more than the computed safety= " 
1359             << fPreviousSafety/mm << " mm  at that point." << G4endl;
1360      G4cerr << "     This difference is " 
1361             << diffShiftSaf/mm << " mm." << G4endl
1362             << "     The tolerated accuracy is "
1363             << fAccuracyForException/mm << " mm." << G4endl;
1364
1365      static G4int warnNow = 0;
1366      if( ((++warnNow % 100) == 1) )
1367      {
1368        G4cerr << "  This problem can be due to either " << G4endl;
1369        G4cerr << "    - a process that has proposed a displacement"
1370               << " larger than the current safety , or" << G4endl;
1371        G4cerr << "    - inaccuracy in the computation of the safety"
1372                   << G4endl;
1373        G4cerr << "  We suggest that you " << G4endl
1374               << "   - find i) what particle is being tracked, and "
1375               << " ii) through what part of your geometry " << G4endl
1376               << "      for example by re-running this event with "
1377               << G4endl
1378               << "         /tracking/verbose 1 "  << G4endl
1379               << "    - check which processes you declare for"
1380               << " this particle (and look at non-standard ones)"
1381               << G4endl
1382               << "   - in case, create a detailed logfile"
1383               << " of this event using:" << G4endl
1384               << "         /tracking/verbose 6 "
1385               << G4endl;
1386      }
1387      G4cout.precision(oldcoutPrec);
1388      G4cerr.precision(oldcerrPrec);
1389    }
1390#ifdef G4DEBUG_NAVIGATION
1391    else
1392    {
1393      G4cerr << "WARNING - G4Navigator::ComputeStep()" << G4endl
1394             << "          The Step's starting point has moved "
1395             << std::sqrt(moveLenSq) << "," << G4endl
1396             << "          which has taken it to the limit of"
1397             << " the current safety. " << G4endl;
1398    }
1399#endif
1400  }
1401  G4double safetyPlus = fPreviousSafety + fAccuracyForException;
1402  if ( shiftOriginSafSq > sqr(safetyPlus) )
1403  {
1404    G4cerr << "ERROR - G4Navigator::ComputeStep()" << G4endl
1405           << "        Position has shifted considerably without"
1406           << " notifying the navigator !" << G4endl
1407           << "        Tolerated safety: " << safetyPlus << G4endl
1408           << "        Computed shift  : " << shiftOriginSafSq << G4endl;
1409    G4Exception("G4Navigator::ComputeStep()",
1410                "SignificantPositionShift", JustWarning,
1411                "May lead to a crash or unreliable results.");
1412  }
1413}
1414
1415// ********************************************************************
1416// Operator <<
1417// ********************************************************************
1418//
1419std::ostream& operator << (std::ostream &os,const G4Navigator &n)
1420{
1421  os << "Current History: " << G4endl << n.fHistory;
1422  return os;
1423}
Note: See TracBrowser for help on using the repository browser.