#include <G4CoupledTransportation.hh>
Inheritance diagram for G4CoupledTransportation:
Definition at line 64 of file G4CoupledTransportation.hh.
G4CoupledTransportation::G4CoupledTransportation | ( | G4int | verbosityLevel = 0 |
) |
Definition at line 64 of file G4CoupledTransportation.cc.
References G4TransportationManager::ActivateNavigator(), COUPLED_TRANSPORTATION, G4cout, G4endl, G4PathFinder::GetInstance(), G4TransportationManager::GetNavigatorForTracking(), G4TransportationManager::GetPropagatorInField(), G4TransportationManager::GetSafetyHelper(), G4TransportationManager::GetTransportationManager(), and G4VProcess::SetProcessSubType().
00065 : G4VProcess( G4String("CoupledTransportation"), fTransportation ), 00066 fParticleIsLooping( false ), 00067 fPreviousSftOrigin( 0.,0.,0. ), 00068 fPreviousMassSafety( 0.0 ), 00069 fPreviousFullSafety( 0.0 ), 00070 fThreshold_Warning_Energy( 100 * MeV ), 00071 fThreshold_Important_Energy( 250 * MeV ), 00072 fThresholdTrials( 10 ), 00073 fNoLooperTrials( 0 ), 00074 fSumEnergyKilled( 0.0 ), fMaxEnergyKilled( 0.0 ), 00075 fUseMagneticMoment( false ), 00076 fVerboseLevel( verbosity ) 00077 { 00078 // set Process Sub Type 00079 SetProcessSubType(static_cast<G4int>(COUPLED_TRANSPORTATION)); 00080 00081 G4TransportationManager* transportMgr ; 00082 00083 transportMgr = G4TransportationManager::GetTransportationManager() ; 00084 00085 fMassNavigator = transportMgr->GetNavigatorForTracking() ; 00086 fFieldPropagator = transportMgr->GetPropagatorInField() ; 00087 // fGlobalFieldMgr = transportMgr->GetFieldManager() ; 00088 fNavigatorId= transportMgr->ActivateNavigator( fMassNavigator ); 00089 if( fVerboseLevel > 0 ) 00090 { 00091 G4cout << " G4CoupledTransportation constructor: ----- " << G4endl; 00092 G4cout << " Verbose level is " << fVerboseLevel << G4endl; 00093 G4cout << " Navigator Id obtained in G4CoupledTransportation constructor " 00094 << fNavigatorId << G4endl; 00095 } 00096 fPathFinder= G4PathFinder::GetInstance(); 00097 fpSafetyHelper = transportMgr->GetSafetyHelper(); // New 00098 00099 // Following assignment is to fix small memory leak from simple use of 'new' 00100 static G4TouchableHandle nullTouchableHandle; // Points to (G4VTouchable*) 0 00101 fCurrentTouchableHandle = nullTouchableHandle; 00102 00103 fEndGlobalTimeComputed = false; 00104 fCandidateEndGlobalTime = 0; 00105 }
G4CoupledTransportation::~G4CoupledTransportation | ( | ) |
Definition at line 109 of file G4CoupledTransportation.cc.
References G4cout, and G4endl.
00110 { 00111 // fCurrentTouchableHandle is a data member - no deletion required 00112 00113 if( (fVerboseLevel > 0) || (fSumEnergyKilled > 0.0 ) ) 00114 { 00115 G4cout << " G4CoupledTransportation: Statistics for looping particles " << G4endl; 00116 G4cout << " Sum of energy of loopers killed: " << fSumEnergyKilled << G4endl; 00117 G4cout << " Max energy of loopers killed: " << fMaxEnergyKilled << G4endl; 00118 } 00119 }
G4VParticleChange * G4CoupledTransportation::AlongStepDoIt | ( | const G4Track & | track, | |
const G4Step & | stepData | |||
) | [virtual] |
Implements G4VProcess.
Definition at line 499 of file G4CoupledTransportation.cc.
References DBL_MAX, fStopAndKill, G4cout, G4endl, G4Track::GetCurrentStepNumber(), G4Track::GetDynamicParticle(), G4Track::GetGlobalTime(), G4Track::GetKineticEnergy(), G4Track::GetLocalTime(), G4DynamicParticle::GetMass(), G4Step::GetPreStepPoint(), G4Track::GetProperTime(), G4Track::GetStepLength(), G4Track::GetTotalEnergy(), G4StepPoint::GetVelocity(), G4Track::GetVelocity(), G4ParticleChangeForTransport::Initialize(), G4ParticleChange::ProposeEnergy(), G4ParticleChange::ProposeGlobalTime(), G4ParticleChange::ProposeLocalTime(), G4ParticleChange::ProposeMomentumDirection(), G4ParticleChange::ProposePolarization(), G4ParticleChange::ProposePosition(), G4ParticleChange::ProposeProperTime(), G4VParticleChange::ProposeTrackStatus(), and G4ParticleChangeForTransport::SetMomentumChanged().
00501 { 00502 static G4int noCalls=0; 00503 noCalls++; 00504 00505 fParticleChange.Initialize(track) ; 00506 // sets all its members to the value of corresponding members in G4Track 00507 00508 // Code specific for Transport 00509 // 00510 fParticleChange.ProposePosition(fTransportEndPosition) ; 00511 // G4cout << " G4CoupledTransportation::AlongStepDoIt" 00512 // << " proposes position = " << fTransportEndPosition 00513 // << " and end momentum direction = " << fTransportEndMomentumDir << G4endl; 00514 fParticleChange.ProposeMomentumDirection(fTransportEndMomentumDir) ; 00515 fParticleChange.ProposeEnergy(fTransportEndKineticEnergy) ; 00516 fParticleChange.SetMomentumChanged(fMomentumChanged) ; 00517 00518 fParticleChange.ProposePolarization(fTransportEndSpin); 00519 00520 G4double deltaTime = 0.0 ; 00521 00522 // Calculate Lab Time of Flight (ONLY if field Equations used it!) 00523 // G4double endTime = fCandidateEndGlobalTime; 00524 // G4double delta_time = endTime - startTime; 00525 00526 G4double startTime = track.GetGlobalTime() ; 00527 00528 if (!fEndGlobalTimeComputed) 00529 { 00530 G4double finalInverseVel= DBL_MAX, initialInverseVel=DBL_MAX; 00531 00532 // The time was not integrated .. make the best estimate possible 00533 // 00534 G4double finalVelocity = track.GetVelocity() ; 00535 if( finalVelocity > 0.0 ) { finalInverseVel= 1.0 / finalVelocity; } 00536 G4double initialVelocity = stepData.GetPreStepPoint()->GetVelocity() ; 00537 if( initialVelocity > 0.0 ) { initialInverseVel= 1.0 / initialVelocity; } 00538 G4double stepLength = track.GetStepLength() ; 00539 00540 if (finalVelocity > 0.0) 00541 { 00542 // deltaTime = stepLength/finalVelocity ; 00543 G4double meanInverseVelocity = 0.5 * ( initialInverseVel + finalInverseVel ); 00544 deltaTime = stepLength * meanInverseVelocity ; 00545 // G4cout << " dt = s * mean(1/v) , with " << " s = " << stepLength 00546 // << " mean(1/v)= " << meanInverseVelocity << G4endl; 00547 } 00548 else 00549 { 00550 deltaTime = stepLength * initialInverseVel ; 00551 // G4cout << " dt = s / initV " << " s = " << stepLength 00552 // << " 1 / initV= " << initialInverseVel << G4endl; 00553 } // Could do with better estimate for final step (finalVelocity = 0) ? 00554 00555 fCandidateEndGlobalTime = startTime + deltaTime ; 00556 fParticleChange.ProposeLocalTime( track.GetLocalTime() + deltaTime) ; 00557 00558 // G4cout << " Calculated global time from start = " << startTime << " and " 00559 // << " delta time = " << deltaTime << G4endl; 00560 } 00561 else 00562 { 00563 deltaTime = fCandidateEndGlobalTime - startTime ; 00564 fParticleChange.ProposeGlobalTime( fCandidateEndGlobalTime ) ; 00565 // G4cout << " Calculated global time from candidate end time = " 00566 // << fCandidateEndGlobalTime << " and start time = " << startTime << G4endl; 00567 } 00568 00569 // G4cout << " G4CoupledTransportation::AlongStepDoIt " 00570 // << " flag whether computed time = " << fEndGlobalTimeComputed << " and " 00571 // << " is proposes end time " << fCandidateEndGlobalTime << G4endl; 00572 00573 // Now Correct by Lorentz factor to get "proper" deltaTime 00574 00575 G4double restMass = track.GetDynamicParticle()->GetMass() ; 00576 G4double deltaProperTime = deltaTime*( restMass/track.GetTotalEnergy() ) ; 00577 00578 fParticleChange.ProposeProperTime(track.GetProperTime() + deltaProperTime) ; 00579 //fParticleChange. ProposeTrueStepLength( track.GetStepLength() ) ; 00580 00581 // If the particle is caught looping or is stuck (in very difficult 00582 // boundaries) in a magnetic field (doing many steps) 00583 // THEN this kills it ... 00584 // 00585 if ( fParticleIsLooping ) 00586 { 00587 G4double endEnergy= fTransportEndKineticEnergy; 00588 00589 if( (endEnergy < fThreshold_Important_Energy) 00590 || (fNoLooperTrials >= fThresholdTrials ) ) 00591 { 00592 // Kill the looping particle 00593 // 00594 fParticleChange.ProposeTrackStatus( fStopAndKill ) ; 00595 00596 // 'Bare' statistics 00597 fSumEnergyKilled += endEnergy; 00598 if( endEnergy > fMaxEnergyKilled) { fMaxEnergyKilled= endEnergy; } 00599 00600 #ifdef G4VERBOSE 00601 if((fVerboseLevel > 1) && ( endEnergy > fThreshold_Warning_Energy )) 00602 { 00603 G4cout << " G4CoupledTransportation is killing track that is looping or stuck " << G4endl 00604 << " This track has " << track.GetKineticEnergy() / MeV 00605 << " MeV energy." << G4endl; 00606 } 00607 if( fVerboseLevel > 0 ) 00608 { 00609 G4cout << " Steps by this track: " << track.GetCurrentStepNumber() << G4endl; 00610 } 00611 #endif 00612 fNoLooperTrials=0; 00613 } 00614 else 00615 { 00616 fNoLooperTrials ++; 00617 #ifdef G4VERBOSE 00618 if( (fVerboseLevel > 2) ) 00619 { 00620 G4cout << " ** G4CoupledTransportation::AlongStepDoIt(): Particle looping - " << G4endl 00621 << " Number of consecutive problem step (this track) = " << fNoLooperTrials << G4endl 00622 << " Steps by this track: " << track.GetCurrentStepNumber() << G4endl 00623 << " Total no of calls to this method (all tracks) = " << noCalls << G4endl; 00624 } 00625 #endif 00626 } 00627 } 00628 else 00629 { 00630 fNoLooperTrials=0; 00631 } 00632 00633 // Another (sometimes better way) is to use a user-limit maximum Step size 00634 // to alleviate this problem .. 00635 00636 // Add smooth curved trajectories to particle-change 00637 // 00638 // fParticleChange.SetPointerToVectorOfAuxiliaryPoints 00639 // (fFieldPropagator->GimmeTrajectoryVectorAndForgetIt() ); 00640 00641 return &fParticleChange ; 00642 }
G4double G4CoupledTransportation::AlongStepGetPhysicalInteractionLength | ( | const G4Track & | track, | |
G4double | previousStepSize, | |||
G4double | currentMinimumStep, | |||
G4double & | currentSafety, | |||
G4GPILSelection * | selection | |||
) | [virtual] |
Implements G4VProcess.
Definition at line 129 of file G4CoupledTransportation.cc.
References CandidateForSelection, G4PathFinder::ComputeSafety(), G4PathFinder::ComputeStep(), G4FieldManager::ConfigureForTrack(), G4FieldManager::DoesFieldChangeEnergy(), FatalException, G4PropagatorInField::FindAndSetFieldManager(), G4cerr, G4cout, G4endl, G4Exception(), G4DynamicParticle::GetCharge(), G4PropagatorInField::GetCurrentFieldManager(), G4PathFinder::GetCurrentSafety(), G4Track::GetCurrentStepNumber(), G4DynamicParticle::GetDefinition(), G4FieldManager::GetDetectorField(), G4Track::GetDynamicParticle(), G4Track::GetGlobalTime(), G4FieldTrack::GetKineticEnergy(), G4Track::GetKineticEnergy(), G4FieldTrack::GetLabTimeOfFlight(), G4DynamicParticle::GetMagneticMoment(), G4FieldTrack::GetMomentumDir(), G4Track::GetMomentumDirection(), G4DynamicParticle::GetMomentumDirection(), G4VPhysicalVolume::GetName(), G4PathFinder::GetNumberGeometriesLimitingStep(), G4ParticleDefinition::GetPDGMass(), G4Track::GetPolarization(), G4FieldTrack::GetPosition(), G4Track::GetPosition(), G4Track::GetProperTime(), G4FieldTrack::GetSpin(), G4DynamicParticle::GetTotalMomentum(), G4Track::GetVolume(), G4Field::IsGravityActive(), G4PropagatorInField::IsParticleLooping(), kSharedTransport, kUnique, G4PathFinder::ObtainSafety(), G4VParticleChange::ProposeTrueStepLength(), ReportInexactEnergy(), G4FieldTrack::SetChargeAndMoments(), G4PropagatorInField::SetChargeMomentumMass(), G4SafetyHelper::SetCurrentSafety(), and sqr().
00134 { 00135 G4double geometryStepLength; 00136 G4double startMassSafety= 0.0; // estimated safety for start point (mass geometry) 00137 G4double startFullSafety= 0.0; // estimated safety for start point (all geometries) 00138 G4double safetyProposal= -1.0; // local copy of proposal 00139 00140 G4ThreeVector EndUnitMomentum ; 00141 G4double lengthAlongCurve=0.0 ; 00142 00143 fParticleIsLooping = false ; 00144 00145 // Initial actions moved to StartTrack() 00146 // -------------------------------------- 00147 // Note: in case another process changes touchable handle 00148 // it will be necessary to add here (for all steps) 00149 // fCurrentTouchableHandle = aTrack->GetTouchableHandle(); 00150 00151 // GPILSelection is set to defaule value of CandidateForSelection 00152 // It is a return value 00153 // 00154 *selection = CandidateForSelection ; 00155 00156 // Get initial Energy/Momentum of the track 00157 // 00158 const G4DynamicParticle* pParticle = track.GetDynamicParticle() ; 00159 const G4ParticleDefinition* pParticleDef = pParticle->GetDefinition() ; 00160 G4ThreeVector startMomentumDir = pParticle->GetMomentumDirection() ; 00161 G4ThreeVector startPosition = track.GetPosition() ; 00162 G4VPhysicalVolume* currentVolume= track.GetVolume(); 00163 00164 #ifdef G4DEBUG_TRANSPORT 00165 if( fVerboseLevel > 1 ) 00166 { 00167 G4cout << "G4CoupledTransportation::AlongStepGPIL> called in volume " 00168 << currentVolume->GetName() << G4endl; 00169 } 00170 #endif 00171 // G4double theTime = track.GetGlobalTime() ; 00172 00173 // The Step Point safety can be limited by other geometries and/or the 00174 // assumptions of any process - it's not always the geometrical safety. 00175 // We calculate the starting point's isotropic safety here. 00176 // 00177 G4ThreeVector OriginShift = startPosition - fPreviousSftOrigin ; 00178 G4double MagSqShift = OriginShift.mag2() ; 00179 startMassSafety = 0.0; 00180 startFullSafety= 0.0; 00181 00182 // Recall that FullSafety <= MassSafety 00183 // Original: if( MagSqShift < sqr(fPreviousMassSafety) ) { 00184 if( MagSqShift < sqr(fPreviousFullSafety) ) // Revision proposed by Alex H, 2 Oct 07 00185 { 00186 G4double mag_shift= std::sqrt(MagSqShift); 00187 startMassSafety = std::max( (fPreviousMassSafety - mag_shift), 0.0); 00188 startFullSafety = std::max( (fPreviousFullSafety - mag_shift), 0.0); 00189 // Need to be consistent between full safety with Mass safety 00190 // in order reproduce results in simple case --> use same calculation method 00191 00192 // Only compute full safety if massSafety > 0. Else it remains 0 00193 // startFullSafety = fPathFinder->ComputeSafety( startPosition ); 00194 } 00195 00196 // Is the particle charged or has it a magnetic moment? 00197 // 00198 G4double particleCharge = pParticle->GetCharge() ; 00199 G4double magneticMoment = pParticle->GetMagneticMoment() ; 00200 G4double restMass = pParticleDef->GetPDGMass() ; 00201 00202 fMassGeometryLimitedStep = false ; // Set default - alex 00203 fAnyGeometryLimitedStep = false; 00204 00205 // fEndGlobalTimeComputed = false ; 00206 00207 // There is no need to locate the current volume. It is Done elsewhere: 00208 // On track construction 00209 // By the tracking, after all AlongStepDoIts, in "Relocation" 00210 00211 // Check if the particle has a force, EM or gravitational, exerted on it 00212 // 00213 G4FieldManager* fieldMgr=0; 00214 G4bool fieldExertsForce = false ; 00215 00216 G4bool gravityOn = false; 00217 const G4Field* ptrField= 0; 00218 00219 fieldMgr = fFieldPropagator->FindAndSetFieldManager( track.GetVolume() ); 00220 if( fieldMgr != 0 ) 00221 { 00222 // Message the field Manager, to configure it for this track 00223 fieldMgr->ConfigureForTrack( &track ); 00224 // Here it can transition from a null field-ptr to a finite field 00225 00226 // If the field manager has no field ptr, the field is zero 00227 // by definition ( = there is no field ! ) 00228 ptrField= fieldMgr->GetDetectorField(); 00229 00230 if( ptrField != 0) 00231 { 00232 gravityOn= ptrField->IsGravityActive(); 00233 if( (particleCharge != 0.0) 00234 || (fUseMagneticMoment && (magneticMoment != 0.0) ) 00235 || (gravityOn && (restMass != 0.0)) 00236 ) 00237 { 00238 fieldExertsForce = true; 00239 } 00240 } 00241 } 00242 G4double momentumMagnitude = pParticle->GetTotalMomentum() ; 00243 00244 fFieldPropagator->SetChargeMomentumMass( particleCharge, // in e+ units 00245 momentumMagnitude, // in Mev/c 00246 restMass ) ; 00247 // This should be obsolete - the call to SetChargeAndMoments below should do the work 00248 00249 G4ThreeVector spin = track.GetPolarization() ; 00250 G4FieldTrack theFieldTrack = G4FieldTrack( startPosition, 00251 track.GetMomentumDirection(), 00252 0.0, 00253 track.GetKineticEnergy(), 00254 restMass, 00255 0.0, // UNUSED: track.GetVelocity(), 00256 track.GetGlobalTime(), // Lab. 00257 track.GetProperTime(), // Part. 00258 &spin ) ; 00259 theFieldTrack.SetChargeAndMoments( particleCharge ); // EM moments -- future extension 00260 00261 G4int stepNo= track.GetCurrentStepNumber(); 00262 00263 ELimited limitedStep; 00264 G4FieldTrack endTrackState('a'); // Default values 00265 00266 fMassGeometryLimitedStep = false ; // default 00267 fAnyGeometryLimitedStep = false ; 00268 if( currentMinimumStep > 0 ) 00269 { 00270 G4double newMassSafety= 0.0; // temp. for recalculation 00271 00272 // Do the Transport in the field (non recti-linear) 00273 // 00274 lengthAlongCurve = fPathFinder->ComputeStep( theFieldTrack, 00275 currentMinimumStep, 00276 fNavigatorId, 00277 stepNo, 00278 newMassSafety, 00279 limitedStep, 00280 endTrackState, 00281 currentVolume ) ; 00282 // G4cout << " PathFinder ComputeStep returns " << lengthAlongCurve << G4endl; 00283 00284 G4double newFullSafety= fPathFinder->GetCurrentSafety(); 00285 // this was estimated already in step above 00286 // G4double newFullStep= fPathFinder->GetMinimumStep(); 00287 00288 if( limitedStep == kUnique || limitedStep == kSharedTransport ) 00289 { 00290 fMassGeometryLimitedStep = true ; 00291 } 00292 00293 fAnyGeometryLimitedStep = (fPathFinder->GetNumberGeometriesLimitingStep() != 0); 00294 00295 //#ifdef G4DEBUG_TRANSPORT 00296 if( fMassGeometryLimitedStep && !fAnyGeometryLimitedStep ) 00297 { 00298 G4cerr << " Error in determining geometries limiting the step" << G4endl; 00299 G4cerr << " Limiting: mass=" << fMassGeometryLimitedStep 00300 << " any= " << fAnyGeometryLimitedStep << G4endl; 00301 G4Exception("G4CoupledTransportation::AlongStepGetPhysicalInteractionLength()", 00302 "PathFinderConfused", FatalException, 00303 "Incompatible conditions - was limited by a geometry?"); 00304 } 00305 //#endif 00306 00307 // Other potential 00308 // fAnyGeometryLimitedStep = newFullStep < currentMinimumStep; 00309 // ^^^ Not good enough; 00310 // Must compare with maximum requested step size 00311 // (eg in case another process requested bigger, got this!) 00312 00313 geometryStepLength = std::min( lengthAlongCurve, currentMinimumStep); 00314 00315 // Momentum: Magnitude and direction can be changed too now ... 00316 // 00317 fMomentumChanged = true ; 00318 fTransportEndMomentumDir = endTrackState.GetMomentumDir() ; 00319 00320 // Remember last safety origin & value. 00321 fPreviousSftOrigin = startPosition ; 00322 fPreviousMassSafety = newMassSafety ; 00323 fPreviousFullSafety = newFullSafety ; 00324 // fpSafetyHelper->SetCurrentSafety( newFullSafety, startPosition); 00325 00326 #ifdef G4DEBUG_TRANSPORT 00327 if( fVerboseLevel > 1 ) 00328 { 00329 G4cout << "G4Transport:CompStep> " 00330 << " called the pathfinder for a new step at " << startPosition 00331 << " and obtained step = " << lengthAlongCurve << G4endl; 00332 G4cout << " New safety (preStep) = " << newMassSafety 00333 << " versus precalculated = " << startMassSafety << G4endl; 00334 } 00335 #endif 00336 00337 // Store as best estimate value 00338 startMassSafety = newMassSafety ; 00339 startFullSafety = newFullSafety ; 00340 00341 // Get the End-Position and End-Momentum (Dir-ection) 00342 fTransportEndPosition = endTrackState.GetPosition() ; 00343 fTransportEndKineticEnergy = endTrackState.GetKineticEnergy() ; 00344 } 00345 else 00346 { 00347 geometryStepLength = lengthAlongCurve= 0.0 ; 00348 fMomentumChanged = false ; 00349 // fMassGeometryLimitedStep = false ; // --- ??? 00350 // fAnyGeometryLimitedStep = true; 00351 fTransportEndMomentumDir = track.GetMomentumDirection(); 00352 fTransportEndKineticEnergy = track.GetKineticEnergy(); 00353 00354 fTransportEndPosition = startPosition; 00355 // If the step length requested is 0, and we are on a boundary 00356 // then a boundary will also limit the step. 00357 if( startMassSafety == 0.0 ) 00358 { 00359 fMassGeometryLimitedStep = true ; 00360 fAnyGeometryLimitedStep = true; 00361 } 00362 // TODO: Add explicit logical status for being at a boundary 00363 } 00364 // G4FieldTrack aTrackState(endTrackState); 00365 00366 if( !fieldExertsForce ) 00367 { 00368 fParticleIsLooping = false ; 00369 fMomentumChanged = false ; 00370 fEndGlobalTimeComputed = false ; 00371 // G4cout << " global time is false " << G4endl; 00372 } 00373 else 00374 { 00375 00376 #ifdef G4DEBUG_TRANSPORT 00377 if( fVerboseLevel > 1 ) 00378 { 00379 G4cout << " G4CT::CS End Position = " << fTransportEndPosition << G4endl; 00380 G4cout << " G4CT::CS End Direction = " << fTransportEndMomentumDir << G4endl; 00381 } 00382 #endif 00383 if( fFieldPropagator->GetCurrentFieldManager()->DoesFieldChangeEnergy() ) 00384 { 00385 // If the field can change energy, then the time must be integrated 00386 // - so this should have been updated 00387 // 00388 fCandidateEndGlobalTime = endTrackState.GetLabTimeOfFlight(); 00389 fEndGlobalTimeComputed = true; 00390 00391 // was ( fCandidateEndGlobalTime != track.GetGlobalTime() ); 00392 // a cleaner way is to have FieldTrack knowing whether time is updated. 00393 } 00394 else 00395 { 00396 // The energy should be unchanged by field transport, 00397 // - so the time changed will be calculated elsewhere 00398 // 00399 fEndGlobalTimeComputed = false; 00400 00401 // Check that the integration preserved the energy 00402 // - and if not correct this! 00403 G4double startEnergy= track.GetKineticEnergy(); 00404 G4double endEnergy= fTransportEndKineticEnergy; 00405 00406 static G4int no_inexact_steps=0; // , no_large_ediff; 00407 G4double absEdiff = std::fabs(startEnergy- endEnergy); 00408 if( absEdiff > perMillion * endEnergy ) 00409 { 00410 no_inexact_steps++; 00411 // Possible statistics keeping here ... 00412 } 00413 #ifdef G4VERBOSE 00414 if( (fVerboseLevel > 1) && ( absEdiff > perThousand * endEnergy) ) 00415 { 00416 ReportInexactEnergy(startEnergy, endEnergy); 00417 } // end of if (fVerboseLevel) 00418 #endif 00419 // Correct the energy for fields that conserve it 00420 // This - hides the integration error 00421 // - but gives a better physical answer 00422 fTransportEndKineticEnergy= track.GetKineticEnergy(); 00423 } 00424 } 00425 00426 endpointDistance = (fTransportEndPosition - startPosition).mag() ; 00427 fParticleIsLooping = fFieldPropagator->IsParticleLooping() ; 00428 00429 fTransportEndSpin = endTrackState.GetSpin(); 00430 00431 // Calculate the safety 00432 safetyProposal= startFullSafety; // used to be startMassSafety 00433 // Changed to accomodate processes that cannot update the safety -- JA 22 Nov 06 00434 00435 // Update safety for the end-point, if becomes negative at the end-point. 00436 00437 if( (startFullSafety < endpointDistance ) 00438 && ( particleCharge != 0.0 ) ) // Only needed to prepare for Mult Scat. 00439 // && !fAnyGeometryLimitedStep ) // To-Try: No safety update if at a boundary 00440 { 00441 G4double endFullSafety = 00442 fPathFinder->ComputeSafety( fTransportEndPosition); 00443 // Expected mission -- only mass geometry's safety 00444 // fMassNavigator->ComputeSafety( fTransportEndPosition) ; 00445 // Yet discrete processes only have poststep -- and this cannot 00446 // currently revise the safety 00447 // ==> so we use the all-geometry safety as a precaution 00448 00449 fpSafetyHelper->SetCurrentSafety( endFullSafety, fTransportEndPosition); 00450 // Pushing safety to Helper avoids recalculation at this point 00451 00452 G4ThreeVector centerPt= G4ThreeVector(0.0, 0.0, 0.0); // Used for return value 00453 G4double endMassSafety= fPathFinder->ObtainSafety( fNavigatorId, centerPt); 00454 // Retrieves the mass value from PathFinder (it calculated it) 00455 00456 fPreviousMassSafety = endMassSafety ; 00457 fPreviousFullSafety = endFullSafety; 00458 fPreviousSftOrigin = fTransportEndPosition ; 00459 00460 // The convention (Stepping Manager's) is safety from the start point 00461 // 00462 safetyProposal = endFullSafety + endpointDistance; 00463 // --> was endMassSafety 00464 // Changed to accomodate processes that cannot update the safety -- JA 22 Nov 06 00465 00466 // #define G4DEBUG_TRANSPORT 1 00467 00468 #ifdef G4DEBUG_TRANSPORT 00469 G4int prec= G4cout.precision(12) ; 00470 G4cout << "***Transportation::AlongStepGPIL ** " << G4endl ; 00471 G4cout << " Revised Safety at endpoint " << fTransportEndPosition 00472 << " give safety values: Mass= " << endMassSafety 00473 << " All= " << endFullSafety << G4endl ; 00474 G4cout << " Adding endpoint distance " << endpointDistance 00475 << " to obtain pseudo-safety= " << safetyProposal << G4endl ; 00476 G4cout.precision(prec); 00477 } 00478 else 00479 { 00480 G4int prec= G4cout.precision(12) ; 00481 G4cout << "***Transportation::AlongStepGPIL ** " << G4endl ; 00482 G4cout << " Quick Safety estimate at endpoint " << fTransportEndPosition 00483 << " gives safety endpoint value = " << startFullSafety - endpointDistance 00484 << " using start-point value " << startFullSafety 00485 << " and endpointDistance " << endpointDistance << G4endl; 00486 G4cout.precision(prec); 00487 #endif 00488 } 00489 00490 proposedSafetyForStart= safetyProposal; 00491 fParticleChange.ProposeTrueStepLength(geometryStepLength) ; 00492 00493 return geometryStepLength ; 00494 }
G4VParticleChange* G4CoupledTransportation::AtRestDoIt | ( | const G4Track & | , | |
const G4Step & | ||||
) | [inline, virtual] |
G4double G4CoupledTransportation::AtRestGetPhysicalInteractionLength | ( | const G4Track & | , | |
G4ForceCondition * | ||||
) | [inline, virtual] |
G4bool G4CoupledTransportation::DoesGlobalFieldExist | ( | ) | [inline, protected] |
Definition at line 52 of file G4CoupledTransportation.icc.
References G4FieldManager::DoesFieldExist(), G4TransportationManager::GetFieldManager(), and G4TransportationManager::GetTransportationManager().
Referenced by StartTracking().
00053 { 00054 G4TransportationManager* transportMgr; 00055 transportMgr= G4TransportationManager::GetTransportationManager(); 00056 00057 // fFieldExists= transportMgr->GetFieldManager()->DoesFieldExist(); 00058 // return fFieldExists; 00059 return transportMgr->GetFieldManager()->DoesFieldExist(); 00060 }
Definition at line 131 of file G4CoupledTransportation.icc.
00132 { 00133 G4bool lastValue= fUseMagneticMoment; 00134 fUseMagneticMoment= useMoment; 00135 return lastValue; 00136 }
void G4CoupledTransportation::EndTracking | ( | ) | [virtual] |
Reimplemented from G4VProcess.
Definition at line 900 of file G4CoupledTransportation.cc.
References G4TransportationManager::GetTransportationManager(), and G4TransportationManager::InactivateAll().
00901 { 00902 G4TransportationManager::GetTransportationManager()->InactivateAll(); 00903 }
G4double G4CoupledTransportation::GetMaxEnergyKilled | ( | ) | const [inline] |
G4PropagatorInField * G4CoupledTransportation::GetPropagatorInField | ( | ) | [inline] |
G4double G4CoupledTransportation::GetSumEnergyKilled | ( | ) | const [inline] |
G4double G4CoupledTransportation::GetThresholdImportantEnergy | ( | ) | const [inline] |
G4int G4CoupledTransportation::GetThresholdTrials | ( | ) | const [inline] |
G4double G4CoupledTransportation::GetThresholdWarningEnergy | ( | ) | const [inline] |
G4int G4CoupledTransportation::GetVerboseLevel | ( | ) | const [inline] |
G4VParticleChange * G4CoupledTransportation::PostStepDoIt | ( | const G4Track & | track, | |
const G4Step & | stepData | |||
) | [virtual] |
Implements G4VProcess.
Definition at line 677 of file G4CoupledTransportation.cc.
References G4PathFinder::CreateTouchableHandle(), fStopAndKill, G4cerr, G4cout, G4endl, G4VPhysicalVolume::GetLogicalVolume(), G4MaterialCutsCouple::GetMaterial(), G4LogicalVolume::GetMaterial(), G4ProductionCutsTable::GetMaterialCutsCouple(), G4LogicalVolume::GetMaterialCutsCouple(), G4Track::GetMomentumDirection(), G4Track::GetPosition(), G4MaterialCutsCouple::GetProductionCuts(), G4ProductionCutsTable::GetProductionCutsTable(), G4LogicalVolume::GetSensitiveDetector(), G4Track::GetTouchableHandle(), G4Track::GetTrackStatus(), G4PathFinder::Locate(), G4VParticleChange::ProposeLastStepInVolume(), G4VParticleChange::ProposeTrackStatus(), G4PathFinder::ReLocate(), ReportMove(), G4ParticleChangeForTransport::SetMaterialCutsCoupleInTouchable(), G4ParticleChangeForTransport::SetMaterialInTouchable(), G4ParticleChangeForTransport::SetSensitiveDetectorInTouchable(), and G4ParticleChangeForTransport::SetTouchableHandle().
00679 { 00680 G4TouchableHandle retCurrentTouchable ; // The one to return 00681 00682 // Initialize ParticleChange (by setting all its members equal 00683 // to corresponding members in G4Track) 00684 // fParticleChange.Initialize(track) ; // To initialise TouchableChange 00685 00686 fParticleChange.ProposeTrackStatus(track.GetTrackStatus()) ; 00687 00688 // Check that the end position and direction are preserved 00689 // since call to AlongStepDoIt 00690 00691 #ifdef G4DEBUG_TRANSPORT 00692 if( ( fVerboseLevel > 0 ) && ((fTransportEndPosition - track.GetPosition()).mag2() >= 1.0e-16) ) 00693 { 00694 ReportMove( track.GetPosition(), fTransportEndPosition, "End of Step Position" ); 00695 G4cerr << " Problem in G4CoupledTransportation::PostStepDoIt " << G4endl; 00696 } 00697 00698 // If the Step was determined by the volume boundary, relocate the particle 00699 // The pathFinder will know that the geometry limited the step (!?) 00700 00701 if( fVerboseLevel > 0 ) 00702 { 00703 G4cout << " Calling PathFinder::Locate() from " 00704 << " G4CoupledTransportation::PostStepDoIt " << G4endl; 00705 G4cout << " fAnyGeometryLimitedStep is " << fAnyGeometryLimitedStep << G4endl; 00706 00707 } 00708 #endif 00709 00710 if(fAnyGeometryLimitedStep) 00711 { 00712 fPathFinder->Locate( track.GetPosition(), 00713 track.GetMomentumDirection(), 00714 true); 00715 00716 // fCurrentTouchable will now become the previous touchable, 00717 // and what was the previous will be freed. 00718 // (Needed because the preStepPoint can point to the previous touchable) 00719 00720 fCurrentTouchableHandle= 00721 fPathFinder->CreateTouchableHandle( fNavigatorId ); 00722 00723 #ifdef G4DEBUG_TRANSPORT 00724 if( fVerboseLevel > 0 ) 00725 { 00726 G4cout << "G4CoupledTransportation::PostStepDoIt --- fNavigatorId = " 00727 << fNavigatorId << G4endl; 00728 } 00729 if( fVerboseLevel > 1 ) 00730 { 00731 G4VPhysicalVolume* vol= fCurrentTouchableHandle->GetVolume(); 00732 G4cout << "CHECK !!!!!!!!!!! fCurrentTouchableHandle->GetVolume() = " << vol; 00733 if( vol ) { G4cout << "Name=" << vol->GetName(); } 00734 G4cout << G4endl; 00735 } 00736 #endif 00737 00738 // Check whether the particle is out of the world volume 00739 // If so it has exited and must be killed. 00740 // 00741 if( fCurrentTouchableHandle->GetVolume() == 0 ) 00742 { 00743 fParticleChange.ProposeTrackStatus( fStopAndKill ) ; 00744 } 00745 retCurrentTouchable = fCurrentTouchableHandle ; 00746 // fParticleChange.SetTouchableHandle( fCurrentTouchableHandle ) ; 00747 00748 // Notify particle change that this is last step in volume 00749 fParticleChange.ProposeLastStepInVolume(true); 00750 } 00751 else // fAnyGeometryLimitedStep is false 00752 { 00753 #ifdef G4DEBUG_TRANSPORT 00754 if( fVerboseLevel > 1 ) 00755 { 00756 G4cout << "G4CoupledTransportation::PostStepDoIt -- " 00757 << " fAnyGeometryLimitedStep = " << fAnyGeometryLimitedStep 00758 << " must be false " << G4endl; 00759 } 00760 #endif 00761 // This serves only to move each of the Navigator's location 00762 // 00763 // fLinearNavigator->LocateGlobalPointWithinVolume( track.GetPosition() ) ; 00764 00765 // G4cout << "G4CoupledTransportation calling PathFinder::ReLocate() " << G4endl; 00766 fPathFinder->ReLocate( track.GetPosition() ); 00767 // track.GetMomentumDirection() ); 00768 00769 // Keep the value of the track's current Touchable is retained, 00770 // and use it to overwrite the (unset) one in particle change. 00771 // Expect this must be fCurrentTouchable too 00772 // - could it be different, eg at the start of a step ? 00773 // 00774 retCurrentTouchable = track.GetTouchableHandle() ; 00775 // fParticleChange.SetTouchableHandle( track.GetTouchableHandle() ) ; 00776 00777 // Have not reached a boundary 00778 fParticleChange.ProposeLastStepInVolume(false); 00779 } // endif ( fAnyGeometryLimitedStep ) 00780 00781 const G4VPhysicalVolume* pNewVol = retCurrentTouchable->GetVolume() ; 00782 const G4Material* pNewMaterial = 0 ; 00783 const G4VSensitiveDetector* pNewSensitiveDetector = 0 ; 00784 00785 if( pNewVol != 0 ) 00786 { 00787 pNewMaterial= pNewVol->GetLogicalVolume()->GetMaterial(); 00788 pNewSensitiveDetector= pNewVol->GetLogicalVolume()->GetSensitiveDetector(); 00789 } 00790 00791 // ( const_cast<G4Material *> pNewMaterial ) ; 00792 // ( const_cast<G4VSensitiveDetetor *> pNewSensitiveDetector) ; 00793 00794 fParticleChange.SetMaterialInTouchable( (G4Material *) pNewMaterial ) ; 00795 fParticleChange.SetSensitiveDetectorInTouchable( (G4VSensitiveDetector *) pNewSensitiveDetector ) ; 00796 // "temporarily" until Get/Set Material of ParticleChange, 00797 // and StepPoint can be made const. 00798 00799 const G4MaterialCutsCouple* pNewMaterialCutsCouple = 0; 00800 if( pNewVol != 0 ) 00801 { 00802 pNewMaterialCutsCouple=pNewVol->GetLogicalVolume()->GetMaterialCutsCouple(); 00803 if( pNewMaterialCutsCouple!=0 00804 && pNewMaterialCutsCouple->GetMaterial()!=pNewMaterial ) 00805 { 00806 // for parametrized volume 00807 // 00808 pNewMaterialCutsCouple = 00809 G4ProductionCutsTable::GetProductionCutsTable() 00810 ->GetMaterialCutsCouple(pNewMaterial, 00811 pNewMaterialCutsCouple->GetProductionCuts()); 00812 } 00813 } 00814 fParticleChange.SetMaterialCutsCoupleInTouchable( pNewMaterialCutsCouple ); 00815 00816 // Must always set the touchable in ParticleChange, whether relocated or not 00817 fParticleChange.SetTouchableHandle(retCurrentTouchable) ; 00818 00819 return &fParticleChange ; 00820 }
G4double G4CoupledTransportation::PostStepGetPhysicalInteractionLength | ( | const G4Track & | , | |
G4double | previousStepSize, | |||
G4ForceCondition * | pForceCond | |||
) | [virtual] |
Implements G4VProcess.
Definition at line 651 of file G4CoupledTransportation.cc.
References DBL_MAX, and Forced.
00654 { 00655 // Must act as PostStep action -- to relocate particle 00656 *pForceCond = Forced ; 00657 return DBL_MAX ; 00658 }
void G4CoupledTransportation::ReportInexactEnergy | ( | G4double | startEnergy, | |
G4double | endEnergy | |||
) | [protected] |
Definition at line 907 of file G4CoupledTransportation.cc.
References G4cerr, G4cout, and G4endl.
Referenced by AlongStepGetPhysicalInteractionLength().
00908 { 00909 static G4int no_warnings= 0, warnModulo=1, moduloFactor= 10, no_large_ediff= 0; 00910 00911 if( std::fabs(startEnergy- endEnergy) > perThousand * endEnergy ) 00912 { 00913 no_large_ediff ++; 00914 if( (no_large_ediff% warnModulo) == 0 ) 00915 { 00916 no_warnings++; 00917 G4cout << "WARNING - G4CoupledTransportation::AlongStepGetPIL() " 00918 << " Energy change in Step is above 1^-3 relative value. " << G4endl 00919 << " Relative change in 'tracking' step = " 00920 << std::setw(15) << (endEnergy-startEnergy)/startEnergy << G4endl 00921 << " Starting E= " << std::setw(12) << startEnergy / MeV << " MeV " << G4endl 00922 << " Ending E= " << std::setw(12) << endEnergy / MeV << " MeV " << G4endl; 00923 G4cout << " Energy has been corrected -- however, review" 00924 << " field propagation parameters for accuracy." << G4endl; 00925 if( (fVerboseLevel > 2 ) || (no_warnings<4) || (no_large_ediff == warnModulo * moduloFactor) ) 00926 { 00927 G4cout << " These include EpsilonStepMax(/Min) in G4FieldManager " 00928 << " which determine fractional error per step for integrated quantities. " << G4endl 00929 << " Note also the influence of the permitted number of integration steps." 00930 << G4endl; 00931 } 00932 G4cerr << "ERROR - G4CoupledTransportation::AlongStepGetPIL()" << G4endl 00933 << " Bad 'endpoint'. Energy change detected" 00934 << " and corrected. " 00935 << " Has occurred already " 00936 << no_large_ediff << " times." << G4endl; 00937 if( no_large_ediff == warnModulo * moduloFactor ) 00938 { 00939 warnModulo *= moduloFactor; 00940 } 00941 } 00942 } 00943 }
void G4CoupledTransportation::ReportMove | ( | G4ThreeVector | OldVector, | |
G4ThreeVector | NewVector, | |||
const G4String & | Quantity | |||
) | [protected] |
Definition at line 661 of file G4CoupledTransportation.cc.
References G4cerr, and G4endl.
Referenced by PostStepDoIt().
00662 { 00663 G4ThreeVector moveVec = ( NewVector - OldVector ); 00664 00665 G4cerr << G4endl 00666 << "**************************************************************" << G4endl; 00667 G4cerr << "Endpoint has moved between value expected from TransportEndPosition " 00668 << " and value from Track in PostStepDoIt. " << G4endl 00669 << "Change of " << Quantity << " is " << moveVec.mag() / mm << " mm long, " 00670 << " and its vector is " << (1.0/mm) * moveVec << " mm " << G4endl 00671 << "Endpoint of ComputeStep was " << OldVector 00672 << " and current position to locate is " << NewVector << G4endl; 00673 }
void G4CoupledTransportation::ResetKilledStatistics | ( | G4int | report = 1 |
) | [inline] |
Definition at line 117 of file G4CoupledTransportation.icc.
References G4cout, and G4endl.
00118 { 00119 if( report ) { 00120 G4cout << " G4CoupledTransportation: Statistics for looping particles " << G4endl; 00121 G4cout << " Sum of energy of loopers killed: " << fSumEnergyKilled << G4endl; 00122 G4cout << " Max energy of loopers killed: " << fMaxEnergyKilled << G4endl; 00123 } 00124 00125 fSumEnergyKilled= 0; 00126 fMaxEnergyKilled= -1.0*CLHEP::GeV; 00127 }
void G4CoupledTransportation::SetPropagatorInField | ( | G4PropagatorInField * | pFieldPropagator | ) | [inline] |
void G4CoupledTransportation::SetThresholdImportantEnergy | ( | G4double | newEnImp | ) | [inline] |
void G4CoupledTransportation::SetThresholdTrials | ( | G4int | newMaxTrials | ) | [inline] |
void G4CoupledTransportation::SetThresholdWarningEnergy | ( | G4double | newEnWarn | ) | [inline] |
void G4CoupledTransportation::SetVerboseLevel | ( | G4int | verboseLevel | ) | [inline] |
void G4CoupledTransportation::StartTracking | ( | G4Track * | aTrack | ) | [virtual] |
Reimplemented from G4VProcess.
Definition at line 828 of file G4CoupledTransportation.cc.
References G4TransportationManager::ActivateNavigator(), G4FieldManagerStore::ClearAllChordFindersState(), G4PropagatorInField::ClearPropagatorState(), DoesGlobalFieldExist(), G4cout, G4endl, G4PropagatorInField::GetChordFinder(), G4FieldManagerStore::GetInstance(), G4Track::GetMomentumDirection(), G4TransportationManager::GetNavigatorForTracking(), G4Track::GetPosition(), G4Track::GetTouchableHandle(), G4TransportationManager::GetTransportationManager(), position, G4PathFinder::PrepareNewTrack(), and G4ChordFinder::ResetStepEstimate().
00829 { 00830 00831 G4TransportationManager* transportMgr = 00832 G4TransportationManager::GetTransportationManager(); 00833 00834 // G4VProcess::StartTracking(aTrack); 00835 00836 // The 'initialising' actions 00837 // once taken in AlongStepGPIL -- if ( track.GetCurrentStepNumber()==1 ) 00838 00839 // fStartedNewTrack= true; 00840 00841 fMassNavigator = transportMgr->GetNavigatorForTracking() ; 00842 fNavigatorId= transportMgr->ActivateNavigator( fMassNavigator ); // Confirm it! 00843 00844 // if( fVerboseLevel > 1 ){ 00845 // G4cout << " Navigator Id obtained in StartTracking " << fNavigatorId << G4endl; 00846 // } 00847 G4ThreeVector position = aTrack->GetPosition(); 00848 G4ThreeVector direction = aTrack->GetMomentumDirection(); 00849 00850 // if( fVerboseLevel > 1 ){ 00851 // G4cout << " Calling PathFinder::PrepareNewTrack from " 00852 // << " G4CoupledTransportation::StartTracking -- which calls Locate()" << G4endl; 00853 // } 00854 fPathFinder->PrepareNewTrack( position, direction); 00855 // This implies a call to fPathFinder->Locate( position, direction ); 00856 00857 // Global field, if any, must exist before tracking is started 00858 fGlobalFieldExists= DoesGlobalFieldExist(); 00859 // reset safety value and center 00860 // 00861 fPreviousMassSafety = 0.0 ; 00862 fPreviousFullSafety = 0.0 ; 00863 fPreviousSftOrigin = G4ThreeVector(0.,0.,0.) ; 00864 00865 // reset looping counter -- for motion in field 00866 fNoLooperTrials= 0; 00867 // Must clear this state .. else it depends on last track's value 00868 // --> a better solution would set this from state of suspended track TODO ? 00869 // Was if( aTrack->GetCurrentStepNumber()==1 ) { .. } 00870 00871 // ChordFinder reset internal state 00872 // 00873 if( fGlobalFieldExists ) 00874 { 00875 fFieldPropagator->ClearPropagatorState(); 00876 // Resets safety values, in case of overlaps. 00877 00878 G4ChordFinder* chordF= fFieldPropagator->GetChordFinder(); 00879 if( chordF ) { chordF->ResetStepEstimate(); } 00880 } 00881 00882 // Clear the chord finders of all fields (ie managers) derived objects 00883 // 00884 G4FieldManagerStore* fieldMgrStore = G4FieldManagerStore::GetInstance(); 00885 fieldMgrStore->ClearAllChordFindersState(); 00886 00887 #ifdef G4DEBUG_TRANSPORT 00888 if( fVerboseLevel > 1 ) 00889 { 00890 G4cout << " Returning touchable handle " << fCurrentTouchableHandle << G4endl; 00891 } 00892 #endif 00893 00894 // Update the current touchable handle (from the track's) 00895 // 00896 fCurrentTouchableHandle = aTrack->GetTouchableHandle(); 00897 }