[819] | 1 | // |
---|
| 2 | // ******************************************************************** |
---|
| 3 | // * License and Disclaimer * |
---|
| 4 | // * * |
---|
| 5 | // * The Geant4 software is copyright of the Copyright Holders of * |
---|
| 6 | // * the Geant4 Collaboration. It is provided under the terms and * |
---|
| 7 | // * conditions of the Geant4 Software License, included in the file * |
---|
| 8 | // * LICENSE and available at http://cern.ch/geant4/license . These * |
---|
| 9 | // * include a list of copyright holders. * |
---|
| 10 | // * * |
---|
| 11 | // * Neither the authors of this software system, nor their employing * |
---|
| 12 | // * institutes,nor the agencies providing financial support for this * |
---|
| 13 | // * work make any representation or warranty, express or implied, * |
---|
| 14 | // * regarding this software system or assume any liability for its * |
---|
| 15 | // * use. Please see the license in the file LICENSE and URL above * |
---|
| 16 | // * for the full disclaimer and the limitation of liability. * |
---|
| 17 | // * * |
---|
| 18 | // * This code implementation is the result of the scientific and * |
---|
| 19 | // * technical work of the GEANT4 collaboration. * |
---|
| 20 | // * By using, copying, modifying or distributing the software (or * |
---|
| 21 | // * any work based on the software) you agree to acknowledge its * |
---|
| 22 | // * use in resulting scientific publications, and indicate your * |
---|
| 23 | // * acceptance of all terms of the Geant4 Software license. * |
---|
| 24 | // ******************************************************************** |
---|
| 25 | // |
---|
| 26 | // |
---|
| 27 | // $Id: G4Transportation.cc,v 1.72.2.1 2007/12/07 17:18:11 japost Exp $ |
---|
| 28 | // GEANT4 tag $Name: geant4-09-01-patch-02 $ |
---|
| 29 | // |
---|
| 30 | // ------------------------------------------------------------ |
---|
| 31 | // GEANT 4 include file implementation |
---|
| 32 | // |
---|
| 33 | // ------------------------------------------------------------ |
---|
| 34 | // |
---|
| 35 | // This class is a process responsible for the transportation of |
---|
| 36 | // a particle, ie the geometrical propagation that encounters the |
---|
| 37 | // geometrical sub-volumes of the detectors. |
---|
| 38 | // |
---|
| 39 | // It is also tasked with part of updating the "safety". |
---|
| 40 | // |
---|
| 41 | // ======================================================================= |
---|
| 42 | // Modified: |
---|
| 43 | // 19 Jan 2006, P.MoraDeFreitas: Fix for suspended tracks (StartTracking) |
---|
| 44 | // 11 Aug 2004, M.Asai: Add G4VSensitiveDetector* for updating stepPoint. |
---|
| 45 | // 21 June 2003, J.Apostolakis: Calling field manager with |
---|
| 46 | // track, to enable it to configure its accuracy |
---|
| 47 | // 13 May 2003, J.Apostolakis: Zero field areas now taken into |
---|
| 48 | // account correclty in all cases (thanks to W Pokorski). |
---|
| 49 | // 29 June 2001, J.Apostolakis, D.Cote-Ahern, P.Gumplinger: |
---|
| 50 | // correction for spin tracking |
---|
| 51 | // 20 Febr 2001, J.Apostolakis: update for new FieldTrack |
---|
| 52 | // 22 Sept 2000, V.Grichine: update of Kinetic Energy |
---|
| 53 | // Created: 19 March 1997, J. Apostolakis |
---|
| 54 | // ======================================================================= |
---|
| 55 | |
---|
| 56 | #include "G4Transportation.hh" |
---|
| 57 | #include "G4ProductionCutsTable.hh" |
---|
| 58 | #include "G4ParticleTable.hh" |
---|
| 59 | #include "G4ChordFinder.hh" |
---|
| 60 | #include "G4FieldManagerStore.hh" |
---|
| 61 | class G4VSensitiveDetector; |
---|
| 62 | |
---|
| 63 | ////////////////////////////////////////////////////////////////////////// |
---|
| 64 | // |
---|
| 65 | // Constructor |
---|
| 66 | |
---|
| 67 | G4Transportation::G4Transportation( G4int verboseLevel ) |
---|
| 68 | : G4VProcess( G4String("Transportation"), fTransportation ), |
---|
| 69 | fParticleIsLooping( false ), |
---|
| 70 | fPreviousSftOrigin (0.,0.,0.), |
---|
| 71 | fPreviousSafety ( 0.0 ), |
---|
| 72 | fThreshold_Warning_Energy( 100 * MeV ), |
---|
| 73 | fThreshold_Important_Energy( 250 * MeV ), |
---|
| 74 | fThresholdTrials( 10 ), |
---|
| 75 | fUnimportant_Energy( 1 * MeV ), |
---|
| 76 | fNoLooperTrials(0), |
---|
| 77 | fSumEnergyKilled( 0.0 ), fMaxEnergyKilled( 0.0 ), |
---|
| 78 | fShortStepOptimisation(false), // Old default: true (=fast short steps) |
---|
| 79 | fVerboseLevel( verboseLevel ) |
---|
| 80 | { |
---|
| 81 | G4TransportationManager* transportMgr ; |
---|
| 82 | |
---|
| 83 | transportMgr = G4TransportationManager::GetTransportationManager() ; |
---|
| 84 | |
---|
| 85 | fLinearNavigator = transportMgr->GetNavigatorForTracking() ; |
---|
| 86 | |
---|
| 87 | // fGlobalFieldMgr = transportMgr->GetFieldManager() ; |
---|
| 88 | |
---|
| 89 | fFieldPropagator = transportMgr->GetPropagatorInField() ; |
---|
| 90 | |
---|
| 91 | // fpSafetyHelper = fpTransportManager->GetSafetyHelper(); // New |
---|
| 92 | |
---|
| 93 | // Cannot determine whether a field exists here, |
---|
| 94 | // because it would only work if the field manager has informed |
---|
| 95 | // about the detector's field before this transportation process |
---|
| 96 | // is constructed. |
---|
| 97 | // Instead later the method DoesGlobalFieldExist() is called |
---|
| 98 | |
---|
| 99 | static G4TouchableHandle nullTouchableHandle; // Points to (G4VTouchable*) 0 |
---|
| 100 | fCurrentTouchableHandle = nullTouchableHandle; |
---|
| 101 | |
---|
| 102 | fEndGlobalTimeComputed = false; |
---|
| 103 | fCandidateEndGlobalTime = 0; |
---|
| 104 | } |
---|
| 105 | |
---|
| 106 | ////////////////////////////////////////////////////////////////////////// |
---|
| 107 | |
---|
| 108 | G4Transportation::~G4Transportation() |
---|
| 109 | { |
---|
| 110 | |
---|
| 111 | // --- Alternative code to delete 'junk data' in touchable handle |
---|
| 112 | // ** Corresponds to the case that the constructor has |
---|
| 113 | // fCurrentTouchableHandle = new G4TouchableHistory(); |
---|
| 114 | // ** Likely incorrect - as at the end of the simulation |
---|
| 115 | // the handle no longer holds the original 'null' history |
---|
| 116 | // G4VTouchable* pTouchable= fCurrentTouchableHandle(); // Incorrect |
---|
| 117 | // delete pTouchable; // Incorrect |
---|
| 118 | |
---|
| 119 | if( (fVerboseLevel > 0) && (fSumEnergyKilled > 0.0 ) ){ |
---|
| 120 | G4cout << " G4Transportation: Statistics for looping particles " << G4endl; |
---|
| 121 | G4cout << " Sum of energy of loopers killed: " << fSumEnergyKilled << G4endl; |
---|
| 122 | G4cout << " Max energy of loopers killed: " << fMaxEnergyKilled << G4endl; |
---|
| 123 | } |
---|
| 124 | } |
---|
| 125 | |
---|
| 126 | ////////////////////////////////////////////////////////////////////////// |
---|
| 127 | // |
---|
| 128 | // Responsibilities: |
---|
| 129 | // Find whether the geometry limits the Step, and to what length |
---|
| 130 | // Calculate the new value of the safety and return it. |
---|
| 131 | // Store the final time, position and momentum. |
---|
| 132 | |
---|
| 133 | G4double G4Transportation:: |
---|
| 134 | AlongStepGetPhysicalInteractionLength( const G4Track& track, |
---|
| 135 | G4double, // previousStepSize |
---|
| 136 | G4double currentMinimumStep, |
---|
| 137 | G4double& currentSafety, |
---|
| 138 | G4GPILSelection* selection ) |
---|
| 139 | { |
---|
| 140 | G4double geometryStepLength, newSafety ; |
---|
| 141 | fParticleIsLooping = false ; |
---|
| 142 | |
---|
| 143 | // Initial actions moved to StartTrack() |
---|
| 144 | // -------------------------------------- |
---|
| 145 | // Note: in case another process changes touchable handle |
---|
| 146 | // it will be necessary to add here (for all steps) |
---|
| 147 | // fCurrentTouchableHandle = aTrack->GetTouchableHandle(); |
---|
| 148 | |
---|
| 149 | // GPILSelection is set to defaule value of CandidateForSelection |
---|
| 150 | // It is a return value |
---|
| 151 | // |
---|
| 152 | *selection = CandidateForSelection ; |
---|
| 153 | |
---|
| 154 | // Get initial Energy/Momentum of the track |
---|
| 155 | // |
---|
| 156 | const G4DynamicParticle* pParticle = track.GetDynamicParticle() ; |
---|
| 157 | const G4ParticleDefinition* pParticleDef = pParticle->GetDefinition() ; |
---|
| 158 | G4ThreeVector startMomentumDir = pParticle->GetMomentumDirection() ; |
---|
| 159 | G4ThreeVector startPosition = track.GetPosition() ; |
---|
| 160 | |
---|
| 161 | // G4double theTime = track.GetGlobalTime() ; |
---|
| 162 | |
---|
| 163 | // The Step Point safety can be limited by other geometries and/or the |
---|
| 164 | // assumptions of any process - it's not always the geometrical safety. |
---|
| 165 | // We calculate the starting point's isotropic safety here. |
---|
| 166 | // |
---|
| 167 | G4ThreeVector OriginShift = startPosition - fPreviousSftOrigin ; |
---|
| 168 | G4double MagSqShift = OriginShift.mag2() ; |
---|
| 169 | if( MagSqShift >= sqr(fPreviousSafety) ) |
---|
| 170 | { |
---|
| 171 | currentSafety = 0.0 ; |
---|
| 172 | } |
---|
| 173 | else |
---|
| 174 | { |
---|
| 175 | currentSafety = fPreviousSafety - std::sqrt(MagSqShift) ; |
---|
| 176 | } |
---|
| 177 | |
---|
| 178 | // Is the particle charged ? |
---|
| 179 | // |
---|
| 180 | G4double particleCharge = pParticle->GetCharge() ; |
---|
| 181 | |
---|
| 182 | fGeometryLimitedStep = false ; |
---|
| 183 | // fEndGlobalTimeComputed = false ; |
---|
| 184 | |
---|
| 185 | // There is no need to locate the current volume. It is Done elsewhere: |
---|
| 186 | // On track construction |
---|
| 187 | // By the tracking, after all AlongStepDoIts, in "Relocation" |
---|
| 188 | |
---|
| 189 | // Check whether the particle have an (EM) field force exerting upon it |
---|
| 190 | // |
---|
| 191 | G4FieldManager* fieldMgr=0; |
---|
| 192 | G4bool fieldExertsForce = false ; |
---|
| 193 | if( (particleCharge != 0.0) ) |
---|
| 194 | { |
---|
| 195 | fieldMgr= fFieldPropagator->FindAndSetFieldManager( track.GetVolume() ); |
---|
| 196 | if (fieldMgr != 0) { |
---|
| 197 | // Message the field Manager, to configure it for this track |
---|
| 198 | fieldMgr->ConfigureForTrack( &track ); |
---|
| 199 | // Moved here, in order to allow a transition |
---|
| 200 | // from a zero-field status (with fieldMgr->(field)0 |
---|
| 201 | // to a finite field status |
---|
| 202 | |
---|
| 203 | // If the field manager has no field, there is no field ! |
---|
| 204 | fieldExertsForce = (fieldMgr->GetDetectorField() != 0); |
---|
| 205 | } |
---|
| 206 | } |
---|
| 207 | |
---|
| 208 | // G4cout << " G4Transport: field exerts force= " << fieldExertsForce |
---|
| 209 | // << " fieldMgr= " << fieldMgr << G4endl; |
---|
| 210 | |
---|
| 211 | // Choose the calculation of the transportation: Field or not |
---|
| 212 | // |
---|
| 213 | if( !fieldExertsForce ) |
---|
| 214 | { |
---|
| 215 | G4double linearStepLength ; |
---|
| 216 | if( fShortStepOptimisation && (currentMinimumStep <= currentSafety) ) |
---|
| 217 | { |
---|
| 218 | // The Step is guaranteed to be taken |
---|
| 219 | // |
---|
| 220 | geometryStepLength = currentMinimumStep ; |
---|
| 221 | fGeometryLimitedStep = false ; |
---|
| 222 | } |
---|
| 223 | else |
---|
| 224 | { |
---|
| 225 | // Find whether the straight path intersects a volume |
---|
| 226 | // |
---|
| 227 | linearStepLength = fLinearNavigator->ComputeStep( startPosition, |
---|
| 228 | startMomentumDir, |
---|
| 229 | currentMinimumStep, |
---|
| 230 | newSafety) ; |
---|
| 231 | // Remember last safety origin & value. |
---|
| 232 | // |
---|
| 233 | fPreviousSftOrigin = startPosition ; |
---|
| 234 | fPreviousSafety = newSafety ; |
---|
| 235 | // fSafetyHelper->SetCurrentSafety( newSafety, startPosition); |
---|
| 236 | |
---|
| 237 | // The safety at the initial point has been re-calculated: |
---|
| 238 | // |
---|
| 239 | currentSafety = newSafety ; |
---|
| 240 | |
---|
| 241 | fGeometryLimitedStep= (linearStepLength <= currentMinimumStep); |
---|
| 242 | if( fGeometryLimitedStep ) |
---|
| 243 | { |
---|
| 244 | // The geometry limits the Step size (an intersection was found.) |
---|
| 245 | geometryStepLength = linearStepLength ; |
---|
| 246 | } |
---|
| 247 | else |
---|
| 248 | { |
---|
| 249 | // The full Step is taken. |
---|
| 250 | geometryStepLength = currentMinimumStep ; |
---|
| 251 | } |
---|
| 252 | } |
---|
| 253 | endpointDistance = geometryStepLength ; |
---|
| 254 | |
---|
| 255 | // Calculate final position |
---|
| 256 | // |
---|
| 257 | fTransportEndPosition = startPosition+geometryStepLength*startMomentumDir ; |
---|
| 258 | |
---|
| 259 | // Momentum direction, energy and polarisation are unchanged by transport |
---|
| 260 | // |
---|
| 261 | fTransportEndMomentumDir = startMomentumDir ; |
---|
| 262 | fTransportEndKineticEnergy = track.GetKineticEnergy() ; |
---|
| 263 | fTransportEndSpin = track.GetPolarization(); |
---|
| 264 | fParticleIsLooping = false ; |
---|
| 265 | fMomentumChanged = false ; |
---|
| 266 | fEndGlobalTimeComputed = false ; |
---|
| 267 | } |
---|
| 268 | else // A field exerts force |
---|
| 269 | { |
---|
| 270 | G4double momentumMagnitude = pParticle->GetTotalMomentum() ; |
---|
| 271 | G4ThreeVector EndUnitMomentum ; |
---|
| 272 | G4double lengthAlongCurve ; |
---|
| 273 | G4double restMass = pParticleDef->GetPDGMass() ; |
---|
| 274 | |
---|
| 275 | fFieldPropagator->SetChargeMomentumMass( particleCharge, // in e+ units |
---|
| 276 | momentumMagnitude, // in Mev/c |
---|
| 277 | restMass ) ; |
---|
| 278 | |
---|
| 279 | G4ThreeVector spin = track.GetPolarization() ; |
---|
| 280 | G4FieldTrack aFieldTrack = G4FieldTrack( startPosition, |
---|
| 281 | track.GetMomentumDirection(), |
---|
| 282 | 0.0, |
---|
| 283 | track.GetKineticEnergy(), |
---|
| 284 | restMass, |
---|
| 285 | track.GetVelocity(), |
---|
| 286 | track.GetGlobalTime(), // Lab. |
---|
| 287 | track.GetProperTime(), // Part. |
---|
| 288 | &spin ) ; |
---|
| 289 | if( currentMinimumStep > 0 ) |
---|
| 290 | { |
---|
| 291 | // Do the Transport in the field (non recti-linear) |
---|
| 292 | // |
---|
| 293 | lengthAlongCurve = fFieldPropagator->ComputeStep( aFieldTrack, |
---|
| 294 | currentMinimumStep, |
---|
| 295 | currentSafety, |
---|
| 296 | track.GetVolume() ) ; |
---|
| 297 | fGeometryLimitedStep= lengthAlongCurve < currentMinimumStep; |
---|
| 298 | if( fGeometryLimitedStep ) { |
---|
| 299 | geometryStepLength = lengthAlongCurve ; |
---|
| 300 | } else { |
---|
| 301 | geometryStepLength = currentMinimumStep ; |
---|
| 302 | } |
---|
| 303 | } |
---|
| 304 | else |
---|
| 305 | { |
---|
| 306 | geometryStepLength = lengthAlongCurve= 0.0 ; |
---|
| 307 | fGeometryLimitedStep = false ; |
---|
| 308 | } |
---|
| 309 | |
---|
| 310 | // Remember last safety origin & value. |
---|
| 311 | // |
---|
| 312 | fPreviousSftOrigin = startPosition ; |
---|
| 313 | fPreviousSafety = currentSafety ; |
---|
| 314 | // fSafetyHelper->SetCurrentSafety( newSafety, startPosition); |
---|
| 315 | |
---|
| 316 | // Get the End-Position and End-Momentum (Dir-ection) |
---|
| 317 | // |
---|
| 318 | fTransportEndPosition = aFieldTrack.GetPosition() ; |
---|
| 319 | |
---|
| 320 | // Momentum: Magnitude and direction can be changed too now ... |
---|
| 321 | // |
---|
| 322 | fMomentumChanged = true ; |
---|
| 323 | fTransportEndMomentumDir = aFieldTrack.GetMomentumDir() ; |
---|
| 324 | |
---|
| 325 | fTransportEndKineticEnergy = aFieldTrack.GetKineticEnergy() ; |
---|
| 326 | |
---|
| 327 | if( fFieldPropagator->GetCurrentFieldManager()->DoesFieldChangeEnergy() ) |
---|
| 328 | { |
---|
| 329 | // If the field can change energy, then the time must be integrated |
---|
| 330 | // - so this should have been updated |
---|
| 331 | // |
---|
| 332 | fCandidateEndGlobalTime = aFieldTrack.GetLabTimeOfFlight(); |
---|
| 333 | fEndGlobalTimeComputed = true; |
---|
| 334 | |
---|
| 335 | // was ( fCandidateEndGlobalTime != track.GetGlobalTime() ); |
---|
| 336 | // a cleaner way is to have FieldTrack knowing whether time is updated. |
---|
| 337 | } |
---|
| 338 | else |
---|
| 339 | { |
---|
| 340 | // The energy should be unchanged by field transport, |
---|
| 341 | // - so the time changed will be calculated elsewhere |
---|
| 342 | // |
---|
| 343 | fEndGlobalTimeComputed = false; |
---|
| 344 | |
---|
| 345 | // Check that the integration preserved the energy |
---|
| 346 | // - and if not correct this! |
---|
| 347 | G4double startEnergy= track.GetKineticEnergy(); |
---|
| 348 | G4double endEnergy= fTransportEndKineticEnergy; |
---|
| 349 | |
---|
| 350 | static G4int no_inexact_steps=0, no_large_ediff; |
---|
| 351 | G4double absEdiff = std::fabs(startEnergy- endEnergy); |
---|
| 352 | if( absEdiff > perMillion * endEnergy ) |
---|
| 353 | { |
---|
| 354 | no_inexact_steps++; |
---|
| 355 | // Possible statistics keeping here ... |
---|
| 356 | } |
---|
| 357 | if( fVerboseLevel > 1 ) |
---|
| 358 | { |
---|
| 359 | if( std::fabs(startEnergy- endEnergy) > perThousand * endEnergy ) |
---|
| 360 | { |
---|
| 361 | static G4int no_warnings= 0, warnModulo=1, moduloFactor= 10; |
---|
| 362 | no_large_ediff ++; |
---|
| 363 | if( (no_large_ediff% warnModulo) == 0 ) |
---|
| 364 | { |
---|
| 365 | no_warnings++; |
---|
| 366 | G4cout << "WARNING - G4Transportation::AlongStepGetPIL() " |
---|
| 367 | << " Energy change in Step is above 1^-3 relative value. " << G4endl |
---|
| 368 | << " Relative change in 'tracking' step = " |
---|
| 369 | << std::setw(15) << (endEnergy-startEnergy)/startEnergy << G4endl |
---|
| 370 | << " Starting E= " << std::setw(12) << startEnergy / MeV << " MeV " << G4endl |
---|
| 371 | << " Ending E= " << std::setw(12) << endEnergy / MeV << " MeV " << G4endl; |
---|
| 372 | G4cout << " Energy has been corrected -- however, review" |
---|
| 373 | << " field propagation parameters for accuracy." << G4endl; |
---|
| 374 | if( (fVerboseLevel > 2 ) || (no_warnings<4) || (no_large_ediff == warnModulo * moduloFactor) ){ |
---|
| 375 | G4cout << " These include EpsilonStepMax(/Min) in G4FieldManager " |
---|
| 376 | << " which determine fractional error per step for integrated quantities. " << G4endl |
---|
| 377 | << " Note also the influence of the permitted number of integration steps." |
---|
| 378 | << G4endl; |
---|
| 379 | } |
---|
| 380 | G4cerr << "ERROR - G4Transportation::AlongStepGetPIL()" << G4endl |
---|
| 381 | << " Bad 'endpoint'. Energy change detected" |
---|
| 382 | << " and corrected. " |
---|
| 383 | << " Has occurred already " |
---|
| 384 | << no_large_ediff << " times." << G4endl; |
---|
| 385 | if( no_large_ediff == warnModulo * moduloFactor ) |
---|
| 386 | { |
---|
| 387 | warnModulo *= moduloFactor; |
---|
| 388 | } |
---|
| 389 | } |
---|
| 390 | } |
---|
| 391 | } // end of if (fVerboseLevel) |
---|
| 392 | |
---|
| 393 | // Correct the energy for fields that conserve it |
---|
| 394 | // This - hides the integration error |
---|
| 395 | // - but gives a better physical answer |
---|
| 396 | fTransportEndKineticEnergy= track.GetKineticEnergy(); |
---|
| 397 | } |
---|
| 398 | |
---|
| 399 | fTransportEndSpin = aFieldTrack.GetSpin(); |
---|
| 400 | fParticleIsLooping = fFieldPropagator->IsParticleLooping() ; |
---|
| 401 | endpointDistance = (fTransportEndPosition - startPosition).mag() ; |
---|
| 402 | } |
---|
| 403 | |
---|
| 404 | // If we are asked to go a step length of 0, and we are on a boundary |
---|
| 405 | // then a boundary will also limit the step -> we must flag this. |
---|
| 406 | // |
---|
| 407 | if( currentMinimumStep == 0.0 ) |
---|
| 408 | { |
---|
| 409 | if( currentSafety == 0.0 ) fGeometryLimitedStep = true ; |
---|
| 410 | } |
---|
| 411 | |
---|
| 412 | // Update the safety starting from the end-point, |
---|
| 413 | // if it will become negative at the end-point. |
---|
| 414 | // |
---|
| 415 | if( currentSafety < endpointDistance ) |
---|
| 416 | { |
---|
| 417 | G4double endSafety = |
---|
| 418 | fLinearNavigator->ComputeSafety( fTransportEndPosition) ; |
---|
| 419 | currentSafety = endSafety ; |
---|
| 420 | fPreviousSftOrigin = fTransportEndPosition ; |
---|
| 421 | fPreviousSafety = currentSafety ; |
---|
| 422 | // fSafetyHelper->SetCurrentSafety( currentSafety, fTransportEndPosition); |
---|
| 423 | |
---|
| 424 | // Because the Stepping Manager assumes it is from the start point, |
---|
| 425 | // add the StepLength |
---|
| 426 | // |
---|
| 427 | currentSafety += endpointDistance ; |
---|
| 428 | |
---|
| 429 | #ifdef G4DEBUG_TRANSPORT |
---|
| 430 | G4cout.precision(12) ; |
---|
| 431 | G4cout << "***G4Transportation::AlongStepGPIL ** " << G4endl ; |
---|
| 432 | G4cout << " Called Navigator->ComputeSafety at " << fTransportEndPosition |
---|
| 433 | << " and it returned safety= " << endSafety << G4endl ; |
---|
| 434 | G4cout << " Adding endpoint distance " << endpointDistance |
---|
| 435 | << " to obtain pseudo-safety= " << currentSafety << G4endl ; |
---|
| 436 | #endif |
---|
| 437 | } |
---|
| 438 | |
---|
| 439 | fParticleChange.ProposeTrueStepLength(geometryStepLength) ; |
---|
| 440 | |
---|
| 441 | return geometryStepLength ; |
---|
| 442 | } |
---|
| 443 | |
---|
| 444 | ////////////////////////////////////////////////////////////////////////// |
---|
| 445 | // |
---|
| 446 | // Initialize ParticleChange (by setting all its members equal |
---|
| 447 | // to corresponding members in G4Track) |
---|
| 448 | |
---|
| 449 | G4VParticleChange* G4Transportation::AlongStepDoIt( const G4Track& track, |
---|
| 450 | const G4Step& stepData ) |
---|
| 451 | { |
---|
| 452 | static G4int noCalls=0; |
---|
| 453 | static const G4ParticleDefinition* fOpticalPhoton = |
---|
| 454 | G4ParticleTable::GetParticleTable()->FindParticle("opticalphoton"); |
---|
| 455 | |
---|
| 456 | noCalls++; |
---|
| 457 | |
---|
| 458 | fParticleChange.Initialize(track) ; |
---|
| 459 | |
---|
| 460 | // Code for specific process |
---|
| 461 | // |
---|
| 462 | fParticleChange.ProposePosition(fTransportEndPosition) ; |
---|
| 463 | fParticleChange.ProposeMomentumDirection(fTransportEndMomentumDir) ; |
---|
| 464 | fParticleChange.ProposeEnergy(fTransportEndKineticEnergy) ; |
---|
| 465 | fParticleChange.SetMomentumChanged(fMomentumChanged) ; |
---|
| 466 | |
---|
| 467 | fParticleChange.ProposePolarization(fTransportEndSpin); |
---|
| 468 | |
---|
| 469 | G4double deltaTime = 0.0 ; |
---|
| 470 | |
---|
| 471 | // Calculate Lab Time of Flight (ONLY if field Equations used it!) |
---|
| 472 | // G4double endTime = fCandidateEndGlobalTime; |
---|
| 473 | // G4double delta_time = endTime - startTime; |
---|
| 474 | |
---|
| 475 | G4double startTime = track.GetGlobalTime() ; |
---|
| 476 | |
---|
| 477 | if (!fEndGlobalTimeComputed) |
---|
| 478 | { |
---|
| 479 | // The time was not integrated .. make the best estimate possible |
---|
| 480 | // |
---|
| 481 | G4double finalVelocity = track.GetVelocity() ; |
---|
| 482 | G4double initialVelocity = stepData.GetPreStepPoint()->GetVelocity() ; |
---|
| 483 | G4double stepLength = track.GetStepLength() ; |
---|
| 484 | |
---|
| 485 | deltaTime= 0.0; // in case initialVelocity = 0 |
---|
| 486 | const G4DynamicParticle* fpDynamicParticle = track.GetDynamicParticle(); |
---|
| 487 | if (fpDynamicParticle->GetDefinition()== fOpticalPhoton) |
---|
| 488 | { |
---|
| 489 | // A photon is in the medium of the final point |
---|
| 490 | // during the step, so it has the final velocity. |
---|
| 491 | deltaTime = stepLength/finalVelocity ; |
---|
| 492 | } |
---|
| 493 | else if (finalVelocity > 0.0) |
---|
| 494 | { |
---|
| 495 | G4double meanInverseVelocity ; |
---|
| 496 | // deltaTime = stepLength/finalVelocity ; |
---|
| 497 | meanInverseVelocity = 0.5 |
---|
| 498 | * ( 1.0 / initialVelocity + 1.0 / finalVelocity ) ; |
---|
| 499 | deltaTime = stepLength * meanInverseVelocity ; |
---|
| 500 | } |
---|
| 501 | else if( initialVelocity > 0.0 ) |
---|
| 502 | { |
---|
| 503 | deltaTime = stepLength/initialVelocity ; |
---|
| 504 | } |
---|
| 505 | fCandidateEndGlobalTime = startTime + deltaTime ; |
---|
| 506 | } |
---|
| 507 | else |
---|
| 508 | { |
---|
| 509 | deltaTime = fCandidateEndGlobalTime - startTime ; |
---|
| 510 | } |
---|
| 511 | |
---|
| 512 | fParticleChange.ProposeGlobalTime( fCandidateEndGlobalTime ) ; |
---|
| 513 | |
---|
| 514 | // Now Correct by Lorentz factor to get "proper" deltaTime |
---|
| 515 | |
---|
| 516 | G4double restMass = track.GetDynamicParticle()->GetMass() ; |
---|
| 517 | G4double deltaProperTime = deltaTime*( restMass/track.GetTotalEnergy() ) ; |
---|
| 518 | |
---|
| 519 | fParticleChange.ProposeProperTime(track.GetProperTime() + deltaProperTime) ; |
---|
| 520 | //fParticleChange. ProposeTrueStepLength( track.GetStepLength() ) ; |
---|
| 521 | |
---|
| 522 | // If the particle is caught looping or is stuck (in very difficult |
---|
| 523 | // boundaries) in a magnetic field (doing many steps) |
---|
| 524 | // THEN this kills it ... |
---|
| 525 | // |
---|
| 526 | if ( fParticleIsLooping ) |
---|
| 527 | { |
---|
| 528 | G4double endEnergy= fTransportEndKineticEnergy; |
---|
| 529 | |
---|
| 530 | if( (endEnergy < fThreshold_Important_Energy) |
---|
| 531 | || (fNoLooperTrials >= fThresholdTrials ) ){ |
---|
| 532 | // Kill the looping particle |
---|
| 533 | // |
---|
| 534 | fParticleChange.ProposeTrackStatus( fStopAndKill ) ; |
---|
| 535 | |
---|
| 536 | // 'Bare' statistics |
---|
| 537 | fSumEnergyKilled += endEnergy; |
---|
| 538 | if( endEnergy > fMaxEnergyKilled) { fMaxEnergyKilled= endEnergy; } |
---|
| 539 | |
---|
| 540 | #ifdef G4VERBOSE |
---|
| 541 | if( (fVerboseLevel > 1) || |
---|
| 542 | ( endEnergy > fThreshold_Warning_Energy ) ) { |
---|
| 543 | G4cout << " G4Transportation is killing track that is looping or stuck " |
---|
| 544 | << G4endl |
---|
| 545 | << " This track has " << track.GetKineticEnergy() / MeV |
---|
| 546 | << " MeV energy." << G4endl; |
---|
| 547 | G4cout << " Number of trials = " << fNoLooperTrials |
---|
| 548 | << " No of calls to AlongStepDoIt = " << noCalls |
---|
| 549 | << G4endl; |
---|
| 550 | } |
---|
| 551 | #endif |
---|
| 552 | fNoLooperTrials=0; |
---|
| 553 | } |
---|
| 554 | else{ |
---|
| 555 | fNoLooperTrials ++; |
---|
| 556 | #ifdef G4VERBOSE |
---|
| 557 | if( (fVerboseLevel > 2) ){ |
---|
| 558 | G4cout << " G4Transportation::AlongStepDoIt(): Particle looping - " |
---|
| 559 | << " Number of trials = " << fNoLooperTrials |
---|
| 560 | << " No of calls to = " << noCalls |
---|
| 561 | << G4endl; |
---|
| 562 | } |
---|
| 563 | #endif |
---|
| 564 | } |
---|
| 565 | }else{ |
---|
| 566 | fNoLooperTrials=0; |
---|
| 567 | } |
---|
| 568 | |
---|
| 569 | // Another (sometimes better way) is to use a user-limit maximum Step size |
---|
| 570 | // to alleviate this problem .. |
---|
| 571 | |
---|
| 572 | // Introduce smooth curved trajectories to particle-change |
---|
| 573 | // |
---|
| 574 | fParticleChange.SetPointerToVectorOfAuxiliaryPoints |
---|
| 575 | (fFieldPropagator->GimmeTrajectoryVectorAndForgetIt() ); |
---|
| 576 | |
---|
| 577 | return &fParticleChange ; |
---|
| 578 | } |
---|
| 579 | |
---|
| 580 | ////////////////////////////////////////////////////////////////////////// |
---|
| 581 | // |
---|
| 582 | // This ensures that the PostStep action is always called, |
---|
| 583 | // so that it can do the relocation if it is needed. |
---|
| 584 | // |
---|
| 585 | |
---|
| 586 | G4double G4Transportation:: |
---|
| 587 | PostStepGetPhysicalInteractionLength( const G4Track&, |
---|
| 588 | G4double, // previousStepSize |
---|
| 589 | G4ForceCondition* pForceCond ) |
---|
| 590 | { |
---|
| 591 | *pForceCond = Forced ; |
---|
| 592 | return DBL_MAX ; // was kInfinity ; but convention now is DBL_MAX |
---|
| 593 | } |
---|
| 594 | |
---|
| 595 | ///////////////////////////////////////////////////////////////////////////// |
---|
| 596 | // |
---|
| 597 | |
---|
| 598 | G4VParticleChange* G4Transportation::PostStepDoIt( const G4Track& track, |
---|
| 599 | const G4Step& ) |
---|
| 600 | { |
---|
| 601 | G4TouchableHandle retCurrentTouchable ; // The one to return |
---|
| 602 | |
---|
| 603 | // Initialize ParticleChange (by setting all its members equal |
---|
| 604 | // to corresponding members in G4Track) |
---|
| 605 | // fParticleChange.Initialize(track) ; // To initialise TouchableChange |
---|
| 606 | |
---|
| 607 | fParticleChange.ProposeTrackStatus(track.GetTrackStatus()) ; |
---|
| 608 | |
---|
| 609 | // If the Step was determined by the volume boundary, |
---|
| 610 | // logically relocate the particle |
---|
| 611 | |
---|
| 612 | if(fGeometryLimitedStep) |
---|
| 613 | { |
---|
| 614 | // fCurrentTouchable will now become the previous touchable, |
---|
| 615 | // and what was the previous will be freed. |
---|
| 616 | // (Needed because the preStepPoint can point to the previous touchable) |
---|
| 617 | |
---|
| 618 | fLinearNavigator->SetGeometricallyLimitedStep() ; |
---|
| 619 | fLinearNavigator-> |
---|
| 620 | LocateGlobalPointAndUpdateTouchableHandle( track.GetPosition(), |
---|
| 621 | track.GetMomentumDirection(), |
---|
| 622 | fCurrentTouchableHandle, |
---|
| 623 | true ) ; |
---|
| 624 | // Check whether the particle is out of the world volume |
---|
| 625 | // If so it has exited and must be killed. |
---|
| 626 | // |
---|
| 627 | if( fCurrentTouchableHandle->GetVolume() == 0 ) |
---|
| 628 | { |
---|
| 629 | fParticleChange.ProposeTrackStatus( fStopAndKill ) ; |
---|
| 630 | } |
---|
| 631 | retCurrentTouchable = fCurrentTouchableHandle ; |
---|
| 632 | fParticleChange.SetTouchableHandle( fCurrentTouchableHandle ) ; |
---|
| 633 | } |
---|
| 634 | else // fGeometryLimitedStep is false |
---|
| 635 | { |
---|
| 636 | // This serves only to move the Navigator's location |
---|
| 637 | // |
---|
| 638 | fLinearNavigator->LocateGlobalPointWithinVolume( track.GetPosition() ) ; |
---|
| 639 | |
---|
| 640 | // The value of the track's current Touchable is retained. |
---|
| 641 | // (and it must be correct because we must use it below to |
---|
| 642 | // overwrite the (unset) one in particle change) |
---|
| 643 | // It must be fCurrentTouchable too ?? |
---|
| 644 | // |
---|
| 645 | fParticleChange.SetTouchableHandle( track.GetTouchableHandle() ) ; |
---|
| 646 | retCurrentTouchable = track.GetTouchableHandle() ; |
---|
| 647 | } // endif ( fGeometryLimitedStep ) |
---|
| 648 | |
---|
| 649 | const G4VPhysicalVolume* pNewVol = retCurrentTouchable->GetVolume() ; |
---|
| 650 | const G4Material* pNewMaterial = 0 ; |
---|
| 651 | const G4VSensitiveDetector* pNewSensitiveDetector = 0 ; |
---|
| 652 | |
---|
| 653 | if( pNewVol != 0 ) |
---|
| 654 | { |
---|
| 655 | pNewMaterial= pNewVol->GetLogicalVolume()->GetMaterial(); |
---|
| 656 | pNewSensitiveDetector= pNewVol->GetLogicalVolume()->GetSensitiveDetector(); |
---|
| 657 | } |
---|
| 658 | |
---|
| 659 | // ( <const_cast> pNewMaterial ) ; |
---|
| 660 | // ( <const_cast> pNewSensitiveDetector) ; |
---|
| 661 | |
---|
| 662 | fParticleChange.SetMaterialInTouchable( (G4Material *) pNewMaterial ) ; |
---|
| 663 | fParticleChange.SetSensitiveDetectorInTouchable( (G4VSensitiveDetector *) pNewSensitiveDetector ) ; |
---|
| 664 | |
---|
| 665 | const G4MaterialCutsCouple* pNewMaterialCutsCouple = 0; |
---|
| 666 | if( pNewVol != 0 ) |
---|
| 667 | { |
---|
| 668 | pNewMaterialCutsCouple=pNewVol->GetLogicalVolume()->GetMaterialCutsCouple(); |
---|
| 669 | } |
---|
| 670 | |
---|
| 671 | if( pNewVol!=0 && pNewMaterialCutsCouple!=0 && pNewMaterialCutsCouple->GetMaterial()!=pNewMaterial ) |
---|
| 672 | { |
---|
| 673 | // for parametrized volume |
---|
| 674 | // |
---|
| 675 | pNewMaterialCutsCouple = |
---|
| 676 | G4ProductionCutsTable::GetProductionCutsTable() |
---|
| 677 | ->GetMaterialCutsCouple(pNewMaterial, |
---|
| 678 | pNewMaterialCutsCouple->GetProductionCuts()); |
---|
| 679 | } |
---|
| 680 | fParticleChange.SetMaterialCutsCoupleInTouchable( pNewMaterialCutsCouple ); |
---|
| 681 | |
---|
| 682 | // temporarily until Get/Set Material of ParticleChange, |
---|
| 683 | // and StepPoint can be made const. |
---|
| 684 | // Set the touchable in ParticleChange |
---|
| 685 | // this must always be done because the particle change always |
---|
| 686 | // uses this value to overwrite the current touchable pointer. |
---|
| 687 | // |
---|
| 688 | fParticleChange.SetTouchableHandle(retCurrentTouchable) ; |
---|
| 689 | |
---|
| 690 | return &fParticleChange ; |
---|
| 691 | } |
---|
| 692 | |
---|
| 693 | // New method takes over the responsibility to reset the state of G4Transportation |
---|
| 694 | // object at the start of a new track or the resumption of a suspended track. |
---|
| 695 | |
---|
| 696 | void |
---|
| 697 | G4Transportation::StartTracking(G4Track* aTrack) |
---|
| 698 | { |
---|
| 699 | G4VProcess::StartTracking(aTrack); |
---|
| 700 | |
---|
| 701 | // The actions here are those that were taken in AlongStepGPIL |
---|
| 702 | // when track.GetCurrentStepNumber()==1 |
---|
| 703 | |
---|
| 704 | // reset safety value and center |
---|
| 705 | // |
---|
| 706 | fPreviousSafety = 0.0 ; |
---|
| 707 | fPreviousSftOrigin = G4ThreeVector(0.,0.,0.) ; |
---|
| 708 | |
---|
| 709 | // reset looping counter -- for motion in field |
---|
| 710 | fNoLooperTrials= 0; |
---|
| 711 | // Must clear this state .. else it depends on last track's value |
---|
| 712 | // --> a better solution would set this from state of suspended track TODO ? |
---|
| 713 | // Was if( aTrack->GetCurrentStepNumber()==1 ) { .. } |
---|
| 714 | |
---|
| 715 | // ChordFinder reset internal state |
---|
| 716 | // |
---|
| 717 | if( DoesGlobalFieldExist() ) { |
---|
| 718 | fFieldPropagator->ClearPropagatorState(); |
---|
| 719 | // Resets all state of field propagator class (ONLY) |
---|
| 720 | // including safety values (in case of overlaps and to wipe for first track). |
---|
| 721 | |
---|
| 722 | // G4ChordFinder* chordF= fFieldPropagator->GetChordFinder(); |
---|
| 723 | // if( chordF ) chordF->ResetStepEstimate(); |
---|
| 724 | } |
---|
| 725 | |
---|
| 726 | // Make sure to clear the chord finders of all fields (ie managers) |
---|
| 727 | static G4FieldManagerStore* fieldMgrStore= G4FieldManagerStore::GetInstance(); |
---|
| 728 | fieldMgrStore->ClearAllChordFindersState(); |
---|
| 729 | |
---|
| 730 | // Update the current touchable handle (from the track's) |
---|
| 731 | // |
---|
| 732 | fCurrentTouchableHandle = aTrack->GetTouchableHandle(); |
---|
| 733 | } |
---|
| 734 | |
---|