Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer << 3 // * DISCLAIMER * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th << 5 // * The following disclaimer summarizes all the specific disclaimers * 6 // * the Geant4 Collaboration. It is provided << 6 // * of contributors to this software. The specific disclaimers,which * 7 // * conditions of the Geant4 Software License << 7 // * govern, are listed with their locations in: * 8 // * LICENSE and available at http://cern.ch/ << 8 // * http://cern.ch/geant4/license * 9 // * include a list of copyright holders. << 10 // * 9 // * * 11 // * Neither the authors of this software syst 10 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 11 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 12 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 13 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file << 14 // * use. * 16 // * for the full disclaimer and the limitatio << 17 // * 15 // * * 18 // * This code implementation is the result << 16 // * This code implementation is the intellectual property of the * 19 // * technical work of the GEANT4 collaboratio << 17 // * GEANT4 collaboration. * 20 // * By using, copying, modifying or distri << 18 // * By copying, distributing or modifying the Program (or any work * 21 // * any work based on the software) you ag << 19 // * based on the Program) you indicate your acceptance of this * 22 // * use in resulting scientific publicati << 20 // * statement, and all its terms. * 23 // * acceptance of all terms of the Geant4 Sof << 24 // ******************************************* 21 // ******************************************************************** 25 // 22 // 26 // 23 // >> 24 // $Id: G4Transportation.cc,v 1.29 2002/05/17 07:10:23 gcosmo Exp $ >> 25 // GEANT4 tag $Name: geant4-04-01 $ 27 // 26 // 28 // ------------------------------------------- 27 // ------------------------------------------------------------ 29 // GEANT 4 include file implementation << 28 // GEANT 4 include file implementation 30 // 29 // 31 // ------------------------------------------- 30 // ------------------------------------------------------------ 32 // 31 // 33 // This class is a process responsible for the 32 // This class is a process responsible for the transportation of 34 // a particle, ie the geometrical propagation 33 // a particle, ie the geometrical propagation that encounters the 35 // geometrical sub-volumes of the detectors. 34 // geometrical sub-volumes of the detectors. 36 // 35 // 37 // It is also tasked with the key role of prop << 36 // It is also tasked with part of updating the "safety". 38 // which will be used to update the post-ste << 39 // 37 // 40 // =========================================== 38 // ======================================================================= >> 39 // Modified: >> 40 // 29 June 2001, J. Apostolakis, D. Cote-Ahern, P. Gumplinger: >> 41 // correction for spin tracking >> 42 // 20 Febr 2001, J. Apostolakis: update for new FieldTrack >> 43 // 22 Sept 2000, V. Grichine: update of Kinetic Energy >> 44 // 9 June 1999, J. Apostolakis & S.Giani: protect full relocation >> 45 // used in DEBUG for track that started on surface >> 46 // and went step < tolerance >> 47 // Also forced fast relocation in all DEBUG cases >> 48 // & changed #if to use DEBUG instead of VERBOSE 41 // Created: 19 March 1997, J. Apostolakis 49 // Created: 19 March 1997, J. Apostolakis 42 // =========================================== 50 // ======================================================================= 43 51 44 #include "G4Transportation.hh" 52 #include "G4Transportation.hh" 45 #include "G4TransportationProcessType.hh" << 46 #include "G4TransportationLogger.hh" << 47 53 48 #include "G4PhysicalConstants.hh" << 54 /////////////////////////////////////////////////////////////////////////////// 49 #include "G4SystemOfUnits.hh" << 50 #include "G4ProductionCutsTable.hh" << 51 #include "G4ParticleTable.hh" << 52 << 53 #include "G4ChargeState.hh" << 54 #include "G4EquationOfMotion.hh" << 55 << 56 #include "G4FieldManagerStore.hh" << 57 << 58 #include "G4Navigator.hh" << 59 #include "G4PropagatorInField.hh" << 60 #include "G4TransportationManager.hh" << 61 << 62 #include "G4TransportationParameters.hh" << 63 << 64 class G4VSensitiveDetector; << 65 << 66 G4bool G4Transportation::fUseMagneticMoment=fa << 67 G4bool G4Transportation::fUseGravity= false; << 68 G4bool G4Transportation::fSilenceLooperWarning << 69 << 70 ////////////////////////////////////////////// << 71 // 55 // 72 // Constructor 56 // Constructor 73 57 74 G4Transportation::G4Transportation( G4int verb << 58 G4Transportation::G4Transportation() : 75 : G4VProcess( aName, fTransportation ), << 59 G4VProcess(G4String("Transportation") ) 76 fFieldExertedForce( false ), << 77 fPreviousSftOrigin( 0.,0.,0. ), << 78 fPreviousSafety( 0.0 ), << 79 fEndPointDistance( -1.0 ), << 80 fShortStepOptimisation( false ) // Old def << 81 { 60 { 82 SetProcessSubType(static_cast<G4int>(TRANSPO << 83 pParticleChange= &fParticleChange; // Requ << 84 SetVerboseLevel(verbosity); << 85 << 86 G4TransportationManager* transportMgr ; 61 G4TransportationManager* transportMgr ; 87 62 88 transportMgr = G4TransportationManager::GetT 63 transportMgr = G4TransportationManager::GetTransportationManager() ; 89 64 90 fLinearNavigator = transportMgr->GetNavigato << 65 fLinearNavigator = transportMgr->GetNavigatorForTracking() ; >> 66 fFieldPropagator = 0 ; 91 67 92 fFieldPropagator = transportMgr->GetPropagat << 68 // fFieldExists= false ; 93 69 94 fpSafetyHelper = transportMgr->GetSafetyHe << 70 fParticleIsLooping = false ; >> 71 >> 72 // fGlobalFieldMgr= transportMgr->GetFieldManager() ; 95 73 96 fpLogger = new G4TransportationLogger("G4Tra << 74 fFieldPropagator= transportMgr->GetPropagatorInField() ; 97 75 98 if( G4TransportationParameters::Exists() ) << 76 // Find out if an electromagnetic field exists 99 { << 77 // 100 auto trParams= G4TransportationParameters: << 78 // fFieldExists= transportMgr->GetFieldManager()->DoesFieldExist() ; 101 << 79 // 102 SetThresholdWarningEnergy( trParams->GetW << 80 // The above code is problematic, because it only works if 103 SetThresholdImportantEnergy( trParams->Get << 81 // the field manager has informed about the detector's field 104 SetThresholdTrials( trParams->GetNumberOfT << 82 // before this transportation process is constructed. 105 G4Transportation::fSilenceLooperWarnings= << 83 // I cannot foresee how the transportation can be informed later. JA 106 } << 84 // The current answer is to ignore this data member and use 107 else { << 85 // the member function DoesGlobalFieldExist() in its place ... 108 SetHighLooperThresholds(); << 86 // John Apostolakis, July 7, 1997 109 // Use the old defaults: Warning = 100 Me << 110 } << 111 87 112 PushThresholdsToLogger(); << 88 fCurrentTouchableHandle = new G4TouchableHistory(); 113 // Should be done by Set methods in SetHighL << 114 89 115 static G4ThreadLocal G4TouchableHandle* pNul << 90 // Initial value for safety and point-of-origin of safety 116 if ( !pNullTouchableHandle) << 117 { << 118 pNullTouchableHandle = new G4TouchableHand << 119 } << 120 fCurrentTouchableHandle = *pNullTouchableHan << 121 // Points to (G4VTouchable*) 0 << 122 91 >> 92 fPreviousSafety = 0.0 ; >> 93 fPreviousSftOrigin = G4ThreeVector(0.,0.,0.) ; 123 94 124 #ifdef G4VERBOSE << 95 fEndGlobalTimeComputed= false; 125 if( verboseLevel > 0) << 96 fCandidateEndGlobalTime= 0; 126 { << 127 G4cout << " G4Transportation constructor> << 128 if ( fShortStepOptimisation ) { G4cout < << 129 else { G4cout < << 130 } << 131 #endif << 132 } 97 } 133 98 134 ////////////////////////////////////////////// << 99 ///////////////////////////////////////////////////////////////////////////// 135 100 136 G4Transportation::~G4Transportation() 101 G4Transportation::~G4Transportation() 137 { 102 { 138 if( fSumEnergyKilled > 0.0 ) << 139 { << 140 PrintStatistics( G4cout ); << 141 } << 142 delete fpLogger; << 143 } << 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 } 103 } 182 104 183 ////////////////////////////////////////////// << 105 ////////////////////////////////////////////////////////////////////////////// 184 // 106 // 185 // Responsibilities: 107 // Responsibilities: 186 // Find whether the geometry limits the Ste 108 // Find whether the geometry limits the Step, and to what length 187 // Calculate the new value of the safety an 109 // Calculate the new value of the safety and return it. 188 // Store the final time, position and momen 110 // Store the final time, position and momentum. 189 111 190 G4double G4Transportation::AlongStepGetPhysica << 112 G4double G4Transportation:: 191 const G4Track& track, << 113 AlongStepGetPhysicalInteractionLength( const G4Track& track, 192 G4double, // previousStepSize << 114 G4double previousStepSize, 193 G4double currentMinimumStep, G4double& curre << 115 G4double currentMinimumStep, 194 G4GPILSelection* selection) << 116 G4double& currentSafety, 195 { << 117 G4GPILSelection* selection ) 196 // Initial actions moved to StartTrack() << 118 { 197 // -------------------------------------- << 119 G4double geometryStepLength, newSafety ; 198 // Note: in case another process changes tou << 120 fParticleIsLooping = false ; 199 // it will be necessary to add here (for << 121 200 // fCurrentTouchableHandle = aTrack->GetTouc << 122 if(track.GetCurrentStepNumber()==1) { >> 123 // reset safety value >> 124 fPreviousSafety = 0.0 ; >> 125 fPreviousSftOrigin = G4ThreeVector(0.,0.,0.) ; >> 126 // ChordFinder reset internal state >> 127 if ( this->DoesGlobalFieldExist() ) >> 128 fFieldPropagator->GetChordFinder()->ResetStepEstimate(); >> 129 // We need to update the current transportation's touchable handle >> 130 // to the track's one >> 131 fCurrentTouchableHandle = track.GetTouchableHandle(); >> 132 } 201 133 202 // GPILSelection is set to defaule value of 134 // GPILSelection is set to defaule value of CandidateForSelection 203 // It is a return value 135 // It is a return value 204 // << 136 205 *selection = CandidateForSelection; << 137 *selection = CandidateForSelection ; 206 138 207 // Get initial Energy/Momentum of the track 139 // Get initial Energy/Momentum of the track 208 // << 140 209 const G4ThreeVector startPosition = track << 141 const G4DynamicParticle* pParticle = track.GetDynamicParticle() ; 210 const G4ThreeVector startMomentumDir = track << 142 G4ThreeVector startMomentumDir = pParticle->GetMomentumDirection() ; 211 << 143 G4ThreeVector startPosition = track.GetPosition() ; 212 // The Step Point safety can be limited by o << 213 // assumptions of any process - it's not alw << 214 // We calculate the starting point's isotrop << 215 { << 216 const G4double MagSqShift = (startPosition << 217 144 218 if(MagSqShift >= sqr(fPreviousSafety)) << 145 // G4double theTime = track.GetGlobalTime() ; 219 currentSafety = 0.0; << 146 220 else << 147 // The Step Point safety is now generalised to mean the limit of assumption 221 currentSafety = fPreviousSafety - std::s << 148 // of all processes, so it is not the previous Step's geometrical safety. >> 149 // We calculate the starting point's safety here. >> 150 >> 151 G4ThreeVector OriginShift = startPosition - fPreviousSftOrigin ; >> 152 G4double MagSqShift = OriginShift.mag2() ; >> 153 if( MagSqShift >= sqr(fPreviousSafety) ) >> 154 { >> 155 currentSafety = 0.0 ; >> 156 } >> 157 else >> 158 { >> 159 currentSafety = fPreviousSafety - sqrt(MagSqShift) ; 222 } 160 } >> 161 // Is the particle charged ? 223 162 224 // Is the particle charged or has it a magne << 163 G4ParticleDefinition* pParticleDef = pParticle->GetDefinition() ; 225 // << 164 G4double particleCharge = pParticleDef->GetPDGCharge() ; 226 const G4DynamicParticle* pParticle = track.G << 227 << 228 const G4double particleMass = pParticle->G << 229 const G4double particleCharge = pParticle->G << 230 const G4double kineticEnergy = pParticle->Ge << 231 165 232 const G4double magneticMoment = pParticle << 166 G4bool fieldExertsForce = false ; 233 const G4ThreeVector particleSpin = pParticle << 167 fGeometryLimitedStep = false ; >> 168 // fEndGlobalTimeComputed = false ; 234 169 235 // There is no need to locate the current vo 170 // There is no need to locate the current volume. It is Done elsewhere: 236 // On track construction << 171 // On track construction 237 // By the tracking, after all AlongStepDoI 172 // By the tracking, after all AlongStepDoIts, in "Relocation" >> 173 // Does the particle have an (EM) field force exerting upon it? >> 174 >> 175 if( (particleCharge != 0.0) ) >> 176 { >> 177 >> 178 fieldExertsForce= this->DoesGlobalFieldExist() ; 238 179 239 // Check if the particle has a force, EM or << 180 // Future: will/can also check whether current volume's field is Zero or 240 // << 181 // set by the user (in the logical volume) to be zero. >> 182 } >> 183 // Choose the calculation of the transportation: Field or not >> 184 >> 185 if( !fieldExertsForce ) >> 186 { >> 187 G4double linearStepLength ; >> 188 >> 189 if( currentMinimumStep <= currentSafety ) >> 190 { >> 191 // The Step is guaranteed to be taken >> 192 >> 193 geometryStepLength = currentMinimumStep ; >> 194 fGeometryLimitedStep = false ; >> 195 } >> 196 else >> 197 { >> 198 // Find whether the straight path intersects a volume >> 199 >> 200 linearStepLength = fLinearNavigator->ComputeStep( startPosition, >> 201 startMomentumDir, >> 202 currentMinimumStep, >> 203 newSafety) ; >> 204 // Remember last safety origin & value. >> 205 >> 206 fPreviousSftOrigin = startPosition ; >> 207 fPreviousSafety = newSafety ; >> 208 >> 209 // The safety at the initial point has been re-calculated: >> 210 >> 211 currentSafety = newSafety ; >> 212 >> 213 if( linearStepLength <= currentMinimumStep) >> 214 { >> 215 // The geometry limits the Step size (an intersection was found.) >> 216 >> 217 geometryStepLength = linearStepLength ; >> 218 fGeometryLimitedStep = true ; >> 219 } >> 220 else >> 221 { >> 222 // The full Step is taken. >> 223 >> 224 geometryStepLength = currentMinimumStep ; >> 225 fGeometryLimitedStep = false ; >> 226 } >> 227 } >> 228 endpointDistance = geometryStepLength ; >> 229 >> 230 // Calculate final position >> 231 >> 232 fTransportEndPosition = startPosition + geometryStepLength*startMomentumDir ; >> 233 >> 234 // Momentum direction, energy and polarisation are unchanged by transport >> 235 >> 236 fTransportEndMomentumDir = startMomentumDir ; >> 237 fTransportEndKineticEnergy = track.GetKineticEnergy() ; >> 238 fTransportEndSpin = track.GetPolarization(); >> 239 fParticleIsLooping = false ; >> 240 fMomentumChanged = false ; >> 241 fEndGlobalTimeComputed = false ; >> 242 } >> 243 else >> 244 { >> 245 G4double momentumMagnitude = pParticle->GetTotalMomentum() ; >> 246 G4ThreeVector EndUnitMomentum ; >> 247 G4double lengthAlongCurve ; >> 248 G4double restMass = pParticleDef->GetPDGMass() ; >> 249 >> 250 fFieldPropagator->SetChargeMomentumMass( particleCharge, // charge in e+ units >> 251 momentumMagnitude, // Momentum in Mev/c >> 252 restMass ) ; >> 253 >> 254 G4ThreeVector spin = track.GetPolarization() ; >> 255 G4FieldTrack aFieldTrack = >> 256 G4FieldTrack( startPosition, >> 257 track.GetMomentumDirection(), >> 258 0.0, >> 259 track.GetKineticEnergy(), >> 260 restMass, >> 261 track.GetVelocity(), >> 262 track.GetGlobalTime(), // Laboratory fr. >> 263 track.GetProperTime(), // Particle rest fr. >> 264 &spin ) ; >> 265 >> 266 if( currentMinimumStep > 0 ) >> 267 { >> 268 // Do the Transport in the field (non recti-linear) >> 269 >> 270 lengthAlongCurve = fFieldPropagator->ComputeStep( aFieldTrack, >> 271 currentMinimumStep, >> 272 currentSafety, >> 273 track.GetVolume() ) ; >> 274 if( lengthAlongCurve < currentMinimumStep) >> 275 { >> 276 geometryStepLength = lengthAlongCurve ; >> 277 fGeometryLimitedStep = true ; >> 278 } >> 279 else >> 280 { >> 281 geometryStepLength = currentMinimumStep ; >> 282 fGeometryLimitedStep = false ; >> 283 } >> 284 } >> 285 else >> 286 { >> 287 geometryStepLength = lengthAlongCurve= 0.0 ; >> 288 fGeometryLimitedStep = false ; >> 289 } >> 290 // Remember last safety origin & value. 241 291 242 G4bool eligibleEM = << 292 fPreviousSftOrigin = startPosition ; 243 (particleCharge != 0.0) || ((magneticMomen << 293 fPreviousSafety = currentSafety ; 244 G4bool eligibleGrav = (particleMass != 0.0) << 294 >> 295 // Get the End-Position and End-Momentum (Dir-ection) 245 296 246 fFieldExertedForce = false; << 297 fTransportEndPosition = aFieldTrack.GetPosition() ; 247 298 248 if(eligibleEM || eligibleGrav) << 299 // Momentum: Magnitude and direction can be changed too now ... 249 { << 250 if(G4FieldManager* fieldMgr = << 251 fFieldPropagator->FindAndSetFieldMana << 252 { << 253 // User can configure the field Manager << 254 fieldMgr->ConfigureForTrack(&track); << 255 // Called here to allow a transition fro << 256 // to finite field (non-zero pointer). << 257 << 258 // If the field manager has no field ptr << 259 // by definition ( = there is no field << 260 if(const G4Field* ptrField = fieldMgr->G << 261 fFieldExertedForce = << 262 eligibleEM || (eligibleGrav && ptrFi << 263 } << 264 } << 265 300 266 G4double geometryStepLength = currentMinimum << 301 fMomentumChanged = true ; >> 302 fTransportEndMomentumDir = aFieldTrack.GetMomentumDir() ; 267 303 268 if(currentMinimumStep == 0.0) << 304 fTransportEndKineticEnergy = aFieldTrack.GetKineticEnergy() ; 269 { << 270 fEndPointDistance = 0.0; << 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 } << 281 else if(!fFieldExertedForce) << 282 { << 283 fGeometryLimitedStep = false; << 284 if(geometryStepLength > currentSafety || ! << 285 { << 286 const G4double linearStepLength = fLinea << 287 startPosition, startMomentumDir, curre << 288 305 289 if(linearStepLength <= currentMinimumSte << 306 // if( (track.GetKineticEnergy() - fTransportEndKineticEnergy) 290 { << 307 // > perMillion * fTransportEndKineticEnergy ){ 291 geometryStepLength = linearStepLength; << 292 fGeometryLimitedStep = true; << 293 } << 294 // Remember last safety origin & value. << 295 // << 296 fPreviousSftOrigin = startPosition; << 297 fPreviousSafety = currentSafety; << 298 fpSafetyHelper->SetCurrentSafety(current << 299 } << 300 308 301 fEndPointDistance = geometryStepLength; << 309 if( fFieldPropagator->GetCurrentFieldManager()->DoesFieldChangeEnergy() ){ 302 310 303 fMomentumChanged = false; << 311 // If the field can changed energy, then the time must be integrated 304 fParticleIsLooping = false; << 312 // - so this should have been updated 305 fEndGlobalTimeComputed = false; << 313 fCandidateEndGlobalTime = aFieldTrack.GetLabTimeOfFlight(); 306 fTransportEndPosition = << 314 fEndGlobalTimeComputed = true; 307 startPosition + geometryStepLength * sta << 315 // was ( fCandidateEndGlobalTime != track.GetGlobalTime() ); 308 fTransportEndMomentumDir = startMomentum << 316 // a cleaner way is to have FieldTrack knowing whether time is updated. 309 fTransportEndKineticEnergy = kineticEnergy << 317 }else{ 310 fTransportEndSpin = particleSpin; << 318 // The energy is unchanged by field transport, 311 } << 319 // so the time changed will be calculated elsewhere 312 else // A field exerts force << 320 fEndGlobalTimeComputed = false; 313 { << 321 } 314 const auto pParticleDef = pParticle->Ge << 315 const auto particlePDGSpin = pParticleDef- << 316 const auto particlePDGMagM = pParticleDef- << 317 << 318 auto equationOfMotion = fFieldPropagator-> << 319 << 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 << 348 fGeometryLimitedStep = fFieldPropagator->I << 349 // << 350 // It is possible that step was reduced in << 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 << 360 fEndGlobalTimeComputed = changesEnergy; << 361 fTransportEndPosition = aFieldTrack.Get << 362 fTransportEndMomentumDir = aFieldTrack.Get << 363 << 364 // G4cout << " G4Transport: End of step pM << 365 << 366 fEndPointDistance = (fTransportEndPositio << 367 << 368 // Ignore change in energy for fields that << 369 // This hides the integration error, but g << 370 fTransportEndKineticEnergy = << 371 changesEnergy ? aFieldTrack.GetKineticEn << 372 fTransportEndSpin = aFieldTrack.GetSpin(); << 373 322 374 if(fEndGlobalTimeComputed) << 323 fTransportEndSpin = aFieldTrack.GetSpin(); 375 { << 376 // If the field can change energy, then << 377 // - so this should have been updated << 378 // << 379 fCandidateEndGlobalTime = aFieldTrack.Ge << 380 324 381 // was ( fCandidateEndGlobalTime != trac << 325 fParticleIsLooping = fFieldPropagator->IsParticleLooping() ; 382 // a cleaner way is to have FieldTrack k << 326 endpointDistance = (fTransportEndPosition - startPosition).mag() ; 383 } << 384 #if defined(G4VERBOSE) || defined(G4DEBUG_TRAN << 385 else << 386 { << 387 // The energy should be unchanged by fie << 388 // - so the time changed will be calc << 389 // << 390 // Check that the integration preserved << 391 // - and if not correct this! << 392 G4double startEnergy = kineticEnergy; << 393 G4double endEnergy = fTransportEndKine << 394 << 395 static G4ThreadLocal G4int no_large_edif << 396 if(verboseLevel > 1) << 397 { << 398 if(std::fabs(startEnergy - endEnergy) << 399 { << 400 static G4ThreadLocal G4int no_warnin << 401 moduloFac << 402 no_large_ediff++; << 403 if((no_large_ediff % warnModulo) == << 404 { << 405 no_warnings++; << 406 std::ostringstream message; << 407 message << "Energy change in Step << 408 << G4endl << " Relativ << 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 { << 435 warnModulo *= moduloFactor; << 436 } << 437 } << 438 } << 439 } // end of if (verboseLevel) << 440 } << 441 #endif << 442 } 327 } >> 328 // If we are asked to go a step length of 0, and we are on a boundary >> 329 // then a boundary will also limit the step -> we must flag this. 443 330 444 // Update the safety starting from the end-p << 331 if (currentMinimumStep == 0.0 ) 445 // if it will become negative at the end-poi << 446 // << 447 if(currentSafety < fEndPointDistance) << 448 { 332 { 449 if(particleCharge != 0.0) << 333 if( currentSafety == 0.0 ) fGeometryLimitedStep = true ; 450 { << 334 } 451 G4double endSafety = << 335 452 fLinearNavigator->ComputeSafety(fTrans << 336 // Update the safety starting from the end-point, if it will become 453 currentSafety = endSafety; << 337 // negative at the end-point. 454 fPreviousSftOrigin = fTransportEndPositi << 338 455 fPreviousSafety = currentSafety; << 339 if( currentSafety < endpointDistance ) 456 fpSafetyHelper->SetCurrentSafety(current << 340 { >> 341 G4double endSafety = fLinearNavigator->ComputeSafety( fTransportEndPosition) ; >> 342 currentSafety = endSafety ; >> 343 fPreviousSftOrigin = fTransportEndPosition ; >> 344 fPreviousSafety = currentSafety ; 457 345 458 // Because the Stepping Manager assumes << 346 // Because the Stepping Manager assumes it is from the start point, 459 // add the StepLength 347 // add the StepLength 460 // << 461 currentSafety += fEndPointDistance; << 462 348 463 #ifdef G4DEBUG_TRANSPORT << 349 currentSafety += endpointDistance ; 464 G4cout.precision(12); << 350 465 G4cout << "***G4Transportation::AlongSte << 351 #ifdef G4DEBUG_TRANSPORT >> 352 G4cout.precision(16) ; >> 353 G4cout << "***Transportation::AlongStepGPIL ** " << G4endl ; 466 G4cout << " Called Navigator->ComputeSa 354 G4cout << " Called Navigator->ComputeSafety at " << fTransportEndPosition 467 << " and it returned safety= << 355 << " and it returned safety= " << endSafety << G4endl ; 468 G4cout << " Adding endpoint distance << 356 G4cout << " Adding endpoint distance " << endpointDistance 469 << " to obtain pseudo-safety= << 357 << " we obtain pseudo-safety= " << currentSafety << G4endl ; 470 } << 471 else << 472 { << 473 G4cout << "***G4Transportation::AlongSte << 474 G4cout << " Avoiding call to ComputeSaf << 475 G4cout << " charge = " << particl << 476 G4cout << " mag moment = " << magneti << 477 #endif 358 #endif 478 } << 359 } 479 } << 480 << 481 fFirstStepInVolume = fNewTrack || fLastStepI << 482 fLastStepInVolume = false; << 483 fNewTrack = false; << 484 360 485 fParticleChange.ProposeFirstStepInVolume(fFi << 361 fParticleChange.SetTrueStepLength(geometryStepLength) ; 486 fParticleChange.ProposeTrueStepLength(geomet << 487 362 488 return geometryStepLength; << 363 return geometryStepLength ; 489 } 364 } 490 365 491 ////////////////////////////////////////////// << 366 ///////////////////////////////////////////////////////////////////////////// 492 // 367 // 493 // Initialize ParticleChange (by setting al 368 // Initialize ParticleChange (by setting all its members equal 494 // to correspond 369 // to corresponding members in G4Track) 495 370 496 G4VParticleChange* G4Transportation::AlongStep 371 G4VParticleChange* G4Transportation::AlongStepDoIt( const G4Track& track, 497 372 const G4Step& stepData ) 498 { 373 { 499 #if defined(G4VERBOSE) || defined(G4DEBUG_TRAN << 500 static G4ThreadLocal G4long noCallsASDI=0; << 501 noCallsASDI++; << 502 #else << 503 #define noCallsASDI 0 << 504 #endif << 505 << 506 if(fGeometryLimitedStep) << 507 { << 508 stepData.GetPostStepPoint()->SetStepStatus << 509 } << 510 << 511 fParticleChange.Initialize(track) ; 374 fParticleChange.Initialize(track) ; 512 375 513 // Code for specific process 376 // Code for specific process 514 // << 377 515 fParticleChange.ProposePosition(fTransportEn << 378 fParticleChange.SetPositionChange(fTransportEndPosition) ; 516 fParticleChange.ProposeMomentumDirection(fTr << 379 fParticleChange.SetMomentumChange(fTransportEndMomentumDir) ; 517 fParticleChange.ProposeEnergy(fTransportEndK << 380 fParticleChange.SetEnergyChange(fTransportEndKineticEnergy) ; 518 fParticleChange.SetMomentumChanged(fMomentum 381 fParticleChange.SetMomentumChanged(fMomentumChanged) ; 519 382 520 fParticleChange.ProposePolarization(fTranspo << 383 fParticleChange.SetPolarizationChange(fTransportEndSpin); 521 << 384 522 G4double deltaTime = 0.0 ; 385 G4double deltaTime = 0.0 ; 523 386 524 // Calculate Lab Time of Flight (ONLY if fi 387 // Calculate Lab Time of Flight (ONLY if field Equations used it!) 525 // G4double endTime = fCandidateEndGlobalT << 388 // G4double endTime = fCandidateEndGlobalTime; 526 // G4double delta_time = endTime - startTime << 389 // G4double delta_time = endTime - startTime; 527 << 390 G4double startTime = track.GetGlobalTime(); 528 G4double startTime = track.GetGlobalTime() ; << 529 391 530 if (!fEndGlobalTimeComputed) << 392 if (!fEndGlobalTimeComputed){ 531 { << 393 // The time was not integrated .. make the best estimate possible 532 // The time was not integrated .. make th << 394 G4double finalVelocity = track.GetVelocity() ; 533 // << 534 G4double initialVelocity = stepData.GetPr 395 G4double initialVelocity = stepData.GetPreStepPoint()->GetVelocity(); 535 G4double stepLength = track.GetStepL << 396 G4double stepLength = track.GetStepLength(); 536 397 537 deltaTime= 0.0; // in case initialVeloci << 398 if (finalVelocity > 0.0) { 538 if ( initialVelocity > 0.0 ) { deltaTime << 399 G4double meanInverseVelocity ; >> 400 // deltaTime = stepLength/finalVelocity ; >> 401 meanInverseVelocity= 0.5 * ( 1.0 / initialVelocity + 1.0 / finalVelocity ); >> 402 deltaTime = stepLength * meanInverseVelocity ; >> 403 }else{ >> 404 deltaTime = stepLength/initialVelocity; >> 405 } >> 406 fCandidateEndGlobalTime = startTime + deltaTime; >> 407 }else >> 408 deltaTime = fCandidateEndGlobalTime - startTime; 539 409 540 fCandidateEndGlobalTime = startTime + d << 410 fParticleChange. SetTimeChange( fCandidateEndGlobalTime ) ; 541 fParticleChange.ProposeLocalTime( track. << 542 } << 543 else << 544 { << 545 deltaTime = fCandidateEndGlobalTime - sta << 546 fParticleChange.ProposeGlobalTime( fCandi << 547 } << 548 411 549 << 412 // Now Correct by Lorentz factor to get "proper" deltaTime 550 // Now Correct by Lorentz factor to get delt << 551 413 552 G4double restMass = track.GetDynamicP 414 G4double restMass = track.GetDynamicParticle()->GetMass() ; 553 G4double deltaProperTime = deltaTime*( restM 415 G4double deltaProperTime = deltaTime*( restMass/track.GetTotalEnergy() ) ; 554 416 555 fParticleChange.ProposeProperTime(track.GetP << 417 fParticleChange.SetProperTimeChange( track.GetProperTime() + deltaProperTime ) ; 556 //fParticleChange.ProposeTrueStepLength( tra << 418 //fParticleChange. SetTrueStepLength( track.GetStepLength() ) ; 557 419 558 // If the particle is caught looping or is s << 420 // If the particle is caught looping or is stuck (in very difficult boundaries) 559 // boundaries) in a magnetic field (doing ma << 421 // in a magnetic field (doing many steps) 560 // << 422 // THEN this kills it ... 561 if ( fParticleIsLooping ) 423 if ( fParticleIsLooping ) 562 { 424 { 563 G4double endEnergy= fTransportEndKinetic << 425 // Kill the looping particle 564 fNoLooperTrials ++; << 426 565 auto particleType= track.GetDynamicParti << 427 fParticleChange.SetStatusChange( fStopAndKill ) ; 566 << 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 { << 575 // Kill the looping particle << 576 // << 577 fParticleChange.ProposeTrackStatus( fS << 578 G4int particlePDG= particleType->GetPD << 579 const G4int electronPDG= 11; // G4Elec << 580 << 581 // Simple statistics << 582 fSumEnergyKilled += endEnergy; << 583 fSumEnerSqKilled += endEnergy * endEne << 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 << 603 if( endEnergy > fThreshold_Warning_Ene << 604 { << 605 fpLogger->ReportLoopingTrack( track, << 606 noCall << 607 } << 608 fNoLooperTrials=0; << 609 } << 610 else << 611 { << 612 fMaxEnergySaved = std::max( endEnergy, << 613 if( fNoLooperTrials == 1 ) { << 614 fSumEnergySaved += endEnergy; << 615 if ( !stable ) << 616 fSumEnergyUnstableSaved += endEne << 617 } << 618 #ifdef G4VERBOSE 428 #ifdef G4VERBOSE 619 if( verboseLevel > 2 && ! fSilenceLoop << 429 G4cout << " G4Transportation is killing track that is looping or stuck " << G4endl 620 { << 430 << " This track has " << track.GetKineticEnergy() << " MeV energy." 621 G4cout << " " << __func__ << 431 << G4endl; 622 << " Particle is looping but << 623 << " Number of trials = " < << 624 << " No of calls to = " << << 625 } << 626 #endif 432 #endif 627 } << 433 // ClearNumberOfInteractionLengthLeft() ; 628 } 434 } 629 else << 630 { << 631 fNoLooperTrials=0; << 632 } << 633 << 634 // Another (sometimes better way) is to use 435 // Another (sometimes better way) is to use a user-limit maximum Step size 635 // to alleviate this problem .. << 436 // to alleviate this problem .. 636 << 637 // Introduce smooth curved trajectories to p << 638 // << 639 fParticleChange.SetPointerToVectorOfAuxiliar << 640 (fFieldPropagator->GimmeTrajectoryVectorAn << 641 437 642 return &fParticleChange ; 438 return &fParticleChange ; >> 439 643 } 440 } 644 441 645 ////////////////////////////////////////////// << 442 //////////////////////////////////////////////////////////////////////////////// 646 // 443 // 647 // This ensures that the PostStep action is a << 444 // This ensures that the PostStep action is always called, 648 // so that it can do the relocation if it is << 445 // so that it can do the relocation if it is needed. 649 // 446 // 650 447 651 G4double G4Transportation:: 448 G4double G4Transportation:: 652 PostStepGetPhysicalInteractionLength( const G4 449 PostStepGetPhysicalInteractionLength( const G4Track&, 653 G4 << 450 G4double previousStepSize, 654 G4 451 G4ForceCondition* pForceCond ) 655 { << 452 { 656 fFieldExertedForce = false; // Not known << 657 *pForceCond = Forced ; 453 *pForceCond = Forced ; 658 return DBL_MAX ; // was kInfinity ; but con << 659 } << 660 << 661 ////////////////////////////////////////////// << 662 454 663 void G4Transportation::SetTouchableInformation << 455 return DBL_MAX ; // was kInfinity ; but convention now is DBL_MAX 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 } 456 } 704 457 705 ////////////////////////////////////////////// 458 ///////////////////////////////////////////////////////////////////////////// 706 // 459 // 707 460 708 G4VParticleChange* G4Transportation::PostStepD 461 G4VParticleChange* G4Transportation::PostStepDoIt( const G4Track& track, 709 << 462 const G4Step& stepData ) 710 { 463 { 711 G4TouchableHandle retCurrentTouchable ; / << 464 G4TouchableHandle retCurrentTouchable ; // The one to return 712 G4bool isLastStep= false; << 713 465 714 // Initialize ParticleChange (by setting al << 466 // Initialize ParticleChange (by setting all its members equal 715 // to correspond << 467 // to corresponding members in G4Track) 716 // fParticleChange.Initialize(track) ; // T 468 // fParticleChange.Initialize(track) ; // To initialise TouchableChange 717 469 718 fParticleChange.ProposeTrackStatus(track.Get << 470 fParticleChange.SetStatusChange(track.GetTrackStatus()) ; 719 471 720 // If the Step was determined by the volume 472 // If the Step was determined by the volume boundary, 721 // logically relocate the particle 473 // logically relocate the particle 722 << 474 723 if(fGeometryLimitedStep) << 475 if( fGeometryLimitedStep ) 724 { 476 { 725 // fCurrentTouchable will now become the p 477 // fCurrentTouchable will now become the previous touchable, 726 // and what was the previous will be freed << 478 // and what was the previous will be freed. 727 // (Needed because the preStepPoint can po 479 // (Needed because the preStepPoint can point to the previous touchable) 728 480 729 fLinearNavigator->SetGeometricallyLimitedS 481 fLinearNavigator->SetGeometricallyLimitedStep() ; 730 fLinearNavigator-> 482 fLinearNavigator-> 731 LocateGlobalPointAndUpdateTouchableHandle( 483 LocateGlobalPointAndUpdateTouchableHandle( track.GetPosition(), 732 484 track.GetMomentumDirection(), 733 485 fCurrentTouchableHandle, 734 486 true ) ; >> 487 735 // Check whether the particle is out of th 488 // Check whether the particle is out of the world volume 736 // If so it has exited and must be killed. << 489 // If so it has exited and must be killed. 737 // << 490 738 if( fCurrentTouchableHandle->GetVolume() = 491 if( fCurrentTouchableHandle->GetVolume() == 0 ) 739 { 492 { 740 fParticleChange.ProposeTrackStatus( fSt << 493 fParticleChange.SetStatusChange( fStopAndKill ) ; 741 } 494 } 742 retCurrentTouchable = fCurrentTouchableHan << 495 retCurrentTouchable = fCurrentTouchableHandle; 743 fParticleChange.SetTouchableHandle( fCurre 496 fParticleChange.SetTouchableHandle( fCurrentTouchableHandle ) ; 744 << 745 // Update the Step flag which identifies t << 746 if( !fFieldExertedForce ) << 747 isLastStep = fLinearNavigator->ExitedM << 748 || fLinearNavigator->Entered << 749 else << 750 isLastStep = fFieldPropagator->IsLastSt << 751 } << 752 else // fGeometryLimitedStep << 753 { << 754 // This serves only to move the Navigator' << 755 // << 756 fLinearNavigator->LocateGlobalPointWithinV << 757 << 758 // The value of the track's current Toucha << 759 // (and it must be correct because we must << 760 // overwrite the (unset) one in particle c << 761 // It must be fCurrentTouchable too ?? << 762 // << 763 fParticleChange.SetTouchableHandle( track. << 764 retCurrentTouchable = track.GetTouchableHa << 765 << 766 isLastStep= false; << 767 } // endif ( fGeometryLimitedStep ) << 768 fLastStepInVolume= isLastStep; << 769 << 770 fParticleChange.ProposeFirstStepInVolume(fFi << 771 fParticleChange.ProposeLastStepInVolume(isLa << 772 << 773 SetTouchableInformation(retCurrentTouchable) << 774 << 775 return &fParticleChange ; << 776 } << 777 << 778 ////////////////////////////////////////////// << 779 // New method takes over the responsibility to << 780 // G4Transportation object at the start of a n << 781 // of a suspended track. << 782 // << 783 << 784 void << 785 G4Transportation::StartTracking(G4Track* aTrac << 786 { << 787 G4VProcess::StartTracking(aTrack); << 788 fNewTrack= true; << 789 fFirstStepInVolume= true; << 790 fLastStepInVolume= false; << 791 << 792 // The actions here are those that were take << 793 // when track.GetCurrentStepNumber()==1 << 794 << 795 // Whether field exists should be determined << 796 G4FieldManagerStore* fieldMgrStore= G4FieldM << 797 fAnyFieldExists = fieldMgrStore->size() > 0; << 798 << 799 // reset safety value and center << 800 // << 801 fPreviousSafety = 0.0 ; << 802 fPreviousSftOrigin = G4ThreeVector(0.,0.,0.) << 803 << 804 // reset looping counter -- for motion in fi << 805 fNoLooperTrials= 0; << 806 // Must clear this state .. else it depends << 807 // --> a better solution would set this fro << 808 // Was if( aTrack->GetCurrentStepNumber()==1 << 809 << 810 // ChordFinder reset internal state << 811 // << 812 if( fFieldPropagator && fAnyFieldExists ) << 813 { << 814 fFieldPropagator->ClearPropagatorState(); << 815 // Resets all state of field propagator << 816 // values (in case of overlaps and to w << 817 } 497 } >> 498 else >> 499 { // fGeometryLimitedStep is false >> 500 #ifdef G4DEBUG_POSTSTEP_TRANSPORT >> 501 // Although the location is changed, we know that the physical >> 502 // volume remains constant. >> 503 // In order to help in checking the user geometry >> 504 // we perform a full-relocation and check its result >> 505 // *except* if we have made a very small step from a boundary >> 506 // (ie remaining inside the tolerance >> 507 >> 508 G4bool startAtSurface_And_MoveEpsilon ; >> 509 startAtSurface_And_MoveEpsilon = >> 510 (stepData.GetPreStepPoint()->GetSafety() == 0.0) && >> 511 (stepData.GetStepLength() < kCarTolerance ) ; 818 512 819 // Make sure to clear the chord finders of a << 513 if( startAtSurface_And_MoveEpsilon) 820 // << 514 { 821 fieldMgrStore->ClearAllChordFindersState(); << 822 << 823 // Update the current touchable handle (fro << 824 // << 825 fCurrentTouchableHandle = aTrack->GetTouchab << 826 << 827 // Inform field propagator of new track << 828 // << 829 fFieldPropagator->PrepareNewTrack(); << 830 } << 831 515 832 ////////////////////////////////////////////// << 516 fLinearNavigator-> 833 // << 517 LocateGlobalPointAndUpdateTouchableHandle( track.GetPosition(), >> 518 track.GetMomentumDirection(), >> 519 fCurrentTouchableHandle, >> 520 true ); >> 521 if( fCurrentTouchableHandle->GetVolume() != track.GetVolume() ) >> 522 { >> 523 G4cerr << " ERROR: A relocation within safety has caused a volume change! " << G4endl ; >> 524 G4cerr << " The old volume is called " >> 525 << track.GetVolume()->GetName() << G4endl ; >> 526 G4cerr << " The new volume is called " ; >> 527 >> 528 if ( fCurrentTouchableHandle->GetVolume() != 0 ) >> 529 { >> 530 G4cerr << fCurrentTouchableHandle->GetVolume()->GetName() << G4endl ; >> 531 } >> 532 else >> 533 { >> 534 G4cerr << "Out of World" << G4endl ; >> 535 } >> 536 G4cerr.precision(7) ; >> 537 G4cerr << " The position is " << track.GetPosition() << G4endl ; >> 538 >> 539 // Let us relocate again, for debuging >> 540 >> 541 fLinearNavigator-> >> 542 LocateGlobalPointAndUpdateTouchableHandle( track.GetPosition(), >> 543 track.GetMomentumDirection(), >> 544 fCurrentTouchableHandle, >> 545 true ) ; >> 546 G4cerr << " The newer volume is called " ; >> 547 >> 548 if ( fCurrentTouchableHandle->GetVolume() != 0 ) >> 549 { >> 550 G4cerr << fCurrentTouchableHandle->GetVolume()->GetName() << G4endl ; >> 551 } >> 552 else >> 553 { >> 554 G4cerr << "Out of World" << G4endl ; >> 555 } >> 556 } >> 557 >> 558 assert( fCurrentTouchableHandle->GetVolume()->GetName() == >> 559 track.GetVolume()->GetName() ) ; >> 560 >> 561 retCurrentTouchable = fCurrentTouchableHandle ; >> 562 fParticleChange.SetTouchableHandle( fCurrentTouchableHandle ) ; >> 563 >> 564 } >> 565 else >> 566 { >> 567 retCurrentTouchable = track.GetTouchableHandle() ; >> 568 fParticleChange.SetTouchableHandle( track.GetTouchableHandle() ) ; >> 569 } >> 570 // This must be done in the above if ( AtSur ) fails >> 571 // We also do it for if (true) in order to get debug/opt to >> 572 // behave as exactly the same way as possible. 834 573 835 G4bool G4Transportation::EnableMagneticMoment( << 574 fLinearNavigator->LocateGlobalPointWithinVolume( track.GetPosition()) ; 836 { << 575 #else 837 G4bool lastValue= fUseMagneticMoment; << 576 // ie #ifndef G4DEBUG_POSTSTEP_TRANSPORT does a quick relocation 838 fUseMagneticMoment= useMoment; << 577 839 return lastValue; << 578 // The serves only to move the Navigator's location 840 } << 841 579 842 ////////////////////////////////////////////// << 580 fLinearNavigator->LocateGlobalPointWithinVolume( track.GetPosition()) ; 843 // << 844 581 845 G4bool G4Transportation::EnableGravity(G4bool << 582 // The value of the track's current Touchable is retained. 846 { << 583 // (and it must be correct because we must use it below to 847 G4bool lastValue= fUseGravity; << 584 // overwrite the (unset) one in particle change) 848 fUseGravity= useGravity; << 585 // Although in general this is fCurrentTouchable, at the start of 849 return lastValue; << 586 // a step it could be different ... ?? 850 } << 851 587 852 ////////////////////////////////////////////// << 588 fParticleChange.SetTouchableHandle( track.GetTouchableHandle() ) ; 853 // << 589 retCurrentTouchable = track.GetTouchableHandle() ; 854 // Supress (or not) warnings about 'looping' << 590 #endif 855 591 856 void G4Transportation::SetSilenceLooperWarning << 592 } // endif ( fGeometryLimitedStep ) 857 { << 858 fSilenceLooperWarnings= val; // Flag to *Su << 859 } << 860 593 861 ////////////////////////////////////////////// << 594 const G4VPhysicalVolume* pNewVol = retCurrentTouchable->GetVolume() ; 862 // << 595 const G4Material* pNewMaterial = 0 ; 863 G4bool G4Transportation::GetSilenceLooperWarni << 596 864 { << 597 if( pNewVol != 0 ) pNewMaterial= pNewVol->GetLogicalVolume()->GetMaterial() ; 865 return fSilenceLooperWarnings; << 866 } << 867 598 >> 599 // ( <const_cast> pNewMaterial ) ; 868 600 869 ////////////////////////////////////////////// << 601 fParticleChange.SetMaterialChange( (G4Material *) pNewMaterial ) ; 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 602 879 G4int maxTrials = 10; << 603 // temporarily until Get/Set Material of ParticleChange, 880 SetThresholdTrials( maxTrials ); << 604 // and StepPoint can be made const. >> 605 // Set the touchable in ParticleChange >> 606 // this must always be done because the particle change always >> 607 // uses this value to overwrite the current touchable pointer. 881 608 882 PushThresholdsToLogger(); // Again, to be s << 609 fParticleChange.SetTouchableHandle(retCurrentTouchable) ; 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 610 893 G4int maxTrials = 30; // A new value - was << 611 return &fParticleChange ; 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 } 612 } 919 613 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 } << 938 614