11 #include <simulation/kernel/ExtEnergyLoss.h>
12 #include <simulation/kernel/EnergyLossForExtrapolator.h>
14 #include <G4ErrorPropagatorData.hh>
16 #include <framework/logging/Logger.h>
20 using namespace Belle2::Simulation;
22 ExtEnergyLoss::ExtEnergyLoss(
const G4String& processName, G4ProcessType type)
23 : G4VContinuousProcess(processName, type), m_energyLossForExtrapolator(NULL)
25 B2DEBUG(200,
"ExtEnergyLoss is created");
48 aParticleChange.Initialize(aTrack);
50 G4ErrorPropagatorData* g4edata = G4ErrorPropagatorData::GetErrorPropagatorData();
52 G4double kinEnergyStart = aTrack.GetKineticEnergy();
53 G4double step_length = aStep.GetStepLength();
55 const G4Material* aMaterial = aTrack.GetMaterial();
56 const G4ParticleDefinition* aParticleDef = aTrack.GetDynamicParticle()->GetDefinition();
57 G4double kinEnergyEnd = kinEnergyStart;
59 if (g4edata->GetMode() == G4ErrorMode(G4ErrorMode_PropBackwards)) {
64 G4double kinEnergyHalfStep = (kinEnergyStart + kinEnergyEnd) * 0.5;
66 B2DEBUG(200,
"ExtEnergyLoss::AlongStepDoIt() BWD end " << kinEnergyEnd <<
" halfstep " << kinEnergyHalfStep);
73 kinEnergyEnd = kinEnergyStart - (kinEnergyHalfStep - kinEnergyEnd);
74 }
else if (g4edata->GetMode() == G4ErrorMode(G4ErrorMode_PropForwards)) {
80 G4double kinEnergyHalfStep = (kinEnergyStart + kinEnergyEnd) * 0.5;
81 B2DEBUG(200,
"ExtEnergyLoss::AlongStepDoIt() FWD end " << kinEnergyEnd <<
" halfstep " << kinEnergyHalfStep);
88 kinEnergyEnd = kinEnergyStart - (kinEnergyHalfStep - kinEnergyEnd);
91 G4double edepo = kinEnergyEnd - kinEnergyStart;
93 B2DEBUG(300,
"ExtEnergyLoss::AlongStepDoIt() Estart= " << kinEnergyStart <<
" Eend " << kinEnergyEnd
94 <<
" Ediff " << -edepo <<
" step= " << step_length <<
" mate= " << aMaterial->GetName()
95 <<
" particle= " << aParticleDef->GetParticleName());
97 aParticleChange.ClearDebugFlag();
98 aParticleChange.ProposeLocalEnergyDeposit(edepo);
99 aParticleChange.SetNumberOfSecondaries(0);
101 aParticleChange.ProposeEnergy(kinEnergyEnd);
103 return &aParticleChange;
108 G4double currentMinimumStep,
111 G4double step = DBL_MAX;
113 G4double kinEnergyStart = aTrack.GetKineticEnergy();
114 G4double kinEnergyLoss = kinEnergyStart;
115 const G4Material* aMaterial = aTrack.GetMaterial();
116 const G4ParticleDefinition* aParticleDef = aTrack.GetDynamicParticle()->GetDefinition();
117 G4ErrorPropagatorData* g4edata = G4ErrorPropagatorData::GetErrorPropagatorData();
118 if (g4edata->GetMode() == G4ErrorMode(G4ErrorMode_PropBackwards)) {
119 kinEnergyLoss = - kinEnergyStart +
121 }
else if (g4edata->GetMode() == G4ErrorMode(G4ErrorMode_PropForwards)) {
122 kinEnergyLoss = kinEnergyStart -
125 B2DEBUG(300,
"ExtEnergyLoss::GetContinuousStepLimit() currentMinimumStep " << currentMinimumStep
126 <<
" kinEnergyLoss " << kinEnergyLoss <<
" kinEnergyStart " << kinEnergyStart);
127 if (kinEnergyLoss / kinEnergyStart >
m_StepLimit) {
128 step =
m_StepLimit / (kinEnergyLoss / kinEnergyStart) * currentMinimumStep;
129 B2DEBUG(300,
"ExtEnergyLoss::GetContinuousStepLimit() limiting Step " << step
130 <<
" energy loss fraction " << kinEnergyLoss / kinEnergyStart <<
" > " <<
m_StepLimit);