Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************* 24 // ******************************************************************** >> 25 // >> 26 // >> 27 // $Id: G4PathFinder.cc 90009 2015-05-08 07:42:39Z gcosmo $ >> 28 // GEANT4 tag $ Name: $ 25 // 29 // 26 // G4PathFinder Implementation << 30 // class G4PathFinder Implementation 27 // 31 // 28 // Original author: John Apostolakis, April 20 << 32 // Original author: John Apostolakis, April 2006 >> 33 // 29 // ------------------------------------------- 34 // -------------------------------------------------------------------- 30 35 31 #include <iomanip> 36 #include <iomanip> 32 37 33 #include "G4PathFinder.hh" 38 #include "G4PathFinder.hh" 34 39 35 #include "G4SystemOfUnits.hh" 40 #include "G4SystemOfUnits.hh" 36 #include "G4GeometryTolerance.hh" 41 #include "G4GeometryTolerance.hh" 37 #include "G4Navigator.hh" 42 #include "G4Navigator.hh" 38 #include "G4PropagatorInField.hh" 43 #include "G4PropagatorInField.hh" 39 #include "G4TransportationManager.hh" 44 #include "G4TransportationManager.hh" 40 #include "G4MultiNavigator.hh" 45 #include "G4MultiNavigator.hh" 41 #include "G4SafetyHelper.hh" 46 #include "G4SafetyHelper.hh" 42 47 43 // Initialise the static instance of the singl 48 // Initialise the static instance of the singleton 44 // 49 // 45 G4ThreadLocal G4PathFinder* G4PathFinder::fpPa << 50 G4ThreadLocal G4PathFinder* G4PathFinder::fpPathFinder=0; 46 51 47 // ------------------------------------------- 52 // ---------------------------------------------------------------------------- 48 // GetInstance() 53 // GetInstance() 49 // 54 // 50 // Retrieve the static instance of the singlet << 55 // Retrieve the static instance of the singleton 51 // 56 // 52 G4PathFinder* G4PathFinder::GetInstance() 57 G4PathFinder* G4PathFinder::GetInstance() 53 { 58 { 54 if( fpPathFinder == nullptr ) << 59 if( ! fpPathFinder ) 55 { 60 { 56 fpPathFinder = new G4PathFinder; 61 fpPathFinder = new G4PathFinder; 57 } 62 } 58 return fpPathFinder; 63 return fpPathFinder; 59 } 64 } 60 65 61 // ------------------------------------------- 66 // ---------------------------------------------------------------------------- 62 // GetInstanceIfExist() << 63 // << 64 // Retrieve the static instance pointer of the << 65 // << 66 G4PathFinder* G4PathFinder::GetInstanceIfExist << 67 { << 68 return fpPathFinder; << 69 } << 70 << 71 // ------------------------------------------- << 72 // Constructor 67 // Constructor 73 // 68 // 74 G4PathFinder::G4PathFinder() 69 G4PathFinder::G4PathFinder() 75 : fEndState( G4ThreeVector(), G4ThreeVector( << 70 : fEndState( G4ThreeVector(), G4ThreeVector(), 0., 0., 0., 0., 0.), >> 71 fFieldExertedForce(false), >> 72 fRelocatedPoint(true), >> 73 fLastStepNo(-1), fCurrentStepNo(-1), >> 74 fVerboseLevel(0) 76 { 75 { 77 fpMultiNavigator = new G4MultiNavigator(); << 76 fpMultiNavigator= new G4MultiNavigator(); 78 77 79 fpTransportManager= G4TransportationManager 78 fpTransportManager= G4TransportationManager::GetTransportationManager(); 80 fpFieldPropagator = fpTransportManager->Get 79 fpFieldPropagator = fpTransportManager->GetPropagatorInField(); 81 80 82 kCarTolerance = G4GeometryTolerance::GetIns 81 kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance(); 83 82 >> 83 fNoActiveNavigators= 0; 84 G4ThreeVector Big3Vector( kInfinity, kInfi 84 G4ThreeVector Big3Vector( kInfinity, kInfinity, kInfinity ); 85 fLastLocatedPosition= Big3Vector; 85 fLastLocatedPosition= Big3Vector; 86 fSafetyLocation= Big3Vector; 86 fSafetyLocation= Big3Vector; 87 fPreSafetyLocation= Big3Vector; 87 fPreSafetyLocation= Big3Vector; 88 fPreStepLocation= Big3Vector; 88 fPreStepLocation= Big3Vector; 89 89 90 for( auto num=0; num<fMaxNav; ++num ) << 90 fPreSafetyMinValue= -1.0; >> 91 fMinSafety_PreStepPt= -1.0; >> 92 fMinSafety_atSafLocation= -1.0; >> 93 fMinStep= -1.0; >> 94 fTrueMinStep= -1.0; >> 95 fPreStepCenterRenewed= false; >> 96 fNewTrack= false; >> 97 fNoGeometriesLimiting= 0; >> 98 >> 99 for( G4int num=0; num< fMaxNav; ++num ) 91 { 100 { 92 fpNavigator[num] = nullptr; << 101 fpNavigator[num] = 0; 93 fLimitTruth[num] = false; 102 fLimitTruth[num] = false; 94 fLimitedStep[num] = kUndefLimited; 103 fLimitedStep[num] = kUndefLimited; 95 fCurrentStepSize[num] = -1.0; 104 fCurrentStepSize[num] = -1.0; 96 fLocatedVolume[num] = nullptr; << 105 fLocatedVolume[num] = 0; 97 fPreSafetyValues[num]= -1.0; 106 fPreSafetyValues[num]= -1.0; 98 fCurrentPreStepSafety[num] = -1.0; 107 fCurrentPreStepSafety[num] = -1.0; 99 fNewSafetyComputed[num]= -1.0; 108 fNewSafetyComputed[num]= -1.0; 100 } 109 } 101 } 110 } 102 111 103 // ------------------------------------------- 112 // ---------------------------------------------------------------------------- 104 // Destructor 113 // Destructor 105 // 114 // 106 G4PathFinder::~G4PathFinder() 115 G4PathFinder::~G4PathFinder() 107 { 116 { 108 delete fpMultiNavigator; 117 delete fpMultiNavigator; 109 fpPathFinder = nullptr; << 118 if (fpPathFinder) { delete fpPathFinder; fpPathFinder=0; } 110 } 119 } 111 120 112 // ------------------------------------------- 121 // ---------------------------------------------------------------------------- 113 // 122 // 114 void 123 void 115 G4PathFinder::EnableParallelNavigation(G4bool 124 G4PathFinder::EnableParallelNavigation(G4bool enableChoice) 116 { 125 { 117 G4Navigator *navigatorForPropagation = null << 126 G4Navigator *navigatorForPropagation=0, *massNavigator=0; 118 127 119 massNavigator = fpTransportManager->GetNavi << 128 massNavigator= fpTransportManager->GetNavigatorForTracking(); 120 if( enableChoice ) 129 if( enableChoice ) 121 { 130 { 122 navigatorForPropagation = fpMultiNavigat << 131 navigatorForPropagation= fpMultiNavigator; 123 132 124 // Enable SafetyHelper to use PF 133 // Enable SafetyHelper to use PF 125 // 134 // 126 fpTransportManager->GetSafetyHelper()->E 135 fpTransportManager->GetSafetyHelper()->EnableParallelNavigation(true); 127 } 136 } 128 else 137 else 129 { 138 { 130 navigatorForPropagation = massNavigator; << 139 navigatorForPropagation= massNavigator; 131 140 132 // Disable SafetyHelper to use PF 141 // Disable SafetyHelper to use PF 133 // 142 // 134 fpTransportManager->GetSafetyHelper()->E 143 fpTransportManager->GetSafetyHelper()->EnableParallelNavigation(false); 135 } 144 } 136 fpFieldPropagator->SetNavigatorForPropagati 145 fpFieldPropagator->SetNavigatorForPropagating(navigatorForPropagation); 137 } 146 } 138 147 139 // ------------------------------------------- 148 // ---------------------------------------------------------------------------- 140 // 149 // 141 G4double 150 G4double 142 G4PathFinder::ComputeStep( const G4FieldTrack& << 151 G4PathFinder::ComputeStep( const G4FieldTrack &InitialFieldTrack, 143 G4double << 152 G4double proposedStepLength, 144 G4int << 153 G4int navigatorNo, 145 G4int << 154 G4int stepNo, // find next step 146 G4double& << 155 G4double &pNewSafety, // for this geom 147 ELimited& << 156 ELimited &limitedStep, 148 G4FieldTrack& << 157 G4FieldTrack &EndState, 149 G4VPhysicalVo 158 G4VPhysicalVolume* currentVolume) 150 { 159 { 151 G4double possibleStep = -1.0; << 160 G4double possibleStep= -1.0; 152 161 153 #ifdef G4DEBUG_PATHFINDER 162 #ifdef G4DEBUG_PATHFINDER 154 if( fVerboseLevel > 2 ) 163 if( fVerboseLevel > 2 ) 155 { 164 { 156 G4cout << " -------------------------" << 165 G4cout << " -------------------------" << G4endl; 157 G4cout << " G4PathFinder::ComputeStep - en 166 G4cout << " G4PathFinder::ComputeStep - entered " << G4endl; 158 G4cout << " - stepNo = " << std::setw(4 167 G4cout << " - stepNo = " << std::setw(4) << stepNo << " " 159 << " navigatorId = " << std::setw(2 168 << " navigatorId = " << std::setw(2) << navigatorNo << " " 160 << " proposed step len = " << propo 169 << " proposed step len = " << proposedStepLength << " " << G4endl; 161 G4cout << " PF::ComputeStep requested step 170 G4cout << " PF::ComputeStep requested step " 162 << " from " << InitialFieldTrack.Ge 171 << " from " << InitialFieldTrack.GetPosition() 163 << " dir " << InitialFieldTrack.Ge 172 << " dir " << InitialFieldTrack.GetMomentumDirection() << G4endl; 164 } 173 } 165 #endif 174 #endif 166 #ifdef G4VERBOSE 175 #ifdef G4VERBOSE 167 if( navigatorNo >= fNoActiveNavigators ) 176 if( navigatorNo >= fNoActiveNavigators ) 168 { 177 { 169 std::ostringstream message; 178 std::ostringstream message; 170 message << "Bad Navigator ID !" << G4endl 179 message << "Bad Navigator ID !" << G4endl 171 << " Requested Navigator ID 180 << " Requested Navigator ID = " << navigatorNo << G4endl 172 << " Number of active navig 181 << " Number of active navigators = " << fNoActiveNavigators; 173 G4Exception("G4PathFinder::ComputeStep()", 182 G4Exception("G4PathFinder::ComputeStep()", "GeomNav0002", 174 FatalException, message); 183 FatalException, message); 175 } 184 } 176 #endif 185 #endif 177 186 178 if( fNewTrack || (stepNo != fLastStepNo) ) 187 if( fNewTrack || (stepNo != fLastStepNo) ) 179 { 188 { 180 // This is a new track or a new step, so w 189 // This is a new track or a new step, so we must make the step 181 // ( else we can simply retrieve its resul 190 // ( else we can simply retrieve its results for this Navigator Id ) 182 191 183 G4FieldTrack currentState = InitialFieldTr << 192 G4FieldTrack currentState= InitialFieldTrack; 184 193 185 fCurrentStepNo = stepNo; 194 fCurrentStepNo = stepNo; 186 195 187 // Check whether a process shifted the pos 196 // Check whether a process shifted the position 188 // since the last step -- by physics proce 197 // since the last step -- by physics processes 189 // 198 // 190 G4ThreeVector newPosition = InitialFieldTr 199 G4ThreeVector newPosition = InitialFieldTrack.GetPosition(); 191 G4ThreeVector moveVector = newPosition - f << 200 G4ThreeVector moveVector= newPosition - fLastLocatedPosition; 192 G4double moveLenSq = moveVector.mag2(); << 201 G4double moveLenSq= moveVector.mag2(); 193 if( moveLenSq > sqr(kCarTolerance) ) << 202 if( moveLenSq > kCarTolerance * kCarTolerance ) 194 { 203 { 195 G4ThreeVector newDirection = InitialFie 204 G4ThreeVector newDirection = InitialFieldTrack.GetMomentumDirection(); 196 #ifdef G4DEBUG_PATHFINDER 205 #ifdef G4DEBUG_PATHFINDER 197 if( fVerboseLevel > 2 ) 206 if( fVerboseLevel > 2 ) 198 { 207 { 199 G4double moveLen= std::sqrt( moveLen 208 G4double moveLen= std::sqrt( moveLenSq ); 200 G4cout << " G4PathFinder::ComputeSte 209 G4cout << " G4PathFinder::ComputeStep : Point moved since last step " 201 << " -- at step # = " << step 210 << " -- at step # = " << stepNo << G4endl 202 << " by " << moveLen << " to 211 << " by " << moveLen << " to " << newPosition << G4endl; 203 } 212 } 204 #endif 213 #endif 205 MovePoint(); // Unintentional changed 214 MovePoint(); // Unintentional changed -- ???? 206 215 207 // Relocate to cope with this move -- e 216 // Relocate to cope with this move -- else could abort !? 208 // 217 // 209 Locate( newPosition, newDirection ); 218 Locate( newPosition, newDirection ); 210 } 219 } 211 220 212 // Check whether the particle have an (EM) 221 // Check whether the particle have an (EM) field force exerting upon it 213 // 222 // 214 G4double particleCharge = currentState.Get << 223 G4double particleCharge= currentState.GetCharge(); 215 224 216 G4FieldManager* fieldMgr = nullptr; << 225 G4FieldManager* fieldMgr=0; 217 G4bool fieldExertsForce = false ; 226 G4bool fieldExertsForce = false ; 218 if( particleCharge != 0.0 ) << 227 if( (particleCharge != 0.0) ) 219 { 228 { 220 fieldMgr = fpFieldPropagator->FindAndS << 229 fieldMgr= fpFieldPropagator->FindAndSetFieldManager( currentVolume ); 221 230 222 // Protect for case where field manage 231 // Protect for case where field manager has no field (= field is zero) 223 // 232 // 224 fieldExertsForce = (fieldMgr != nullpt << 233 fieldExertsForce = (fieldMgr != 0) 225 && (fieldMgr->GetDetec << 234 && (fieldMgr->GetDetectorField() != 0); 226 } 235 } 227 fFieldExertedForce = fieldExertsForce; // 236 fFieldExertedForce = fieldExertsForce; // Store for use in later calls 228 // 237 // referring to this 'step'. 229 238 230 fNoGeometriesLimiting = -1; // At start o << 239 fNoGeometriesLimiting= -1; // At start of track, no process limited step 231 if( fieldExertsForce ) 240 if( fieldExertsForce ) 232 { 241 { 233 DoNextCurvedStep( currentState, propose 242 DoNextCurvedStep( currentState, proposedStepLength, currentVolume ); 234 //-------------- 243 //-------------- 235 }else{ 244 }else{ 236 DoNextLinearStep( currentState, propose 245 DoNextLinearStep( currentState, proposedStepLength ); 237 //-------------- 246 //-------------- 238 } 247 } 239 fLastStepNo = stepNo; << 248 fLastStepNo= stepNo; 240 fRelocatedPoint = false; << 241 249 242 #ifdef G4DEBUG_PATHFINDER 250 #ifdef G4DEBUG_PATHFINDER 243 if ( (fNoGeometriesLimiting < 0) 251 if ( (fNoGeometriesLimiting < 0) 244 || (fNoGeometriesLimiting > fNoActiveNav 252 || (fNoGeometriesLimiting > fNoActiveNavigators) ) 245 { 253 { 246 std::ostringstream message; 254 std::ostringstream message; 247 message << "Number of geometries limitin 255 message << "Number of geometries limiting the step not set." << G4endl 248 << " Number of geometries 256 << " Number of geometries limiting step = " 249 << fNoGeometriesLimiting; 257 << fNoGeometriesLimiting; 250 G4Exception("G4PathFinder::ComputeStep() 258 G4Exception("G4PathFinder::ComputeStep()", 251 "GeomNav0002", FatalExceptio 259 "GeomNav0002", FatalException, message); 252 } 260 } 253 #endif 261 #endif 254 } 262 } 255 #ifdef G4DEBUG_PATHFINDER 263 #ifdef G4DEBUG_PATHFINDER 256 else 264 else 257 { 265 { 258 const G4double checkTolerance = 1.0e-9; 266 const G4double checkTolerance = 1.0e-9; 259 if( proposedStepLength < fTrueMinStep * ( << 267 if( proposedStepLength < fTrueMinStep * ( 1.0 + checkTolerance) ) // For 2nd+ geometry 260 { << 268 { 261 std::ostringstream message; 269 std::ostringstream message; 262 message.precision( 12 ); 270 message.precision( 12 ); 263 message << "Problem in step size reques 271 message << "Problem in step size request." << G4endl 264 << " Being requested to 272 << " Being requested to make a step which is shorter" 265 << " than the minimum Step " << 273 << " than the minimum Step " << G4endl 266 << " already computed fo 274 << " already computed for any Navigator/geometry during" 267 << " this tracking-step: " << G 275 << " this tracking-step: " << G4endl 268 << " This could happen d 276 << " This could happen due to an error in process ordering." 269 << G4endl 277 << G4endl 270 << " Check that all phys 278 << " Check that all physics processes are registered" 271 << " before all processe 279 << " before all processes with a navigator/geometry." 272 << G4endl 280 << G4endl 273 << " If using pre-packag 281 << " If using pre-packaged physics list and/or" 274 << " functionality, plea 282 << " functionality, please report this error." 275 << G4endl << G4endl 283 << G4endl << G4endl 276 << " Additional informat 284 << " Additional information for problem: " << G4endl 277 << " Steps request/propo 285 << " Steps request/proposed = " << proposedStepLength 278 << G4endl 286 << G4endl 279 << " MinimumStep (true) 287 << " MinimumStep (true) = " << fTrueMinStep 280 << G4endl 288 << G4endl 281 << " MinimumStep (navraw 289 << " MinimumStep (navraw) = " << fMinStep 282 << G4endl 290 << G4endl 283 << " Navigator raw retur 291 << " Navigator raw return value" << G4endl 284 << " Requested step now 292 << " Requested step now = " << proposedStepLength 285 << G4endl 293 << G4endl 286 << " Difference min-req 294 << " Difference min-req (absolute) = " 287 << fTrueMinStep-proposedStepLen 295 << fTrueMinStep-proposedStepLength << G4endl 288 << " Relative (to max of 296 << " Relative (to max of two) = " 289 << (fTrueMinStep-proposedStepLe 297 << (fTrueMinStep-proposedStepLength) 290 / std::max(proposedStepLengt 298 / std::max(proposedStepLength, fTrueMinStep) << G4endl 291 << " -- Step info> stepNo= 299 << " -- Step info> stepNo= " << stepNo 292 << " last= " << fLastStepNo 300 << " last= " << fLastStepNo 293 << " newTr= " << fNewTrack << G 301 << " newTr= " << fNewTrack << G4endl; 294 G4Exception("G4PathFinder::ComputeStep 302 G4Exception("G4PathFinder::ComputeStep()", 295 "GeomNav0003", FatalExcept 303 "GeomNav0003", FatalException, message); 296 } 304 } 297 else 305 else 298 { 306 { 299 // This is neither a new track nor a n 307 // This is neither a new track nor a new step -- just another 300 // client accessing information for th 308 // client accessing information for the current track, step 301 // We will simply retrieve the results 309 // We will simply retrieve the results of the synchronous 302 // stepping for this Navigator Id belo 310 // stepping for this Navigator Id below. 303 // 311 // 304 if( fVerboseLevel > 1 ) 312 if( fVerboseLevel > 1 ) 305 { 313 { 306 G4cout << " G4P::CS -> Not calling 314 G4cout << " G4P::CS -> Not calling DoNextLinearStep: " 307 << " stepNo= " << stepNo << 315 << " stepNo= " << stepNo << " last= " << fLastStepNo 308 << " new= " << fNewTrack << 316 << " new= " << fNewTrack << " Step already done" << G4endl; 309 } 317 } 310 } 318 } 311 } 319 } 312 #endif 320 #endif 313 321 314 fNewTrack = false; << 322 fNewTrack= false; 315 323 316 // Prepare the information to return 324 // Prepare the information to return 317 325 318 pNewSafety = fCurrentPreStepSafety[ navigat 326 pNewSafety = fCurrentPreStepSafety[ navigatorNo ]; 319 limitedStep = fLimitedStep[ navigatorNo ]; 327 limitedStep = fLimitedStep[ navigatorNo ]; >> 328 fRelocatedPoint= false; 320 329 321 possibleStep = std::min(proposedStepLength, << 330 possibleStep= std::min(proposedStepLength, fCurrentStepSize[ navigatorNo ]); 322 EndState = fEndState; // now corrected for 331 EndState = fEndState; // now corrected for smaller step, if needed 323 332 324 #ifdef G4DEBUG_PATHFINDER 333 #ifdef G4DEBUG_PATHFINDER 325 if( fVerboseLevel > 0 ) 334 if( fVerboseLevel > 0 ) 326 { 335 { 327 G4cout << " G4PathFinder::ComputeStep retu 336 G4cout << " G4PathFinder::ComputeStep returns " 328 << fCurrentStepSize[ navigatorNo ] 337 << fCurrentStepSize[ navigatorNo ] 329 << " for Navigator " << navigatorNo 338 << " for Navigator " << navigatorNo 330 << " Limited step = " << limitedSte 339 << " Limited step = " << limitedStep 331 << " Safety(mm) = " << pNewSafety / 340 << " Safety(mm) = " << pNewSafety / mm 332 << G4endl; 341 << G4endl; 333 } 342 } 334 #endif 343 #endif 335 344 336 return possibleStep; 345 return possibleStep; 337 } 346 } 338 347 339 // ------------------------------------------- 348 // ---------------------------------------------------------------------- 340 349 341 void 350 void 342 G4PathFinder::PrepareNewTrack( const G4ThreeVe 351 G4PathFinder::PrepareNewTrack( const G4ThreeVector& position, 343 const G4ThreeVe 352 const G4ThreeVector& direction, 344 G4VPhysicalVolu << 353 G4VPhysicalVolume* massStartVol) 345 { 354 { 346 // Key purposes: 355 // Key purposes: 347 // - Check and cache set of active navigat 356 // - Check and cache set of active navigators 348 // - Reset state for new track 357 // - Reset state for new track 349 358 350 G4int num=0; 359 G4int num=0; 351 360 352 EnableParallelNavigation(true); 361 EnableParallelNavigation(true); 353 // Switch PropagatorInField to use MultiNa 362 // Switch PropagatorInField to use MultiNavigator 354 363 355 fpTransportManager->GetSafetyHelper()->Initi 364 fpTransportManager->GetSafetyHelper()->InitialiseHelper(); 356 // Reinitialise state of safety helper -- 365 // Reinitialise state of safety helper -- avoid problems with overlaps 357 366 358 fNewTrack = true; << 367 fNewTrack= true; 359 this->MovePoint(); // Signal further that 368 this->MovePoint(); // Signal further that the last status is wiped 360 369 361 fpFieldPropagator->PrepareNewTrack(); // Inf << 370 fpFieldPropagator->PrepareNewTrack(); // Inform field propagator of new track 362 371 363 // Message the G4NavigatorPanel / Dispatcher 372 // Message the G4NavigatorPanel / Dispatcher to find active navigators 364 // 373 // 365 std::vector<G4Navigator*>::iterator pNavigat 374 std::vector<G4Navigator*>::iterator pNavigatorIter; 366 375 367 fNoActiveNavigators = (G4int)fpTransportMana << 376 fNoActiveNavigators= fpTransportManager-> GetNoActiveNavigators(); 368 if( fNoActiveNavigators > fMaxNav ) 377 if( fNoActiveNavigators > fMaxNav ) 369 { 378 { 370 std::ostringstream message; 379 std::ostringstream message; 371 message << "Too many active Navigators / w 380 message << "Too many active Navigators / worlds." << G4endl 372 << " Transportation Manager 381 << " Transportation Manager has " 373 << fNoActiveNavigators << " active 382 << fNoActiveNavigators << " active navigators." << G4endl 374 << " This is more than the 383 << " This is more than the number allowed = " 375 << fMaxNav << " !"; 384 << fMaxNav << " !"; 376 G4Exception("G4PathFinder::PrepareNewTrack 385 G4Exception("G4PathFinder::PrepareNewTrack()", "GeomNav0002", 377 FatalException, message); 386 FatalException, message); 378 } 387 } 379 388 380 fpMultiNavigator->PrepareNavigators(); 389 fpMultiNavigator->PrepareNavigators(); 381 //------------------------------------ 390 //------------------------------------ 382 391 383 pNavigatorIter = fpTransportManager->GetActi << 392 pNavigatorIter= fpTransportManager->GetActiveNavigatorsIterator(); 384 for( num=0; num< fNoActiveNavigators; ++pNav 393 for( num=0; num< fNoActiveNavigators; ++pNavigatorIter,++num ) 385 { 394 { 386 // Keep information in C-array ... for cr 395 // Keep information in C-array ... for creating touchables - at least 387 // << 396 388 fpNavigator[num] = *pNavigatorIter; << 397 fpNavigator[num] = *pNavigatorIter; 389 fLimitTruth[num] = false; 398 fLimitTruth[num] = false; 390 fLimitedStep[num] = kDoNot; 399 fLimitedStep[num] = kDoNot; 391 fCurrentStepSize[num] = 0.0; 400 fCurrentStepSize[num] = 0.0; 392 fLocatedVolume[num] = nullptr; << 401 fLocatedVolume[num] = 0; 393 } 402 } 394 fNoGeometriesLimiting = 0; // At start of t << 403 fNoGeometriesLimiting= 0; // At start of track, no process limited step 395 404 396 // In case of one geometry, the tracking wil 405 // In case of one geometry, the tracking will have done the locating!! 397 406 398 if( fNoActiveNavigators > 1 ) 407 if( fNoActiveNavigators > 1 ) 399 { 408 { 400 Locate( position, direction, false ); 409 Locate( position, direction, false ); 401 } 410 } 402 else 411 else 403 { 412 { 404 // Update state -- depending on the track 413 // Update state -- depending on the tracking's call to Mass Navigator 405 414 406 fLastLocatedPosition = position; << 415 fLastLocatedPosition= position; 407 fLocatedVolume[0] = massStartVol; // This << 416 fLocatedVolume[0]= massStartVol; // This information must be given 408 // by t << 417 // by transportation 409 fLimitedStep[0] = kDoNot; << 418 fLimitedStep[0] = kDoNot; 410 fCurrentStepSize[0] = 0.0; 419 fCurrentStepSize[0] = 0.0; 411 } 420 } 412 421 413 // Reset Safety Information -- as in case of 422 // Reset Safety Information -- as in case of overlaps this can cause 414 // inconsistencies ... 423 // inconsistencies ... 415 // 424 // 416 fMinSafety_PreStepPt = fPreSafetyMinValue = << 425 fMinSafety_PreStepPt= fPreSafetyMinValue= fMinSafety_atSafLocation= 0.0; 417 426 418 for( num=0; num<fNoActiveNavigators; ++num ) << 427 for( num=0; num< fNoActiveNavigators; ++num ) 419 { 428 { 420 fPreSafetyValues[num] = 0.0; << 429 fPreSafetyValues[num]= 0.0; 421 fNewSafetyComputed[num] = 0.0; << 430 fNewSafetyComputed[num]= 0.0; 422 fCurrentPreStepSafety[num] = 0.0; 431 fCurrentPreStepSafety[num] = 0.0; 423 } 432 } 424 433 425 // The first location for each Navigator mus 434 // The first location for each Navigator must be non-relative 426 // or else call ResetStackAndState() for eac 435 // or else call ResetStackAndState() for each Navigator 427 436 428 fRelocatedPoint = false; << 437 fRelocatedPoint= false; 429 } 438 } 430 439 431 440 432 void G4PathFinder::EndTrack() 441 void G4PathFinder::EndTrack() 433 // Signal end of tracking of current track. << 442 // Signal end of tracking of current track. 434 // Reset TransportationManager to use 'ordin << 443 // Reset TransportationManager to use 'ordinary' Navigator 435 // Reset internal state, if needed << 444 // Reset internal state, if needed 436 { 445 { 437 EnableParallelNavigation(false); // Else it 446 EnableParallelNavigation(false); // Else it will be continue to be used 438 } 447 } 439 448 440 void G4PathFinder::ReportMove( const G4ThreeVe 449 void G4PathFinder::ReportMove( const G4ThreeVector& OldVector, 441 const G4ThreeVe 450 const G4ThreeVector& NewVector, 442 const G4String& 451 const G4String& Quantity ) const 443 { 452 { 444 G4ThreeVector moveVec = ( NewVector - OldV 453 G4ThreeVector moveVec = ( NewVector - OldVector ); 445 454 >> 455 G4int prc= G4cerr.precision(12); 446 std::ostringstream message; 456 std::ostringstream message; 447 message.precision(16); << 448 message << "Endpoint moved between value r 457 message << "Endpoint moved between value returned by ComputeStep()" 449 << " and call to Locate(). " << G4 458 << " and call to Locate(). " << G4endl 450 << " Change of " << Quant 459 << " Change of " << Quantity << " is " 451 << moveVec.mag() / mm << " mm long 460 << moveVec.mag() / mm << " mm long" << G4endl 452 << " and its vector is " 461 << " and its vector is " 453 << (1.0/mm) * moveVec << " mm " << 462 << (1.0/mm) * moveVec << " mm " << G4endl 454 << " Endpoint of ComputeS << 463 << " Endpoint of ComputeStep() was " << OldVector << G4endl 455 << G4endl << 456 << " and current position 464 << " and current position to locate is " << NewVector; 457 G4Exception("G4PathFinder::ReportMove()", 465 G4Exception("G4PathFinder::ReportMove()", "GeomNav1002", 458 JustWarning, message); 466 JustWarning, message); >> 467 G4cerr.precision(prc); 459 } 468 } 460 469 461 void G4PathFinder::Locate( const G4ThreeVector << 470 void 462 const G4ThreeVector << 471 G4PathFinder::Locate( const G4ThreeVector& position, 463 G4bool relati << 472 const G4ThreeVector& direction, >> 473 G4bool relative) 464 { 474 { 465 // Locate the point in each geometry 475 // Locate the point in each geometry 466 476 467 auto pNavIter = fpTransportManager->GetActiv << 477 std::vector<G4Navigator*>::iterator pNavIter= 468 << 478 fpTransportManager->GetActiveNavigatorsIterator(); 469 G4ThreeVector lastEndPosition = fRelocatedPo << 470 ? fLastLocated << 471 : fEndState.Ge << 472 fLastLocatedPosition = position; << 473 << 474 #ifdef G4DEBUG_PATHFINDER << 475 static const G4double movLenTol = 10*sqr(kCa << 476 479 477 G4ThreeVector moveVec = ( position - lastEnd << 480 G4ThreeVector lastEndPosition= fEndState.GetPosition(); 478 G4double moveLenSq = moveVec.mag2(); << 481 G4ThreeVector moveVec = (position - lastEndPosition ); 479 if( (!fNewTrack) && ( moveLenSq > movLenTol << 482 G4double moveLenSq= moveVec.mag2(); >> 483 if( (!fNewTrack) && (!fRelocatedPoint) >> 484 && ( moveLenSq> 10*kCarTolerance*kCarTolerance ) ) 480 { 485 { 481 ReportMove( lastEndPosition, position, << 486 ReportMove( position, lastEndPosition, "Position" ); 482 " (End) Position / G4PathFind << 483 } 487 } >> 488 fLastLocatedPosition= position; 484 489 >> 490 #ifdef G4DEBUG_PATHFINDER 485 if( fVerboseLevel > 2 ) 491 if( fVerboseLevel > 2 ) 486 { 492 { 487 G4cout << G4endl; 493 G4cout << G4endl; 488 G4cout << " G4PathFinder::Locate : entered 494 G4cout << " G4PathFinder::Locate : entered " << G4endl; 489 G4cout << " -------------------- ------- 495 G4cout << " -------------------- -------" << G4endl; 490 G4cout << " Locating at position " << po 496 G4cout << " Locating at position " << position 491 << " with direction " << direction 497 << " with direction " << direction 492 << " relative= " << relative << G4 498 << " relative= " << relative << G4endl; 493 if ( (fVerboseLevel > 1) || ( moveLenSq > 499 if ( (fVerboseLevel > 1) || ( moveLenSq > 0.0) ) 494 { 500 { 495 G4cout << " lastEndPosition = " << las 501 G4cout << " lastEndPosition = " << lastEndPosition 496 << " moveVec = " << moveVec 502 << " moveVec = " << moveVec 497 << " newTr = " << fNewTrack 503 << " newTr = " << fNewTrack 498 << " relocated = " << fRelocate 504 << " relocated = " << fRelocatedPoint << G4endl; 499 } 505 } 500 506 501 G4cout << " Located at " << position ; 507 G4cout << " Located at " << position ; 502 if( fNoActiveNavigators > 1 ) { G4cout << 508 if( fNoActiveNavigators > 1 ) { G4cout << G4endl; } 503 } 509 } 504 #endif 510 #endif 505 511 506 for ( auto num=0; num<fNoActiveNavigators ; << 512 for ( G4int num=0; num< fNoActiveNavigators ; ++pNavIter,++num ) 507 { 513 { 508 // ... who limited the step .... 514 // ... who limited the step .... 509 515 510 if( fLimitTruth[num] ) { (*pNavIter)->Set 516 if( fLimitTruth[num] ) { (*pNavIter)->SetGeometricallyLimitedStep(); } 511 517 512 G4VPhysicalVolume *pLocated= 518 G4VPhysicalVolume *pLocated= 513 (*pNavIter)->LocateGlobalPointAndSetup( p 519 (*pNavIter)->LocateGlobalPointAndSetup( position, &direction, 514 r 520 relative, 515 f 521 false); 516 // Set the state related to the location 522 // Set the state related to the location 517 // 523 // 518 fLocatedVolume[num] = pLocated; 524 fLocatedVolume[num] = pLocated; 519 525 520 // Clear state related to the step 526 // Clear state related to the step 521 // 527 // 522 fLimitedStep[num] = kDoNot; << 528 fLimitedStep[num] = kDoNot; 523 fCurrentStepSize[num] = 0.0; 529 fCurrentStepSize[num] = 0.0; 524 530 525 #ifdef G4DEBUG_PATHFINDER 531 #ifdef G4DEBUG_PATHFINDER 526 if( fVerboseLevel > 2 ) 532 if( fVerboseLevel > 2 ) 527 { 533 { 528 G4cout << " - In world " << num << " ge 534 G4cout << " - In world " << num << " geomLimStep= " << fLimitTruth[num] 529 << " gives volume= " << pLocate 535 << " gives volume= " << pLocated ; 530 if( pLocated ) 536 if( pLocated ) 531 { 537 { 532 G4cout << " name = '" << pLocated->G 538 G4cout << " name = '" << pLocated->GetName() << "'"; 533 G4cout << " - CopyNo= " << pLocated-> 539 G4cout << " - CopyNo= " << pLocated->GetCopyNo(); 534 } 540 } 535 G4cout << G4endl; 541 G4cout << G4endl; 536 } 542 } 537 #endif 543 #endif 538 } 544 } 539 545 540 fRelocatedPoint = false; << 546 fRelocatedPoint= false; 541 } 547 } 542 548 543 void G4PathFinder::ReLocate( const G4ThreeVect 549 void G4PathFinder::ReLocate( const G4ThreeVector& position ) 544 { 550 { 545 // Locate the point in each geometry 551 // Locate the point in each geometry 546 552 547 auto pNavIter = fpTransportManager->GetActiv << 553 std::vector<G4Navigator*>::iterator pNavIter= >> 554 fpTransportManager->GetActiveNavigatorsIterator(); 548 555 549 #ifdef G4DEBUG_PATHFINDER 556 #ifdef G4DEBUG_PATHFINDER 550 557 551 // Check that this relocation does not viola 558 // Check that this relocation does not violate safety 552 // - at endpoint (computed from start poin 559 // - at endpoint (computed from start point) AND 553 // - at last safety location (likely just 560 // - at last safety location (likely just called) 554 561 555 G4ThreeVector lastEndPosition = fEndState.Ge << 562 G4ThreeVector lastEndPosition= fEndState.GetPosition(); 556 563 557 // Calculate end-point safety ... 564 // Calculate end-point safety ... 558 // 565 // 559 G4double DistanceStartEnd = (lastEndPosition << 566 G4double DistanceStartEnd= (lastEndPosition - fPreStepLocation).mag(); 560 G4double endPointSafety_raw = fMinSafety_Pre << 567 G4double endPointSafety_raw = fMinSafety_PreStepPt - DistanceStartEnd; 561 G4double endPointSafety_Est1 = std::max( 0.0 << 568 G4double endPointSafety_Est1 = std::max( 0.0, endPointSafety_raw ); 562 569 563 // ... and check move from endpoint against 570 // ... and check move from endpoint against this endpoint safety 564 // 571 // 565 G4ThreeVector moveVecEndPos = position - la 572 G4ThreeVector moveVecEndPos = position - lastEndPosition; 566 G4double moveLenEndPosSq = moveVecEndPo 573 G4double moveLenEndPosSq = moveVecEndPos.mag2(); 567 574 568 // Check that move from endpoint of last ste 575 // Check that move from endpoint of last step is within safety 569 // -- or check against last location or relo 576 // -- or check against last location or relocation ?? 570 // 577 // 571 G4ThreeVector moveVecSafety = position - fSa << 578 G4ThreeVector moveVecSafety= position - fSafetyLocation; 572 G4double moveLenSafSq = moveVecSafety.m << 579 G4double moveLenSafSq= moveVecSafety.mag2(); 573 580 574 G4double distCheckEnd_sq = ( moveLenEndPosSq << 581 G4double distCheckEnd_sq= ( moveLenEndPosSq - endPointSafety_Est1 575 << 582 *endPointSafety_Est1 ); 576 G4double distCheckSaf_sq = ( moveLenSafSq - << 583 G4double distCheckSaf_sq= ( moveLenSafSq - fMinSafety_atSafLocation 577 * << 584 *fMinSafety_atSafLocation ); 578 585 579 G4bool longMoveEnd = distCheckEnd_sq > 0.0; 586 G4bool longMoveEnd = distCheckEnd_sq > 0.0; 580 G4bool longMoveSaf = distCheckSaf_sq > 0.0; 587 G4bool longMoveSaf = distCheckSaf_sq > 0.0; 581 588 582 G4double revisedSafety = 0.0; << 589 G4double revisedSafety= 0.0; 583 590 584 if( (!fNewTrack) && ( longMoveEnd && longMov 591 if( (!fNewTrack) && ( longMoveEnd && longMoveSaf ) ) 585 { 592 { 586 // Recompute ComputeSafety for end positi 593 // Recompute ComputeSafety for end position 587 // 594 // 588 revisedSafety = ComputeSafety(lastEndPosi << 595 revisedSafety= ComputeSafety(lastEndPosition); 589 596 590 const G4double kRadTolerance = 597 const G4double kRadTolerance = 591 G4GeometryTolerance::GetInstance()- << 598 G4GeometryTolerance::GetInstance()->GetRadialTolerance(); 592 const G4double cErrorTolerance = 1e-12; << 599 const G4double cErrorTolerance=1e-12; 593 // Maximum relative error from roundoff 600 // Maximum relative error from roundoff of arithmetic 594 601 595 G4double distCheckRevisedEnd = moveLenEnd << 602 G4double distCheckRevisedEnd= moveLenEndPosSq-revisedSafety*revisedSafety; 596 603 597 G4bool longMoveRevisedEnd = ( distCheckR << 604 G4bool longMoveRevisedEnd= ( distCheckRevisedEnd > 0. ) ; 598 605 599 G4double moveMinusSafety = 0.0; << 606 G4double moveMinusSafety= 0.0; 600 G4double moveLenEndPosition = std::sqrt( << 607 G4double moveLenEndPosition= std::sqrt( moveLenEndPosSq ); 601 moveMinusSafety = moveLenEndPosition - re 608 moveMinusSafety = moveLenEndPosition - revisedSafety; 602 609 603 if ( longMoveRevisedEnd && ( moveMinusSaf << 610 if ( longMoveRevisedEnd && (moveMinusSafety > 0.0 ) 604 && ( revisedSafet << 611 && ( revisedSafety > 0.0 ) ) 605 { 612 { 606 // Take into account possibility of ro 613 // Take into account possibility of roundoff error causing 607 // this apparent move further than saf 614 // this apparent move further than safety 608 615 609 if( fVerboseLevel > 0 ) 616 if( fVerboseLevel > 0 ) 610 { 617 { 611 G4cout << " G4PF:Relocate> Ratio to 618 G4cout << " G4PF:Relocate> Ratio to revised safety is " 612 << std::fabs(moveMinusSafety 619 << std::fabs(moveMinusSafety)/revisedSafety << G4endl; 613 } 620 } 614 621 615 G4double absMoveMinusSafety = std::fab << 622 G4double absMoveMinusSafety= std::fabs(moveMinusSafety); 616 G4bool smallRatio = absMoveMinusSafety << 623 G4bool smallRatio= absMoveMinusSafety < kRadTolerance * revisedSafety ; 617 G4double maxCoordPos = std::max( 624 G4double maxCoordPos = std::max( 618 std::max 625 std::max( std::fabs(position.x()), 619 626 std::fabs(position.y())), 620 std::fab 627 std::fabs(position.z()) ); 621 G4bool smallValue= absMoveMinusSafety 628 G4bool smallValue= absMoveMinusSafety < cErrorTolerance * maxCoordPos; 622 if( !(smallRatio || smallValue) ) << 629 if( ! (smallRatio || smallValue) ) 623 { 630 { 624 G4cout << " G4PF:Relocate> Ratio to 631 G4cout << " G4PF:Relocate> Ratio to revised safety is " 625 << std::fabs(moveMinusSafety 632 << std::fabs(moveMinusSafety)/revisedSafety << G4endl; 626 G4cout << " Difference of move and 633 G4cout << " Difference of move and safety is not very small." 627 << G4endl; 634 << G4endl; 628 } 635 } 629 else 636 else 630 { 637 { 631 moveMinusSafety = 0.0; 638 moveMinusSafety = 0.0; 632 longMoveRevisedEnd = false; // Num 639 longMoveRevisedEnd = false; // Numerical issue -- not too long! 633 640 634 G4cout << " Difference of move & saf 641 G4cout << " Difference of move & safety is very small in magnitude, " 635 << absMoveMinusSafety << G4en 642 << absMoveMinusSafety << G4endl; 636 if( smallRatio ) 643 if( smallRatio ) 637 { 644 { 638 G4cout << " ratio to safety " << r 645 G4cout << " ratio to safety " << revisedSafety 639 << " is " << absMoveMinusS 646 << " is " << absMoveMinusSafety / revisedSafety 640 << "smaller than " << kRadT 647 << "smaller than " << kRadTolerance << " of safety "; 641 } 648 } 642 else 649 else 643 { 650 { 644 G4cout << " as fraction " << absMo 651 G4cout << " as fraction " << absMoveMinusSafety / maxCoordPos 645 << " of position vector max 652 << " of position vector max-coord " << maxCoordPos 646 << " smaller than " << cErr 653 << " smaller than " << cErrorTolerance ; 647 } 654 } 648 G4cout << " -- reset moveMinusSafety 655 G4cout << " -- reset moveMinusSafety to " 649 << moveMinusSafety << G4endl; 656 << moveMinusSafety << G4endl; 650 } 657 } 651 } 658 } 652 659 653 if ( longMoveEnd && longMoveSaf 660 if ( longMoveEnd && longMoveSaf 654 && longMoveRevisedEnd && (moveMinusSafe 661 && longMoveRevisedEnd && (moveMinusSafety>0.0) ) 655 { 662 { >> 663 G4int oldPrec= G4cout.precision(9); 656 std::ostringstream message; 664 std::ostringstream message; 657 message.precision(9); << 658 message << "ReLocation is further than 665 message << "ReLocation is further than end-safety value." << G4endl 659 << " Moved from last endpoint 666 << " Moved from last endpoint by " << moveLenEndPosition 660 << " compared to end safety (f 667 << " compared to end safety (from preStep point) = " 661 << endPointSafety_Est1 << G4en 668 << endPointSafety_Est1 << G4endl 662 << " --> last PreSafety Locat 669 << " --> last PreSafety Location was " << fPreSafetyLocation 663 << G4endl 670 << G4endl 664 << " safety value = " < 671 << " safety value = " << fPreSafetyMinValue << G4endl 665 << " --> last PreStep Locatio 672 << " --> last PreStep Location was " << fPreStepLocation 666 << G4endl 673 << G4endl 667 << " safety value = " < 674 << " safety value = " << fMinSafety_PreStepPt << G4endl 668 << " --> last EndStep Locatio 675 << " --> last EndStep Location was " << lastEndPosition 669 << G4endl 676 << G4endl 670 << " safety value = " < 677 << " safety value = " << endPointSafety_Est1 671 << " raw-value = " << endPoint 678 << " raw-value = " << endPointSafety_raw << G4endl 672 << " --> Calling again at thi 679 << " --> Calling again at this endpoint, we get " 673 << revisedSafety << " as safe 680 << revisedSafety << " as safety value." << G4endl 674 << " --> last position for sa 681 << " --> last position for safety " << fSafetyLocation 675 << G4endl 682 << G4endl 676 << " its safety value = 683 << " its safety value = " << fMinSafety_atSafLocation 677 << G4endl 684 << G4endl 678 << " move from safety lo 685 << " move from safety location = " 679 << std::sqrt(moveLenSafSq) << 686 << std::sqrt(moveLenSafSq) << G4endl 680 << " again= " << moveV 687 << " again= " << moveVecSafety.mag() << G4endl 681 << " safety - Move-from- 688 << " safety - Move-from-end= " 682 << revisedSafety - moveLenEndP 689 << revisedSafety - moveLenEndPosition 683 << " (negative is Bad.)" << G4 690 << " (negative is Bad.)" << G4endl 684 << " Debug: distCheckRevisedE 691 << " Debug: distCheckRevisedEnd = " 685 << distCheckRevisedEnd; 692 << distCheckRevisedEnd; 686 ReportMove( lastEndPosition, position, 693 ReportMove( lastEndPosition, position, "Position" ); 687 G4Exception("G4PathFinder::ReLocate", 694 G4Exception("G4PathFinder::ReLocate", "GeomNav0003", 688 FatalException, message); 695 FatalException, message); >> 696 G4cout.precision(oldPrec); 689 } 697 } 690 } 698 } 691 699 692 if( fVerboseLevel > 2 ) 700 if( fVerboseLevel > 2 ) 693 { 701 { 694 G4cout << G4endl; 702 G4cout << G4endl; 695 G4cout << " G4PathFinder::ReLocate : enter 703 G4cout << " G4PathFinder::ReLocate : entered " << G4endl; 696 G4cout << " ---------------------- ----- 704 G4cout << " ---------------------- -------" << G4endl; 697 G4cout << " *Re*Locating at position " << 705 G4cout << " *Re*Locating at position " << position << G4endl; >> 706 // << " with direction " << direction >> 707 // << " relative= " << relative << G4endl; 698 if ( (fVerboseLevel > -1) || ( moveLenEndP 708 if ( (fVerboseLevel > -1) || ( moveLenEndPosSq > 0.0) ) 699 { 709 { 700 G4cout << " lastEndPosition = " << las 710 G4cout << " lastEndPosition = " << lastEndPosition 701 << " moveVec from step-end = " 711 << " moveVec from step-end = " << moveVecEndPos 702 << " is new Track = " << fNewTr 712 << " is new Track = " << fNewTrack 703 << " relocated = " << fRelocate 713 << " relocated = " << fRelocatedPoint << G4endl; 704 } 714 } 705 } 715 } 706 #endif // G4DEBUG_PATHFINDER 716 #endif // G4DEBUG_PATHFINDER 707 717 708 for ( auto num=0; num< fNoActiveNavigators ; << 718 for ( G4int num=0; num< fNoActiveNavigators ; ++pNavIter,++num ) 709 { 719 { 710 // ... none limited the step 720 // ... none limited the step 711 721 712 (*pNavIter)->LocateGlobalPointWithinVolum 722 (*pNavIter)->LocateGlobalPointWithinVolume( position ); 713 723 714 // Clear state related to the step 724 // Clear state related to the step 715 // 725 // 716 fLimitedStep[num] = kDoNot; << 726 fLimitedStep[num] = kDoNot; 717 fCurrentStepSize[num] = 0.0; 727 fCurrentStepSize[num] = 0.0; 718 fLimitTruth[num] = false; 728 fLimitTruth[num] = false; 719 } 729 } 720 730 721 fLastLocatedPosition = position; << 731 fLastLocatedPosition= position; 722 fRelocatedPoint = true; << 732 fRelocatedPoint= false; 723 733 724 #ifdef G4DEBUG_PATHFINDER 734 #ifdef G4DEBUG_PATHFINDER 725 if( fVerboseLevel > 2 ) 735 if( fVerboseLevel > 2 ) 726 { 736 { 727 G4cout << " G4PathFinder::ReLocate : exiti 737 G4cout << " G4PathFinder::ReLocate : exiting " 728 << " at position " << fLastLocated 738 << " at position " << fLastLocatedPosition << G4endl << G4endl; 729 } 739 } 730 #endif 740 #endif 731 } 741 } 732 742 733 // ------------------------------------------- 743 // ----------------------------------------------------------------------------- 734 744 735 G4double G4PathFinder::ComputeSafety( const G 745 G4double G4PathFinder::ComputeSafety( const G4ThreeVector& position ) 736 { 746 { 737 // Recompute safety for the relevant point 747 // Recompute safety for the relevant point 738 748 739 G4double minSafety = kInfinity; << 749 G4double minSafety= kInfinity; 740 750 741 std::vector<G4Navigator*>::iterator pNaviga 751 std::vector<G4Navigator*>::iterator pNavigatorIter; 742 pNavigatorIter = fpTransportManager->GetAct << 752 pNavigatorIter= fpTransportManager->GetActiveNavigatorsIterator(); 743 753 744 for( auto num=0; num<fNoActiveNavigators; + << 754 for( G4int num=0; num<fNoActiveNavigators; ++pNavigatorIter,++num ) 745 { 755 { 746 G4double safety = (*pNavigatorIter)->Com << 756 G4double safety = (*pNavigatorIter)->ComputeSafety( position, DBL_MAX, true ); 747 if( safety < minSafety ) { minSafety = s 757 if( safety < minSafety ) { minSafety = safety; } 748 fNewSafetyComputed[num] = safety; << 758 fNewSafetyComputed[num]= safety; 749 } 759 } 750 760 751 fSafetyLocation = position; << 761 fSafetyLocation= position; 752 fMinSafety_atSafLocation = minSafety; 762 fMinSafety_atSafLocation = minSafety; 753 763 754 #ifdef G4DEBUG_PATHFINDER 764 #ifdef G4DEBUG_PATHFINDER 755 if( fVerboseLevel > 1 ) 765 if( fVerboseLevel > 1 ) 756 { 766 { 757 G4cout << " G4PathFinder::ComputeSafety - 767 G4cout << " G4PathFinder::ComputeSafety - returns " 758 << minSafety << " at location " << 768 << minSafety << " at location " << position << G4endl; 759 } 769 } 760 #endif 770 #endif 761 return minSafety; 771 return minSafety; 762 } 772 } 763 773 764 774 765 // ------------------------------------------- 775 // ----------------------------------------------------------------------------- 766 776 767 G4TouchableHandle 777 G4TouchableHandle 768 G4PathFinder::CreateTouchableHandle( G4int nav 778 G4PathFinder::CreateTouchableHandle( G4int navId ) const 769 { 779 { 770 #ifdef G4DEBUG_PATHFINDER 780 #ifdef G4DEBUG_PATHFINDER 771 if( fVerboseLevel > 2 ) 781 if( fVerboseLevel > 2 ) 772 { 782 { 773 G4cout << "G4PathFinder::CreateTouchableHa 783 G4cout << "G4PathFinder::CreateTouchableHandle : navId = " 774 << navId << " -- " << GetNavigator( 784 << navId << " -- " << GetNavigator(navId) << G4endl; 775 } 785 } 776 #endif 786 #endif 777 787 778 G4TouchableHistory* touchHist; 788 G4TouchableHistory* touchHist; 779 touchHist = GetNavigator(navId)->CreateTouch << 789 touchHist= GetNavigator(navId) -> CreateTouchableHistory(); 780 790 781 G4VPhysicalVolume* locatedVolume = fLocatedV << 791 G4VPhysicalVolume* locatedVolume= fLocatedVolume[navId]; 782 if( locatedVolume == nullptr ) << 792 if( locatedVolume == 0 ) 783 { 793 { 784 // Workaround to ensure that the touchabl 794 // Workaround to ensure that the touchable is fixed !! // TODO: fix 785 795 786 touchHist->UpdateYourself( locatedVolume, 796 touchHist->UpdateYourself( locatedVolume, touchHist->GetHistory() ); 787 } 797 } 788 798 789 #ifdef G4DEBUG_PATHFINDER 799 #ifdef G4DEBUG_PATHFINDER 790 if( fVerboseLevel > 2 ) 800 if( fVerboseLevel > 2 ) 791 { 801 { 792 G4String VolumeName("None"); 802 G4String VolumeName("None"); 793 if( locatedVolume ) { VolumeName = located << 803 if( locatedVolume ) { VolumeName= locatedVolume->GetName(); } 794 G4cout << " Touchable History created at a 804 G4cout << " Touchable History created at address " << touchHist 795 << "; volume = " << locatedVolume < << 805 << " volume = " << locatedVolume << " name= " << VolumeName 796 << G4endl; 806 << G4endl; 797 } 807 } 798 #endif 808 #endif 799 809 800 return {touchHist}; << 810 return G4TouchableHandle(touchHist); 801 } 811 } 802 812 803 G4double 813 G4double 804 G4PathFinder::DoNextLinearStep( const G4FieldT << 814 G4PathFinder::DoNextLinearStep( const G4FieldTrack &initialState, 805 G4double 815 G4double proposedStepLength ) 806 { 816 { 807 std::vector<G4Navigator*>::iterator pNavigat 817 std::vector<G4Navigator*>::iterator pNavigatorIter; 808 G4double safety = 0.0, step =0.0; << 818 G4double safety= 0.0, step=0.0; 809 G4double minSafety = kInfinity, minStep = kI << 819 G4double minSafety= kInfinity, minStep; 810 << 820 811 const G4int IdTransport = 0; // Id of Mass << 821 const G4int IdTransport= 0; // Id of Mass Navigator !! 812 G4int num = 0; << 822 G4int num=0; 813 823 814 #ifdef G4DEBUG_PATHFINDER 824 #ifdef G4DEBUG_PATHFINDER 815 if( fVerboseLevel > 2 ) 825 if( fVerboseLevel > 2 ) 816 { 826 { 817 G4cout << " G4PathFinder::DoNextLinearStep 827 G4cout << " G4PathFinder::DoNextLinearStep : entered " << G4endl; 818 G4cout << " Input field track= " << init 828 G4cout << " Input field track= " << initialState << G4endl; 819 G4cout << " Requested step= " << propose 829 G4cout << " Requested step= " << proposedStepLength << G4endl; 820 } 830 } 821 #endif 831 #endif 822 832 823 G4ThreeVector initialPosition = initialState << 833 G4ThreeVector initialPosition= initialState.GetPosition(); 824 G4ThreeVector initialDirection = initialStat << 834 G4ThreeVector initialDirection= initialState.GetMomentumDirection(); 825 835 826 G4ThreeVector OriginShift = initialPosition 836 G4ThreeVector OriginShift = initialPosition - fPreSafetyLocation; 827 G4double MagSqShift = OriginShift.mag2 837 G4double MagSqShift = OriginShift.mag2() ; 828 G4double MagShift; // Only given value 838 G4double MagShift; // Only given value if it larger than minimum safety 829 839 830 // Potential optimisation using Maximum Valu 840 // Potential optimisation using Maximum Value of safety! 831 // if( MagSqShift >= sqr(fPreSafetyMaxValue 841 // if( MagSqShift >= sqr(fPreSafetyMaxValue ) ){ 832 // MagShift= kInfinity; // Not a useful 842 // MagShift= kInfinity; // Not a useful value -- all will not use/ignore 833 // else 843 // else 834 // MagShift= std::sqrt(MagSqShift) ; 844 // MagShift= std::sqrt(MagSqShift) ; 835 845 836 MagShift= std::sqrt(MagSqShift) ; 846 MagShift= std::sqrt(MagSqShift) ; 837 847 838 #ifdef G4PATHFINDER_OPTIMISATION 848 #ifdef G4PATHFINDER_OPTIMISATION 839 849 840 G4double fullSafety; // For all geometries, 850 G4double fullSafety; // For all geometries, for prestep point 841 851 842 if( MagSqShift >= sqr(fPreSafetyMinValue) ) << 852 if( MagSqShift >= sqr(fPreSafetyMinValue ) ) 843 { 853 { 844 fullSafety = 0.0 ; 854 fullSafety = 0.0 ; 845 } 855 } 846 else 856 else 847 { 857 { 848 fullSafety = fPreSafetyMinValue - MagShif 858 fullSafety = fPreSafetyMinValue - MagShift; 849 } 859 } 850 if( proposedStepLength < fullSafety ) 860 if( proposedStepLength < fullSafety ) 851 { 861 { 852 // Move is smaller than all safeties 862 // Move is smaller than all safeties 853 // -> so we do not have to move the safe 863 // -> so we do not have to move the safety center 854 864 855 fPreStepCenterRenewed = false; << 865 fPreStepCenterRenewed= false; 856 866 857 for( num=0; num< fNoActiveNavigators; ++n 867 for( num=0; num< fNoActiveNavigators; ++num ) 858 { 868 { 859 fCurrentStepSize[num] = kInfinity; << 869 fCurrentStepSize[num]= kInfinity; 860 safety = std::max( 0.0, fPreSafetyVal 870 safety = std::max( 0.0, fPreSafetyValues[num] - MagShift); 861 minSafety= std::min( safety, minSafety 871 minSafety= std::min( safety, minSafety ); 862 fCurrentPreStepSafety[num] = safety; << 872 fCurrentPreStepSafety[num]= safety; 863 } 873 } 864 minStep = kInfinity; << 874 minStep= kInfinity; 865 875 866 #ifdef G4DEBUG_PATHFINDER 876 #ifdef G4DEBUG_PATHFINDER 867 if( fVerboseLevel > 2 ) 877 if( fVerboseLevel > 2 ) 868 { 878 { 869 G4cout << "G4PathFinder::DoNextLinearSt 879 G4cout << "G4PathFinder::DoNextLinearStep : Quick Stepping. " << G4endl 870 << " proposedStepLength " << p 880 << " proposedStepLength " << proposedStepLength 871 << " < (full) safety = " << ful 881 << " < (full) safety = " << fullSafety 872 << " at " << initialPosition 882 << " at " << initialPosition 873 << G4endl; 883 << G4endl; 874 } 884 } 875 #endif 885 #endif 876 } 886 } 877 else 887 else 878 #endif // End of G4PATHFINDER_OPTIMISATION 1 888 #endif // End of G4PATHFINDER_OPTIMISATION 1 879 { 889 { 880 // Move is larger than at least one of th 890 // Move is larger than at least one of the safeties 881 // -> so we must move the safety center! 891 // -> so we must move the safety center! 882 892 883 fPreStepCenterRenewed = true; << 893 fPreStepCenterRenewed= true; 884 pNavigatorIter = fpTransportManager-> Get << 894 pNavigatorIter= fpTransportManager-> GetActiveNavigatorsIterator(); 885 895 886 minStep = kInfinity; // Not proposedStep << 896 minStep= kInfinity; // Not proposedStepLength; 887 897 888 for( num=0; num< fNoActiveNavigators; ++p 898 for( num=0; num< fNoActiveNavigators; ++pNavigatorIter,++num ) 889 { 899 { 890 safety = std::max( 0.0, fPreSafetyVal 900 safety = std::max( 0.0, fPreSafetyValues[num] - MagShift); 891 901 892 #ifdef G4PATHFINDER_OPTIMISATION 902 #ifdef G4PATHFINDER_OPTIMISATION 893 if( proposedStepLength <= safety ) // 903 if( proposedStepLength <= safety ) // Should be just < safety ? 894 { 904 { 895 // The Step is guaranteed to be tak 905 // The Step is guaranteed to be taken 896 906 897 step = kInfinity; // ComputeSte << 907 step= kInfinity; // ComputeStep Would return this 898 908 899 #ifdef G4DEBUG_PATHFINDER 909 #ifdef G4DEBUG_PATHFINDER 900 G4cout.precision(8); 910 G4cout.precision(8); 901 G4cout << "PathFinder::ComputeStep> 911 G4cout << "PathFinder::ComputeStep> small proposed step = " 902 << proposedStepLength 912 << proposedStepLength 903 << " <= safety = " << safet 913 << " <= safety = " << safety << " for nav " << num 904 << " Step fully taken. " << 914 << " Step fully taken. " << G4endl; 905 #endif 915 #endif 906 } 916 } 907 else 917 else 908 #endif // End of G4PATHFINDER_OPTIMISATION 2 918 #endif // End of G4PATHFINDER_OPTIMISATION 2 909 { 919 { 910 #ifdef G4DEBUG_PATHFINDER 920 #ifdef G4DEBUG_PATHFINDER 911 G4double previousSafety = safety; << 921 G4double previousSafety= safety; 912 #endif 922 #endif 913 step = (*pNavigatorIter)->ComputeSt << 923 step= (*pNavigatorIter)->ComputeStep( initialPosition, 914 << 924 initialDirection, 915 << 925 proposedStepLength, 916 << 926 safety ); 917 minStep = std::min(step, minStep); << 927 minStep = std::min( step, minStep); 918 << 928 919 << 929 // TODO: consider whether/how to reduce the proposed step >> 930 // to the latest minStep value - to reduce calculations >> 931 920 #ifdef G4DEBUG_PATHFINDER 932 #ifdef G4DEBUG_PATHFINDER 921 if( fVerboseLevel > 0) 933 if( fVerboseLevel > 0) 922 { 934 { 923 G4cout.precision(8); 935 G4cout.precision(8); 924 G4cout << "PathFinder::ComputeSte 936 G4cout << "PathFinder::ComputeStep> long proposed step = " 925 << proposedStepLength 937 << proposedStepLength 926 << " > safety = " << pre 938 << " > safety = " << previousSafety 927 << " for nav " << num 939 << " for nav " << num 928 << " . New safety = " << 940 << " . New safety = " << safety << " step= " << step 929 << G4endl; 941 << G4endl; 930 } 942 } 931 #endif 943 #endif 932 } 944 } 933 fCurrentStepSize[num] = step; // Raw << 945 fCurrentStepSize[num] = step; 934 946 935 // TODO: consider whether/how to r << 936 // to the latest minStep val << 937 // ( If so, much change 1st << 938 << 939 // Save safety value, must be done for 947 // Save safety value, must be done for all geometries "together" 940 // (even if not recomputed using call 948 // (even if not recomputed using call to ComputeStep) 941 // since they share the fPreSafetyLoca 949 // since they share the fPreSafetyLocation 942 950 943 fPreSafetyValues[num] = safety; << 951 fPreSafetyValues[num]= safety; 944 fCurrentPreStepSafety[num] = safety; << 952 fCurrentPreStepSafety[num]= safety; 945 953 946 minSafety = std::min( safety, minSafet << 954 minSafety= std::min( safety, minSafety ); 947 955 948 #ifdef G4DEBUG_PATHFINDER 956 #ifdef G4DEBUG_PATHFINDER 949 if( fVerboseLevel > 2 ) 957 if( fVerboseLevel > 2 ) 950 { 958 { 951 G4cout << "G4PathFinder::DoNextLinea 959 G4cout << "G4PathFinder::DoNextLinearStep : Navigator [" 952 << num << "] -- step size " < 960 << num << "] -- step size " << step << G4endl; 953 } 961 } 954 #endif 962 #endif 955 } 963 } 956 964 957 // Only change these when safety is recal 965 // Only change these when safety is recalculated 958 // it is good/relevant only for safety ca 966 // it is good/relevant only for safety calculations 959 967 960 fPreSafetyLocation = initialPosition; << 968 fPreSafetyLocation= initialPosition; 961 fPreSafetyMinValue = minSafety; << 969 fPreSafetyMinValue= minSafety; 962 } // end of else for if( proposedStepLength 970 } // end of else for if( proposedStepLength <= fullSafety) 963 971 964 // For use in Relocation, need PreStep point 972 // For use in Relocation, need PreStep point location, min-safety 965 // 973 // 966 fPreStepLocation = initialPosition; << 974 fPreStepLocation= initialPosition; 967 fMinSafety_PreStepPt = minSafety; << 975 fMinSafety_PreStepPt= minSafety; 968 976 969 fMinStep = minStep; << 977 fMinStep= minStep; 970 978 971 if( fMinStep == kInfinity ) 979 if( fMinStep == kInfinity ) 972 { 980 { 973 minStep = proposedStepLength; // Use t 981 minStep = proposedStepLength; // Use this below for endpoint !! 974 } 982 } 975 fTrueMinStep = minStep; 983 fTrueMinStep = minStep; 976 984 977 // Set the EndState 985 // Set the EndState 978 986 979 G4ThreeVector endPosition; 987 G4ThreeVector endPosition; 980 988 981 fEndState = initialState; << 989 fEndState= initialState; 982 endPosition = initialPosition + minStep * in << 990 endPosition= initialPosition + minStep * initialDirection ; 983 991 984 #ifdef G4DEBUG_PATHFINDER 992 #ifdef G4DEBUG_PATHFINDER 985 if( fVerboseLevel > 1 ) 993 if( fVerboseLevel > 1 ) 986 { 994 { 987 G4int oldPrec= G4cout.precision(14); << 988 G4cout << "G4PathFinder::DoNextLinearStep 995 G4cout << "G4PathFinder::DoNextLinearStep : " 989 << " initialPosition = " << initial 996 << " initialPosition = " << initialPosition 990 << " and endPosition = " << endPosi 997 << " and endPosition = " << endPosition<< G4endl; 991 G4cout.precision(oldPrec); << 992 } 998 } 993 #endif 999 #endif 994 1000 995 fEndState.SetPosition( endPosition ); 1001 fEndState.SetPosition( endPosition ); 996 fEndState.SetProperTimeOfFlight( -1.000 ); 1002 fEndState.SetProperTimeOfFlight( -1.000 ); // Not defined YET 997 1003 998 if( fNoActiveNavigators == 1 ) 1004 if( fNoActiveNavigators == 1 ) 999 { 1005 { 1000 G4bool transportLimited = (fMinStep!= kI 1006 G4bool transportLimited = (fMinStep!= kInfinity); 1001 fLimitTruth[IdTransport] = transportLimi 1007 fLimitTruth[IdTransport] = transportLimited; 1002 fLimitedStep[IdTransport] = transportLim 1008 fLimitedStep[IdTransport] = transportLimited ? kUnique : kDoNot; 1003 1009 1004 // Set fNoGeometriesLimiting - as WhichL 1010 // Set fNoGeometriesLimiting - as WhichLimited does 1005 fNoGeometriesLimiting = transportLimited 1011 fNoGeometriesLimiting = transportLimited ? 1 : 0; 1006 } 1012 } 1007 else 1013 else 1008 { 1014 { 1009 WhichLimited(); 1015 WhichLimited(); 1010 } 1016 } 1011 1017 1012 #ifdef G4DEBUG_PATHFINDER 1018 #ifdef G4DEBUG_PATHFINDER 1013 if( fVerboseLevel > 2 ) 1019 if( fVerboseLevel > 2 ) 1014 { 1020 { 1015 G4cout << " G4PathFinder::DoNextLinearSte 1021 G4cout << " G4PathFinder::DoNextLinearStep : exits returning " 1016 << minStep << G4endl; 1022 << minStep << G4endl; 1017 G4cout << " - Endpoint values = " << fEnd << 1023 G4cout << " Endpoint values = " << fEndState << G4endl; 1018 G4cout << G4endl; 1024 G4cout << G4endl; 1019 } 1025 } 1020 #endif 1026 #endif 1021 1027 1022 return minStep; 1028 return minStep; 1023 } 1029 } 1024 1030 1025 void G4PathFinder::WhichLimited() 1031 void G4PathFinder::WhichLimited() 1026 { 1032 { 1027 // Flag which processes limited the step 1033 // Flag which processes limited the step 1028 1034 1029 G4int num = -1, last = -1; << 1035 G4int num=-1, last=-1; 1030 G4int noLimited = 0; << 1036 G4int noLimited=0; 1031 ELimited shared = kSharedOther; << 1037 ELimited shared= kSharedOther; 1032 1038 1033 const G4int IdTransport = 0; // Id of Mass << 1039 const G4int IdTransport= 0; // Id of Mass Navigator !! 1034 1040 1035 // Assume that [IdTransport] is Mass / Tran 1041 // Assume that [IdTransport] is Mass / Transport 1036 // 1042 // 1037 G4bool transportLimited = (fCurrentStepSize 1043 G4bool transportLimited = (fCurrentStepSize[IdTransport] == fMinStep) 1038 && (fMinStep != kInf << 1044 && ( fMinStep!= kInfinity) ; 1039 1045 1040 if( transportLimited ) << 1046 if( transportLimited ) { 1041 { << 1042 shared= kSharedTransport; 1047 shared= kSharedTransport; 1043 } 1048 } 1044 1049 1045 for ( num = 0; num < fNoActiveNavigators; + << 1050 for ( num= 0; num < fNoActiveNavigators; num++ ) 1046 { 1051 { 1047 G4bool limitedStep; 1052 G4bool limitedStep; 1048 1053 1049 G4double step = fCurrentStepSize[num]; << 1054 G4double step= fCurrentStepSize[num]; 1050 1055 1051 limitedStep = ( std::fabs(step - fMinStep 1056 limitedStep = ( std::fabs(step - fMinStep) < kCarTolerance ) 1052 && ( step != kInfinity); << 1057 && ( step != kInfinity); 1053 1058 1054 fLimitTruth[ num ] = limitedStep; 1059 fLimitTruth[ num ] = limitedStep; 1055 if( limitedStep ) 1060 if( limitedStep ) 1056 { 1061 { 1057 ++noLimited; << 1062 noLimited++; 1058 fLimitedStep[num] = shared; 1063 fLimitedStep[num] = shared; 1059 last= num; 1064 last= num; 1060 } 1065 } 1061 else 1066 else 1062 { 1067 { 1063 fLimitedStep[num] = kDoNot; 1068 fLimitedStep[num] = kDoNot; 1064 } 1069 } 1065 } 1070 } 1066 fNoGeometriesLimiting= noLimited; // Save 1071 fNoGeometriesLimiting= noLimited; // Save # processes limiting step 1067 1072 1068 if( (last > -1) && (noLimited == 1 ) ) 1073 if( (last > -1) && (noLimited == 1 ) ) 1069 { 1074 { 1070 fLimitedStep[ last ] = kUnique; 1075 fLimitedStep[ last ] = kUnique; 1071 } 1076 } 1072 1077 1073 #ifdef G4DEBUG_PATHFINDER 1078 #ifdef G4DEBUG_PATHFINDER 1074 if( fVerboseLevel > 1 ) 1079 if( fVerboseLevel > 1 ) 1075 { 1080 { 1076 PrintLimited(); // --> for tracing 1081 PrintLimited(); // --> for tracing 1077 if( fVerboseLevel > 4 ) << 1082 if( fVerboseLevel > 4 ) { 1078 { << 1079 G4cout << " G4PathFinder::WhichLimited 1083 G4cout << " G4PathFinder::WhichLimited - exiting. " << G4endl; 1080 } 1084 } 1081 } 1085 } 1082 #endif 1086 #endif 1083 } 1087 } 1084 1088 1085 void G4PathFinder::PrintLimited() 1089 void G4PathFinder::PrintLimited() 1086 { 1090 { 1087 // Report results -- for checking 1091 // Report results -- for checking 1088 1092 1089 G4cout << "G4PathFinder::PrintLimited repor 1093 G4cout << "G4PathFinder::PrintLimited reports: " ; 1090 G4cout << " Minimum step (true)= " << fTru 1094 G4cout << " Minimum step (true)= " << fTrueMinStep 1091 << " reported min = " << fMinStep 1095 << " reported min = " << fMinStep 1092 << G4endl; 1096 << G4endl; 1093 if( (fCurrentStepNo <= 2) || (fVerboseLeve 1097 if( (fCurrentStepNo <= 2) || (fVerboseLevel>=2) ) 1094 { 1098 { 1095 G4cout << std::setw(5) << " Step#" << " 1099 G4cout << std::setw(5) << " Step#" << " " 1096 << std::setw(5) << " NavId" << " 1100 << std::setw(5) << " NavId" << " " 1097 << std::setw(12) << " step-size " 1101 << std::setw(12) << " step-size " << " " 1098 << std::setw(12) << " raw-size " 1102 << std::setw(12) << " raw-size " << " " 1099 << std::setw(12) << " pre-safety " 1103 << std::setw(12) << " pre-safety " << " " 1100 << std::setw(15) << " Limited / fl 1104 << std::setw(15) << " Limited / flag" << " " 1101 << std::setw(15) << " World " << 1105 << std::setw(15) << " World " << " " 1102 << G4endl; 1106 << G4endl; 1103 } 1107 } 1104 for ( auto num = 0; num < fNoActiveNavigato << 1108 G4int num; >> 1109 for ( num= 0; num < fNoActiveNavigators; num++ ) 1105 { 1110 { 1106 G4double rawStep = fCurrentStepSize[num]; 1111 G4double rawStep = fCurrentStepSize[num]; 1107 G4double stepLen = fCurrentStepSize[num]; 1112 G4double stepLen = fCurrentStepSize[num]; 1108 if( stepLen > fTrueMinStep ) 1113 if( stepLen > fTrueMinStep ) 1109 { 1114 { 1110 stepLen = fTrueMinStep; // did not 1115 stepLen = fTrueMinStep; // did not limit (went as far as asked) 1111 } 1116 } 1112 G4long oldPrec = G4cout.precision(9); << 1117 G4int oldPrec= G4cout.precision(9); 1113 1118 1114 G4cout << std::setw(5) << fCurrentStepNo 1119 G4cout << std::setw(5) << fCurrentStepNo << " " 1115 << std::setw(5) << num << " " 1120 << std::setw(5) << num << " " 1116 << std::setw(12) << stepLen << " " 1121 << std::setw(12) << stepLen << " " 1117 << std::setw(12) << rawStep << " " 1122 << std::setw(12) << rawStep << " " 1118 << std::setw(12) << fCurrentPreSte 1123 << std::setw(12) << fCurrentPreStepSafety[num] << " " 1119 << std::setw(5) << (fLimitTruth[nu 1124 << std::setw(5) << (fLimitTruth[num] ? "YES" : " NO") << " "; 1120 G4String limitedStr= LimitedString(fLimit 1125 G4String limitedStr= LimitedString(fLimitedStep[num]); 1121 G4cout << " " << std::setw(15) << limited 1126 G4cout << " " << std::setw(15) << limitedStr << " "; 1122 G4cout.precision(oldPrec); 1127 G4cout.precision(oldPrec); 1123 1128 1124 G4Navigator* pNav = GetNavigator( num ); << 1129 G4Navigator *pNav= GetNavigator( num ); 1125 G4String WorldName( "Not-Set" ); 1130 G4String WorldName( "Not-Set" ); 1126 if (pNav != nullptr) << 1131 if (pNav) 1127 { 1132 { 1128 G4VPhysicalVolume *pWorld = pNav->GetW << 1133 G4VPhysicalVolume *pWorld= pNav->GetWorldVolume(); 1129 if( pWorld != nullptr ) << 1134 if( pWorld ) 1130 { 1135 { 1131 WorldName = pWorld->GetName(); << 1136 WorldName = pWorld->GetName(); 1132 } 1137 } 1133 } 1138 } 1134 G4cout << " " << WorldName ; 1139 G4cout << " " << WorldName ; 1135 G4cout << G4endl; 1140 G4cout << G4endl; 1136 } 1141 } 1137 1142 1138 if( fVerboseLevel > 4 ) 1143 if( fVerboseLevel > 4 ) 1139 { 1144 { 1140 G4cout << " G4PathFinder::PrintLimited - 1145 G4cout << " G4PathFinder::PrintLimited - exiting. " << G4endl; 1141 } 1146 } 1142 } 1147 } 1143 1148 1144 G4double 1149 G4double 1145 G4PathFinder::DoNextCurvedStep( const G4Field 1150 G4PathFinder::DoNextCurvedStep( const G4FieldTrack &initialState, 1146 G4doubl << 1151 G4double proposedStepLength, 1147 G4VPhysicalVo 1152 G4VPhysicalVolume* pCurrentPhysicalVolume ) 1148 { 1153 { 1149 const G4double toleratedRelativeError = 1.0 << 1154 const G4double toleratedRelativeError= 1.0e-10; 1150 G4double minStep= kInfinity, newSafety = 0. << 1155 G4double minStep= kInfinity, newSafety=0.0; 1151 G4int numNav; 1156 G4int numNav; 1152 G4FieldTrack fieldTrack = initialState; << 1157 G4FieldTrack fieldTrack= initialState; 1153 G4ThreeVector startPoint = initialState.Get << 1158 G4ThreeVector startPoint= initialState.GetPosition(); 1154 1159 1155 1160 1156 G4EquationOfMotion* equationOfMotion = 1161 G4EquationOfMotion* equationOfMotion = 1157 fpFieldPropagator->GetCurrentEquationOfM << 1162 (fpFieldPropagator->GetChordFinder()->GetIntegrationDriver()->GetStepper()) >> 1163 ->GetEquationOfMotion(); 1158 1164 1159 equationOfMotion->SetChargeMomentumMass( *( 1165 equationOfMotion->SetChargeMomentumMass( *(initialState.GetChargeState()), 1160 in << 1166 initialState.GetMomentum().mag2(), 1161 in 1167 initialState.GetRestMass() ); 1162 1168 1163 #ifdef G4DEBUG_PATHFINDER 1169 #ifdef G4DEBUG_PATHFINDER 1164 G4int prc = G4cout.precision(9); << 1170 G4int prc= G4cout.precision(9); 1165 if( fVerboseLevel > 2 ) 1171 if( fVerboseLevel > 2 ) 1166 { 1172 { 1167 G4cout << " G4PathFinder::DoNextCurvedSte 1173 G4cout << " G4PathFinder::DoNextCurvedStep ****** " << G4endl; 1168 G4cout << " Initial value of field track 1174 G4cout << " Initial value of field track is " << fieldTrack 1169 << " and proposed step= " << propo 1175 << " and proposed step= " << proposedStepLength << G4endl; 1170 } 1176 } 1171 #endif 1177 #endif 1172 1178 1173 fPreStepCenterRenewed = true; // Always upd << 1179 fPreStepCenterRenewed= true; // Always update PreSafety with PreStep point 1174 1180 1175 if( fNoActiveNavigators > 1 ) 1181 if( fNoActiveNavigators > 1 ) 1176 { 1182 { 1177 // Calculate the safety values before ma 1183 // Calculate the safety values before making the step 1178 1184 1179 G4double minSafety= kInfinity, safety; 1185 G4double minSafety= kInfinity, safety; 1180 for( numNav=0; numNav < fNoActiveNavigat 1186 for( numNav=0; numNav < fNoActiveNavigators; ++numNav ) 1181 { 1187 { 1182 safety= fpNavigator[numNav]->ComputeS << 1188 safety= fpNavigator[numNav]->ComputeSafety( startPoint, DBL_MAX, false ); 1183 fPreSafetyValues[numNav] = safety; << 1189 fPreSafetyValues[numNav]= safety; 1184 fCurrentPreStepSafety[numNav] = safet << 1190 fCurrentPreStepSafety[numNav]= safety; 1185 minSafety = std::min( safety, minSafe 1191 minSafety = std::min( safety, minSafety ); 1186 } 1192 } 1187 1193 1188 // Save safety value, related position 1194 // Save safety value, related position 1189 1195 1190 fPreSafetyLocation = startPoint; << 1196 fPreSafetyLocation= startPoint; 1191 fPreSafetyMinValue = minSafety; << 1197 fPreSafetyMinValue= minSafety; 1192 fPreStepLocation = startPoint; << 1198 fPreStepLocation= startPoint; 1193 fMinSafety_PreStepPt = minSafety; << 1199 fMinSafety_PreStepPt= minSafety; 1194 } 1200 } 1195 1201 1196 // Allow Propagator In Field to do the hard 1202 // Allow Propagator In Field to do the hard work, calling G4MultiNavigator 1197 // 1203 // 1198 minStep = fpFieldPropagator->ComputeStep( f << 1204 minStep= fpFieldPropagator->ComputeStep( fieldTrack, 1199 p 1205 proposedStepLength, 1200 n 1206 newSafety, 1201 p << 1207 pCurrentPhysicalVolume ); 1202 f << 1203 1208 1204 // fieldTrack now contains the endpoint inf 1209 // fieldTrack now contains the endpoint information 1205 // 1210 // 1206 fEndState = fieldTrack; << 1211 fEndState= fieldTrack; 1207 fMinStep = minStep; << 1212 fMinStep= minStep; 1208 fTrueMinStep = std::min( minStep, proposedS 1213 fTrueMinStep = std::min( minStep, proposedStepLength ); 1209 1214 1210 if( fNoActiveNavigators == 1 ) << 1215 if( fNoActiveNavigators== 1 ) 1211 { 1216 { 1212 // Update the 'PreSafety' sphere - as an 1217 // Update the 'PreSafety' sphere - as any ComputeStep was called 1213 // (must be done anyway in field) 1218 // (must be done anyway in field) 1214 1219 1215 fPreSafetyValues[0] = newSafety; << 1220 fPreSafetyValues[0]= newSafety; 1216 fPreSafetyLocation = startPoint; << 1221 fPreSafetyLocation= startPoint; 1217 fPreSafetyMinValue = newSafety; << 1222 fPreSafetyMinValue= newSafety; 1218 1223 1219 // Update the current 'PreStep' point's 1224 // Update the current 'PreStep' point's values - mandatory 1220 // 1225 // 1221 fCurrentPreStepSafety[0] = newSafety; << 1226 fCurrentPreStepSafety[0]= newSafety; 1222 fPreStepLocation = startPoint; << 1227 fPreStepLocation= startPoint; 1223 fMinSafety_PreStepPt= newSafety; 1228 fMinSafety_PreStepPt= newSafety; 1224 } 1229 } 1225 1230 1226 #ifdef G4DEBUG_PATHFINDER 1231 #ifdef G4DEBUG_PATHFINDER 1227 if( fVerboseLevel > 2 ) 1232 if( fVerboseLevel > 2 ) 1228 { 1233 { 1229 G4cout << "G4PathFinder::DoNextCurvedStep 1234 G4cout << "G4PathFinder::DoNextCurvedStep : " << G4endl 1230 << " initialState = " << initialSt 1235 << " initialState = " << initialState << G4endl 1231 << " and endState = " << fEndState 1236 << " and endState = " << fEndState << G4endl; 1232 G4cout << "G4PathFinder::DoNextCurvedStep 1237 G4cout << "G4PathFinder::DoNextCurvedStep : " 1233 << " minStep = " << minStep 1238 << " minStep = " << minStep 1234 << " proposedStepLength " << propo 1239 << " proposedStepLength " << proposedStepLength 1235 << " safety = " << newSafety << G4 1240 << " safety = " << newSafety << G4endl; 1236 } 1241 } 1237 #endif 1242 #endif 1238 G4double currentStepSize; // = 0.0; 1243 G4double currentStepSize; // = 0.0; 1239 if( minStep < proposedStepLength ) // if == 1244 if( minStep < proposedStepLength ) // if == , then a boundary found at end ?? 1240 { 1245 { 1241 // Recover the remaining information from 1246 // Recover the remaining information from MultiNavigator 1242 // especially regarding which Navigator l 1247 // especially regarding which Navigator limited the step 1243 1248 1244 G4int noLimited = 0; // No geometries << 1249 G4int noLimited= 0; // No geometries limiting step 1245 for( numNav=0; numNav < fNoActiveNavigato 1250 for( numNav=0; numNav < fNoActiveNavigators; ++numNav ) 1246 { 1251 { 1247 G4double finalStep, lastPreSafety = 0.0 << 1252 G4double finalStep, lastPreSafety=0.0, minStepLast; 1248 ELimited didLimit; 1253 ELimited didLimit; 1249 G4bool limited; 1254 G4bool limited; 1250 1255 1251 finalStep= fpMultiNavigator->ObtainFin 1256 finalStep= fpMultiNavigator->ObtainFinalStep( numNav, lastPreSafety, 1252 1257 minStepLast, didLimit ); 1253 1258 1254 // Calculate the step for this geometry 1259 // Calculate the step for this geometry, using the 1255 // final step (the only one which can d 1260 // final step (the only one which can differ.) 1256 1261 1257 currentStepSize = fTrueMinStep; 1262 currentStepSize = fTrueMinStep; 1258 G4double diffStep = 0.0; << 1263 G4double diffStep= 0.0; 1259 if( (minStepLast != kInfinity) ) 1264 if( (minStepLast != kInfinity) ) 1260 { 1265 { 1261 diffStep = (finalStep-minStepLast); 1266 diffStep = (finalStep-minStepLast); 1262 if ( std::abs(diffStep) <= toleratedR 1267 if ( std::abs(diffStep) <= toleratedRelativeError * finalStep ) 1263 { 1268 { 1264 diffStep = 0.0; 1269 diffStep = 0.0; 1265 } 1270 } 1266 currentStepSize += diffStep; 1271 currentStepSize += diffStep; 1267 } 1272 } 1268 fCurrentStepSize[numNav] = currentStepS 1273 fCurrentStepSize[numNav] = currentStepSize; 1269 1274 1270 // TODO: could refine the way to obtain 1275 // TODO: could refine the way to obtain safeties for > 1 geometries 1271 // - for pre step safety 1276 // - for pre step safety 1272 // notify MultiNavigator about n 1277 // notify MultiNavigator about new set of sub-steps 1273 // allow it to return this value 1278 // allow it to return this value in ObtainFinalStep 1274 // instead of lastPreSafety (or 1279 // instead of lastPreSafety (or as well?) 1275 // - for final step start (availabl 1280 // - for final step start (available) 1276 // get final Step start from Mul 1281 // get final Step start from MultiNavigator 1277 // and corresponding safety valu 1282 // and corresponding safety values 1278 // and/or ALSO calculate ComputeSafety 1283 // and/or ALSO calculate ComputeSafety at endpoint 1279 // endSafety= fpNavigator[numNav]-> 1284 // endSafety= fpNavigator[numNav]->ComputeSafety( endPoint ); 1280 1285 1281 fLimitedStep[numNav] = didLimit; 1286 fLimitedStep[numNav] = didLimit; 1282 fLimitTruth[numNav] = limited = (didLim 1287 fLimitTruth[numNav] = limited = (didLimit != kDoNot ); 1283 if( limited ) { ++noLimited; } << 1288 if( limited ) { noLimited++; } 1284 1289 1285 #ifdef G4DEBUG_PATHFINDER 1290 #ifdef G4DEBUG_PATHFINDER 1286 G4bool StepError = (currentStepSize < 0 << 1291 G4bool StepError= (currentStepSize < 0) 1287 || ( (minStepLast != kI << 1292 || ( (minStepLast != kInfinity) && (diffStep < 0) ) ; 1288 if( StepError || (fVerboseLevel > 2) ) 1293 if( StepError || (fVerboseLevel > 2) ) 1289 { 1294 { 1290 G4String limitedString = LimitedStr << 1295 G4String limitedString= LimitedString( fLimitedStep[numNav] ); 1291 1296 1292 G4cout << " G4PathFinder::ComputeStep 1297 G4cout << " G4PathFinder::ComputeStep. Geometry " << numNav 1293 << " step= " << fCurrentStepS 1298 << " step= " << fCurrentStepSize[numNav] 1294 << " from final-step= " << fin 1299 << " from final-step= " << finalStep 1295 << " fTrueMinStep= " << fTrueM 1300 << " fTrueMinStep= " << fTrueMinStep 1296 << " minStepLast= " << minSte 1301 << " minStepLast= " << minStepLast 1297 << " limited = " << (fLimitTr 1302 << " limited = " << (fLimitTruth[numNav] ? "YES" : " NO") 1298 << " "; 1303 << " "; 1299 G4cout << " status = " << limitedStr 1304 G4cout << " status = " << limitedString << " #= " << didLimit 1300 << G4endl; 1305 << G4endl; 1301 1306 1302 if( StepError ) 1307 if( StepError ) 1303 { 1308 { 1304 std::ostringstream message; 1309 std::ostringstream message; 1305 message << "Incorrect calculation o 1310 message << "Incorrect calculation of step size for one navigator" 1306 << G4endl 1311 << G4endl 1307 << " currentStepSize 1312 << " currentStepSize = " << currentStepSize 1308 << ", diffStep= " << diffSt 1313 << ", diffStep= " << diffStep << G4endl 1309 << "ERROR in computing step 1314 << "ERROR in computing step size for this navigator."; 1310 G4Exception("G4PathFinder::DoNextCu 1315 G4Exception("G4PathFinder::DoNextCurvedStep", 1311 "GeomNav0003", FatalExc 1316 "GeomNav0003", FatalException, message); 1312 } 1317 } 1313 } 1318 } 1314 #endif 1319 #endif 1315 } // for num Navigators 1320 } // for num Navigators 1316 1321 1317 fNoGeometriesLimiting = noLimited; // Sa << 1322 fNoGeometriesLimiting= noLimited; // Save # processes limiting step 1318 } 1323 } 1319 else if ( (minStep == proposedStepLength) 1324 else if ( (minStep == proposedStepLength) 1320 || (minStep == kInfinity) 1325 || (minStep == kInfinity) 1321 || ( std::abs(minStep-proposedSte 1326 || ( std::abs(minStep-proposedStepLength) 1322 < toleratedRelativeError * pro 1327 < toleratedRelativeError * proposedStepLength ) ) 1323 { 1328 { 1324 // In case the step was not limited, use 1329 // In case the step was not limited, use default responses 1325 // --> all Navigators 1330 // --> all Navigators 1326 // Also avoid problems in case of PathFin 1331 // Also avoid problems in case of PathFinder using safety to optimise 1327 // - it is possible that the Navigators 1332 // - it is possible that the Navigators were not called 1328 // if the safety was already satisfact 1333 // if the safety was already satisfactory. 1329 // (In that case calling ObtainFinalSt 1334 // (In that case calling ObtainFinalStep gives invalid results.) 1330 1335 1331 currentStepSize = minStep; << 1336 currentStepSize= minStep; 1332 for( numNav=0; numNav < fNoActiveNavigato 1337 for( numNav=0; numNav < fNoActiveNavigators; ++numNav ) 1333 { 1338 { 1334 fCurrentStepSize[numNav] = minStep; 1339 fCurrentStepSize[numNav] = minStep; 1335 // Safety for endpoint ?? // Can event 1340 // Safety for endpoint ?? // Can eventuall improve it -- see TODO above 1336 fLimitedStep[numNav] = kDoNot; 1341 fLimitedStep[numNav] = kDoNot; 1337 fLimitTruth[numNav] = false; 1342 fLimitTruth[numNav] = false; 1338 } 1343 } 1339 fNoGeometriesLimiting = 0; // Save # pro << 1344 fNoGeometriesLimiting= 0; // Save # processes limiting step 1340 } 1345 } 1341 else // (minStep > proposedStepLength) 1346 else // (minStep > proposedStepLength) and not (minStep == kInfinity) 1342 { 1347 { 1343 std::ostringstream message; 1348 std::ostringstream message; 1344 message << "Incorrect calculation of step 1349 message << "Incorrect calculation of step size for one navigator." << G4endl 1345 << " currentStepSize = " < 1350 << " currentStepSize = " << minStep << " is larger than " 1346 << " proposed StepSize = " << pro 1351 << " proposed StepSize = " << proposedStepLength << "."; 1347 G4Exception("G4PathFinder::DoNextCurvedSt 1352 G4Exception("G4PathFinder::DoNextCurvedStep()", 1348 "GeomNav0003", FatalException 1353 "GeomNav0003", FatalException, message); 1349 } 1354 } 1350 1355 1351 #ifdef G4DEBUG_PATHFINDER 1356 #ifdef G4DEBUG_PATHFINDER 1352 if( fVerboseLevel > 2 ) 1357 if( fVerboseLevel > 2 ) 1353 { 1358 { 1354 G4cout << " Exiting G4PathFinder::DoNextC 1359 G4cout << " Exiting G4PathFinder::DoNextCurvedStep " << G4endl; 1355 PrintLimited(); 1360 PrintLimited(); 1356 } 1361 } 1357 G4cout.precision(prc); 1362 G4cout.precision(prc); 1358 #endif 1363 #endif 1359 1364 1360 return minStep; 1365 return minStep; 1361 } 1366 } 1362 1367 >> 1368 >> 1369 G4bool G4PathFinder::RecheckDistanceToCurrentBoundary( >> 1370 const G4ThreeVector &pGlobalPoint, >> 1371 const G4ThreeVector &pDirection, >> 1372 const G4double aProposedMove, >> 1373 G4double *prDistance, >> 1374 G4double *prNewSafety)const >> 1375 { >> 1376 G4bool retval= true; >> 1377 >> 1378 if( fNoActiveNavigators > 0 ) >> 1379 { >> 1380 // Calculate the safety values before making the step >> 1381 >> 1382 G4double minSafety= kInfinity; >> 1383 G4double minMove= kInfinity; >> 1384 int numNav; >> 1385 for( numNav=0; numNav < fNoActiveNavigators; ++numNav ) >> 1386 { >> 1387 G4double distance, safety; >> 1388 G4bool moveIsOK; >> 1389 moveIsOK= fpNavigator[numNav]->RecheckDistanceToCurrentBoundary( >> 1390 pGlobalPoint, >> 1391 pDirection, >> 1392 aProposedMove, >> 1393 &distance, >> 1394 &safety); >> 1395 minSafety = std::min( safety, minSafety ); >> 1396 minMove = std::min( distance, minMove ); >> 1397 // The first surface encountered will determine it >> 1398 // - even if it is at a negative distance. >> 1399 retval &= moveIsOK; >> 1400 } >> 1401 >> 1402 *prDistance= minMove; >> 1403 if( prNewSafety ) *prNewSafety= minSafety; >> 1404 >> 1405 }else{ >> 1406 retval= false; >> 1407 } >> 1408 >> 1409 return retval; >> 1410 } >> 1411 >> 1412 >> 1413 1363 G4String& G4PathFinder::LimitedString( ELimit 1414 G4String& G4PathFinder::LimitedString( ELimited lim ) 1364 { 1415 { 1365 static G4String StrDoNot("DoNot"), 1416 static G4String StrDoNot("DoNot"), 1366 StrUnique("Unique"), 1417 StrUnique("Unique"), 1367 StrUndefined("Undefined"), 1418 StrUndefined("Undefined"), 1368 StrSharedTransport("SharedT 1419 StrSharedTransport("SharedTransport"), 1369 StrSharedOther("SharedOther 1420 StrSharedOther("SharedOther"); 1370 1421 1371 G4String* limitedStr; 1422 G4String* limitedStr; 1372 switch ( lim ) 1423 switch ( lim ) 1373 { 1424 { 1374 case kDoNot: limitedStr = &StrDoNot; br << 1425 case kDoNot: limitedStr= &StrDoNot; break; 1375 case kUnique: limitedStr = &StrUnique; b 1426 case kUnique: limitedStr = &StrUnique; break; 1376 case kSharedTransport: limitedStr = &St << 1427 case kSharedTransport: limitedStr= &StrSharedTransport; break; 1377 case kSharedOther: limitedStr = &StrShar 1428 case kSharedOther: limitedStr = &StrSharedOther; break; 1378 default: limitedStr = &StrUndefined; bre 1429 default: limitedStr = &StrUndefined; break; 1379 } 1430 } 1380 return *limitedStr; 1431 return *limitedStr; 1381 } 1432 } 1382 1433 1383 void G4PathFinder::PushPostSafetyToPreSafety( 1434 void G4PathFinder::PushPostSafetyToPreSafety() 1384 { 1435 { 1385 fPreSafetyLocation = fSafetyLocation; << 1436 fPreSafetyLocation= fSafetyLocation; 1386 fPreSafetyMinValue = fMinSafety_atSafLocati << 1437 fPreSafetyMinValue= fMinSafety_atSafLocation; 1387 for( auto nav=0; nav < fNoActiveNavigators; << 1438 for( G4int nav=0; nav < fNoActiveNavigators; ++nav ) 1388 { 1439 { 1389 fPreSafetyValues[nav] = fNewSafetyComput << 1440 fPreSafetyValues[nav]= fNewSafetyComputed[nav]; 1390 } 1441 } 1391 } 1442 } 1392 1443