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