G4PathFinder.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 //
00027 // $Id$
00028 // GEANT4 tag $ Name:  $
00029 // 
00030 // class G4PathFinder Implementation
00031 //
00032 // Original author:  John Apostolakis,  April 2006
00033 // 
00034 // --------------------------------------------------------------------
00035 
00036 #include <iomanip>
00037 
00038 #include "G4PathFinder.hh"
00039 
00040 #include "G4SystemOfUnits.hh"
00041 #include "G4GeometryTolerance.hh"
00042 #include "G4Navigator.hh"
00043 #include "G4PropagatorInField.hh"
00044 #include "G4TransportationManager.hh"
00045 #include "G4MultiNavigator.hh"
00046 #include "G4SafetyHelper.hh"
00047 
00048 // Initialise the static instance of the singleton
00049 //
00050 G4PathFinder* G4PathFinder::fpPathFinder=0;
00051 
00052 // ----------------------------------------------------------------------------
00053 // GetInstance()
00054 //
00055 // Retrieve the static instance of the singleton
00056 //
00057 G4PathFinder* G4PathFinder::GetInstance()
00058 {
00059    static G4PathFinder theInstance; 
00060    if( ! fpPathFinder )
00061    {
00062      fpPathFinder = &theInstance; 
00063    }
00064    return fpPathFinder;
00065 }
00066 
00067 // ----------------------------------------------------------------------------
00068 // Constructor
00069 //
00070 G4PathFinder::G4PathFinder() 
00071   : fEndState( G4ThreeVector(), G4ThreeVector(), 0., 0., 0., 0., 0.),
00072     fFieldExertedForce(false),
00073     fRelocatedPoint(true),
00074     fLastStepNo(-1), fCurrentStepNo(-1),
00075     fVerboseLevel(0)
00076 {
00077    fpMultiNavigator= new G4MultiNavigator(); 
00078 
00079    fpTransportManager= G4TransportationManager::GetTransportationManager();
00080    fpFieldPropagator = fpTransportManager->GetPropagatorInField();
00081 
00082    kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
00083 
00084    fNoActiveNavigators= 0; 
00085    G4ThreeVector  Big3Vector( kInfinity, kInfinity, kInfinity );
00086    fLastLocatedPosition= Big3Vector;
00087    fSafetyLocation= Big3Vector; 
00088    fPreSafetyLocation= Big3Vector;
00089    fPreStepLocation= Big3Vector;
00090 
00091    fPreSafetyMinValue=  -1.0; 
00092    fMinSafety_PreStepPt= -1.0; 
00093    fMinSafety_atSafLocation= -1.0; 
00094    fMinStep=   -1.0;
00095    fTrueMinStep= -1.0;
00096    fPreStepCenterRenewed= false;
00097    fNewTrack= false; 
00098    fNoGeometriesLimiting= 0; 
00099 
00100    for( register int num=0; num< fMaxNav; ++num )
00101    {
00102       fpNavigator[num] =  0;   
00103       fLimitTruth[num] = false;
00104       fLimitedStep[num] = kUndefLimited;
00105       fCurrentStepSize[num] = -1.0; 
00106       fLocatedVolume[num] = 0;
00107       fPreSafetyValues[num]= -1.0; 
00108       fCurrentPreStepSafety[num] = -1.0;
00109       fNewSafetyComputed[num]= -1.0; 
00110    }
00111 }
00112 
00113 // ----------------------------------------------------------------------------
00114 // Destructor
00115 //
00116 G4PathFinder::~G4PathFinder() 
00117 {
00118    delete fpMultiNavigator; 
00119 }
00120 
00121 // ----------------------------------------------------------------------------
00122 //
00123 void
00124 G4PathFinder::EnableParallelNavigation(G4bool enableChoice)
00125 {
00126    G4Navigator *navigatorForPropagation=0, *massNavigator=0;
00127 
00128    massNavigator= fpTransportManager->GetNavigatorForTracking(); 
00129    if( enableChoice )
00130    {
00131       navigatorForPropagation= fpMultiNavigator;
00132 
00133       // Enable SafetyHelper to use PF
00134       //
00135       fpTransportManager->GetSafetyHelper()->EnableParallelNavigation(true);
00136    }
00137    else
00138    {
00139       navigatorForPropagation= massNavigator;
00140        
00141       // Disable SafetyHelper to use PF
00142       //
00143       fpTransportManager->GetSafetyHelper()->EnableParallelNavigation(false);
00144    }
00145    fpFieldPropagator->SetNavigatorForPropagating(navigatorForPropagation);
00146 }
00147 
00148 // ----------------------------------------------------------------------------
00149 //
00150 G4double 
00151 G4PathFinder::ComputeStep( const G4FieldTrack &InitialFieldTrack, 
00152                                  G4double     proposedStepLength,
00153                                  G4int        navigatorNo, 
00154                                  G4int        stepNo,       // find next step 
00155                                  G4double     &pNewSafety,  // for this geom 
00156                                  ELimited     &limitedStep, 
00157                                  G4FieldTrack &EndState,
00158                                  G4VPhysicalVolume* currentVolume)
00159 {
00160   G4double  possibleStep= -1.0; 
00161 
00162 #ifdef G4DEBUG_PATHFINDER
00163   if( fVerboseLevel > 2 )
00164   { 
00165     G4cout << " -------------------------" <<  G4endl;
00166     G4cout << " G4PathFinder::ComputeStep - entered " << G4endl;
00167     G4cout << "   - stepNo = "  << std::setw(4) << stepNo  << " "
00168            << " navigatorId = " << std::setw(2) << navigatorNo  << " "
00169            << " proposed step len = " << proposedStepLength << " " << G4endl;
00170     G4cout << " PF::ComputeStep requested step " 
00171            << " from " << InitialFieldTrack.GetPosition()
00172            << " dir  " << InitialFieldTrack.GetMomentumDirection() << G4endl;
00173   }
00174 #endif
00175 #ifdef G4VERBOSE
00176   if( navigatorNo >= fNoActiveNavigators )
00177   {
00178     std::ostringstream message;
00179     message << "Bad Navigator ID !" << G4endl
00180             << "        Requested Navigator ID = " << navigatorNo << G4endl
00181             << "        Number of active navigators = " << fNoActiveNavigators;
00182     G4Exception("G4PathFinder::ComputeStep()", "GeomNav0002",
00183                 FatalException, message); 
00184   }
00185 #endif
00186 
00187   if( fNewTrack || (stepNo != fLastStepNo)  )
00188   {
00189     // This is a new track or a new step, so we must make the step
00190     // ( else we can simply retrieve its results for this Navigator Id )    
00191 
00192     G4FieldTrack currentState= InitialFieldTrack;
00193 
00194     fCurrentStepNo = stepNo; 
00195 
00196     // Check whether a process shifted the position 
00197     // since the last step -- by physics processes
00198     //
00199     G4ThreeVector newPosition = InitialFieldTrack.GetPosition();   
00200     G4ThreeVector moveVector= newPosition - fLastLocatedPosition; 
00201     G4double moveLenSq= moveVector.mag2(); 
00202     if( moveLenSq > kCarTolerance * kCarTolerance )
00203     { 
00204        G4ThreeVector newDirection = InitialFieldTrack.GetMomentumDirection();   
00205 #ifdef G4DEBUG_PATHFINDER
00206        if( fVerboseLevel > 2 )
00207        { 
00208           G4double moveLen= std::sqrt( moveLenSq ); 
00209           G4cout << " G4PathFinder::ComputeStep : Point moved since last step " 
00210                  << " -- at step # = " << stepNo << G4endl
00211                  << " by " << moveLen  << " to " << newPosition << G4endl;      
00212        } 
00213 #endif
00214        MovePoint();  // Unintentional changed -- ????
00215 
00216        // Relocate to cope with this move -- else could abort !?
00217        //
00218        Locate( newPosition, newDirection ); 
00219     }
00220 
00221     // Check whether the particle have an (EM) field force exerting upon it
00222     //
00223     G4double particleCharge=  currentState.GetCharge(); 
00224 
00225     G4FieldManager* fieldMgr=0;
00226     G4bool          fieldExertsForce = false ;
00227     if( (particleCharge != 0.0) )
00228     {
00229         fieldMgr= fpFieldPropagator->FindAndSetFieldManager( currentVolume );
00230 
00231         // Protect for case where field manager has no field (= field is zero)
00232         //
00233         fieldExertsForce = (fieldMgr != 0) 
00234                         && (fieldMgr->GetDetectorField() != 0);
00235     }
00236     fFieldExertedForce = fieldExertsForce;  // Store for use in later calls
00237                                             // referring to this 'step'.
00238 
00239     fNoGeometriesLimiting= -1;  // At start of track, no process limited step
00240     if( fieldExertsForce )
00241     {
00242        DoNextCurvedStep( currentState, proposedStepLength, currentVolume ); 
00243        //--------------
00244     }else{
00245        DoNextLinearStep( currentState, proposedStepLength ); 
00246        //--------------
00247     }
00248     fLastStepNo= stepNo; 
00249 
00250 #ifdef  G4DEBUG_PATHFINDER
00251     if ( (fNoGeometriesLimiting < 0)
00252       || (fNoGeometriesLimiting > fNoActiveNavigators) )
00253     {
00254       std::ostringstream message;
00255       message << "Number of geometries limiting the step not set." << G4endl
00256               << "        Number of geometries limiting step = "
00257               << fNoGeometriesLimiting;
00258       G4Exception("G4PathFinder::ComputeStep()", 
00259                   "GeomNav0002", FatalException, message); 
00260     }
00261 #endif
00262   }
00263 #ifdef G4DEBUG_PATHFINDER      
00264   else
00265   {
00266      if( proposedStepLength < fTrueMinStep )  // For 2nd+ geometry 
00267      { 
00268        std::ostringstream message;
00269        message << "Problem in step size request." << G4endl
00270                << "        Error can be caused by incorrect process ordering."
00271                << "        Being requested to make a step which is shorter"
00272                << " than the minimum Step " << G4endl
00273                << "        already computed for any Navigator/geometry during"
00274                << " this tracking-step: " << G4endl
00275                << "        This can happen due to an error in process ordering."
00276                << G4endl
00277                << "        Check that all physics processes are registered"
00278                << G4endl
00279                << "        before all processes with a navigator/geometry."
00280                << G4endl
00281                << "        If using pre-packaged physics list and/or"
00282                << G4endl
00283                << "        functionality, please report this error."
00284                << G4endl << G4endl
00285                << "        Additional information for problem: "  << G4endl
00286                << "        Steps request/proposed = " << proposedStepLength
00287                << G4endl
00288                << "        MinimumStep (true) = " << fTrueMinStep
00289                << G4endl
00290                << "        MinimumStep (navraw)  = " << fMinStep
00291                << G4endl
00292                << "        Navigator raw return value" << G4endl
00293                << "        Requested step now = " << proposedStepLength
00294                << G4endl
00295                << "        Difference min-req = "
00296                << fTrueMinStep-proposedStepLength << G4endl
00297                << "     -- Step info> stepNo= " << stepNo
00298                << " last= " << fLastStepNo 
00299                << " newTr= " << fNewTrack << G4endl;
00300         G4Exception("G4PathFinder::ComputeStep()", 
00301                     "GeomNav0003", FatalException, message);
00302      }
00303      else
00304      { 
00305         // This is neither a new track nor a new step -- just another 
00306         // client accessing information for the current track, step 
00307         // We will simply retrieve the results of the synchronous
00308         // stepping for this Navigator Id below.
00309         //
00310         if( fVerboseLevel > 1 )
00311         { 
00312            G4cout << " G4P::CS -> Not calling DoNextLinearStep: " 
00313                   << " stepNo= " << stepNo << " last= " << fLastStepNo 
00314                   << " new= " << fNewTrack << " Step already done" << G4endl; 
00315         }
00316      } 
00317   }
00318 #endif
00319 
00320   fNewTrack= false; 
00321 
00322   // Prepare the information to return
00323 
00324   pNewSafety  = fCurrentPreStepSafety[ navigatorNo ]; 
00325   limitedStep = fLimitedStep[ navigatorNo ];
00326   fRelocatedPoint= false;
00327 
00328   possibleStep= std::min(proposedStepLength, fCurrentStepSize[ navigatorNo ]);
00329   EndState = fEndState;  //  now corrected for smaller step, if needed
00330 
00331 #ifdef G4DEBUG_PATHFINDER
00332   if( fVerboseLevel > 0 )
00333   { 
00334     G4cout << " G4PathFinder::ComputeStep returns "
00335            << fCurrentStepSize[ navigatorNo ]
00336            << " for Navigator " << navigatorNo 
00337            << " Limited step = " << limitedStep 
00338            << " Safety(mm) = " << pNewSafety / mm 
00339            << G4endl; 
00340   }
00341 #endif
00342 
00343   return possibleStep;
00344 }
00345 
00346 // ----------------------------------------------------------------------
00347 
00348 void
00349 G4PathFinder::PrepareNewTrack( const G4ThreeVector& position, 
00350                                const G4ThreeVector& direction,
00351                                G4VPhysicalVolume*  massStartVol)
00352 {
00353   // Key purposes:
00354   //   - Check and cache set of active navigators
00355   //   - Reset state for new track
00356 
00357   G4int num=0; 
00358 
00359   EnableParallelNavigation(true); 
00360     // Switch PropagatorInField to use MultiNavigator
00361 
00362   fpTransportManager->GetSafetyHelper()->InitialiseHelper(); 
00363     // Reinitialise state of safety helper -- avoid problems with overlaps
00364 
00365   fNewTrack= true; 
00366   this->MovePoint();   // Signal further that the last status is wiped
00367 
00368   // Message the G4NavigatorPanel / Dispatcher to find active navigators
00369   //
00370   std::vector<G4Navigator*>::iterator pNavigatorIter; 
00371 
00372   fNoActiveNavigators=  fpTransportManager-> GetNoActiveNavigators();
00373   if( fNoActiveNavigators > fMaxNav )
00374   {
00375     std::ostringstream message;
00376     message << "Too many active Navigators / worlds." << G4endl
00377             << "        Transportation Manager has "
00378             << fNoActiveNavigators << " active navigators." << G4endl
00379             << "        This is more than the number allowed = "
00380             << fMaxNav << " !";
00381     G4Exception("G4PathFinder::PrepareNewTrack()", "GeomNav0002",  
00382                 FatalException, message); 
00383   }
00384 
00385   fpMultiNavigator->PrepareNavigators(); 
00386   //------------------------------------
00387 
00388   pNavigatorIter= fpTransportManager->GetActiveNavigatorsIterator();
00389   for( num=0; num< fNoActiveNavigators; ++pNavigatorIter,++num )
00390   {
00391      // Keep information in C-array ... for creating touchables - at least
00392 
00393      fpNavigator[num] =  *pNavigatorIter;   
00394      fLimitTruth[num] = false;
00395      fLimitedStep[num] = kDoNot;
00396      fCurrentStepSize[num] = 0.0; 
00397      fLocatedVolume[num] = 0; 
00398   }
00399   fNoGeometriesLimiting= 0;  // At start of track, no process limited step
00400 
00401   // In case of one geometry, the tracking will have done the locating!!
00402 
00403   if( fNoActiveNavigators > 1 )
00404   {
00405      Locate( position, direction, false );   
00406   }
00407   else
00408   {
00409      // Update state -- depending on the tracking's call to Mass Navigator
00410 
00411      fLastLocatedPosition= position; 
00412      fLocatedVolume[0]= massStartVol; // This information must be given
00413                                       // by transportation
00414      fLimitedStep[0]   = kDoNot; 
00415      fCurrentStepSize[0] = 0.0;
00416   }
00417 
00418   // Reset Safety Information -- as in case of overlaps this can cause
00419   // inconsistencies ...
00420   //
00421   fMinSafety_PreStepPt= fPreSafetyMinValue= fMinSafety_atSafLocation= 0.0; 
00422  
00423   for( num=0; num< fNoActiveNavigators; ++num )
00424   {
00425      fPreSafetyValues[num]= 0.0; 
00426      fNewSafetyComputed[num]= 0.0; 
00427      fCurrentPreStepSafety[num] = 0.0;
00428   }
00429 
00430   // The first location for each Navigator must be non-relative
00431   // or else call ResetStackAndState() for each Navigator
00432 
00433   fRelocatedPoint= false; 
00434 }
00435 
00436 void G4PathFinder::ReportMove( const G4ThreeVector& OldVector, 
00437                                const G4ThreeVector& NewVector, 
00438                                const G4String& Quantity ) const
00439 {
00440     G4ThreeVector moveVec = ( NewVector - OldVector );
00441 
00442     G4int prc= G4cerr.precision(12); 
00443     std::ostringstream message;
00444     message << "Endpoint moved between value returned by ComputeStep()"
00445             << " and call to Locate(). " << G4endl
00446             << "          Change of " << Quantity << " is "
00447             << moveVec.mag() / mm << " mm long" << G4endl
00448             << "          and its vector is "
00449             << (1.0/mm) * moveVec << " mm " << G4endl
00450             << "          Endpoint of ComputeStep() was " << OldVector << G4endl
00451             << "          and current position to locate is " << NewVector;
00452     G4Exception("G4PathFinder::ReportMove()", "GeomNav1002",  
00453                 JustWarning, message); 
00454     G4cerr.precision(prc); 
00455 }
00456 
00457 void
00458 G4PathFinder::Locate( const   G4ThreeVector& position, 
00459                       const   G4ThreeVector& direction,
00460                       G4bool  relative)
00461 {
00462   // Locate the point in each geometry
00463 
00464   std::vector<G4Navigator*>::iterator pNavIter=
00465      fpTransportManager->GetActiveNavigatorsIterator(); 
00466 
00467   G4ThreeVector lastEndPosition= fEndState.GetPosition(); 
00468   G4ThreeVector moveVec = (position - lastEndPosition );
00469   G4double      moveLenSq= moveVec.mag2();
00470   if( (!fNewTrack) && (!fRelocatedPoint)
00471    && ( moveLenSq> kCarTolerance*kCarTolerance ) )
00472   {
00473      ReportMove( position, lastEndPosition, "Position" ); 
00474   }
00475   fLastLocatedPosition= position; 
00476 
00477 #ifdef G4DEBUG_PATHFINDER
00478   if( fVerboseLevel > 2 )
00479   {
00480     G4cout << G4endl; 
00481     G4cout << " G4PathFinder::Locate : entered " << G4endl;
00482     G4cout << " --------------------   -------" <<  G4endl;
00483     G4cout << "   Locating at position " << position
00484            << "  with direction " << direction 
00485            << "  relative= " << relative << G4endl;
00486     if ( (fVerboseLevel > 1) || ( moveLenSq > 0.0) )
00487     { 
00488        G4cout << "  lastEndPosition = " << lastEndPosition
00489               << "  moveVec = " << moveVec
00490               << "  newTr = " << fNewTrack 
00491               << "  relocated = " << fRelocatedPoint << G4endl;
00492     }
00493 
00494     G4cout << " Located at " << position ; 
00495     if( fNoActiveNavigators > 1 )  { G4cout << G4endl; }
00496   }
00497 #endif
00498 
00499   for ( register G4int num=0; num< fNoActiveNavigators ; ++pNavIter,++num )
00500   {
00501      //  ... who limited the step ....
00502 
00503      if( fLimitTruth[num] ) { (*pNavIter)->SetGeometricallyLimitedStep(); }
00504 
00505      G4VPhysicalVolume *pLocated= 
00506      (*pNavIter)->LocateGlobalPointAndSetup( position, &direction,
00507                                              relative,  
00508                                              false);   
00509      // Set the state related to the location
00510      //
00511      fLocatedVolume[num] = pLocated; 
00512 
00513      // Clear state related to the step
00514      //
00515      fLimitedStep[num]   = kDoNot; 
00516      fCurrentStepSize[num] = 0.0;      
00517     
00518 #ifdef G4DEBUG_PATHFINDER
00519      if( fVerboseLevel > 2 )
00520      {
00521        G4cout << " - In world " << num << " geomLimStep= " << fLimitTruth[num]
00522               << "  gives volume= " << pLocated ; 
00523        if( pLocated )
00524        { 
00525          G4cout << "  name = '" << pLocated->GetName() << "'"; 
00526          G4cout << " - CopyNo= " << pLocated->GetCopyNo(); 
00527        } 
00528        G4cout  << G4endl; 
00529      }
00530 #endif
00531   }
00532 
00533   fRelocatedPoint= false;
00534 }
00535 
00536 void G4PathFinder::ReLocate( const G4ThreeVector& position )
00537 {
00538   // Locate the point in each geometry
00539 
00540   std::vector<G4Navigator*>::iterator pNavIter=
00541     fpTransportManager->GetActiveNavigatorsIterator(); 
00542 
00543 #ifdef G4DEBUG_PATHFINDER
00544 
00545   // Check that this relocation does not violate safety
00546   //   - at endpoint (computed from start point) AND
00547   //   - at last safety location  (likely just called)
00548 
00549   G4ThreeVector lastEndPosition= fEndState.GetPosition();
00550 
00551   // Calculate end-point safety ...
00552   //
00553   G4double      DistanceStartEnd= (lastEndPosition - fPreStepLocation).mag();
00554   G4double      endPointSafety_raw = fMinSafety_PreStepPt - DistanceStartEnd; 
00555   G4double      endPointSafety_Est1 = std::max( 0.0, endPointSafety_raw ); 
00556 
00557   // ... and check move from endpoint against this endpoint safety
00558   //
00559   G4ThreeVector moveVecEndPos  = position - lastEndPosition;
00560   G4double      moveLenEndPosSq = moveVecEndPos.mag2(); 
00561 
00562   // Check that move from endpoint of last step is within safety
00563   // -- or check against last location or relocation ?? 
00564   //
00565   G4ThreeVector moveVecSafety=  position - fSafetyLocation; 
00566   G4double      moveLenSafSq=   moveVecSafety.mag2();
00567 
00568   G4double distCheckEnd_sq= ( moveLenEndPosSq - endPointSafety_Est1 
00569                                                *endPointSafety_Est1 ); 
00570   G4double distCheckSaf_sq=   ( moveLenSafSq -  fMinSafety_atSafLocation
00571                                                *fMinSafety_atSafLocation ); 
00572 
00573   G4bool longMoveEnd = distCheckEnd_sq > 0.0; 
00574   G4bool longMoveSaf = distCheckSaf_sq > 0.0; 
00575 
00576   G4double revisedSafety= 0.0;
00577 
00578   if( (!fNewTrack) && ( longMoveEnd && longMoveSaf ) )
00579   {  
00580      // Recompute ComputeSafety for end position
00581      //
00582      revisedSafety= ComputeSafety(lastEndPosition); 
00583 
00584      const G4double kRadTolerance =
00585        G4GeometryTolerance::GetInstance()->GetRadialTolerance();
00586      const G4double cErrorTolerance=1e-12;   
00587        // Maximum relative error from roundoff of arithmetic 
00588 
00589      G4double  distCheckRevisedEnd= moveLenEndPosSq-revisedSafety*revisedSafety;
00590 
00591      G4bool  longMoveRevisedEnd=  ( distCheckRevisedEnd > 0. ) ; 
00592 
00593      G4double  moveMinusSafety= 0.0; 
00594      G4double  moveLenEndPosition= std::sqrt( moveLenEndPosSq );
00595      moveMinusSafety = moveLenEndPosition - revisedSafety; 
00596 
00597      if ( longMoveRevisedEnd && (moveMinusSafety > 0.0 )
00598        && ( revisedSafety > 0.0 ) )
00599      {
00600         // Take into account possibility of roundoff error causing
00601         // this apparent move further than safety
00602 
00603         if( fVerboseLevel > 0 )
00604         {
00605            G4cout << " G4PF:Relocate> Ratio to revised safety is " 
00606                   << std::fabs(moveMinusSafety)/revisedSafety << G4endl;
00607         }
00608 
00609         G4double  absMoveMinusSafety= std::fabs(moveMinusSafety);
00610         G4bool smallRatio= absMoveMinusSafety < kRadTolerance * revisedSafety ; 
00611         G4double maxCoordPos = std::max( 
00612                                       std::max( std::fabs(position.x()), 
00613                                                 std::fabs(position.y())), 
00614                                       std::fabs(position.z()) );
00615         G4bool smallValue= absMoveMinusSafety < cErrorTolerance * maxCoordPos;
00616         if( ! (smallRatio || smallValue) )
00617         {
00618            G4cout << " G4PF:Relocate> Ratio to revised safety is " 
00619                   << std::fabs(moveMinusSafety)/revisedSafety << G4endl;
00620            G4cout << " Difference of move and safety is not very small."
00621                   << G4endl;
00622         }
00623         else
00624         {
00625           moveMinusSafety = 0.0; 
00626           longMoveRevisedEnd = false;   // Numerical issue -- not too long!
00627 
00628           G4cout << " Difference of move & safety is very small in magnitude, "
00629                  << absMoveMinusSafety << G4endl;
00630           if( smallRatio )
00631           {
00632             G4cout << " ratio to safety " << revisedSafety 
00633                    << " is " <<  absMoveMinusSafety / revisedSafety
00634                    << "smaller than " << kRadTolerance << " of safety ";
00635           }
00636           else
00637           {
00638             G4cout << " as fraction " << absMoveMinusSafety / maxCoordPos 
00639                    << " of position vector max-coord " << maxCoordPos
00640                    << " smaller than " << cErrorTolerance ;
00641           }
00642           G4cout << " -- reset moveMinusSafety to "
00643                  << moveMinusSafety << G4endl;
00644         }
00645      }
00646 
00647      if ( longMoveEnd && longMoveSaf
00648        && longMoveRevisedEnd && (moveMinusSafety>0.0) )
00649      { 
00650         G4int oldPrec= G4cout.precision(9); 
00651         std::ostringstream message;
00652         message << "ReLocation is further than end-safety value." << G4endl
00653                 << " Moved from last endpoint by " << moveLenEndPosition 
00654                 << " compared to end safety (from preStep point) = " 
00655                 << endPointSafety_Est1 << G4endl
00656                 << "  --> last PreSafety Location was " << fPreSafetyLocation
00657                 << G4endl
00658                 << "       safety value =  " << fPreSafetyMinValue << G4endl
00659                 << "  --> last PreStep Location was " << fPreStepLocation
00660                 << G4endl
00661                 << "       safety value =  " << fMinSafety_PreStepPt << G4endl
00662                 << "  --> last EndStep Location was " << lastEndPosition
00663                 << G4endl
00664                 << "       safety value =  " << endPointSafety_Est1 
00665                 << " raw-value = " << endPointSafety_raw << G4endl
00666                 << "  --> Calling again at this endpoint, we get "
00667                 <<  revisedSafety << " as safety value."  << G4endl
00668                 << "  --> last position for safety " << fSafetyLocation
00669                 << G4endl
00670                 << "       its safety value =  " << fMinSafety_atSafLocation
00671                 << G4endl
00672                 << "       move from safety location = "
00673                 << std::sqrt(moveLenSafSq) << G4endl
00674                 << "         again= " << moveVecSafety.mag() << G4endl
00675                 << "       safety - Move-from-end= " 
00676                 << revisedSafety - moveLenEndPosition
00677                 << " (negative is Bad.)" << G4endl
00678                 << " Debug:  distCheckRevisedEnd = "
00679                 << distCheckRevisedEnd;
00680         ReportMove( lastEndPosition, position, "Position" ); 
00681         G4Exception("G4PathFinder::ReLocate", "GeomNav0003", 
00682                     FatalException, message); 
00683         G4cout.precision(oldPrec); 
00684     }
00685   }
00686 
00687   if( fVerboseLevel > 2 )
00688   {
00689     G4cout << G4endl; 
00690     G4cout << " G4PathFinder::ReLocate : entered " << G4endl;
00691     G4cout << " ----------------------   -------" <<  G4endl;
00692     G4cout << "  *Re*Locating at position " << position  << G4endl; 
00693       // << "  with direction " << direction 
00694       // << "  relative= " << relative << G4endl;
00695     if ( (fVerboseLevel > -1) || ( moveLenEndPosSq > 0.0) )
00696     {
00697        G4cout << "  lastEndPosition = " << lastEndPosition
00698               << "  moveVec from step-end = " << moveVecEndPos
00699               << "  is new Track = " << fNewTrack 
00700               << "  relocated = " << fRelocatedPoint << G4endl;
00701     }
00702   }
00703 #endif // G4DEBUG_PATHFINDER
00704 
00705   for ( register G4int num=0; num< fNoActiveNavigators ; ++pNavIter,++num )
00706   {
00707      //  ... none limited the step
00708 
00709      (*pNavIter)->LocateGlobalPointWithinVolume( position ); 
00710 
00711      // Clear state related to the step
00712      //
00713      fLimitedStep[num]   = kDoNot; 
00714      fCurrentStepSize[num] = 0.0;      
00715      fLimitTruth[num] = false;   
00716   }
00717 
00718   fLastLocatedPosition= position; 
00719   fRelocatedPoint= false;
00720 
00721 #ifdef G4DEBUG_PATHFINDER
00722   if( fVerboseLevel > 2 )
00723   {
00724     G4cout << " G4PathFinder::ReLocate : exiting " 
00725            << "  at position " << fLastLocatedPosition << G4endl << G4endl;
00726   }
00727 #endif
00728 }
00729 
00730 // -----------------------------------------------------------------------------
00731 
00732 G4double  G4PathFinder::ComputeSafety( const G4ThreeVector& position )
00733 {
00734     // Recompute safety for the relevant point
00735 
00736    G4double minSafety= kInfinity; 
00737   
00738    std::vector<G4Navigator*>::iterator pNavigatorIter;
00739    pNavigatorIter= fpTransportManager->GetActiveNavigatorsIterator();
00740 
00741    for( register G4int num=0; num<fNoActiveNavigators; ++pNavigatorIter,++num )
00742    {
00743       G4double safety = (*pNavigatorIter)->ComputeSafety( position,true );
00744       if( safety < minSafety ) { minSafety = safety; } 
00745       fNewSafetyComputed[num]= safety;
00746    } 
00747 
00748    fSafetyLocation= position;
00749    fMinSafety_atSafLocation = minSafety;
00750 
00751 #ifdef G4DEBUG_PATHFINDER
00752    if( fVerboseLevel > 1 )
00753    { 
00754      G4cout << " G4PathFinder::ComputeSafety - returns " 
00755             << minSafety << " at location " << position << G4endl;
00756    }
00757 #endif
00758    return minSafety; 
00759 }
00760 
00761 
00762 // -----------------------------------------------------------------------------
00763 
00764 G4TouchableHandle 
00765 G4PathFinder::CreateTouchableHandle( G4int navId ) const
00766 {
00767 #ifdef G4DEBUG_PATHFINDER
00768   if( fVerboseLevel > 2 )
00769   {
00770     G4cout << "G4PathFinder::CreateTouchableHandle : navId = "
00771            << navId << " -- " << GetNavigator(navId) << G4endl;
00772   }
00773 #endif
00774 
00775   G4TouchableHistory* touchHist;
00776   touchHist= GetNavigator(navId) -> CreateTouchableHistory(); 
00777 
00778   G4VPhysicalVolume* locatedVolume= fLocatedVolume[navId]; 
00779   if( locatedVolume == 0 )
00780   {
00781      // Workaround to ensure that the touchable is fixed !! // TODO: fix
00782 
00783      touchHist->UpdateYourself( locatedVolume, touchHist->GetHistory() );
00784   }
00785  
00786 #ifdef G4DEBUG_PATHFINDER
00787   if( fVerboseLevel > 2 )
00788   {   
00789     G4String VolumeName("None"); 
00790     if( locatedVolume ) { VolumeName= locatedVolume->GetName(); }
00791     G4cout << " Touchable History created at address " << touchHist
00792            << "  volume = " << locatedVolume << " name= " << VolumeName
00793            << G4endl;
00794   }
00795 #endif
00796 
00797   return G4TouchableHandle(touchHist); 
00798 }
00799 
00800 G4double
00801 G4PathFinder::DoNextLinearStep( const G4FieldTrack &initialState,
00802                                       G4double      proposedStepLength )
00803 {
00804   std::vector<G4Navigator*>::iterator pNavigatorIter;
00805   G4double safety= 0.0, step=0.0;
00806   G4double minSafety= kInfinity, minStep;
00807 
00808   const G4int IdTransport= 0;  // Id of Mass Navigator !!
00809   register G4int num=0; 
00810 
00811 #ifdef G4DEBUG_PATHFINDER
00812   if( fVerboseLevel > 2 )
00813   {
00814     G4cout << " G4PathFinder::DoNextLinearStep : entered " << G4endl;
00815     G4cout << "   Input field track= " << initialState << G4endl;
00816     G4cout << "   Requested step= " << proposedStepLength << G4endl;
00817   }
00818 #endif
00819 
00820   G4ThreeVector initialPosition= initialState.GetPosition(); 
00821   G4ThreeVector initialDirection= initialState.GetMomentumDirection();
00822   
00823   G4ThreeVector OriginShift = initialPosition - fPreSafetyLocation;
00824   G4double      MagSqShift  = OriginShift.mag2() ;
00825   G4double      MagShift;  // Only given value if it larger than minimum safety
00826 
00827   // Potential optimisation using Maximum Value of safety!
00828   // if( MagSqShift >= sqr(fPreSafetyMaxValue ) ){ 
00829   //   MagShift= kInfinity;   // Not a useful value -- all will not use/ignore
00830   // else
00831   //  MagShift= std::sqrt(MagSqShift) ;
00832 
00833   MagShift= std::sqrt(MagSqShift) ;
00834 
00835 #ifdef G4PATHFINDER_OPTIMISATION
00836 
00837   G4double fullSafety;  // For all geometries, for prestep point
00838 
00839   if( MagSqShift >= sqr(fPreSafetyMinValue ) )
00840   {
00841      fullSafety = 0.0 ;     
00842   }
00843   else
00844   {
00845      fullSafety = fPreSafetyMinValue - MagShift;
00846   }
00847   if( proposedStepLength < fullSafety ) 
00848   {
00849      // Move is smaller than all safeties
00850      //  -> so we do not have to move the safety center
00851 
00852      fPreStepCenterRenewed= false;
00853 
00854      for( num=0; num< fNoActiveNavigators; ++num )
00855      {
00856         fCurrentStepSize[num]= kInfinity; 
00857         safety = std::max( 0.0,  fPreSafetyValues[num] - MagShift); 
00858         minSafety= std::min( safety, minSafety ); 
00859         fCurrentPreStepSafety[num]= safety; 
00860      }
00861      minStep= kInfinity;
00862 
00863 #ifdef G4DEBUG_PATHFINDER
00864      if( fVerboseLevel > 2 )
00865      {
00866        G4cout << "G4PathFinder::DoNextLinearStep : Quick Stepping. " << G4endl
00867                << " proposedStepLength " <<  proposedStepLength
00868                << " < (full) safety = " << fullSafety 
00869                << " at " << initialPosition 
00870                << G4endl;
00871      }
00872 #endif
00873   }
00874   else
00875 #endif   // End of G4PATHFINDER_OPTIMISATION 1
00876   {
00877      // Move is larger than at least one of the safeties
00878      //  -> so we must move the safety center!
00879 
00880      fPreStepCenterRenewed= true;
00881      pNavigatorIter= fpTransportManager-> GetActiveNavigatorsIterator();
00882 
00883      minStep= kInfinity;  // Not proposedStepLength; 
00884 
00885      for( num=0; num< fNoActiveNavigators; ++pNavigatorIter,++num ) 
00886      {
00887         safety = std::max( 0.0,  fPreSafetyValues[num] - MagShift); 
00888 
00889 #ifdef G4PATHFINDER_OPTIMISATION
00890         if( proposedStepLength <= safety )  // Should be just < safety ?
00891         {
00892            // The Step is guaranteed to be taken
00893 
00894            step= kInfinity;    //  ComputeStep Would return this
00895 
00896 #ifdef G4DEBUG_PATHFINDER
00897            G4cout.precision(8); 
00898            G4cout << "PathFinder::ComputeStep> small proposed step = "
00899                   << proposedStepLength
00900                   << " <=  safety = " << safety << " for nav " << num 
00901                   << " Step fully taken. " << G4endl;
00902 #endif
00903         }
00904         else
00905 #endif   // End of G4PATHFINDER_OPTIMISATION 2
00906         {
00907 #ifdef G4DEBUG_PATHFINDER
00908            G4double previousSafety= safety; 
00909 #endif
00910            step= (*pNavigatorIter)->ComputeStep( initialPosition, 
00911                                                  initialDirection,
00912                                                  proposedStepLength,
00913                                                  safety ); 
00914            minStep  = std::min( step,  minStep);
00915 
00916            //  TODO: consider whether/how to reduce the proposed step 
00917            //        to the latest minStep value - to reduce calculations
00918 
00919 #ifdef G4DEBUG_PATHFINDER
00920            if( fVerboseLevel > 0)
00921            {
00922              G4cout.precision(8); 
00923              G4cout << "PathFinder::ComputeStep> long  proposed step = "
00924                     << proposedStepLength
00925                     << "  >  safety = " << previousSafety
00926                     << " for nav " << num 
00927                     << " .  New safety = " << safety << " step= " << step
00928                     << G4endl;      
00929            } 
00930 #endif
00931         }
00932         fCurrentStepSize[num] = step; 
00933 
00934         // Save safety value, must be done for all geometries "together"
00935         // (even if not recomputed using call to ComputeStep)
00936         // since they share the fPreSafetyLocation
00937 
00938         fPreSafetyValues[num]= safety; 
00939         fCurrentPreStepSafety[num]= safety; 
00940 
00941         minSafety= std::min( safety, minSafety ); 
00942            
00943 #ifdef G4DEBUG_PATHFINDER
00944         if( fVerboseLevel > 2 )
00945         {
00946           G4cout << "G4PathFinder::DoNextLinearStep : Navigator ["
00947                  << num << "] -- step size " << step << G4endl;
00948         }
00949 #endif
00950      }
00951 
00952      // Only change these when safety is recalculated
00953      // it is good/relevant only for safety calculations
00954 
00955      fPreSafetyLocation=  initialPosition; 
00956      fPreSafetyMinValue=  minSafety;
00957   } // end of else for  if( proposedStepLength <= fullSafety)
00958 
00959   // For use in Relocation, need PreStep point location, min-safety
00960   //
00961   fPreStepLocation= initialPosition; 
00962   fMinSafety_PreStepPt= minSafety; 
00963 
00964   fMinStep=   minStep; 
00965 
00966   if( fMinStep == kInfinity )
00967   {
00968      minStep = proposedStepLength;   //  Use this below for endpoint !!
00969   }
00970   fTrueMinStep = minStep;
00971 
00972   // Set the EndState
00973 
00974   G4ThreeVector endPosition;
00975 
00976   fEndState= initialState; 
00977   endPosition= initialPosition + minStep * initialDirection ; 
00978 
00979 #ifdef G4DEBUG_PATHFINDER
00980   if( fVerboseLevel > 1 )
00981   {
00982     G4cout << "G4PathFinder::DoNextLinearStep : "
00983            << " initialPosition = " << initialPosition 
00984            << " and endPosition = " << endPosition<< G4endl;
00985   }
00986 #endif
00987 
00988   fEndState.SetPosition( endPosition ); 
00989   fEndState.SetProperTimeOfFlight( -1.000 );   // Not defined YET
00990 
00991   if( fNoActiveNavigators == 1 )
00992   { 
00993      G4bool transportLimited = (fMinStep!= kInfinity); 
00994      fLimitTruth[IdTransport] = transportLimited; 
00995      fLimitedStep[IdTransport] = transportLimited ? kUnique : kDoNot;
00996 
00997      // Set fNoGeometriesLimiting - as WhichLimited does
00998      fNoGeometriesLimiting = transportLimited ? 1 : 0;  
00999   }
01000   else
01001   {
01002      WhichLimited(); 
01003   }
01004 
01005 #ifdef G4DEBUG_PATHFINDER
01006   if( fVerboseLevel > 2 )
01007   {
01008     G4cout << " G4PathFinder::DoNextLinearStep : exits returning "
01009            << minStep << G4endl;
01010     G4cout << "   Endpoint values = " << fEndState << G4endl;
01011     G4cout << G4endl;
01012   }
01013 #endif
01014 
01015   return minStep;
01016 }
01017 
01018 void G4PathFinder::WhichLimited()
01019 {
01020   // Flag which processes limited the step
01021 
01022   G4int num=-1, last=-1; 
01023   G4int noLimited=0; 
01024   ELimited shared= kSharedOther; 
01025 
01026   const G4int IdTransport= 0;  // Id of Mass Navigator !!
01027 
01028   // Assume that [IdTransport] is Mass / Transport
01029   //
01030   G4bool transportLimited = (fCurrentStepSize[IdTransport] == fMinStep)
01031                            && ( fMinStep!= kInfinity) ; 
01032 
01033   if( transportLimited )  { 
01034      shared= kSharedTransport;
01035   }
01036 
01037   for ( num= 0; num < fNoActiveNavigators; num++ )
01038   { 
01039     G4bool limitedStep;
01040 
01041     G4double step= fCurrentStepSize[num]; 
01042 
01043     limitedStep = ( std::fabs(step - fMinStep) < kCarTolerance ) 
01044                  && ( step != kInfinity); 
01045 
01046     fLimitTruth[ num ] = limitedStep; 
01047     if( limitedStep )
01048     {
01049       noLimited++;  
01050       fLimitedStep[num] = shared;
01051       last= num; 
01052     }
01053     else
01054     {
01055       fLimitedStep[num] = kDoNot;
01056     }
01057   }
01058   fNoGeometriesLimiting= noLimited;  // Save # processes limiting step
01059 
01060   if( (last > -1) && (noLimited == 1 ) )
01061   {
01062     fLimitedStep[ last ] = kUnique; 
01063   }
01064 
01065 #ifdef G4DEBUG_PATHFINDER
01066   if( fVerboseLevel > 1 )
01067   {
01068     PrintLimited();   // --> for tracing 
01069     if( fVerboseLevel > 4 ) {
01070       G4cout << " G4PathFinder::WhichLimited - exiting. " << G4endl;
01071     }
01072   }
01073 #endif
01074 }
01075 
01076 void G4PathFinder::PrintLimited()
01077 {
01078   // Report results -- for checking   
01079 
01080   G4cout << "G4PathFinder::PrintLimited reports: " ; 
01081   G4cout << "  Minimum step (true)= " << fTrueMinStep 
01082          << "  reported min = " << fMinStep 
01083          << G4endl; 
01084   if(  (fCurrentStepNo <= 2) || (fVerboseLevel>=2) )
01085   {
01086     G4cout << std::setw(5) << " Step#"  << " "
01087            << std::setw(5) << " NavId"  << " "
01088            << std::setw(12) << " step-size " << " "
01089            << std::setw(12) << " raw-size "  << " "
01090            << std::setw(12) << " pre-safety " << " " 
01091            << std::setw(15) << " Limited / flag"  << " "
01092            << std::setw(15) << "  World "  << " "
01093            << G4endl;  
01094   }
01095   G4int num;
01096   for ( num= 0; num < fNoActiveNavigators; num++ )
01097   { 
01098     G4double rawStep = fCurrentStepSize[num]; 
01099     G4double stepLen = fCurrentStepSize[num]; 
01100     if( stepLen > fTrueMinStep )
01101     { 
01102       stepLen = fTrueMinStep;     // did not limit (went as far as asked)
01103     }
01104     G4int oldPrec= G4cout.precision(9); 
01105 
01106     G4cout << std::setw(5) << fCurrentStepNo  << " " 
01107            << std::setw(5) << num  << " "
01108            << std::setw(12) << stepLen << " "
01109            << std::setw(12) << rawStep << " "
01110            << std::setw(12) << fCurrentPreStepSafety[num] << " "
01111            << std::setw(5) << (fLimitTruth[num] ? "YES" : " NO") << " ";
01112     G4String limitedStr= LimitedString(fLimitedStep[num]); 
01113     G4cout << " " << std::setw(15) << limitedStr << " ";  
01114     G4cout.precision(oldPrec); 
01115 
01116     G4Navigator *pNav= GetNavigator( num ); 
01117     G4String  WorldName( "Not-Set" ); 
01118     if (pNav)
01119     {
01120        G4VPhysicalVolume *pWorld= pNav->GetWorldVolume(); 
01121        if( pWorld )
01122        {
01123            WorldName = pWorld->GetName(); 
01124        }
01125     }
01126     G4cout << " " << WorldName ; 
01127     G4cout << G4endl;
01128   }
01129 
01130   if( fVerboseLevel > 4 )
01131   {
01132     G4cout << " G4PathFinder::PrintLimited - exiting. " << G4endl;
01133   }
01134 }
01135 
01136 G4double
01137 G4PathFinder::DoNextCurvedStep( const G4FieldTrack &initialState,
01138                                 G4double      proposedStepLength,
01139                                 G4VPhysicalVolume* pCurrentPhysicalVolume )
01140 {
01141   const G4double toleratedRelativeError= 1.0e-10; 
01142   G4double minStep= kInfinity, newSafety=0.0;
01143   G4int numNav; 
01144   G4FieldTrack  fieldTrack= initialState;
01145   G4ThreeVector startPoint= initialState.GetPosition(); 
01146 
01147 #ifdef G4DEBUG_PATHFINDER
01148   G4int prc= G4cout.precision(9);
01149   if( fVerboseLevel > 2 )
01150   {
01151     G4cout << " G4PathFinder::DoNextCurvedStep ****** " << G4endl;
01152     G4cout << " Initial value of field track is " << fieldTrack 
01153            << " and proposed step= " << proposedStepLength  << G4endl;
01154   }
01155 #endif
01156 
01157   fPreStepCenterRenewed= true; // Always update PreSafety with PreStep point
01158 
01159   if( fNoActiveNavigators > 1 )
01160   { 
01161      // Calculate the safety values before making the step
01162 
01163      G4double minSafety= kInfinity, safety; 
01164      for( numNav=0; numNav < fNoActiveNavigators; ++numNav )
01165      {
01166         safety= fpNavigator[numNav]->ComputeSafety( startPoint, false );
01167         fPreSafetyValues[numNav]= safety; 
01168         fCurrentPreStepSafety[numNav]= safety; 
01169         minSafety = std::min( safety, minSafety ); 
01170      }
01171 
01172      // Save safety value, related position
01173 
01174      fPreSafetyLocation=  startPoint;   
01175      fPreSafetyMinValue=  minSafety;
01176      fPreStepLocation=    startPoint;
01177      fMinSafety_PreStepPt= minSafety;
01178   }
01179 
01180   // Allow Propagator In Field to do the hard work, calling G4MultiNavigator
01181   //
01182   minStep=  fpFieldPropagator->ComputeStep( fieldTrack,
01183                                             proposedStepLength,
01184                                             newSafety, 
01185                                             pCurrentPhysicalVolume );
01186 
01187   // fieldTrack now contains the endpoint information
01188   //
01189   fEndState= fieldTrack; 
01190   fMinStep=   minStep; 
01191   fTrueMinStep = std::min( minStep, proposedStepLength );
01192 
01193   if( fNoActiveNavigators== 1 )
01194   { 
01195      // Update the 'PreSafety' sphere - as any ComputeStep was called 
01196      // (must be done anyway in field)
01197 
01198      fPreSafetyValues[0]=   newSafety;
01199      fPreSafetyLocation= startPoint;   
01200      fPreSafetyMinValue= newSafety;
01201 
01202      // Update the current 'PreStep' point's values - mandatory
01203      //
01204      fCurrentPreStepSafety[0]= newSafety; 
01205      fPreStepLocation=  startPoint;
01206      fMinSafety_PreStepPt= newSafety;
01207   }
01208 
01209 #ifdef G4DEBUG_PATHFINDER
01210   if( fVerboseLevel > 2 )
01211   {
01212     G4cout << "G4PathFinder::DoNextCurvedStep : " << G4endl
01213            << " initialState = " << initialState << G4endl
01214            << " and endState = " << fEndState << G4endl;
01215     G4cout << "G4PathFinder::DoNextCurvedStep : " 
01216            << " minStep = " << minStep 
01217            << " proposedStepLength " << proposedStepLength 
01218            << " safety = " << newSafety << G4endl;
01219   }
01220 #endif
01221   G4double currentStepSize;   // = 0.0; 
01222   if( minStep < proposedStepLength ) // if == , then a boundary found at end ??
01223   {   
01224     // Recover the remaining information from MultiNavigator
01225     // especially regarding which Navigator limited the step
01226 
01227     G4int noLimited= 0;  //   No geometries limiting step
01228     for( numNav=0; numNav < fNoActiveNavigators; ++numNav )
01229     {
01230       G4double finalStep, lastPreSafety=0.0, minStepLast;
01231       ELimited didLimit; 
01232       G4bool limited; 
01233 
01234       finalStep=  fpMultiNavigator->ObtainFinalStep( numNav, lastPreSafety, 
01235                                                      minStepLast, didLimit );
01236 
01237       // Calculate the step for this geometry, using the 
01238       // final step (the only one which can differ.)
01239 
01240       currentStepSize = fTrueMinStep;  
01241       G4double diffStep= 0.0; 
01242       if( (minStepLast != kInfinity) )
01243       { 
01244         diffStep = (finalStep-minStepLast);
01245         if ( std::abs(diffStep) <= toleratedRelativeError * finalStep ) 
01246         {
01247           diffStep = 0.0;
01248         }
01249         currentStepSize += diffStep; 
01250       }
01251       fCurrentStepSize[numNav] = currentStepSize;  
01252       
01253       // TODO: could refine the way to obtain safeties for > 1 geometries
01254       //     - for pre step safety
01255       //        notify MultiNavigator about new set of sub-steps
01256       //        allow it to return this value in ObtainFinalStep 
01257       //        instead of lastPreSafety (or as well?)
01258       //     - for final step start (available)
01259       //        get final Step start from MultiNavigator
01260       //        and corresponding safety values
01261       // and/or ALSO calculate ComputeSafety at endpoint
01262       //     endSafety= fpNavigator[numNav]->ComputeSafety( endPoint ); 
01263 
01264       fLimitedStep[numNav] = didLimit; 
01265       fLimitTruth[numNav] = limited = (didLimit != kDoNot ); 
01266       if( limited ) { noLimited++; }
01267 
01268 #ifdef G4DEBUG_PATHFINDER
01269       G4bool StepError= (currentStepSize < 0) 
01270                    || ( (minStepLast != kInfinity) && (diffStep < 0) ) ; 
01271       if( StepError || (fVerboseLevel > 2) )
01272       {
01273         G4String  limitedString=  LimitedString( fLimitedStep[numNav] ); 
01274         
01275         G4cout << " G4PathFinder::ComputeStep. Geometry " << numNav
01276                << "  step= " << fCurrentStepSize[numNav] 
01277                << " from final-step= " << finalStep 
01278                << " fTrueMinStep= " << fTrueMinStep 
01279                << " minStepLast= "  << minStepLast 
01280                << "  limited = " << (fLimitTruth[numNav] ? "YES" : " NO")
01281                << " ";
01282         G4cout << "  status = " << limitedString << " #= " << didLimit
01283                << G4endl;
01284         
01285         if( StepError )
01286         { 
01287           std::ostringstream message;
01288           message << "Incorrect calculation of step size for one navigator"
01289                   << G4endl
01290                   << "        currentStepSize = " << currentStepSize 
01291                   << ", diffStep= " << diffStep << G4endl
01292                   << "ERROR in computing step size for this navigator.";
01293           G4Exception("G4PathFinder::DoNextCurvedStep",
01294                       "GeomNav0003", FatalException, message);
01295         }
01296       }
01297 #endif
01298     } // for num Navigators
01299 
01300     fNoGeometriesLimiting= noLimited;  // Save # processes limiting step
01301   } 
01302   else if ( (minStep == proposedStepLength)  
01303             || (minStep == kInfinity)  
01304             || ( std::abs(minStep-proposedStepLength)
01305                < toleratedRelativeError * proposedStepLength ) )
01306   { 
01307     // In case the step was not limited, use default responses
01308     //  --> all Navigators 
01309     // Also avoid problems in case of PathFinder using safety to optimise
01310     //  - it is possible that the Navigators were not called
01311     //    if the safety was already satisfactory.
01312     //    (In that case calling ObtainFinalStep gives invalid results.)
01313 
01314     currentStepSize= minStep;
01315     for( numNav=0; numNav < fNoActiveNavigators; ++numNav )
01316     {
01317       fCurrentStepSize[numNav] = minStep; 
01318       // Safety for endpoint ??  // Can eventuall improve it -- see TODO above
01319       fLimitedStep[numNav] = kDoNot; 
01320       fLimitTruth[numNav] = false; 
01321     }
01322     fNoGeometriesLimiting= 0;  // Save # processes limiting step
01323   } 
01324   else    //  (minStep > proposedStepLength) and not (minStep == kInfinity)
01325   {
01326     std::ostringstream message;
01327     message << "Incorrect calculation of step size for one navigator." << G4endl
01328             << "        currentStepSize = " << minStep << " is larger than "
01329             << " proposed StepSize = " << proposedStepLength << ".";
01330     G4Exception("G4PathFinder::DoNextCurvedStep()",
01331                 "GeomNav0003", FatalException, message); 
01332   }
01333 
01334 #ifdef G4DEBUG_PATHFINDER
01335   if( fVerboseLevel > 2 )
01336   {
01337     G4cout << " Exiting G4PathFinder::DoNextCurvedStep " << G4endl;
01338     PrintLimited(); 
01339   }
01340   G4cout.precision(prc); 
01341 #endif
01342 
01343   return minStep; 
01344 }
01345 
01346 G4String& G4PathFinder::LimitedString( ELimited lim )
01347 {
01348   static G4String StrDoNot("DoNot"),
01349                   StrUnique("Unique"),
01350                   StrUndefined("Undefined"),
01351                   StrSharedTransport("SharedTransport"),  
01352                   StrSharedOther("SharedOther");
01353 
01354   G4String* limitedStr;
01355   switch ( lim )
01356   {
01357      case kDoNot:  limitedStr= &StrDoNot; break;
01358      case kUnique: limitedStr = &StrUnique; break; 
01359      case kSharedTransport:  limitedStr= &StrSharedTransport; break; 
01360      case kSharedOther: limitedStr = &StrSharedOther; break;
01361      default: limitedStr = &StrUndefined; break;
01362   }
01363   return *limitedStr;
01364 }
01365 
01366 void G4PathFinder::PushPostSafetyToPreSafety()
01367 {
01368   fPreSafetyLocation= fSafetyLocation;
01369   fPreSafetyMinValue= fMinSafety_atSafLocation;
01370   for( register G4int nav=0; nav < fNoActiveNavigators; ++nav )
01371   {
01372      fPreSafetyValues[nav]= fNewSafetyComputed[nav];
01373   }
01374 }

Generated on Mon May 27 17:49:17 2013 for Geant4 by  doxygen 1.4.7