Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************* 24 // ******************************************************************** 25 // 25 // 26 // 26 // >> 27 // $Id: G4Transportation.cc 2011/06/10 16:19:46 japost Exp japost $ 27 // 28 // 28 // ------------------------------------------- 29 // ------------------------------------------------------------ 29 // GEANT 4 include file implementation 30 // GEANT 4 include file implementation 30 // 31 // 31 // ------------------------------------------- 32 // ------------------------------------------------------------ 32 // 33 // 33 // This class is a process responsible for the 34 // This class is a process responsible for the transportation of 34 // a particle, ie the geometrical propagation 35 // a particle, ie the geometrical propagation that encounters the 35 // geometrical sub-volumes of the detectors. 36 // geometrical sub-volumes of the detectors. 36 // 37 // 37 // It is also tasked with the key role of prop 38 // It is also tasked with the key role of proposing the "isotropic safety", 38 // which will be used to update the post-ste 39 // which will be used to update the post-step point's safety. 39 // 40 // 40 // =========================================== 41 // ======================================================================= >> 42 // Modified: >> 43 // 10 Jan 2015, M.Kelsey: Use G4DynamicParticle mass, NOT PDGMass >> 44 // 28 Oct 2011, P.Gumpl./J.Ap: Detect gravity field, use magnetic moment >> 45 // 20 Nov 2008, J.Apostolakis: Push safety to helper - after ComputeSafety >> 46 // 9 Nov 2007, J.Apostolakis: Flag for short steps, push safety to helper >> 47 // 19 Jan 2006, P.MoraDeFreitas: Fix for suspended tracks (StartTracking) >> 48 // 11 Aug 2004, M.Asai: Add G4VSensitiveDetector* for updating stepPoint. >> 49 // 21 June 2003, J.Apostolakis: Calling field manager with >> 50 // track, to enable it to configure its accuracy >> 51 // 13 May 2003, J.Apostolakis: Zero field areas now taken into >> 52 // account correclty in all cases (thanks to W Pokorski). >> 53 // 29 June 2001, J.Apostolakis, D.Cote-Ahern, P.Gumplinger: >> 54 // correction for spin tracking >> 55 // 20 Febr 2001, J.Apostolakis: update for new FieldTrack >> 56 // 22 Sept 2000, V.Grichine: update of Kinetic Energy 41 // Created: 19 March 1997, J. Apostolakis 57 // Created: 19 March 1997, J. Apostolakis 42 // =========================================== 58 // ======================================================================= 43 59 44 #include "G4Transportation.hh" 60 #include "G4Transportation.hh" 45 #include "G4TransportationProcessType.hh" 61 #include "G4TransportationProcessType.hh" 46 #include "G4TransportationLogger.hh" << 47 62 48 #include "G4PhysicalConstants.hh" 63 #include "G4PhysicalConstants.hh" 49 #include "G4SystemOfUnits.hh" 64 #include "G4SystemOfUnits.hh" 50 #include "G4ProductionCutsTable.hh" 65 #include "G4ProductionCutsTable.hh" 51 #include "G4ParticleTable.hh" 66 #include "G4ParticleTable.hh" 52 67 53 #include "G4ChargeState.hh" 68 #include "G4ChargeState.hh" 54 #include "G4EquationOfMotion.hh" 69 #include "G4EquationOfMotion.hh" 55 70 56 #include "G4FieldManagerStore.hh" 71 #include "G4FieldManagerStore.hh" 57 72 58 #include "G4Navigator.hh" << 59 #include "G4PropagatorInField.hh" << 60 #include "G4TransportationManager.hh" << 61 << 62 #include "G4TransportationParameters.hh" << 63 << 64 class G4VSensitiveDetector; 73 class G4VSensitiveDetector; 65 74 66 G4bool G4Transportation::fUseMagneticMoment=fa 75 G4bool G4Transportation::fUseMagneticMoment=false; 67 G4bool G4Transportation::fUseGravity= false; << 76 68 G4bool G4Transportation::fSilenceLooperWarning << 77 // #define G4DEBUG_TRANSPORT 1 69 78 70 ////////////////////////////////////////////// 79 ////////////////////////////////////////////////////////////////////////// 71 // 80 // 72 // Constructor 81 // Constructor 73 82 74 G4Transportation::G4Transportation( G4int verb << 83 G4Transportation::G4Transportation( G4int verbosity ) 75 : G4VProcess( aName, fTransportation ), << 84 : G4VProcess( G4String("Transportation"), fTransportation ), >> 85 fTransportEndPosition( 0.0, 0.0, 0.0 ), >> 86 fTransportEndMomentumDir( 0.0, 0.0, 0.0 ), >> 87 fTransportEndKineticEnergy( 0.0 ), >> 88 fTransportEndSpin( 0.0, 0.0, 0.0 ), >> 89 fMomentumChanged(true), >> 90 fEndGlobalTimeComputed(false), >> 91 fCandidateEndGlobalTime(0.0), >> 92 fParticleIsLooping( false ), >> 93 fNewTrack( true ), >> 94 fFirstStepInVolume( true ), >> 95 fLastStepInVolume( false ), >> 96 fGeometryLimitedStep(true), 76 fFieldExertedForce( false ), 97 fFieldExertedForce( false ), 77 fPreviousSftOrigin( 0.,0.,0. ), 98 fPreviousSftOrigin( 0.,0.,0. ), 78 fPreviousSafety( 0.0 ), 99 fPreviousSafety( 0.0 ), >> 100 // fParticleChange(), 79 fEndPointDistance( -1.0 ), 101 fEndPointDistance( -1.0 ), >> 102 fThreshold_Warning_Energy( 100 * MeV ), >> 103 fThreshold_Important_Energy( 250 * MeV ), >> 104 fThresholdTrials( 10 ), >> 105 fNoLooperTrials( 0 ), >> 106 fSumEnergyKilled( 0.0 ), fMaxEnergyKilled( 0.0 ), 80 fShortStepOptimisation( false ) // Old def 107 fShortStepOptimisation( false ) // Old default: true (=fast short steps) 81 { 108 { >> 109 // set Process Sub Type 82 SetProcessSubType(static_cast<G4int>(TRANSPO 110 SetProcessSubType(static_cast<G4int>(TRANSPORTATION)); 83 pParticleChange= &fParticleChange; // Requ 111 pParticleChange= &fParticleChange; // Required to conform to G4VProcess 84 SetVerboseLevel(verbosity); 112 SetVerboseLevel(verbosity); 85 113 86 G4TransportationManager* transportMgr ; 114 G4TransportationManager* transportMgr ; 87 115 88 transportMgr = G4TransportationManager::GetT 116 transportMgr = G4TransportationManager::GetTransportationManager() ; 89 117 90 fLinearNavigator = transportMgr->GetNavigato 118 fLinearNavigator = transportMgr->GetNavigatorForTracking() ; 91 119 92 fFieldPropagator = transportMgr->GetPropagat 120 fFieldPropagator = transportMgr->GetPropagatorInField() ; 93 121 94 fpSafetyHelper = transportMgr->GetSafetyHe 122 fpSafetyHelper = transportMgr->GetSafetyHelper(); // New 95 123 96 fpLogger = new G4TransportationLogger("G4Tra << 124 // Cannot determine whether a field exists here, as it would 97 << 125 // depend on the relative order of creating the detector's 98 if( G4TransportationParameters::Exists() ) << 126 // field and this process. That order is not guaranted. 99 { << 127 // Instead later the method DoesGlobalFieldExist() is called 100 auto trParams= G4TransportationParameters: << 101 << 102 SetThresholdWarningEnergy( trParams->GetW << 103 SetThresholdImportantEnergy( trParams->Get << 104 SetThresholdTrials( trParams->GetNumberOfT << 105 G4Transportation::fSilenceLooperWarnings= << 106 } << 107 else { << 108 SetHighLooperThresholds(); << 109 // Use the old defaults: Warning = 100 Me << 110 } << 111 128 112 PushThresholdsToLogger(); << 113 // Should be done by Set methods in SetHighL << 114 << 115 static G4ThreadLocal G4TouchableHandle* pNul 129 static G4ThreadLocal G4TouchableHandle* pNullTouchableHandle = 0; 116 if ( !pNullTouchableHandle) << 130 if ( !pNullTouchableHandle) { pNullTouchableHandle = new G4TouchableHandle; } 117 { << 118 pNullTouchableHandle = new G4TouchableHand << 119 } << 120 fCurrentTouchableHandle = *pNullTouchableHan 131 fCurrentTouchableHandle = *pNullTouchableHandle; 121 // Points to (G4VTouchable*) 0 132 // Points to (G4VTouchable*) 0 122 133 123 << 124 #ifdef G4VERBOSE 134 #ifdef G4VERBOSE 125 if( verboseLevel > 0) 135 if( verboseLevel > 0) 126 { 136 { 127 G4cout << " G4Transportation constructor> 137 G4cout << " G4Transportation constructor> set fShortStepOptimisation to "; 128 if ( fShortStepOptimisation ) { G4cout < << 138 if ( fShortStepOptimisation ) G4cout << "true" << G4endl; 129 else { G4cout < << 139 else G4cout << "false" << G4endl; 130 } 140 } 131 #endif 141 #endif 132 } 142 } 133 143 134 ////////////////////////////////////////////// 144 ////////////////////////////////////////////////////////////////////////// 135 145 136 G4Transportation::~G4Transportation() 146 G4Transportation::~G4Transportation() 137 { 147 { 138 if( fSumEnergyKilled > 0.0 ) << 148 if( (verboseLevel > 0) && (fSumEnergyKilled > 0.0 ) ) 139 { << 149 { 140 PrintStatistics( G4cout ); << 150 G4cout << " G4Transportation: Statistics for looping particles " << G4endl; 141 } << 151 G4cout << " Sum of energy of loopers killed: " << fSumEnergyKilled << G4endl; 142 delete fpLogger; << 152 G4cout << " Max energy of loopers killed: " << fMaxEnergyKilled << G4endl; 143 } << 153 } 144 << 145 ////////////////////////////////////////////// << 146 << 147 void << 148 G4Transportation::PrintStatistics( std::ostrea << 149 { << 150 outStr << " G4Transportation: Statistics fo << 151 if( fSumEnergyKilled > 0.0 || fNumLoopersKi << 152 { << 153 outStr << " Sum of energy of looping t << 154 << fSumEnergyKilled / CLHEP::MeV << 155 << " from " << fNumLoopersKilled << 156 << " Sum of energy of non-elect << 157 << fSumEnergyKilled_NonElectron / << 158 << " from " << fNumLoopersKilled << 159 << G4endl; << 160 outStr << " Max energy of *any type* << 161 << " its PDG was " << fMaxEner << 162 if( fMaxEnergyKilled_NonElectron > 0.0 ) << 163 { << 164 outStr << " Max energy of non-elect << 165 << fMaxEnergyKilled_NonElectro << 166 << " its PDG was " << fMaxE << 167 } << 168 if( fMaxEnergySaved > 0.0 ) << 169 { << 170 outStr << " Max energy of loopers ' << 171 outStr << " Sum of energy of looper << 172 << fSumEnergySaved << G4endl; << 173 outStr << " Sum of energy of unstab << 174 << fSumEnergyUnstableSaved << << 175 } << 176 } << 177 else << 178 { << 179 outStr << " No looping tracks found or k << 180 } << 181 } 154 } 182 155 183 ////////////////////////////////////////////// 156 ////////////////////////////////////////////////////////////////////////// 184 // 157 // 185 // Responsibilities: 158 // Responsibilities: 186 // Find whether the geometry limits the Ste 159 // Find whether the geometry limits the Step, and to what length 187 // Calculate the new value of the safety an 160 // Calculate the new value of the safety and return it. 188 // Store the final time, position and momen 161 // Store the final time, position and momentum. 189 162 190 G4double G4Transportation::AlongStepGetPhysica << 163 G4double G4Transportation:: 191 const G4Track& track, << 164 AlongStepGetPhysicalInteractionLength( const G4Track& track, 192 G4double, // previousStepSize << 165 G4double, // previousStepSize 193 G4double currentMinimumStep, G4double& curre << 166 G4double currentMinimumStep, 194 G4GPILSelection* selection) << 167 G4double& currentSafety, >> 168 G4GPILSelection* selection ) 195 { 169 { 196 // Initial actions moved to StartTrack() << 170 G4double geometryStepLength= -1.0, newSafety= -1.0; >> 171 fParticleIsLooping = false ; >> 172 >> 173 // Initial actions moved to StartTrack() 197 // -------------------------------------- 174 // -------------------------------------- 198 // Note: in case another process changes tou 175 // Note: in case another process changes touchable handle 199 // it will be necessary to add here (for << 176 // it will be necessary to add here (for all steps) 200 // fCurrentTouchableHandle = aTrack->GetTouc 177 // fCurrentTouchableHandle = aTrack->GetTouchableHandle(); 201 178 202 // GPILSelection is set to defaule value of 179 // GPILSelection is set to defaule value of CandidateForSelection 203 // It is a return value 180 // It is a return value 204 // 181 // 205 *selection = CandidateForSelection; << 182 *selection = CandidateForSelection ; >> 183 >> 184 fFirstStepInVolume= fNewTrack || fLastStepInVolume; >> 185 fLastStepInVolume= false; >> 186 fNewTrack = false; 206 187 207 // Get initial Energy/Momentum of the track 188 // Get initial Energy/Momentum of the track 208 // 189 // 209 const G4ThreeVector startPosition = track << 190 const G4DynamicParticle* pParticle = track.GetDynamicParticle() ; 210 const G4ThreeVector startMomentumDir = track << 191 const G4ParticleDefinition* pParticleDef = pParticle->GetDefinition() ; >> 192 G4ThreeVector startMomentumDir = pParticle->GetMomentumDirection() ; >> 193 G4ThreeVector startPosition = track.GetPosition() ; 211 194 212 // The Step Point safety can be limited by o << 195 // G4double theTime = track.GetGlobalTime() ; >> 196 >> 197 // The Step Point safety can be limited by other geometries and/or the 213 // assumptions of any process - it's not alw 198 // assumptions of any process - it's not always the geometrical safety. 214 // We calculate the starting point's isotrop 199 // We calculate the starting point's isotropic safety here. >> 200 // >> 201 G4ThreeVector OriginShift = startPosition - fPreviousSftOrigin ; >> 202 G4double MagSqShift = OriginShift.mag2() ; >> 203 if( MagSqShift >= sqr(fPreviousSafety) ) 215 { 204 { 216 const G4double MagSqShift = (startPosition << 205 currentSafety = 0.0 ; 217 << 206 } 218 if(MagSqShift >= sqr(fPreviousSafety)) << 207 else 219 currentSafety = 0.0; << 208 { 220 else << 209 currentSafety = fPreviousSafety - std::sqrt(MagSqShift) ; 221 currentSafety = fPreviousSafety - std::s << 222 } 210 } 223 211 224 // Is the particle charged or has it a magne 212 // Is the particle charged or has it a magnetic moment? 225 // 213 // 226 const G4DynamicParticle* pParticle = track.G << 214 G4double particleCharge = pParticle->GetCharge() ; >> 215 G4double magneticMoment = pParticle->GetMagneticMoment() ; >> 216 G4double restMass = pParticle->GetMass() ; 227 217 228 const G4double particleMass = pParticle->G << 218 fGeometryLimitedStep = false ; 229 const G4double particleCharge = pParticle->G << 219 // fEndGlobalTimeComputed = false ; 230 const G4double kineticEnergy = pParticle->Ge << 231 << 232 const G4double magneticMoment = pParticle << 233 const G4ThreeVector particleSpin = pParticle << 234 220 235 // There is no need to locate the current vo 221 // There is no need to locate the current volume. It is Done elsewhere: 236 // On track construction << 222 // On track construction 237 // By the tracking, after all AlongStepDoI 223 // By the tracking, after all AlongStepDoIts, in "Relocation" 238 224 239 // Check if the particle has a force, EM or 225 // Check if the particle has a force, EM or gravitational, exerted on it 240 // 226 // >> 227 G4FieldManager* fieldMgr=0; >> 228 G4bool fieldExertsForce = false ; 241 229 242 G4bool eligibleEM = << 230 G4bool gravityOn = false; 243 (particleCharge != 0.0) || ((magneticMomen << 231 G4bool fieldExists= false; // Field is not 0 (null pointer) 244 G4bool eligibleGrav = (particleMass != 0.0) << 245 << 246 fFieldExertedForce = false; << 247 232 248 if(eligibleEM || eligibleGrav) << 233 fieldMgr = fFieldPropagator->FindAndSetFieldManager( track.GetVolume() ); >> 234 if( fieldMgr != 0 ) 249 { 235 { 250 if(G4FieldManager* fieldMgr = << 236 // Message the field Manager, to configure it for this track 251 fFieldPropagator->FindAndSetFieldMana << 237 fieldMgr->ConfigureForTrack( &track ); 252 { << 238 // Is here to allow a transition from no-field pointer 253 // User can configure the field Manager << 239 // to finite field (non-zero pointer). 254 fieldMgr->ConfigureForTrack(&track); << 240 255 // Called here to allow a transition fro << 241 // If the field manager has no field ptr, the field is zero 256 // to finite field (non-zero pointer). << 242 // by definition ( = there is no field ! ) 257 << 243 const G4Field* ptrField= fieldMgr->GetDetectorField(); 258 // If the field manager has no field ptr << 244 fieldExists = (ptrField!=0) ; 259 // by definition ( = there is no field << 245 if( fieldExists ) 260 if(const G4Field* ptrField = fieldMgr->G << 246 { 261 fFieldExertedForce = << 247 gravityOn= ptrField->IsGravityActive(); 262 eligibleEM || (eligibleGrav && ptrFi << 248 263 } << 249 if( (particleCharge != 0.0) 264 } << 250 || (fUseMagneticMoment && (magneticMoment != 0.0) ) 265 << 251 || (gravityOn && (restMass != 0.0) ) 266 G4double geometryStepLength = currentMinimum << 252 ) 267 << 253 { 268 if(currentMinimumStep == 0.0) << 254 fieldExertsForce = fieldExists; 269 { << 255 } 270 fEndPointDistance = 0.0; << 256 } 271 fGeometryLimitedStep = false; // Old cod << 272 // Changed to avoid problems when setting << 273 fMomentumChanged = false; << 274 fParticleIsLooping = false; << 275 fEndGlobalTimeComputed = false; << 276 fTransportEndPosition = startPosition << 277 fTransportEndMomentumDir = startMomentum << 278 fTransportEndKineticEnergy = kineticEnergy << 279 fTransportEndSpin = particleSpin; << 280 } 257 } 281 else if(!fFieldExertedForce) << 258 // G4cout << " G4Transport: field exerts force= " << fieldExertsForce >> 259 // << " fieldMgr= " << fieldMgr << G4endl; >> 260 fFieldExertedForce = fieldExertsForce; >> 261 >> 262 if( !fieldExertsForce ) 282 { 263 { 283 fGeometryLimitedStep = false; << 264 G4double linearStepLength ; 284 if(geometryStepLength > currentSafety || ! << 265 if( fShortStepOptimisation && (currentMinimumStep <= currentSafety) ) 285 { << 266 { 286 const G4double linearStepLength = fLinea << 267 // The Step is guaranteed to be taken 287 startPosition, startMomentumDir, curre << 268 // 288 << 269 geometryStepLength = currentMinimumStep ; 289 if(linearStepLength <= currentMinimumSte << 270 fGeometryLimitedStep = false ; 290 { << 271 } 291 geometryStepLength = linearStepLength; << 272 else 292 fGeometryLimitedStep = true; << 273 { 293 } << 274 // Find whether the straight path intersects a volume 294 // Remember last safety origin & value. << 275 // 295 // << 276 linearStepLength = fLinearNavigator->ComputeStep( startPosition, 296 fPreviousSftOrigin = startPosition; << 277 startMomentumDir, 297 fPreviousSafety = currentSafety; << 278 currentMinimumStep, 298 fpSafetyHelper->SetCurrentSafety(current << 279 newSafety) ; 299 } << 280 // Remember last safety origin & value. 300 << 281 // 301 fEndPointDistance = geometryStepLength; << 282 fPreviousSftOrigin = startPosition ; 302 << 283 fPreviousSafety = newSafety ; 303 fMomentumChanged = false; << 284 fpSafetyHelper->SetCurrentSafety( newSafety, startPosition); 304 fParticleIsLooping = false; << 285 305 fEndGlobalTimeComputed = false; << 286 currentSafety = newSafety ; 306 fTransportEndPosition = << 287 307 startPosition + geometryStepLength * sta << 288 fGeometryLimitedStep= (linearStepLength <= currentMinimumStep); 308 fTransportEndMomentumDir = startMomentum << 289 if( fGeometryLimitedStep ) 309 fTransportEndKineticEnergy = kineticEnergy << 290 { 310 fTransportEndSpin = particleSpin; << 291 // The geometry limits the Step size (an intersection was found.) 311 } << 292 geometryStepLength = linearStepLength ; 312 else // A field exerts force << 293 } 313 { << 294 else 314 const auto pParticleDef = pParticle->Ge << 295 { 315 const auto particlePDGSpin = pParticleDef- << 296 // The full Step is taken. 316 const auto particlePDGMagM = pParticleDef- << 297 geometryStepLength = currentMinimumStep ; 317 << 298 } 318 auto equationOfMotion = fFieldPropagator-> << 299 } 319 << 300 fEndPointDistance = geometryStepLength ; 320 // The charge can change (dynamic), theref << 321 // << 322 equationOfMotion->SetChargeMomentumMass( << 323 G4ChargeState(particleCharge, magneticMo << 324 pParticle->GetTotalMomentum(), particleM << 325 << 326 G4FieldTrack aFieldTrack(startPosition, << 327 track.GetGlobalTi << 328 startMomentumDir, << 329 particleCharge, p << 330 0.0, // Length a << 331 particlePDGSpin); << 332 << 333 // Do the Transport in the field (non rect << 334 // << 335 const G4double lengthAlongCurve = fFieldPr << 336 aFieldTrack, currentMinimumStep, current << 337 kineticEnergy < fThreshold_Important_Ene << 338 << 339 if(lengthAlongCurve < geometryStepLength) << 340 geometryStepLength = lengthAlongCurve; << 341 << 342 // Remember last safety origin & value. << 343 // << 344 fPreviousSftOrigin = startPosition; << 345 fPreviousSafety = currentSafety; << 346 fpSafetyHelper->SetCurrentSafety(currentSa << 347 301 348 fGeometryLimitedStep = fFieldPropagator->I << 302 // Calculate final position 349 // << 303 // 350 // It is possible that step was reduced in << 304 fTransportEndPosition = startPosition+geometryStepLength*startMomentumDir ; 351 // previous zero steps. To cope with case << 352 // in full, we must rely on PiF to obtain << 353 << 354 G4bool changesEnergy = << 355 fFieldPropagator->GetCurrentFieldManager << 356 << 357 fMomentumChanged = true; << 358 fParticleIsLooping = fFieldPropagator->IsP << 359 305 360 fEndGlobalTimeComputed = changesEnergy; << 306 // Momentum direction, energy and polarisation are unchanged by transport 361 fTransportEndPosition = aFieldTrack.Get << 307 // 362 fTransportEndMomentumDir = aFieldTrack.Get << 308 fTransportEndMomentumDir = startMomentumDir ; >> 309 fTransportEndKineticEnergy = track.GetKineticEnergy() ; >> 310 fTransportEndSpin = track.GetPolarization(); >> 311 fParticleIsLooping = false ; >> 312 fMomentumChanged = false ; >> 313 fEndGlobalTimeComputed = false ; >> 314 } >> 315 else // A field exerts force >> 316 { >> 317 G4double momentumMagnitude = pParticle->GetTotalMomentum() ; >> 318 G4ThreeVector EndUnitMomentum ; >> 319 G4double lengthAlongCurve ; >> 320 >> 321 G4ChargeState chargeState(particleCharge, // The charge can change (dynamic) >> 322 magneticMoment, >> 323 pParticleDef->GetPDGSpin() ); >> 324 // For insurance, could set it again >> 325 // chargeState.SetPDGSpin(pParticleDef->GetPDGSpin() ); // Provisionally in same object >> 326 >> 327 G4EquationOfMotion* equationOfMotion = fFieldPropagator->GetCurrentEquationOfMotion(); >> 328 >> 329 equationOfMotion->SetChargeMomentumMass( chargeState, >> 330 momentumMagnitude, >> 331 restMass); >> 332 >> 333 G4FieldTrack aFieldTrack = G4FieldTrack( startPosition, >> 334 track.GetGlobalTime(), // Lab. >> 335 // track.GetProperTime(), // Particle rest frame >> 336 track.GetMomentumDirection(), >> 337 track.GetKineticEnergy(), >> 338 restMass, >> 339 particleCharge, >> 340 track.GetPolarization(), >> 341 pParticleDef->GetPDGMagneticMoment(), >> 342 0.0, // Length along track >> 343 pParticleDef->GetPDGSpin() >> 344 ) ; >> 345 >> 346 if( currentMinimumStep > 0 ) >> 347 { >> 348 // Do the Transport in the field (non recti-linear) >> 349 // >> 350 lengthAlongCurve = fFieldPropagator->ComputeStep( aFieldTrack, >> 351 currentMinimumStep, >> 352 currentSafety, >> 353 track.GetVolume() ) ; >> 354 >> 355 fGeometryLimitedStep= fFieldPropagator->IsLastStepInVolume(); >> 356 // It is possible that step was reduced in PropagatorInField due to previous zero steps >> 357 // To cope with case that reduced step is taken in full, we must rely on PiF to obtain this >> 358 // value. 363 359 364 // G4cout << " G4Transport: End of step pM << 360 geometryStepLength = std::min( lengthAlongCurve, currentMinimumStep ); >> 361 >> 362 // Remember last safety origin & value. >> 363 // >> 364 fPreviousSftOrigin = startPosition ; >> 365 fPreviousSafety = currentSafety ; >> 366 fpSafetyHelper->SetCurrentSafety( currentSafety, startPosition); >> 367 } >> 368 else >> 369 { >> 370 geometryStepLength = lengthAlongCurve= 0.0 ; >> 371 fGeometryLimitedStep = false ; >> 372 } >> 373 >> 374 // Get the End-Position and End-Momentum (Dir-ection) >> 375 // >> 376 fTransportEndPosition = aFieldTrack.GetPosition() ; 365 377 366 fEndPointDistance = (fTransportEndPositio << 378 // Momentum: Magnitude and direction can be changed too now ... >> 379 // >> 380 fMomentumChanged = true ; >> 381 fTransportEndMomentumDir = aFieldTrack.GetMomentumDir() ; 367 382 368 // Ignore change in energy for fields that << 383 fTransportEndKineticEnergy = aFieldTrack.GetKineticEnergy() ; 369 // This hides the integration error, but g << 370 fTransportEndKineticEnergy = << 371 changesEnergy ? aFieldTrack.GetKineticEn << 372 fTransportEndSpin = aFieldTrack.GetSpin(); << 373 384 374 if(fEndGlobalTimeComputed) << 385 if( fFieldPropagator->GetCurrentFieldManager()->DoesFieldChangeEnergy() ) 375 { << 386 { 376 // If the field can change energy, then << 387 // If the field can change energy, then the time must be integrated 377 // - so this should have been updated << 388 // - so this should have been updated 378 // << 389 // 379 fCandidateEndGlobalTime = aFieldTrack.Ge << 390 fCandidateEndGlobalTime = aFieldTrack.GetLabTimeOfFlight(); >> 391 fEndGlobalTimeComputed = true; 380 392 381 // was ( fCandidateEndGlobalTime != trac << 393 // was ( fCandidateEndGlobalTime != track.GetGlobalTime() ); 382 // a cleaner way is to have FieldTrack k << 394 // a cleaner way is to have FieldTrack knowing whether time is updated. 383 } << 395 } 384 #if defined(G4VERBOSE) || defined(G4DEBUG_TRAN << 396 else 385 else << 397 { 386 { << 398 // The energy should be unchanged by field transport, 387 // The energy should be unchanged by fie << 399 // - so the time changed will be calculated elsewhere 388 // - so the time changed will be calc << 400 // 389 // << 401 fEndGlobalTimeComputed = false; 390 // Check that the integration preserved << 391 // - and if not correct this! << 392 G4double startEnergy = kineticEnergy; << 393 G4double endEnergy = fTransportEndKine << 394 402 395 static G4ThreadLocal G4int no_large_edif << 403 // Check that the integration preserved the energy 396 if(verboseLevel > 1) << 404 // - and if not correct this! 397 { << 405 G4double startEnergy= track.GetKineticEnergy(); 398 if(std::fabs(startEnergy - endEnergy) << 406 G4double endEnergy= fTransportEndKineticEnergy; >> 407 >> 408 static G4ThreadLocal G4int no_inexact_steps=0, no_large_ediff; >> 409 G4double absEdiff = std::fabs(startEnergy- endEnergy); >> 410 if( absEdiff > perMillion * endEnergy ) 399 { 411 { 400 static G4ThreadLocal G4int no_warnin << 412 no_inexact_steps++; 401 moduloFac << 413 // Possible statistics keeping here ... 402 no_large_ediff++; << 414 } 403 if((no_large_ediff % warnModulo) == << 415 if( verboseLevel > 1 ) >> 416 { >> 417 if( std::fabs(startEnergy- endEnergy) > perThousand * endEnergy ) 404 { 418 { 405 no_warnings++; << 419 static G4ThreadLocal G4int no_warnings= 0, warnModulo=1, 406 std::ostringstream message; << 420 moduloFactor= 10; 407 message << "Energy change in Step << 421 no_large_ediff ++; 408 << G4endl << " Relativ << 422 if( (no_large_ediff% warnModulo) == 0 ) 409 << std::setw(15) << (endEn << 410 << G4endl << " Startin << 411 << startEnergy / MeV << " << 412 << " Ending E= " << << 413 << " MeV " << G4endl << 414 << "Energy has been correc << 415 << " field propagation par << 416 if((verboseLevel > 2) || (no_warni << 417 (no_large_ediff == warnModulo * << 418 { << 419 message << "These include Epsilo << 420 << G4endl << 421 << "which determine frac << 422 "integrated quantitie << 423 << G4endl << 424 << "Note also the influe << 425 "integration steps." << 426 << G4endl; << 427 } << 428 message << "Bad 'endpoint'. Energy << 429 << G4endl << "Has occurred << 430 << " times."; << 431 G4Exception("G4Transportation::Alo << 432 JustWarning, message); << 433 if(no_large_ediff == warnModulo * << 434 { 423 { 435 warnModulo *= moduloFactor; << 424 no_warnings++; >> 425 G4cout << "WARNING - G4Transportation::AlongStepGetPIL() " >> 426 << " Energy change in Step is above 1^-3 relative value. " << G4endl >> 427 << " Relative change in 'tracking' step = " >> 428 << std::setw(15) << (endEnergy-startEnergy)/startEnergy << G4endl >> 429 << " Starting E= " << std::setw(12) << startEnergy / MeV << " MeV " << G4endl >> 430 << " Ending E= " << std::setw(12) << endEnergy / MeV << " MeV " << G4endl; >> 431 G4cout << " Energy has been corrected -- however, review" >> 432 << " field propagation parameters for accuracy." << G4endl; >> 433 if( (verboseLevel > 2 ) || (no_warnings<4) || (no_large_ediff == warnModulo * moduloFactor) ) >> 434 { >> 435 G4cout << " These include EpsilonStepMax(/Min) in G4FieldManager " >> 436 << " which determine fractional error per step for integrated quantities. " << G4endl >> 437 << " Note also the influence of the permitted number of integration steps." >> 438 << G4endl; >> 439 } >> 440 G4cerr << "ERROR - G4Transportation::AlongStepGetPIL()" << G4endl >> 441 << " Bad 'endpoint'. Energy change detected" >> 442 << " and corrected. " >> 443 << " Has occurred already " >> 444 << no_large_ediff << " times." << G4endl; >> 445 if( no_large_ediff == warnModulo * moduloFactor ) >> 446 { >> 447 warnModulo *= moduloFactor; >> 448 } 436 } 449 } 437 } 450 } 438 } << 451 } // end of if (verboseLevel) 439 } // end of if (verboseLevel) << 452 440 } << 453 // Correct the energy for fields that conserve it 441 #endif << 454 // This - hides the integration error >> 455 // - but gives a better physical answer >> 456 fTransportEndKineticEnergy= track.GetKineticEnergy(); >> 457 } >> 458 >> 459 fTransportEndSpin = aFieldTrack.GetSpin(); >> 460 fParticleIsLooping = fFieldPropagator->IsParticleLooping() ; >> 461 fEndPointDistance = (fTransportEndPosition - startPosition).mag() ; >> 462 } >> 463 >> 464 // If we are asked to go a step length of 0, and we are on a boundary >> 465 // then a boundary will also limit the step -> we must flag this. >> 466 // >> 467 if( currentMinimumStep == 0.0 ) >> 468 { >> 469 if( currentSafety == 0.0 ) { fGeometryLimitedStep = true; } 442 } 470 } 443 471 444 // Update the safety starting from the end-p 472 // Update the safety starting from the end-point, 445 // if it will become negative at the end-poi 473 // if it will become negative at the end-point. 446 // 474 // 447 if(currentSafety < fEndPointDistance) << 475 if( currentSafety < fEndPointDistance ) 448 { 476 { 449 if(particleCharge != 0.0) << 477 if( particleCharge != 0.0 ) 450 { << 478 { 451 G4double endSafety = << 479 G4double endSafety = 452 fLinearNavigator->ComputeSafety(fTrans << 480 fLinearNavigator->ComputeSafety( fTransportEndPosition) ; 453 currentSafety = endSafety; << 481 currentSafety = endSafety ; 454 fPreviousSftOrigin = fTransportEndPositi << 482 fPreviousSftOrigin = fTransportEndPosition ; 455 fPreviousSafety = currentSafety; << 483 fPreviousSafety = currentSafety ; 456 fpSafetyHelper->SetCurrentSafety(current << 484 fpSafetyHelper->SetCurrentSafety( currentSafety, fTransportEndPosition); 457 << 485 458 // Because the Stepping Manager assumes << 486 // Because the Stepping Manager assumes it is from the start point, 459 // add the StepLength << 487 // add the StepLength 460 // << 488 // 461 currentSafety += fEndPointDistance; << 489 currentSafety += fEndPointDistance ; 462 << 490 463 #ifdef G4DEBUG_TRANSPORT << 491 #ifdef G4DEBUG_TRANSPORT 464 G4cout.precision(12); << 492 G4cout.precision(12) ; 465 G4cout << "***G4Transportation::AlongSte << 493 G4cout << "***G4Transportation::AlongStepGPIL ** " << G4endl ; 466 G4cout << " Called Navigator->ComputeSa << 494 G4cout << " Called Navigator->ComputeSafety at " << fTransportEndPosition 467 << " and it returned safety= << 495 << " and it returned safety= " << endSafety << G4endl ; 468 G4cout << " Adding endpoint distance << 496 G4cout << " Adding endpoint distance " << fEndPointDistance 469 << " to obtain pseudo-safety= << 497 << " to obtain pseudo-safety= " << currentSafety << G4endl ; 470 } << 498 } 471 else << 499 else 472 { << 500 { 473 G4cout << "***G4Transportation::AlongSte << 501 G4cout << "***G4Transportation::AlongStepGPIL ** " << G4endl ; 474 G4cout << " Avoiding call to ComputeSaf << 502 G4cout << " Avoiding call to ComputeSafety : " << G4endl; 475 G4cout << " charge = " << particl << 503 G4cout << " charge = " << particleCharge << G4endl; 476 G4cout << " mag moment = " << magneti << 504 G4cout << " mag moment = " << magneticMoment << G4endl; 477 #endif 505 #endif 478 } << 506 } 479 } << 507 } 480 << 481 fFirstStepInVolume = fNewTrack || fLastStepI << 482 fLastStepInVolume = false; << 483 fNewTrack = false; << 484 508 485 fParticleChange.ProposeFirstStepInVolume(fFi << 509 fParticleChange.ProposeTrueStepLength(geometryStepLength) ; 486 fParticleChange.ProposeTrueStepLength(geomet << 487 510 488 return geometryStepLength; << 511 return geometryStepLength ; 489 } 512 } 490 513 491 ////////////////////////////////////////////// 514 ////////////////////////////////////////////////////////////////////////// 492 // 515 // 493 // Initialize ParticleChange (by setting al 516 // Initialize ParticleChange (by setting all its members equal 494 // to correspond 517 // to corresponding members in G4Track) 495 518 496 G4VParticleChange* G4Transportation::AlongStep 519 G4VParticleChange* G4Transportation::AlongStepDoIt( const G4Track& track, 497 520 const G4Step& stepData ) 498 { 521 { 499 #if defined(G4VERBOSE) || defined(G4DEBUG_TRAN << 522 static G4ThreadLocal G4int noCalls=0; 500 static G4ThreadLocal G4long noCallsASDI=0; << 523 noCalls++; 501 noCallsASDI++; << 502 #else << 503 #define noCallsASDI 0 << 504 #endif << 505 << 506 if(fGeometryLimitedStep) << 507 { << 508 stepData.GetPostStepPoint()->SetStepStatus << 509 } << 510 524 511 fParticleChange.Initialize(track) ; 525 fParticleChange.Initialize(track) ; 512 526 513 // Code for specific process 527 // Code for specific process 514 // 528 // 515 fParticleChange.ProposePosition(fTransportEn 529 fParticleChange.ProposePosition(fTransportEndPosition) ; 516 fParticleChange.ProposeMomentumDirection(fTr 530 fParticleChange.ProposeMomentumDirection(fTransportEndMomentumDir) ; 517 fParticleChange.ProposeEnergy(fTransportEndK 531 fParticleChange.ProposeEnergy(fTransportEndKineticEnergy) ; 518 fParticleChange.SetMomentumChanged(fMomentum 532 fParticleChange.SetMomentumChanged(fMomentumChanged) ; 519 533 520 fParticleChange.ProposePolarization(fTranspo 534 fParticleChange.ProposePolarization(fTransportEndSpin); 521 535 522 G4double deltaTime = 0.0 ; 536 G4double deltaTime = 0.0 ; 523 537 524 // Calculate Lab Time of Flight (ONLY if fi 538 // Calculate Lab Time of Flight (ONLY if field Equations used it!) 525 // G4double endTime = fCandidateEndGlobalT 539 // G4double endTime = fCandidateEndGlobalTime; 526 // G4double delta_time = endTime - startTime 540 // G4double delta_time = endTime - startTime; 527 541 528 G4double startTime = track.GetGlobalTime() ; 542 G4double startTime = track.GetGlobalTime() ; 529 543 530 if (!fEndGlobalTimeComputed) 544 if (!fEndGlobalTimeComputed) 531 { 545 { 532 // The time was not integrated .. make th 546 // The time was not integrated .. make the best estimate possible 533 // 547 // 534 G4double initialVelocity = stepData.GetPr 548 G4double initialVelocity = stepData.GetPreStepPoint()->GetVelocity(); 535 G4double stepLength = track.GetStepL 549 G4double stepLength = track.GetStepLength(); 536 550 537 deltaTime= 0.0; // in case initialVeloci 551 deltaTime= 0.0; // in case initialVelocity = 0 538 if ( initialVelocity > 0.0 ) { deltaTime 552 if ( initialVelocity > 0.0 ) { deltaTime = stepLength/initialVelocity; } 539 553 540 fCandidateEndGlobalTime = startTime + d 554 fCandidateEndGlobalTime = startTime + deltaTime ; 541 fParticleChange.ProposeLocalTime( track. 555 fParticleChange.ProposeLocalTime( track.GetLocalTime() + deltaTime) ; 542 } 556 } 543 else 557 else 544 { 558 { 545 deltaTime = fCandidateEndGlobalTime - sta 559 deltaTime = fCandidateEndGlobalTime - startTime ; 546 fParticleChange.ProposeGlobalTime( fCandi 560 fParticleChange.ProposeGlobalTime( fCandidateEndGlobalTime ) ; 547 } 561 } 548 562 549 563 550 // Now Correct by Lorentz factor to get delt 564 // Now Correct by Lorentz factor to get delta "proper" Time 551 565 552 G4double restMass = track.GetDynamicP 566 G4double restMass = track.GetDynamicParticle()->GetMass() ; 553 G4double deltaProperTime = deltaTime*( restM 567 G4double deltaProperTime = deltaTime*( restMass/track.GetTotalEnergy() ) ; 554 568 555 fParticleChange.ProposeProperTime(track.GetP 569 fParticleChange.ProposeProperTime(track.GetProperTime() + deltaProperTime) ; 556 //fParticleChange.ProposeTrueStepLength( tra << 570 //fParticleChange. ProposeTrueStepLength( track.GetStepLength() ) ; 557 571 558 // If the particle is caught looping or is s 572 // If the particle is caught looping or is stuck (in very difficult 559 // boundaries) in a magnetic field (doing ma << 573 // boundaries) in a magnetic field (doing many steps) >> 574 // THEN this kills it ... 560 // 575 // 561 if ( fParticleIsLooping ) 576 if ( fParticleIsLooping ) 562 { 577 { 563 G4double endEnergy= fTransportEndKinetic 578 G4double endEnergy= fTransportEndKineticEnergy; 564 fNoLooperTrials ++; << 579 565 auto particleType= track.GetDynamicParti << 580 if( (endEnergy < fThreshold_Important_Energy) 566 << 581 || (fNoLooperTrials >= fThresholdTrials ) ) 567 G4bool stable = particleType->GetPDGStab << 568 G4bool candidateForEnd = (endEnergy < fT << 569 || (fNoLooperTrial << 570 G4bool unstableAndKillable = !stable && << 571 G4bool unstableForEnd = (endEnergy < fTh << 572 && (fNoLooperTri << 573 if( (candidateForEnd && stable) || (unst << 574 { 582 { 575 // Kill the looping particle 583 // Kill the looping particle 576 // 584 // 577 fParticleChange.ProposeTrackStatus( fS 585 fParticleChange.ProposeTrackStatus( fStopAndKill ) ; 578 G4int particlePDG= particleType->GetPD << 579 const G4int electronPDG= 11; // G4Elec << 580 586 581 // Simple statistics << 587 // 'Bare' statistics 582 fSumEnergyKilled += endEnergy; << 588 fSumEnergyKilled += endEnergy; 583 fSumEnerSqKilled += endEnergy * endEne << 589 if( endEnergy > fMaxEnergyKilled) { fMaxEnergyKilled= endEnergy; } 584 fNumLoopersKilled++; << 585 << 586 if( endEnergy > fMaxEnergyKilled ) { << 587 fMaxEnergyKilled = endEnergy; << 588 fMaxEnergyKilledPDG = particlePDG; << 589 } << 590 if( particleType->GetPDGEncoding() != << 591 { << 592 fSumEnergyKilled_NonElectron += end << 593 fSumEnerSqKilled_NonElectron += end << 594 fNumLoopersKilled_NonElectron++; << 595 << 596 if( endEnergy > fMaxEnergyKilled_No << 597 { << 598 fMaxEnergyKilled_NonElectron = e << 599 fMaxEnergyKilled_NonElecPDG = p << 600 } << 601 } << 602 590 603 if( endEnergy > fThreshold_Warning_Ene << 591 #ifdef G4VERBOSE 604 { << 592 if( (verboseLevel > 1) && 605 fpLogger->ReportLoopingTrack( track, << 593 ( endEnergy > fThreshold_Warning_Energy ) ) 606 noCall << 594 { >> 595 G4cout << " G4Transportation is killing track that is looping or stuck " >> 596 << G4endl >> 597 << " This track has " << track.GetKineticEnergy() / MeV >> 598 << " MeV energy." << G4endl; >> 599 G4cout << " Number of trials = " << fNoLooperTrials >> 600 << " No of calls to AlongStepDoIt = " << noCalls >> 601 << G4endl; 607 } 602 } >> 603 #endif 608 fNoLooperTrials=0; 604 fNoLooperTrials=0; 609 } 605 } 610 else 606 else 611 { 607 { 612 fMaxEnergySaved = std::max( endEnergy, << 608 fNoLooperTrials ++; 613 if( fNoLooperTrials == 1 ) { << 614 fSumEnergySaved += endEnergy; << 615 if ( !stable ) << 616 fSumEnergyUnstableSaved += endEne << 617 } << 618 #ifdef G4VERBOSE 609 #ifdef G4VERBOSE 619 if( verboseLevel > 2 && ! fSilenceLoop << 610 if( (verboseLevel > 2) ) 620 { 611 { 621 G4cout << " " << __func__ << 612 G4cout << " G4Transportation::AlongStepDoIt(): Particle looping - " 622 << " Particle is looping but << 613 << " Number of trials = " << fNoLooperTrials 623 << " Number of trials = " < << 614 << " No of calls to = " << noCalls 624 << " No of calls to = " << << 615 << G4endl; 625 } 616 } 626 #endif 617 #endif 627 } 618 } 628 } 619 } 629 else 620 else 630 { 621 { 631 fNoLooperTrials=0; 622 fNoLooperTrials=0; 632 } 623 } 633 624 634 // Another (sometimes better way) is to use 625 // Another (sometimes better way) is to use a user-limit maximum Step size 635 // to alleviate this problem .. 626 // to alleviate this problem .. 636 627 637 // Introduce smooth curved trajectories to p 628 // Introduce smooth curved trajectories to particle-change 638 // 629 // 639 fParticleChange.SetPointerToVectorOfAuxiliar 630 fParticleChange.SetPointerToVectorOfAuxiliaryPoints 640 (fFieldPropagator->GimmeTrajectoryVectorAn 631 (fFieldPropagator->GimmeTrajectoryVectorAndForgetIt() ); 641 632 642 return &fParticleChange ; 633 return &fParticleChange ; 643 } 634 } 644 635 645 ////////////////////////////////////////////// 636 ////////////////////////////////////////////////////////////////////////// 646 // 637 // 647 // This ensures that the PostStep action is a 638 // This ensures that the PostStep action is always called, 648 // so that it can do the relocation if it is 639 // so that it can do the relocation if it is needed. 649 // 640 // 650 641 651 G4double G4Transportation:: 642 G4double G4Transportation:: 652 PostStepGetPhysicalInteractionLength( const G4 643 PostStepGetPhysicalInteractionLength( const G4Track&, 653 G4 644 G4double, // previousStepSize 654 G4 645 G4ForceCondition* pForceCond ) 655 { 646 { 656 fFieldExertedForce = false; // Not known 647 fFieldExertedForce = false; // Not known 657 *pForceCond = Forced ; 648 *pForceCond = Forced ; 658 return DBL_MAX ; // was kInfinity ; but con 649 return DBL_MAX ; // was kInfinity ; but convention now is DBL_MAX 659 } 650 } 660 651 661 ////////////////////////////////////////////// 652 ///////////////////////////////////////////////////////////////////////////// 662 << 663 void G4Transportation::SetTouchableInformation << 664 { << 665 const G4VPhysicalVolume* pNewVol = touchable << 666 const G4Material* pNewMaterial = 0 ; << 667 G4VSensitiveDetector* pNewSensitiveDetector << 668 << 669 if( pNewVol != 0 ) << 670 { << 671 pNewMaterial= pNewVol->GetLogicalVolume()- << 672 pNewSensitiveDetector= pNewVol->GetLogical << 673 } << 674 << 675 fParticleChange.SetMaterialInTouchable( (G4M << 676 fParticleChange.SetSensitiveDetectorInToucha << 677 // temporarily until Get/Set Material of Par << 678 // and StepPoint can be made const. << 679 << 680 const G4MaterialCutsCouple* pNewMaterialCuts << 681 if( pNewVol != 0 ) << 682 { << 683 pNewMaterialCutsCouple=pNewVol->GetLogical << 684 } << 685 << 686 if ( pNewVol!=0 && pNewMaterialCutsCouple!=0 << 687 && pNewMaterialCutsCouple->GetMaterial()!= << 688 { << 689 // for parametrized volume << 690 // << 691 pNewMaterialCutsCouple = << 692 G4ProductionCutsTable::GetProductionCuts << 693 ->GetMaterialCuts << 694 pNewMaterialCut << 695 } << 696 fParticleChange.SetMaterialCutsCoupleInTouch << 697 << 698 // Set the touchable in ParticleChange << 699 // this must always be done because the part << 700 // uses this value to overwrite the current << 701 // << 702 fParticleChange.SetTouchableHandle(touchable << 703 } << 704 << 705 ////////////////////////////////////////////// << 706 // 653 // 707 654 708 G4VParticleChange* G4Transportation::PostStepD 655 G4VParticleChange* G4Transportation::PostStepDoIt( const G4Track& track, 709 656 const G4Step& ) 710 { 657 { 711 G4TouchableHandle retCurrentTouchable ; / 658 G4TouchableHandle retCurrentTouchable ; // The one to return 712 G4bool isLastStep= false; 659 G4bool isLastStep= false; 713 660 714 // Initialize ParticleChange (by setting al 661 // Initialize ParticleChange (by setting all its members equal 715 // to correspond 662 // to corresponding members in G4Track) 716 // fParticleChange.Initialize(track) ; // T 663 // fParticleChange.Initialize(track) ; // To initialise TouchableChange 717 664 718 fParticleChange.ProposeTrackStatus(track.Get 665 fParticleChange.ProposeTrackStatus(track.GetTrackStatus()) ; 719 666 720 // If the Step was determined by the volume 667 // If the Step was determined by the volume boundary, 721 // logically relocate the particle 668 // logically relocate the particle 722 669 723 if(fGeometryLimitedStep) 670 if(fGeometryLimitedStep) 724 { 671 { 725 // fCurrentTouchable will now become the p 672 // fCurrentTouchable will now become the previous touchable, 726 // and what was the previous will be freed 673 // and what was the previous will be freed. 727 // (Needed because the preStepPoint can po 674 // (Needed because the preStepPoint can point to the previous touchable) 728 675 729 fLinearNavigator->SetGeometricallyLimitedS 676 fLinearNavigator->SetGeometricallyLimitedStep() ; 730 fLinearNavigator-> 677 fLinearNavigator-> 731 LocateGlobalPointAndUpdateTouchableHandle( 678 LocateGlobalPointAndUpdateTouchableHandle( track.GetPosition(), 732 679 track.GetMomentumDirection(), 733 680 fCurrentTouchableHandle, 734 681 true ) ; 735 // Check whether the particle is out of th 682 // Check whether the particle is out of the world volume 736 // If so it has exited and must be killed. 683 // If so it has exited and must be killed. 737 // 684 // 738 if( fCurrentTouchableHandle->GetVolume() = 685 if( fCurrentTouchableHandle->GetVolume() == 0 ) 739 { 686 { 740 fParticleChange.ProposeTrackStatus( fSt 687 fParticleChange.ProposeTrackStatus( fStopAndKill ) ; 741 } 688 } 742 retCurrentTouchable = fCurrentTouchableHan 689 retCurrentTouchable = fCurrentTouchableHandle ; 743 fParticleChange.SetTouchableHandle( fCurre 690 fParticleChange.SetTouchableHandle( fCurrentTouchableHandle ) ; 744 691 745 // Update the Step flag which identifies t 692 // Update the Step flag which identifies the Last Step in a volume 746 if( !fFieldExertedForce ) 693 if( !fFieldExertedForce ) 747 isLastStep = fLinearNavigator->ExitedM 694 isLastStep = fLinearNavigator->ExitedMotherVolume() 748 || fLinearNavigator->Entered << 695 | fLinearNavigator->EnteredDaughterVolume() ; 749 else 696 else 750 isLastStep = fFieldPropagator->IsLastSt 697 isLastStep = fFieldPropagator->IsLastStepInVolume(); 751 } 698 } 752 else // fGeometryLimitedStep 699 else // fGeometryLimitedStep is false 753 { 700 { 754 // This serves only to move the Navigator' 701 // This serves only to move the Navigator's location 755 // 702 // 756 fLinearNavigator->LocateGlobalPointWithinV 703 fLinearNavigator->LocateGlobalPointWithinVolume( track.GetPosition() ) ; 757 704 758 // The value of the track's current Toucha 705 // The value of the track's current Touchable is retained. 759 // (and it must be correct because we must 706 // (and it must be correct because we must use it below to 760 // overwrite the (unset) one in particle c 707 // overwrite the (unset) one in particle change) 761 // It must be fCurrentTouchable too ?? 708 // It must be fCurrentTouchable too ?? 762 // 709 // 763 fParticleChange.SetTouchableHandle( track. 710 fParticleChange.SetTouchableHandle( track.GetTouchableHandle() ) ; 764 retCurrentTouchable = track.GetTouchableHa 711 retCurrentTouchable = track.GetTouchableHandle() ; 765 712 766 isLastStep= false; 713 isLastStep= false; 767 } // endif ( fGeometryLimitedStep ) 714 } // endif ( fGeometryLimitedStep ) 768 fLastStepInVolume= isLastStep; 715 fLastStepInVolume= isLastStep; 769 716 770 fParticleChange.ProposeFirstStepInVolume(fFi 717 fParticleChange.ProposeFirstStepInVolume(fFirstStepInVolume); 771 fParticleChange.ProposeLastStepInVolume(isLa 718 fParticleChange.ProposeLastStepInVolume(isLastStep); 772 719 773 SetTouchableInformation(retCurrentTouchable) << 720 const G4VPhysicalVolume* pNewVol = retCurrentTouchable->GetVolume() ; >> 721 const G4Material* pNewMaterial = 0 ; >> 722 const G4VSensitiveDetector* pNewSensitiveDetector = 0 ; >> 723 >> 724 if( pNewVol != 0 ) >> 725 { >> 726 pNewMaterial= pNewVol->GetLogicalVolume()->GetMaterial(); >> 727 pNewSensitiveDetector= pNewVol->GetLogicalVolume()->GetSensitiveDetector(); >> 728 } >> 729 >> 730 // ( <const_cast> pNewMaterial ) ; >> 731 // ( <const_cast> pNewSensitiveDetector) ; >> 732 >> 733 fParticleChange.SetMaterialInTouchable( (G4Material *) pNewMaterial ) ; >> 734 fParticleChange.SetSensitiveDetectorInTouchable( (G4VSensitiveDetector *) pNewSensitiveDetector ) ; >> 735 >> 736 const G4MaterialCutsCouple* pNewMaterialCutsCouple = 0; >> 737 if( pNewVol != 0 ) >> 738 { >> 739 pNewMaterialCutsCouple=pNewVol->GetLogicalVolume()->GetMaterialCutsCouple(); >> 740 } >> 741 >> 742 if( pNewVol!=0 && pNewMaterialCutsCouple!=0 && pNewMaterialCutsCouple->GetMaterial()!=pNewMaterial ) >> 743 { >> 744 // for parametrized volume >> 745 // >> 746 pNewMaterialCutsCouple = >> 747 G4ProductionCutsTable::GetProductionCutsTable() >> 748 ->GetMaterialCutsCouple(pNewMaterial, >> 749 pNewMaterialCutsCouple->GetProductionCuts()); >> 750 } >> 751 fParticleChange.SetMaterialCutsCoupleInTouchable( pNewMaterialCutsCouple ); >> 752 >> 753 // temporarily until Get/Set Material of ParticleChange, >> 754 // and StepPoint can be made const. >> 755 // Set the touchable in ParticleChange >> 756 // this must always be done because the particle change always >> 757 // uses this value to overwrite the current touchable pointer. >> 758 // >> 759 fParticleChange.SetTouchableHandle(retCurrentTouchable) ; 774 760 775 return &fParticleChange ; 761 return &fParticleChange ; 776 } 762 } 777 763 778 ////////////////////////////////////////////// << 764 // New method takes over the responsibility to reset the state of G4Transportation 779 // New method takes over the responsibility to << 765 // object at the start of a new track or the resumption of a suspended track. 780 // G4Transportation object at the start of a n << 781 // of a suspended track. << 782 // << 783 766 784 void 767 void 785 G4Transportation::StartTracking(G4Track* aTrac 768 G4Transportation::StartTracking(G4Track* aTrack) 786 { 769 { 787 G4VProcess::StartTracking(aTrack); 770 G4VProcess::StartTracking(aTrack); 788 fNewTrack= true; 771 fNewTrack= true; 789 fFirstStepInVolume= true; 772 fFirstStepInVolume= true; 790 fLastStepInVolume= false; 773 fLastStepInVolume= false; 791 774 792 // The actions here are those that were take 775 // The actions here are those that were taken in AlongStepGPIL 793 // when track.GetCurrentStepNumber()==1 << 776 // when track.GetCurrentStepNumber()==1 794 777 795 // Whether field exists should be determined << 796 G4FieldManagerStore* fieldMgrStore= G4FieldM << 797 fAnyFieldExists = fieldMgrStore->size() > 0; << 798 << 799 // reset safety value and center 778 // reset safety value and center 800 // 779 // 801 fPreviousSafety = 0.0 ; 780 fPreviousSafety = 0.0 ; 802 fPreviousSftOrigin = G4ThreeVector(0.,0.,0.) 781 fPreviousSftOrigin = G4ThreeVector(0.,0.,0.) ; 803 782 804 // reset looping counter -- for motion in fi 783 // reset looping counter -- for motion in field 805 fNoLooperTrials= 0; 784 fNoLooperTrials= 0; 806 // Must clear this state .. else it depends 785 // Must clear this state .. else it depends on last track's value 807 // --> a better solution would set this fro 786 // --> a better solution would set this from state of suspended track TODO ? 808 // Was if( aTrack->GetCurrentStepNumber()==1 787 // Was if( aTrack->GetCurrentStepNumber()==1 ) { .. } 809 788 810 // ChordFinder reset internal state 789 // ChordFinder reset internal state 811 // 790 // 812 if( fFieldPropagator && fAnyFieldExists ) << 791 if( DoesGlobalFieldExist() ) 813 { 792 { 814 fFieldPropagator->ClearPropagatorState(); 793 fFieldPropagator->ClearPropagatorState(); 815 // Resets all state of field propagator << 794 // Resets all state of field propagator class (ONLY) 816 // values (in case of overlaps and to w << 795 // including safety values (in case of overlaps and to wipe for first track). >> 796 >> 797 // G4ChordFinder* chordF= fFieldPropagator->GetChordFinder(); >> 798 // if( chordF ) chordF->ResetStepEstimate(); 817 } 799 } 818 800 819 // Make sure to clear the chord finders of a << 801 // Make sure to clear the chord finders of all fields (ie managers) 820 // 802 // >> 803 G4FieldManagerStore* fieldMgrStore = G4FieldManagerStore::GetInstance(); 821 fieldMgrStore->ClearAllChordFindersState(); 804 fieldMgrStore->ClearAllChordFindersState(); 822 805 823 // Update the current touchable handle (fro 806 // Update the current touchable handle (from the track's) 824 // 807 // 825 fCurrentTouchableHandle = aTrack->GetTouchab 808 fCurrentTouchableHandle = aTrack->GetTouchableHandle(); 826 809 827 // Inform field propagator of new track 810 // Inform field propagator of new track 828 // << 829 fFieldPropagator->PrepareNewTrack(); 811 fFieldPropagator->PrepareNewTrack(); 830 } 812 } 831 813 832 ////////////////////////////////////////////// << 814 #include "G4CoupledTransportation.hh" 833 // << 815 G4bool G4Transportation::EnableUseMagneticMoment(G4bool useMoment) 834 << 835 G4bool G4Transportation::EnableMagneticMoment( << 836 { 816 { 837 G4bool lastValue= fUseMagneticMoment; 817 G4bool lastValue= fUseMagneticMoment; 838 fUseMagneticMoment= useMoment; 818 fUseMagneticMoment= useMoment; >> 819 G4CoupledTransportation::fUseMagneticMoment= useMoment; 839 return lastValue; 820 return lastValue; 840 } << 841 << 842 ////////////////////////////////////////////// << 843 // << 844 << 845 G4bool G4Transportation::EnableGravity(G4bool << 846 { << 847 G4bool lastValue= fUseGravity; << 848 fUseGravity= useGravity; << 849 return lastValue; << 850 } << 851 << 852 ////////////////////////////////////////////// << 853 // << 854 // Supress (or not) warnings about 'looping' << 855 << 856 void G4Transportation::SetSilenceLooperWarning << 857 { << 858 fSilenceLooperWarnings= val; // Flag to *Su << 859 } << 860 << 861 ////////////////////////////////////////////// << 862 // << 863 G4bool G4Transportation::GetSilenceLooperWarni << 864 { << 865 return fSilenceLooperWarnings; << 866 } << 867 << 868 << 869 ////////////////////////////////////////////// << 870 // << 871 void G4Transportation::SetHighLooperThresholds << 872 { << 873 // Restores the old high values -- potential << 874 // HEP experiments. << 875 // Caution: All tracks with E < 100 MeV tha << 876 SetThresholdWarningEnergy( 100.0 * CLHEP: << 877 SetThresholdImportantEnergy( 250.0 * CLHEP: << 878 << 879 G4int maxTrials = 10; << 880 SetThresholdTrials( maxTrials ); << 881 << 882 PushThresholdsToLogger(); // Again, to be s << 883 if( verboseLevel ) ReportLooperThresholds() << 884 } << 885 << 886 ////////////////////////////////////////////// << 887 void G4Transportation::SetLowLooperThresholds( << 888 { << 889 // These values were the default in Geant4 1 << 890 SetThresholdWarningEnergy( 1.0 * CLHEP:: << 891 SetThresholdImportantEnergy( 1.0 * CLHEP:: << 892 << 893 G4int maxTrials = 30; // A new value - was << 894 SetThresholdTrials( maxTrials ); << 895 << 896 PushThresholdsToLogger(); // Again, to be s << 897 if( verboseLevel ) ReportLooperThresholds() << 898 } << 899 << 900 ////////////////////////////////////////////// << 901 // << 902 void << 903 G4Transportation::ReportMissingLogger( const c << 904 { << 905 const char* message= "Logger object missing << 906 G4String classAndMethod= G4String("G4Transp << 907 G4Exception(classAndMethod, "Missing Logger << 908 } << 909 << 910 << 911 ////////////////////////////////////////////// << 912 // << 913 void << 914 G4Transportation::ReportLooperThresholds() << 915 { << 916 PushThresholdsToLogger(); // To be absolut << 917 fpLogger->ReportLooperThresholds("G4Transpo << 918 } << 919 << 920 ////////////////////////////////////////////// << 921 // << 922 void G4Transportation::ProcessDescription(std: << 923 << 924 // StreamInfo(std::ostream& out, const G4Parti << 925 << 926 { << 927 G4String indent = " "; // : ""); << 928 G4long oldPrec= outStr.precision(6); << 929 // outStr << std::setprecision(6); << 930 outStr << G4endl << indent << GetProcessName << 931 << 932 outStr << " Parameters for looping particl << 933 << " warning-E = " << fThreshold_ << 934 << " important E = " << fThreshol << 935 << " thresholdTrials " << fThresh << 936 outStr.precision(oldPrec); << 937 } 821 } 938 822