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 // class G4PropagatorInField Implementation 26 // class G4PropagatorInField Implementation 27 // 27 // 28 // This class implements an algorithm to trac 28 // This class implements an algorithm to track a particle in a 29 // non-uniform magnetic field. It utilises an 29 // non-uniform magnetic field. It utilises an ODE solver (with 30 // the Runge - Kutta method) to evolve the pa 30 // the Runge - Kutta method) to evolve the particle, and drives it 31 // until the particle has traveled a set dist 31 // until the particle has traveled a set distance or it enters a new 32 // volume. 32 // volume. 33 // 33 // 34 // 14.10.96 John Apostolakis, design and imple 34 // 14.10.96 John Apostolakis, design and implementation 35 // 17.03.97 John Apostolakis, renaming new set 35 // 17.03.97 John Apostolakis, renaming new set functions being added 36 // ------------------------------------------- 36 // --------------------------------------------------------------------------- 37 37 38 #include <iomanip> 38 #include <iomanip> 39 39 40 #include "G4PropagatorInField.hh" 40 #include "G4PropagatorInField.hh" 41 #include "G4ios.hh" 41 #include "G4ios.hh" 42 #include "G4SystemOfUnits.hh" 42 #include "G4SystemOfUnits.hh" 43 #include "G4ThreeVector.hh" 43 #include "G4ThreeVector.hh" 44 #include "G4Material.hh" 44 #include "G4Material.hh" 45 #include "G4VPhysicalVolume.hh" 45 #include "G4VPhysicalVolume.hh" 46 #include "G4Navigator.hh" 46 #include "G4Navigator.hh" 47 #include "G4GeometryTolerance.hh" 47 #include "G4GeometryTolerance.hh" 48 #include "G4VCurvedTrajectoryFilter.hh" 48 #include "G4VCurvedTrajectoryFilter.hh" 49 #include "G4ChordFinder.hh" 49 #include "G4ChordFinder.hh" 50 #include "G4MultiLevelLocator.hh" 50 #include "G4MultiLevelLocator.hh" 51 51 52 << 53 // ------------------------------------------- 52 // --------------------------------------------------------------------------- 54 // Constructors and destructor 53 // Constructors and destructor 55 // 54 // 56 G4PropagatorInField::G4PropagatorInField( G4Na 55 G4PropagatorInField::G4PropagatorInField( G4Navigator* theNavigator, 57 G4Fi 56 G4FieldManager* detectorFieldMgr, 58 G4VI 57 G4VIntersectionLocator* vLocator ) 59 : fDetectorFieldMgr(detectorFieldMgr), 58 : fDetectorFieldMgr(detectorFieldMgr), 60 fNavigator(theNavigator), 59 fNavigator(theNavigator), 61 fCurrentFieldMgr(detectorFieldMgr), 60 fCurrentFieldMgr(detectorFieldMgr), 62 End_PointAndTangent(G4ThreeVector(0.,0.,0. 61 End_PointAndTangent(G4ThreeVector(0.,0.,0.), 63 G4ThreeVector(0.,0.,0. 62 G4ThreeVector(0.,0.,0.),0.0,0.0,0.0,0.0,0.0) 64 { 63 { 65 fEpsilonStep = (fDetectorFieldMgr != nullptr 64 fEpsilonStep = (fDetectorFieldMgr != nullptr) 66 ? fDetectorFieldMgr->GetMaximum 65 ? fDetectorFieldMgr->GetMaximumEpsilonStep() : 1.0e-5; 67 << 66 fLargestAcceptableStep = 1000.0 * meter; 68 67 69 fPreviousSftOrigin = G4ThreeVector(0.,0.,0.) 68 fPreviousSftOrigin = G4ThreeVector(0.,0.,0.); 70 kCarTolerance = G4GeometryTolerance::GetInst 69 kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance(); 71 fZeroStepThreshold = std::max( 1.0e5 * kCarT 70 fZeroStepThreshold = std::max( 1.0e5 * kCarTolerance, 1.0e-1 * micrometer ); 72 71 73 fLargestAcceptableStep = 100.0 * meter; // << 74 fMaxStepSizeMultiplier= 0.1 ; // 0.1 in << 75 fMinBigDistance= 100. * CLHEP::mm; << 76 #ifdef G4DEBUG_FIELD 72 #ifdef G4DEBUG_FIELD 77 G4cout << " PiF: Zero Step Threshold set to 73 G4cout << " PiF: Zero Step Threshold set to " 78 << fZeroStepThreshold / millimeter 74 << fZeroStepThreshold / millimeter 79 << " mm." << G4endl; 75 << " mm." << G4endl; 80 G4cout << " PiF: Value of kCarTolerance = 76 G4cout << " PiF: Value of kCarTolerance = " 81 << kCarTolerance / millimeter 77 << kCarTolerance / millimeter 82 << " mm. " << G4endl; 78 << " mm. " << G4endl; 83 fVerboseLevel = 2; 79 fVerboseLevel = 2; 84 fVerbTracePiF = true; 80 fVerbTracePiF = true; 85 #endif 81 #endif 86 82 87 // Defining Intersection Locator and his par 83 // Defining Intersection Locator and his parameters 88 if ( vLocator == nullptr ) 84 if ( vLocator == nullptr ) 89 { 85 { 90 fIntersectionLocator = new G4MultiLevelLoc 86 fIntersectionLocator = new G4MultiLevelLocator(theNavigator); 91 fAllocatedLocator = true; 87 fAllocatedLocator = true; 92 } 88 } 93 else 89 else 94 { 90 { 95 fIntersectionLocator = vLocator; 91 fIntersectionLocator = vLocator; 96 fAllocatedLocator = false; 92 fAllocatedLocator = false; 97 } 93 } 98 RefreshIntersectionLocator(); // Copy all 94 RefreshIntersectionLocator(); // Copy all relevant parameters 99 } 95 } 100 96 101 // ------------------------------------------- 97 // --------------------------------------------------------------------------- 102 // 98 // 103 G4PropagatorInField::~G4PropagatorInField() 99 G4PropagatorInField::~G4PropagatorInField() 104 { 100 { 105 if(fAllocatedLocator) { delete fIntersecti 101 if(fAllocatedLocator) { delete fIntersectionLocator; } 106 } 102 } 107 103 108 // ------------------------------------------- 104 // --------------------------------------------------------------------------- 109 // Update the IntersectionLocator with current 105 // Update the IntersectionLocator with current parameters 110 // 106 // 111 void G4PropagatorInField::RefreshIntersectionL 107 void G4PropagatorInField::RefreshIntersectionLocator() 112 { 108 { 113 fIntersectionLocator->SetEpsilonStepFor(fEps 109 fIntersectionLocator->SetEpsilonStepFor(fEpsilonStep); 114 fIntersectionLocator->SetDeltaIntersectionFo 110 fIntersectionLocator->SetDeltaIntersectionFor(fCurrentFieldMgr->GetDeltaIntersection()); 115 fIntersectionLocator->SetChordFinderFor(GetC 111 fIntersectionLocator->SetChordFinderFor(GetChordFinder()); 116 fIntersectionLocator->SetSafetyParametersFor 112 fIntersectionLocator->SetSafetyParametersFor( fUseSafetyForOptimisation); 117 } 113 } 118 114 119 // ------------------------------------------- 115 // --------------------------------------------------------------------------- 120 // Compute the next geometric Step 116 // Compute the next geometric Step 121 // 117 // 122 G4double G4PropagatorInField::ComputeStep( 118 G4double G4PropagatorInField::ComputeStep( 123 G4FieldTrack& pFieldTrack 119 G4FieldTrack& pFieldTrack, 124 G4double CurrentProp 120 G4double CurrentProposedStepLength, 125 G4double& currentSafe 121 G4double& currentSafety, // IN/OUT 126 G4VPhysicalVolume* pPhysVol, 122 G4VPhysicalVolume* pPhysVol, 127 G4bool canRelaxDel 123 G4bool canRelaxDeltaChord) 128 { 124 { 129 GetChordFinder()->OnComputeStep(&pFieldTrack << 125 GetChordFinder()->OnComputeStep(); 130 const G4double deltaChord = GetChordFinder() 126 const G4double deltaChord = GetChordFinder()->GetDeltaChord(); 131 127 132 // If CurrentProposedStepLength is too small 128 // If CurrentProposedStepLength is too small for finding Chords 133 // then return with no action (for now - TOD 129 // then return with no action (for now - TODO: some action) 134 // 130 // 135 const char* methodName = "G4PropagatorInFiel 131 const char* methodName = "G4PropagatorInField::ComputeStep"; 136 if (CurrentProposedStepLength<kCarTolerance) 132 if (CurrentProposedStepLength<kCarTolerance) 137 { 133 { 138 return kInfinity; 134 return kInfinity; 139 } 135 } 140 136 141 // Introducing smooth trajectory display (ja 137 // Introducing smooth trajectory display (jacek 01/11/2002) 142 // 138 // 143 if (fpTrajectoryFilter != nullptr) << 139 if (fpTrajectoryFilter) 144 { 140 { 145 fpTrajectoryFilter->CreateNewTrajectorySeg 141 fpTrajectoryFilter->CreateNewTrajectorySegment(); 146 } 142 } 147 143 148 fFirstStepInVolume = fNewTrack ? true : fLas 144 fFirstStepInVolume = fNewTrack ? true : fLastStepInVolume; 149 fLastStepInVolume = false; 145 fLastStepInVolume = false; 150 fNewTrack = false; 146 fNewTrack = false; 151 147 152 if( fVerboseLevel > 2 ) 148 if( fVerboseLevel > 2 ) 153 { 149 { 154 G4cout << methodName << " called" << G4end 150 G4cout << methodName << " called" << G4endl; 155 G4cout << " Starting FT: " << pFieldTrac 151 G4cout << " Starting FT: " << pFieldTrack; 156 G4cout << " Requested length = " << Curr 152 G4cout << " Requested length = " << CurrentProposedStepLength << G4endl; 157 G4cout << " PhysVol = "; 153 G4cout << " PhysVol = "; 158 if( pPhysVol != nullptr ) << 154 if( pPhysVol ) 159 { << 160 G4cout << pPhysVol->GetName() << G4endl 155 G4cout << pPhysVol->GetName() << G4endl; 161 } << 162 else 156 else 163 { << 164 G4cout << " N/A "; 157 G4cout << " N/A "; 165 } << 166 G4cout << G4endl; 158 G4cout << G4endl; 167 } 159 } 168 160 169 // Parameters for adaptive Runge-Kutta integ 161 // Parameters for adaptive Runge-Kutta integration 170 162 171 G4double h_TrialStepSize; // 1st Step 163 G4double h_TrialStepSize; // 1st Step Size 172 G4double TruePathLength = CurrentProposedSte 164 G4double TruePathLength = CurrentProposedStepLength; 173 G4double StepTaken = 0.0; 165 G4double StepTaken = 0.0; 174 G4double s_length_taken, epsilon; 166 G4double s_length_taken, epsilon; 175 G4bool intersects; 167 G4bool intersects; 176 G4bool first_substep = true; 168 G4bool first_substep = true; 177 169 178 G4double NewSafety; 170 G4double NewSafety; 179 fParticleIsLooping = false; 171 fParticleIsLooping = false; 180 172 181 // If not yet done, 173 // If not yet done, 182 // Set the field manager to the local one 174 // Set the field manager to the local one if the volume has one, 183 // or to the global one 175 // or to the global one if not 184 // 176 // 185 if( !fSetFieldMgr ) 177 if( !fSetFieldMgr ) 186 { 178 { 187 fCurrentFieldMgr = FindAndSetFieldManager( 179 fCurrentFieldMgr = FindAndSetFieldManager( pPhysVol ); 188 } 180 } 189 fSetFieldMgr = false; // For next call, the 181 fSetFieldMgr = false; // For next call, the field manager must be set again 190 182 191 G4FieldTrack CurrentState(pFieldTrack); 183 G4FieldTrack CurrentState(pFieldTrack); 192 G4FieldTrack OriginalState = CurrentState; 184 G4FieldTrack OriginalState = CurrentState; 193 185 194 // If the Step length is "infinite", then an 186 // If the Step length is "infinite", then an approximate-maximum Step 195 // length (used to calculate the relative ac 187 // length (used to calculate the relative accuracy) must be guessed 196 // 188 // 197 if( CurrentProposedStepLength >= fLargestAcc 189 if( CurrentProposedStepLength >= fLargestAcceptableStep ) 198 { 190 { 199 G4ThreeVector StartPointA, VelocityUnit; 191 G4ThreeVector StartPointA, VelocityUnit; 200 StartPointA = pFieldTrack.GetPosition(); 192 StartPointA = pFieldTrack.GetPosition(); 201 VelocityUnit = pFieldTrack.GetMomentumDir( 193 VelocityUnit = pFieldTrack.GetMomentumDir(); 202 194 203 G4double trialProposedStep = fMaxStepSizeM << 195 G4double trialProposedStep = 1.e2 * ( 10.0 * cm + 204 fNavigator->GetWorldVolume()->GetLogical 196 fNavigator->GetWorldVolume()->GetLogicalVolume()-> 205 GetSolid()->DistanceToOut(St 197 GetSolid()->DistanceToOut(StartPointA, VelocityUnit) ); 206 CurrentProposedStepLength = std::min( tria 198 CurrentProposedStepLength = std::min( trialProposedStep, 207 fLar 199 fLargestAcceptableStep ); 208 } 200 } 209 epsilon = fCurrentFieldMgr->GetDeltaOneStep( 201 epsilon = fCurrentFieldMgr->GetDeltaOneStep() / CurrentProposedStepLength; 210 G4double epsilonMin= fCurrentFieldMgr->GetMi 202 G4double epsilonMin= fCurrentFieldMgr->GetMinimumEpsilonStep(); 211 G4double epsilonMax= fCurrentFieldMgr->GetMa 203 G4double epsilonMax= fCurrentFieldMgr->GetMaximumEpsilonStep(); 212 if( epsilon < epsilonMin ) { epsilon = epsi 204 if( epsilon < epsilonMin ) { epsilon = epsilonMin; } 213 if( epsilon > epsilonMax ) { epsilon = epsi 205 if( epsilon > epsilonMax ) { epsilon = epsilonMax; } 214 SetEpsilonStep( epsilon ); 206 SetEpsilonStep( epsilon ); 215 207 216 // Values for Intersection Locator has to be 208 // Values for Intersection Locator has to be updated on each call for the 217 // case that CurrentFieldManager has changed 209 // case that CurrentFieldManager has changed from the one of previous step 218 // 210 // 219 RefreshIntersectionLocator(); 211 RefreshIntersectionLocator(); 220 212 221 // Shorten the proposed step in case of earl 213 // Shorten the proposed step in case of earlier problems (zero steps) 222 // 214 // 223 if( fNoZeroStep > fActionThreshold_NoZeroSte 215 if( fNoZeroStep > fActionThreshold_NoZeroSteps ) 224 { 216 { 225 G4double stepTrial; 217 G4double stepTrial; 226 218 227 stepTrial = fFull_CurveLen_of_LastAttempt; 219 stepTrial = fFull_CurveLen_of_LastAttempt; 228 if( (stepTrial <= 0.0) && (fLast_ProposedS 220 if( (stepTrial <= 0.0) && (fLast_ProposedStepLength > 0.0) ) 229 { 221 { 230 stepTrial = fLast_ProposedStepLength; 222 stepTrial = fLast_ProposedStepLength; 231 } 223 } 232 224 233 G4double decreaseFactor = 0.9; // Unused d 225 G4double decreaseFactor = 0.9; // Unused default 234 if( (fNoZeroStep < fSevereActionThreshol 226 if( (fNoZeroStep < fSevereActionThreshold_NoZeroSteps) 235 && (stepTrial > 100.0*fZeroStepThreshol 227 && (stepTrial > 100.0*fZeroStepThreshold) ) 236 { 228 { 237 // Attempt quick convergence 229 // Attempt quick convergence 238 // 230 // 239 decreaseFactor= 0.25; 231 decreaseFactor= 0.25; 240 } 232 } 241 else 233 else 242 { 234 { 243 // We are in significant difficulties, p 235 // We are in significant difficulties, probably at a boundary that 244 // is either geometrically sharp or betw 236 // is either geometrically sharp or between very different materials. 245 // Careful decreases to cope with tolera 237 // Careful decreases to cope with tolerance are required 246 // 238 // 247 if( stepTrial > 100.0*fZeroStepThreshold << 239 if( stepTrial > 100.0*fZeroStepThreshold ) 248 decreaseFactor = 0.35; // Try decr 240 decreaseFactor = 0.35; // Try decreasing slower 249 } else if( stepTrial > 30.0*fZeroStepThr << 241 else if( stepTrial > 30.0*fZeroStepThreshold ) 250 decreaseFactor= 0.5; // Try yet 242 decreaseFactor= 0.5; // Try yet slower decrease 251 } else if( stepTrial > 10.0*fZeroStepThr << 243 else if( stepTrial > 10.0*fZeroStepThreshold ) 252 decreaseFactor= 0.75; // Try even 244 decreaseFactor= 0.75; // Try even slower decreases 253 } else { << 245 else 254 decreaseFactor= 0.9; // Try very 246 decreaseFactor= 0.9; // Try very slow decreases 255 } << 256 } 247 } 257 stepTrial *= decreaseFactor; 248 stepTrial *= decreaseFactor; 258 249 259 #ifdef G4DEBUG_FIELD 250 #ifdef G4DEBUG_FIELD 260 if( fVerboseLevel > 2 251 if( fVerboseLevel > 2 261 || (fNoZeroStep >= fSevereActionThreshol 252 || (fNoZeroStep >= fSevereActionThreshold_NoZeroSteps) ) 262 { 253 { 263 G4cerr << " " << methodName 254 G4cerr << " " << methodName 264 << " Decreasing step after " < 255 << " Decreasing step after " << fNoZeroStep << " zero steps " 265 << " - in volume " << pPhysVol; 256 << " - in volume " << pPhysVol; 266 if( pPhysVol ) 257 if( pPhysVol ) 267 G4cerr << " with name " << pPhysVol 258 G4cerr << " with name " << pPhysVol->GetName(); 268 else 259 else 269 G4cerr << " i.e. *unknown* volume." 260 G4cerr << " i.e. *unknown* volume."; 270 G4cerr << G4endl; 261 G4cerr << G4endl; 271 PrintStepLengthDiagnostic(CurrentPropo 262 PrintStepLengthDiagnostic(CurrentProposedStepLength, decreaseFactor, 272 stepTrial, p 263 stepTrial, pFieldTrack); 273 } 264 } 274 #endif 265 #endif 275 if( stepTrial == 0.0 ) // Change to mak 266 if( stepTrial == 0.0 ) // Change to make it < 0.1 * kCarTolerance ?? 276 { 267 { 277 std::ostringstream message; 268 std::ostringstream message; 278 message << "Particle abandoned due to l 269 message << "Particle abandoned due to lack of progress in field." 279 << G4endl 270 << G4endl 280 << " Properties : " << pFieldT 271 << " Properties : " << pFieldTrack << G4endl 281 << " Attempting a zero step = 272 << " Attempting a zero step = " << stepTrial << G4endl 282 << " while attempting to progr 273 << " while attempting to progress after " << fNoZeroStep 283 << " trial steps. Will abandon 274 << " trial steps. Will abandon step."; 284 G4Exception(methodName, "GeomNav1002", 275 G4Exception(methodName, "GeomNav1002", JustWarning, message); 285 fParticleIsLooping = true; 276 fParticleIsLooping = true; 286 return 0; // = stepTrial; 277 return 0; // = stepTrial; 287 } 278 } 288 if( stepTrial < CurrentProposedStepLength 279 if( stepTrial < CurrentProposedStepLength ) 289 { 280 { 290 CurrentProposedStepLength = stepTrial; 281 CurrentProposedStepLength = stepTrial; 291 } 282 } 292 } 283 } 293 fLast_ProposedStepLength = CurrentProposedSt 284 fLast_ProposedStepLength = CurrentProposedStepLength; 294 285 295 G4int do_loop_count = 0; 286 G4int do_loop_count = 0; 296 do // Loop checking, 07.10.2016, JA 287 do // Loop checking, 07.10.2016, JA 297 { 288 { 298 G4FieldTrack SubStepStartState = CurrentSt 289 G4FieldTrack SubStepStartState = CurrentState; 299 G4ThreeVector SubStartPoint = CurrentState 290 G4ThreeVector SubStartPoint = CurrentState.GetPosition(); 300 291 301 if(!first_substep) 292 if(!first_substep) 302 { 293 { 303 if( fVerboseLevel > 4 ) 294 if( fVerboseLevel > 4 ) 304 { 295 { 305 G4cout << " PiF: Calling Nav/Locate Gl 296 G4cout << " PiF: Calling Nav/Locate Global Point within-Volume " 306 << G4endl; 297 << G4endl; 307 } 298 } 308 fNavigator->LocateGlobalPointWithinVolum 299 fNavigator->LocateGlobalPointWithinVolume( SubStartPoint ); 309 } 300 } 310 301 311 // How far to attempt to move the particle 302 // How far to attempt to move the particle ! 312 // 303 // 313 h_TrialStepSize = CurrentProposedStepLengt 304 h_TrialStepSize = CurrentProposedStepLength - StepTaken; 314 305 315 if (canRelaxDeltaChord && 306 if (canRelaxDeltaChord && 316 fIncreaseChordDistanceThreshold > 0 & 307 fIncreaseChordDistanceThreshold > 0 && 317 do_loop_count > fIncreaseChordDistance 308 do_loop_count > fIncreaseChordDistanceThreshold && 318 do_loop_count % fIncreaseChordDistance 309 do_loop_count % fIncreaseChordDistanceThreshold == 0) 319 { 310 { 320 GetChordFinder()->SetDeltaChord( 311 GetChordFinder()->SetDeltaChord( 321 GetChordFinder()->GetDeltaChord() * 312 GetChordFinder()->GetDeltaChord() * 2.0 322 ); 313 ); 323 } 314 } 324 315 325 // Integrate as far as "chord miss" rule a 316 // Integrate as far as "chord miss" rule allows. 326 // 317 // 327 s_length_taken = GetChordFinder()->Advance 318 s_length_taken = GetChordFinder()->AdvanceChordLimited( 328 CurrentState, 319 CurrentState, // Position & velocity 329 h_TrialStepSize, 320 h_TrialStepSize, 330 fEpsilonStep, 321 fEpsilonStep, 331 fPreviousSftOrigi 322 fPreviousSftOrigin, 332 fPreviousSafety ) 323 fPreviousSafety ); 333 // CurrentState is now updated with the 324 // CurrentState is now updated with the final position and velocity 334 325 335 fFull_CurveLen_of_LastAttempt = s_length_t 326 fFull_CurveLen_of_LastAttempt = s_length_taken; 336 327 337 G4ThreeVector EndPointB = CurrentState.Get 328 G4ThreeVector EndPointB = CurrentState.GetPosition(); 338 G4ThreeVector InterSectionPointE; 329 G4ThreeVector InterSectionPointE; 339 G4double LinearStepLength; 330 G4double LinearStepLength; 340 331 341 // Intersect chord AB with geometry 332 // Intersect chord AB with geometry 342 // 333 // 343 intersects= IntersectChord( SubStartPoint, 334 intersects= IntersectChord( SubStartPoint, EndPointB, 344 NewSafety, Lin 335 NewSafety, LinearStepLength, 345 InterSectionPo 336 InterSectionPointE ); 346 // E <- Intersection Point of chord AB a 337 // E <- Intersection Point of chord AB and either volume A's surface 347 // or a 338 // or a daughter volume's surface .. 348 339 349 if( first_substep ) 340 if( first_substep ) 350 { 341 { 351 currentSafety = NewSafety; 342 currentSafety = NewSafety; 352 } // Updating safety in other steps is pot 343 } // Updating safety in other steps is potential future extention 353 344 354 if( intersects ) 345 if( intersects ) 355 { 346 { 356 G4FieldTrack IntersectPointVelct_G(Curr 347 G4FieldTrack IntersectPointVelct_G(CurrentState); // FT-Def-Construct 357 348 358 // Find the intersection point of AB tr 349 // Find the intersection point of AB true path with the surface 359 // of vol(A), if it exists. Start wit 350 // of vol(A), if it exists. Start with point E as first "estimate". 360 G4bool recalculatedEndPt = false; 351 G4bool recalculatedEndPt = false; 361 352 362 G4bool found_intersection = fIntersecti 353 G4bool found_intersection = fIntersectionLocator-> 363 EstimateIntersectionPoint( SubStepSta 354 EstimateIntersectionPoint( SubStepStartState, CurrentState, 364 InterSecti 355 InterSectionPointE, IntersectPointVelct_G, 365 recalculat 356 recalculatedEndPt, fPreviousSafety, 366 fPreviousS 357 fPreviousSftOrigin); 367 intersects = found_intersection; 358 intersects = found_intersection; 368 if( found_intersection ) 359 if( found_intersection ) 369 { 360 { 370 End_PointAndTangent= IntersectPointV 361 End_PointAndTangent= IntersectPointVelct_G; // G is our EndPoint ... 371 StepTaken = TruePathLength = Interse 362 StepTaken = TruePathLength = IntersectPointVelct_G.GetCurveLength() 372 - Origina 363 - OriginalState.GetCurveLength(); 373 } 364 } 374 else 365 else 375 { 366 { 376 // Either "minor" chords do not inte 367 // Either "minor" chords do not intersect 377 // or else stopped (due to too many 368 // or else stopped (due to too many steps) 378 // 369 // 379 if( recalculatedEndPt ) 370 if( recalculatedEndPt ) 380 { 371 { 381 G4double endAchieved = IntersectP 372 G4double endAchieved = IntersectPointVelct_G.GetCurveLength(); 382 G4double endExpected = CurrentSta 373 G4double endExpected = CurrentState.GetCurveLength(); 383 374 384 // Detect failure - due to too ma 375 // Detect failure - due to too many steps 385 G4bool shortEnd = endAchieved 376 G4bool shortEnd = endAchieved 386 < (endExpected*(1 377 < (endExpected*(1.0-CLHEP::perMillion)); 387 378 388 G4double stepAchieved = endAchiev 379 G4double stepAchieved = endAchieved 389 - SubStepSt 380 - SubStepStartState.GetCurveLength(); 390 381 391 // Update remaining state - must 382 // Update remaining state - must work for 'full' step or 392 // abandonned intersection 383 // abandonned intersection 393 // 384 // 394 CurrentState = IntersectPointVelc 385 CurrentState = IntersectPointVelct_G; 395 s_length_taken = stepAchieved; 386 s_length_taken = stepAchieved; 396 if( shortEnd ) 387 if( shortEnd ) 397 { 388 { 398 fParticleIsLooping = true; 389 fParticleIsLooping = true; 399 } 390 } 400 } 391 } 401 } 392 } 402 } 393 } 403 if( !intersects ) 394 if( !intersects ) 404 { 395 { 405 StepTaken += s_length_taken; 396 StepTaken += s_length_taken; 406 397 407 if (fpTrajectoryFilter != nullptr) // Fo << 398 if (fpTrajectoryFilter) // For smooth trajectory display (jacek 1/11/2002) 408 { 399 { 409 fpTrajectoryFilter->TakeIntermediatePo 400 fpTrajectoryFilter->TakeIntermediatePoint(CurrentState.GetPosition()); 410 } 401 } 411 } 402 } 412 first_substep = false; 403 first_substep = false; 413 404 414 #ifdef G4DEBUG_FIELD 405 #ifdef G4DEBUG_FIELD 415 if( fNoZeroStep > fActionThreshold_NoZeroS 406 if( fNoZeroStep > fActionThreshold_NoZeroSteps ) 416 { 407 { 417 if( fNoZeroStep > fSevereActionThreshold 408 if( fNoZeroStep > fSevereActionThreshold_NoZeroSteps ) 418 G4cout << " Above 'Severe Action' thre 409 G4cout << " Above 'Severe Action' threshold -- for Zero steps. "; 419 else 410 else 420 G4cout << " Above 'action' threshold - 411 G4cout << " Above 'action' threshold -- for Zero steps. "; 421 G4cout << " Number of zero steps = " << 412 G4cout << " Number of zero steps = " << fNoZeroStep << G4endl; 422 printStatus( SubStepStartState, // or O 413 printStatus( SubStepStartState, // or OriginalState, 423 CurrentState, CurrentPropos 414 CurrentState, CurrentProposedStepLength, 424 NewSafety, do_loop_count, p 415 NewSafety, do_loop_count, pPhysVol ); 425 } 416 } 426 if( (fVerboseLevel > 1) && (do_loop_count 417 if( (fVerboseLevel > 1) && (do_loop_count > fMax_loop_count-10 )) 427 { 418 { 428 if( do_loop_count == fMax_loop_count-9 ) 419 if( do_loop_count == fMax_loop_count-9 ) 429 { 420 { 430 G4cout << " G4PropagatorInField::Compu 421 G4cout << " G4PropagatorInField::ComputeStep(): " << G4endl 431 << " Difficult track - taking 422 << " Difficult track - taking many sub steps." << G4endl; 432 printStatus( SubStepStartState, SubSte 423 printStatus( SubStepStartState, SubStepStartState, CurrentProposedStepLength, 433 NewSafety, 0, pPhysVol ); 424 NewSafety, 0, pPhysVol ); 434 } 425 } 435 printStatus( SubStepStartState, CurrentS 426 printStatus( SubStepStartState, CurrentState, CurrentProposedStepLength, 436 NewSafety, do_loop_count, p 427 NewSafety, do_loop_count, pPhysVol ); 437 } 428 } 438 #endif 429 #endif 439 430 440 ++do_loop_count; 431 ++do_loop_count; 441 432 442 } while( (!intersects ) 433 } while( (!intersects ) 443 && (!fParticleIsLooping) 434 && (!fParticleIsLooping) 444 && (StepTaken + kCarTolerance < Curren 435 && (StepTaken + kCarTolerance < CurrentProposedStepLength) 445 && ( do_loop_count < fMax_loop_count ) 436 && ( do_loop_count < fMax_loop_count ) ); 446 437 447 if( do_loop_count >= fMax_loop_count 438 if( do_loop_count >= fMax_loop_count 448 && (StepTaken + kCarTolerance < CurrentPro 439 && (StepTaken + kCarTolerance < CurrentProposedStepLength) ) 449 { 440 { 450 fParticleIsLooping = true; 441 fParticleIsLooping = true; 451 } 442 } 452 if ( ( fParticleIsLooping ) && (fVerboseLeve 443 if ( ( fParticleIsLooping ) && (fVerboseLevel > 0) ) 453 { 444 { 454 ReportLoopingParticle( do_loop_count, Step 445 ReportLoopingParticle( do_loop_count, StepTaken, 455 CurrentProposedStep 446 CurrentProposedStepLength, methodName, 456 CurrentState.GetMom 447 CurrentState.GetMomentum(), pPhysVol ); 457 } 448 } 458 449 459 if( !intersects ) 450 if( !intersects ) 460 { 451 { 461 // Chord AB or "minor chords" do not inter 452 // Chord AB or "minor chords" do not intersect 462 // B is the endpoint Step of the current S 453 // B is the endpoint Step of the current Step. 463 // 454 // 464 End_PointAndTangent = CurrentState; 455 End_PointAndTangent = CurrentState; 465 TruePathLength = StepTaken; // Original 456 TruePathLength = StepTaken; // Original code 466 457 467 // Tried the following to avoid potential 458 // Tried the following to avoid potential issue with round-off error 468 // - but has issues... Suppressing this ch 459 // - but has issues... Suppressing this change JA 2015/05/02 469 // TruePathLength = CurrentProposedStepLen 460 // TruePathLength = CurrentProposedStepLength; 470 } 461 } 471 fLastStepInVolume = intersects; 462 fLastStepInVolume = intersects; 472 463 473 // Set pFieldTrack to the return value 464 // Set pFieldTrack to the return value 474 // 465 // 475 pFieldTrack = End_PointAndTangent; 466 pFieldTrack = End_PointAndTangent; 476 467 477 #ifdef G4VERBOSE 468 #ifdef G4VERBOSE 478 // Check that "s" is correct 469 // Check that "s" is correct 479 // 470 // 480 if( std::fabs(OriginalState.GetCurveLength() 471 if( std::fabs(OriginalState.GetCurveLength() + TruePathLength 481 - End_PointAndTangent.GetCurveLength()) 472 - End_PointAndTangent.GetCurveLength()) > 3.e-4 * TruePathLength ) 482 { 473 { 483 std::ostringstream message; 474 std::ostringstream message; 484 message << "Curve length mis-match between 475 message << "Curve length mis-match between original state " 485 << "and proposed endpoint of propa 476 << "and proposed endpoint of propagation." << G4endl 486 << " The curve length of the endp 477 << " The curve length of the endpoint should be: " 487 << OriginalState.GetCurveLength() 478 << OriginalState.GetCurveLength() + TruePathLength << G4endl 488 << " and it is instead: " 479 << " and it is instead: " 489 << End_PointAndTangent.GetCurveLen 480 << End_PointAndTangent.GetCurveLength() << "." << G4endl 490 << " A difference of: " 481 << " A difference of: " 491 << OriginalState.GetCurveLength() 482 << OriginalState.GetCurveLength() + TruePathLength 492 - End_PointAndTangent.GetCurveL 483 - End_PointAndTangent.GetCurveLength() << G4endl 493 << " Original state = " << Origin 484 << " Original state = " << OriginalState << G4endl 494 << " Proposed state = " << End_Po 485 << " Proposed state = " << End_PointAndTangent; 495 G4Exception(methodName, "GeomNav0003", Fat 486 G4Exception(methodName, "GeomNav0003", FatalException, message); 496 } 487 } 497 #endif 488 #endif 498 489 499 if( TruePathLength+kCarTolerance >= CurrentP 490 if( TruePathLength+kCarTolerance >= CurrentProposedStepLength ) 500 { 491 { 501 fNoZeroStep = 0; 492 fNoZeroStep = 0; 502 } 493 } 503 else 494 else 504 { 495 { 505 // In particular anomalous cases, we can 496 // In particular anomalous cases, we can get repeated zero steps 506 // We identify these cases and take corre 497 // We identify these cases and take corrective action when they occur. 507 // 498 // 508 if( TruePathLength < std::max( fZeroStepT 499 if( TruePathLength < std::max( fZeroStepThreshold, 0.5*kCarTolerance ) ) 509 { 500 { 510 ++fNoZeroStep; 501 ++fNoZeroStep; 511 } 502 } 512 else 503 else 513 { 504 { 514 fNoZeroStep = 0; 505 fNoZeroStep = 0; 515 } 506 } 516 } 507 } 517 if( fNoZeroStep > fAbandonThreshold_NoZeroSt 508 if( fNoZeroStep > fAbandonThreshold_NoZeroSteps ) 518 { 509 { 519 fParticleIsLooping = true; 510 fParticleIsLooping = true; 520 ReportStuckParticle( fNoZeroStep, Current 511 ReportStuckParticle( fNoZeroStep, CurrentProposedStepLength, 521 fFull_CurveLen_of_La 512 fFull_CurveLen_of_LastAttempt, pPhysVol ); 522 fNoZeroStep = 0; 513 fNoZeroStep = 0; 523 } 514 } 524 515 525 GetChordFinder()->SetDeltaChord(deltaChord); 516 GetChordFinder()->SetDeltaChord(deltaChord); 526 return TruePathLength; 517 return TruePathLength; 527 } 518 } 528 519 529 // ------------------------------------------- 520 // --------------------------------------------------------------------------- 530 // Dumps status of propagator 521 // Dumps status of propagator 531 // 522 // 532 void 523 void 533 G4PropagatorInField::printStatus( const G4Fiel 524 G4PropagatorInField::printStatus( const G4FieldTrack& StartFT, 534 const G4Fiel 525 const G4FieldTrack& CurrentFT, 535 G4doub 526 G4double requestStep, 536 G4doub 527 G4double safety, 537 G4int 528 G4int stepNo, 538 G4VPhy 529 G4VPhysicalVolume* startVolume) 539 { 530 { 540 const G4int verboseLevel = fVerboseLevel; 531 const G4int verboseLevel = fVerboseLevel; 541 const G4ThreeVector StartPosition = St 532 const G4ThreeVector StartPosition = StartFT.GetPosition(); 542 const G4ThreeVector StartUnitVelocity = St 533 const G4ThreeVector StartUnitVelocity = StartFT.GetMomentumDir(); 543 const G4ThreeVector CurrentPosition = Cu 534 const G4ThreeVector CurrentPosition = CurrentFT.GetPosition(); 544 const G4ThreeVector CurrentUnitVelocity = Cu 535 const G4ThreeVector CurrentUnitVelocity = CurrentFT.GetMomentumDir(); 545 536 546 G4double step_len = CurrentFT.GetCurveLength 537 G4double step_len = CurrentFT.GetCurveLength() - StartFT.GetCurveLength(); 547 538 548 G4long oldprec; // cout/cerr precision set 539 G4long oldprec; // cout/cerr precision settings 549 540 550 if( ((stepNo == 0) && (verboseLevel <3)) || 541 if( ((stepNo == 0) && (verboseLevel <3)) || (verboseLevel >= 3) ) 551 { 542 { 552 oldprec = G4cout.precision(4); 543 oldprec = G4cout.precision(4); 553 G4cout << std::setw( 5) << "Step#" 544 G4cout << std::setw( 5) << "Step#" 554 << std::setw(10) << " s " << " " 545 << std::setw(10) << " s " << " " 555 << std::setw(10) << "X(mm)" << " " 546 << std::setw(10) << "X(mm)" << " " 556 << std::setw(10) << "Y(mm)" << " " 547 << std::setw(10) << "Y(mm)" << " " 557 << std::setw(10) << "Z(mm)" << " " 548 << std::setw(10) << "Z(mm)" << " " 558 << std::setw( 7) << " N_x " << " " 549 << std::setw( 7) << " N_x " << " " 559 << std::setw( 7) << " N_y " << " " 550 << std::setw( 7) << " N_y " << " " 560 << std::setw( 7) << " N_z " << " " 551 << std::setw( 7) << " N_z " << " " ; 561 G4cout << std::setw( 7) << " Delta|N|" << 552 G4cout << std::setw( 7) << " Delta|N|" << " " 562 << std::setw( 9) << "StepLen" << " 553 << std::setw( 9) << "StepLen" << " " 563 << std::setw(12) << "StartSafety" < 554 << std::setw(12) << "StartSafety" << " " 564 << std::setw( 9) << "PhsStep" << " 555 << std::setw( 9) << "PhsStep" << " "; 565 if( startVolume != nullptr ) 556 if( startVolume != nullptr ) 566 { G4cout << std::setw(18) << "NextVolume 557 { G4cout << std::setw(18) << "NextVolume" << " "; } 567 G4cout.precision(oldprec); 558 G4cout.precision(oldprec); 568 G4cout << G4endl; 559 G4cout << G4endl; 569 } 560 } 570 if((stepNo == 0) && (verboseLevel <=3)) 561 if((stepNo == 0) && (verboseLevel <=3)) 571 { 562 { 572 // Recurse to print the start values 563 // Recurse to print the start values 573 // 564 // 574 printStatus( StartFT, StartFT, -1.0, safet 565 printStatus( StartFT, StartFT, -1.0, safety, -1, startVolume); 575 } 566 } 576 if( verboseLevel <= 3 ) 567 if( verboseLevel <= 3 ) 577 { 568 { 578 if( stepNo >= 0) 569 if( stepNo >= 0) 579 { G4cout << std::setw( 4) << stepNo << " 570 { G4cout << std::setw( 4) << stepNo << " "; } 580 else 571 else 581 { G4cout << std::setw( 5) << "Start" ; } 572 { G4cout << std::setw( 5) << "Start" ; } 582 oldprec = G4cout.precision(8); 573 oldprec = G4cout.precision(8); 583 G4cout << std::setw(10) << CurrentFT.GetCu 574 G4cout << std::setw(10) << CurrentFT.GetCurveLength() << " "; 584 G4cout.precision(8); 575 G4cout.precision(8); 585 G4cout << std::setw(10) << CurrentPosition 576 G4cout << std::setw(10) << CurrentPosition.x() << " " 586 << std::setw(10) << CurrentPosition 577 << std::setw(10) << CurrentPosition.y() << " " 587 << std::setw(10) << CurrentPosition 578 << std::setw(10) << CurrentPosition.z() << " "; 588 G4cout.precision(4); 579 G4cout.precision(4); 589 G4cout << std::setw( 7) << CurrentUnitVelo 580 G4cout << std::setw( 7) << CurrentUnitVelocity.x() << " " 590 << std::setw( 7) << CurrentUnitVelo 581 << std::setw( 7) << CurrentUnitVelocity.y() << " " 591 << std::setw( 7) << CurrentUnitVelo 582 << std::setw( 7) << CurrentUnitVelocity.z() << " "; 592 G4cout.precision(3); 583 G4cout.precision(3); 593 G4cout << std::setw( 7) 584 G4cout << std::setw( 7) 594 << CurrentFT.GetMomentum().mag()-St 585 << CurrentFT.GetMomentum().mag()-StartFT.GetMomentum().mag() << " "; 595 G4cout << std::setw( 9) << step_len << " " 586 G4cout << std::setw( 9) << step_len << " "; 596 G4cout << std::setw(12) << safety << " "; 587 G4cout << std::setw(12) << safety << " "; 597 if( requestStep != -1.0 ) 588 if( requestStep != -1.0 ) 598 { G4cout << std::setw( 9) << requestStep 589 { G4cout << std::setw( 9) << requestStep << " "; } 599 else 590 else 600 { G4cout << std::setw( 9) << "Init/NotKn 591 { G4cout << std::setw( 9) << "Init/NotKnown" << " "; } 601 if( startVolume != nullptr) << 592 if( startVolume != 0) 602 { G4cout << std::setw(12) << startVolume 593 { G4cout << std::setw(12) << startVolume->GetName() << " "; } 603 G4cout.precision(oldprec); 594 G4cout.precision(oldprec); 604 G4cout << G4endl; 595 G4cout << G4endl; 605 } 596 } 606 else // if( verboseLevel > 3 ) 597 else // if( verboseLevel > 3 ) 607 { 598 { 608 // Multi-line output 599 // Multi-line output 609 600 610 G4cout << "Step taken was " << step_len 601 G4cout << "Step taken was " << step_len 611 << " out of PhysicalStep = " << re 602 << " out of PhysicalStep = " << requestStep << G4endl; 612 G4cout << "Final safety is: " << safety << 603 G4cout << "Final safety is: " << safety << G4endl; 613 G4cout << "Chord length = " << (CurrentPos 604 G4cout << "Chord length = " << (CurrentPosition-StartPosition).mag() 614 << G4endl; 605 << G4endl; 615 G4cout << G4endl; 606 G4cout << G4endl; 616 } 607 } 617 } 608 } 618 609 619 // ------------------------------------------- 610 // --------------------------------------------------------------------------- 620 // Prints Step diagnostics 611 // Prints Step diagnostics 621 // 612 // 622 void 613 void 623 G4PropagatorInField::PrintStepLengthDiagnostic 614 G4PropagatorInField::PrintStepLengthDiagnostic( 624 G4double CurrentProp 615 G4double CurrentProposedStepLength, 625 G4double decreaseFac 616 G4double decreaseFactor, 626 G4double stepTrial, 617 G4double stepTrial, 627 const G4FieldTrack& ) 618 const G4FieldTrack& ) 628 { 619 { 629 G4long iprec= G4cout.precision(8); 620 G4long iprec= G4cout.precision(8); 630 G4cout << " " << std::setw(12) << " PiF: NoZ 621 G4cout << " " << std::setw(12) << " PiF: NoZeroStep " 631 << " " << std::setw(20) << " CurrentP 622 << " " << std::setw(20) << " CurrentProposed len " 632 << " " << std::setw(18) << " Full_cur 623 << " " << std::setw(18) << " Full_curvelen_last" 633 << " " << std::setw(18) << " last pro 624 << " " << std::setw(18) << " last proposed len " 634 << " " << std::setw(18) << " decrease 625 << " " << std::setw(18) << " decrease factor " 635 << " " << std::setw(15) << " step tri 626 << " " << std::setw(15) << " step trial " 636 << G4endl; 627 << G4endl; 637 628 638 G4cout << " " << std::setw(10) << fNoZeroSte 629 G4cout << " " << std::setw(10) << fNoZeroStep << " " 639 << " " << std::setw(20) << CurrentPro 630 << " " << std::setw(20) << CurrentProposedStepLength 640 << " " << std::setw(18) << fFull_Curv 631 << " " << std::setw(18) << fFull_CurveLen_of_LastAttempt 641 << " " << std::setw(18) << fLast_Prop 632 << " " << std::setw(18) << fLast_ProposedStepLength 642 << " " << std::setw(18) << decreaseFa 633 << " " << std::setw(18) << decreaseFactor 643 << " " << std::setw(15) << stepTrial 634 << " " << std::setw(15) << stepTrial 644 << G4endl; 635 << G4endl; 645 G4cout.precision( iprec ); 636 G4cout.precision( iprec ); 646 } 637 } 647 638 648 // Access the points which have passed through 639 // Access the points which have passed through the filter. The 649 // points are stored as ThreeVectors for the i 640 // points are stored as ThreeVectors for the initial impelmentation 650 // only (jacek 30/10/2002) 641 // only (jacek 30/10/2002) 651 // Responsibility for deleting the points lies 642 // Responsibility for deleting the points lies with 652 // SmoothTrajectoryPoint, which is the points' 643 // SmoothTrajectoryPoint, which is the points' final 653 // destination. The points pointer is set to N 644 // destination. The points pointer is set to NULL, to ensure that 654 // the points are not re-used in subsequent st 645 // the points are not re-used in subsequent steps, therefore THIS 655 // METHOD MUST BE CALLED EXACTLY ONCE PER STEP 646 // METHOD MUST BE CALLED EXACTLY ONCE PER STEP. (jacek 08/11/2002) 656 647 657 std::vector<G4ThreeVector>* 648 std::vector<G4ThreeVector>* 658 G4PropagatorInField::GimmeTrajectoryVectorAndF 649 G4PropagatorInField::GimmeTrajectoryVectorAndForgetIt() const 659 { 650 { 660 // NB, GimmeThePointsAndForgetThem really fo 651 // NB, GimmeThePointsAndForgetThem really forgets them, so it can 661 // only be called (exactly) once for each st 652 // only be called (exactly) once for each step. 662 653 663 if (fpTrajectoryFilter != nullptr) 654 if (fpTrajectoryFilter != nullptr) 664 { 655 { 665 return fpTrajectoryFilter->GimmeThePointsA 656 return fpTrajectoryFilter->GimmeThePointsAndForgetThem(); 666 } 657 } 667 return nullptr; << 658 else >> 659 { >> 660 return nullptr; >> 661 } 668 } 662 } 669 663 670 // ------------------------------------------- 664 // --------------------------------------------------------------------------- 671 // 665 // 672 void 666 void 673 G4PropagatorInField::SetTrajectoryFilter(G4VCu 667 G4PropagatorInField::SetTrajectoryFilter(G4VCurvedTrajectoryFilter* filter) 674 { 668 { 675 fpTrajectoryFilter = filter; 669 fpTrajectoryFilter = filter; 676 } 670 } 677 671 678 // ------------------------------------------- 672 // --------------------------------------------------------------------------- 679 // 673 // 680 void G4PropagatorInField::ClearPropagatorState 674 void G4PropagatorInField::ClearPropagatorState() 681 { 675 { 682 // Goal: Clear all memory of previous steps, 676 // Goal: Clear all memory of previous steps, cached information 683 677 684 fParticleIsLooping = false; 678 fParticleIsLooping = false; 685 fNoZeroStep = 0; 679 fNoZeroStep = 0; 686 680 687 fSetFieldMgr = false; // Has field-manager 681 fSetFieldMgr = false; // Has field-manager been set for the current step? 688 fEpsilonStep= 1.0e-5; // Relative accuracy 682 fEpsilonStep= 1.0e-5; // Relative accuracy of current Step 689 683 690 End_PointAndTangent= G4FieldTrack( G4ThreeVe 684 End_PointAndTangent= G4FieldTrack( G4ThreeVector(0.,0.,0.), 691 G4ThreeVe 685 G4ThreeVector(0.,0.,0.), 692 0.0,0.0,0 686 0.0,0.0,0.0,0.0,0.0); 693 fFull_CurveLen_of_LastAttempt = -1; 687 fFull_CurveLen_of_LastAttempt = -1; 694 fLast_ProposedStepLength = -1; 688 fLast_ProposedStepLength = -1; 695 689 696 fPreviousSftOrigin= G4ThreeVector(0.,0.,0.); 690 fPreviousSftOrigin= G4ThreeVector(0.,0.,0.); 697 fPreviousSafety= 0.0; 691 fPreviousSafety= 0.0; 698 692 699 fNewTrack = true; 693 fNewTrack = true; 700 } 694 } 701 695 702 // ------------------------------------------- 696 // --------------------------------------------------------------------------- 703 // 697 // 704 G4FieldManager* G4PropagatorInField:: 698 G4FieldManager* G4PropagatorInField:: 705 FindAndSetFieldManager( G4VPhysicalVolume* pCu 699 FindAndSetFieldManager( G4VPhysicalVolume* pCurrentPhysicalVolume ) 706 { 700 { 707 G4FieldManager* currentFieldMgr; 701 G4FieldManager* currentFieldMgr; 708 702 709 currentFieldMgr = fDetectorFieldMgr; 703 currentFieldMgr = fDetectorFieldMgr; 710 if( pCurrentPhysicalVolume != nullptr ) 704 if( pCurrentPhysicalVolume != nullptr ) 711 { 705 { 712 G4FieldManager *pRegionFieldMgr = nullptr 706 G4FieldManager *pRegionFieldMgr = nullptr, *localFieldMgr = nullptr; 713 G4LogicalVolume* pLogicalVol = pCurrentPh 707 G4LogicalVolume* pLogicalVol = pCurrentPhysicalVolume->GetLogicalVolume(); 714 708 715 if( pLogicalVol != nullptr ) 709 if( pLogicalVol != nullptr ) 716 { 710 { 717 // Value for Region, if any, overrides 711 // Value for Region, if any, overrides 718 // 712 // 719 G4Region* pRegion = pLogicalVol->GetR 713 G4Region* pRegion = pLogicalVol->GetRegion(); 720 if( pRegion != nullptr ) 714 if( pRegion != nullptr ) 721 { 715 { 722 pRegionFieldMgr = pRegion->GetField 716 pRegionFieldMgr = pRegion->GetFieldManager(); 723 if( pRegionFieldMgr != nullptr ) 717 if( pRegionFieldMgr != nullptr ) 724 { 718 { 725 currentFieldMgr= pRegionFieldMgr 719 currentFieldMgr= pRegionFieldMgr; 726 } 720 } 727 } 721 } 728 722 729 // 'Local' Value from logical volume, 723 // 'Local' Value from logical volume, if any, overrides 730 // 724 // 731 localFieldMgr = pLogicalVol->GetFieldM 725 localFieldMgr = pLogicalVol->GetFieldManager(); 732 if ( localFieldMgr != nullptr ) 726 if ( localFieldMgr != nullptr ) 733 { 727 { 734 currentFieldMgr = localFieldMgr; 728 currentFieldMgr = localFieldMgr; 735 } 729 } 736 } 730 } 737 } 731 } 738 fCurrentFieldMgr = currentFieldMgr; 732 fCurrentFieldMgr = currentFieldMgr; 739 733 740 // Flag that field manager has been set 734 // Flag that field manager has been set 741 // 735 // 742 fSetFieldMgr = true; 736 fSetFieldMgr = true; 743 737 744 return currentFieldMgr; 738 return currentFieldMgr; 745 } 739 } 746 740 747 // ------------------------------------------- 741 // --------------------------------------------------------------------------- 748 // 742 // 749 G4int G4PropagatorInField::SetVerboseLevel( G4 743 G4int G4PropagatorInField::SetVerboseLevel( G4int level ) 750 { 744 { 751 G4int oldval = fVerboseLevel; 745 G4int oldval = fVerboseLevel; 752 fVerboseLevel = level; 746 fVerboseLevel = level; 753 747 754 // Forward the verbose level 'reduced' to Ch 748 // Forward the verbose level 'reduced' to ChordFinder, 755 // MagIntegratorDriver ... ? 749 // MagIntegratorDriver ... ? 756 // 750 // 757 auto integrDriver = GetChordFinder()->GetInt 751 auto integrDriver = GetChordFinder()->GetIntegrationDriver(); 758 integrDriver->SetVerboseLevel( fVerboseLevel 752 integrDriver->SetVerboseLevel( fVerboseLevel - 2 ); 759 G4cout << "Set Driver verbosity to " << fVer 753 G4cout << "Set Driver verbosity to " << fVerboseLevel - 2 << G4endl; 760 754 761 return oldval; 755 return oldval; 762 } 756 } 763 757 764 // ------------------------------------------- 758 // --------------------------------------------------------------------------- 765 // 759 // 766 void G4PropagatorInField::ReportLoopingParticl 760 void G4PropagatorInField::ReportLoopingParticle( G4int count, 767 761 G4double StepTaken, 768 762 G4double StepRequested, 769 763 const char* methodName, 770 << 764 G4ThreeVector momentumVec, 771 765 G4VPhysicalVolume* pPhysVol ) 772 { 766 { 773 std::ostringstream message; 767 std::ostringstream message; 774 G4double fraction = StepTaken / StepRequest 768 G4double fraction = StepTaken / StepRequested; 775 message << " Unfinished integration of trac 769 message << " Unfinished integration of track (likely looping particle) " 776 << " of momentum " << momentumVec < 770 << " of momentum " << momentumVec << " ( magnitude = " 777 << momentumVec.mag() << " ) " << G4 771 << momentumVec.mag() << " ) " << G4endl 778 << " after " << count << " field su 772 << " after " << count << " field substeps " 779 << " totaling " << std::setprecisio 773 << " totaling " << std::setprecision(12) << StepTaken / mm << " mm " 780 << " out of requested step " << std 774 << " out of requested step " << std::setprecision(12) 781 << StepRequested / mm << " mm "; 775 << StepRequested / mm << " mm "; 782 message << " a fraction of "; 776 message << " a fraction of "; 783 G4int prec = 4; 777 G4int prec = 4; 784 if( fraction > 0.99 ) 778 if( fraction > 0.99 ) 785 { 779 { 786 prec = 7; 780 prec = 7; 787 } 781 } 788 else 782 else 789 { 783 { 790 if (fraction > 0.97 ) { prec = 5; } 784 if (fraction > 0.97 ) { prec = 5; } 791 } 785 } 792 message << std::setprecision(prec) 786 message << std::setprecision(prec) 793 << 100. * StepTaken / StepRequested 787 << 100. * StepTaken / StepRequested << " % " << G4endl ; 794 if( pPhysVol != nullptr ) << 788 if( pPhysVol ) 795 { 789 { 796 message << " in volume " << pPhysVol->Get << 790 message << " in volume " << pPhysVol->GetName() ; 797 auto material = pPhysVol->GetLogicalVolum << 791 auto material = pPhysVol->GetLogicalVolume()->GetMaterial(); 798 if( material != nullptr ) << 792 if( material != nullptr ) 799 { << 793 message << " with material " << material->GetName() 800 message << " with material " << materia << 794 << " ( density = " 801 << " ( density = " << 795 << material->GetDensity() / ( g/(cm*cm*cm) ) << " g / cm^3 ) "; 802 << material->GetDensity() / ( g << 803 } << 804 } 796 } 805 else 797 else 806 { 798 { 807 message << " in unknown (null) volume. " << 799 message << " in unknown (null) volume. " ; 808 } 800 } 809 G4Exception(methodName, "GeomNav1002", Just 801 G4Exception(methodName, "GeomNav1002", JustWarning, message); 810 } 802 } 811 803 812 // ------------------------------------------- 804 // --------------------------------------------------------------------------- 813 // 805 // 814 void G4PropagatorInField::ReportStuckParticle( 806 void G4PropagatorInField::ReportStuckParticle( G4int noZeroSteps, 815 807 G4double proposedStep, 816 808 G4double lastTriedStep, 817 809 G4VPhysicalVolume* physVol ) 818 { 810 { 819 std::ostringstream message; 811 std::ostringstream message; 820 message << "Particle is stuck; it will be k 812 message << "Particle is stuck; it will be killed." << G4endl 821 << " Zero progress for " << noZero 813 << " Zero progress for " << noZeroSteps << " attempted steps." 822 << G4endl 814 << G4endl 823 << " Proposed Step is " << propose 815 << " Proposed Step is " << proposedStep 824 << " but Step Taken is "<< lastTrie 816 << " but Step Taken is "<< lastTriedStep << G4endl; 825 if( physVol != nullptr ) 817 if( physVol != nullptr ) 826 { << 827 message << " in volume " << physVol->Get 818 message << " in volume " << physVol->GetName() ; 828 } << 829 else 819 else 830 { << 831 message << " in unknown or null volume. 820 message << " in unknown or null volume. " ; 832 } << 833 G4Exception("G4PropagatorInField::ComputeSt 821 G4Exception("G4PropagatorInField::ComputeStep()", 834 "GeomNav1002", JustWarning, mes 822 "GeomNav1002", JustWarning, message); 835 } << 836 << 837 // ------------------------------------------- << 838 << 839 // ------------------------------------------- << 840 // Methods to alter Parameters << 841 // ------------------------------------------- << 842 << 843 // Was a data member (of an object) -- now mov << 844 G4double G4PropagatorInField::GetLargestAccep << 845 { << 846 return fLargestAcceptableStep; << 847 } << 848 << 849 // ------------------------------------------- << 850 // << 851 void G4PropagatorInField::SetLargestAcceptable << 852 { << 853 if( fLargestAcceptableStep>0.0 ) << 854 { << 855 fLargestAcceptableStep = newBigDist; << 856 } << 857 } << 858 << 859 // ------------------------------------------- << 860 << 861 G4double G4PropagatorInField::GetMaxStepSizeMu << 862 { << 863 return fMaxStepSizeMultiplier; << 864 } << 865 << 866 // ------------------------------------------- << 867 << 868 void G4PropagatorInField::SetMaxStepSizeMu << 869 { << 870 fMaxStepSizeMultiplier=vm; << 871 } << 872 << 873 // ------------------------------------------- << 874 << 875 G4double G4PropagatorInField::GetMinBigDistanc << 876 { << 877 return fMinBigDistance; << 878 } << 879 << 880 // ------------------------------------------- << 881 << 882 void G4PropagatorInField::SetMinBigDistanc << 883 { << 884 fMinBigDistance= val; << 885 } 823 } 886 824