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