Belle II Software development
ExtEnergyLoss.cc
1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * *
5 * See git log for contributors and copyright holders. *
6 * This file is licensed under LGPL-3.0, see LICENSE.md. *
7 **************************************************************************/
8
9#include <simulation/kernel/ExtEnergyLoss.h>
10#include <simulation/kernel/EnergyLossForExtrapolator.h>
11
12#include <G4ErrorPropagatorData.hh>
13
14#include <framework/logging/Logger.h>
15
16using namespace std;
17using namespace Belle2;
18using namespace Belle2::Simulation;
19
20ExtEnergyLoss::ExtEnergyLoss(const G4String& processName, G4ProcessType type)
21 : G4VContinuousProcess(processName, type), m_energyLossForExtrapolator(NULL)
22{
23 B2DEBUG(200, "ExtEnergyLoss is created");
24 if (m_energyLossForExtrapolator == NULL) {
26 }
27 m_StepLimit = 1.0; // fraction of kinetic energy that could be lost in one step
28 if (false) {
29 G4Track aTrack;
30 G4Step aStep;
31 double dummy;
32 AlongStepDoIt(aTrack, aStep);
33 GetContinuousStepLimit(aTrack, 0.0, 0.0, dummy);
34 }
35}
36
38{
39 if (m_energyLossForExtrapolator != NULL) {
41 }
42}
43
44G4VParticleChange* ExtEnergyLoss::AlongStepDoIt(const G4Track& aTrack, const G4Step& aStep)
45{
46 aParticleChange.Initialize(aTrack);
47
48 G4ErrorPropagatorData* g4edata = G4ErrorPropagatorData::GetErrorPropagatorData();
49
50 G4double kinEnergyStart = aTrack.GetKineticEnergy();
51 G4double step_length = aStep.GetStepLength();
52
53 const G4Material* aMaterial = aTrack.GetMaterial();
54 const G4ParticleDefinition* aParticleDef = aTrack.GetDynamicParticle()->GetDefinition();
55 G4double kinEnergyEnd = kinEnergyStart;
56
57 if (g4edata->GetMode() == G4ErrorMode(G4ErrorMode_PropBackwards)) {
58 kinEnergyEnd = m_energyLossForExtrapolator->EnergyBeforeStep(kinEnergyStart,
59 step_length,
60 aMaterial,
61 aParticleDef);
62 G4double kinEnergyHalfStep = (kinEnergyStart + kinEnergyEnd) * 0.5;
63
64 B2DEBUG(200, "ExtEnergyLoss::AlongStepDoIt() BWD end " << kinEnergyEnd << " halfstep " << kinEnergyHalfStep);
65
66 // rescale to energy lost at midpoint of step
67 kinEnergyEnd = m_energyLossForExtrapolator->EnergyBeforeStep(kinEnergyHalfStep,
68 step_length,
69 aMaterial,
70 aParticleDef);
71 kinEnergyEnd = kinEnergyStart - (kinEnergyHalfStep - kinEnergyEnd);
72 } else if (g4edata->GetMode() == G4ErrorMode(G4ErrorMode_PropForwards)) {
73
74 kinEnergyEnd = m_energyLossForExtrapolator->EnergyAfterStep(kinEnergyStart,
75 step_length,
76 aMaterial,
77 aParticleDef);
78 G4double kinEnergyHalfStep = (kinEnergyStart + kinEnergyEnd) * 0.5;
79 B2DEBUG(200, "ExtEnergyLoss::AlongStepDoIt() FWD end " << kinEnergyEnd << " halfstep " << kinEnergyHalfStep);
80
81 // rescale to energy lost at midpoint of step
82 kinEnergyEnd = m_energyLossForExtrapolator->EnergyAfterStep(kinEnergyHalfStep,
83 step_length,
84 aMaterial,
85 aParticleDef);
86 kinEnergyEnd = kinEnergyStart - (kinEnergyHalfStep - kinEnergyEnd);
87 }
88
89 G4double edepo = kinEnergyEnd - kinEnergyStart;
90
91 B2DEBUG(300, "ExtEnergyLoss::AlongStepDoIt() Estart= " << kinEnergyStart << " Eend " << kinEnergyEnd
92 << " Ediff " << -edepo << " step= " << step_length << " mate= " << aMaterial->GetName()
93 << " particle= " << aParticleDef->GetParticleName());
94
95 aParticleChange.ClearDebugFlag();
96 aParticleChange.ProposeLocalEnergyDeposit(edepo);
97 aParticleChange.SetNumberOfSecondaries(0);
98
99 aParticleChange.ProposeEnergy(kinEnergyEnd);
100
101 return &aParticleChange;
102}
103
104G4double ExtEnergyLoss::GetContinuousStepLimit(const G4Track& aTrack,
105 G4double,
106 G4double currentMinimumStep,
107 G4double&)
108{
109 G4double step = DBL_MAX;
110 if (m_StepLimit < 1.0) {
111 G4double kinEnergyStart = aTrack.GetKineticEnergy();
112 G4double kinEnergyLoss = kinEnergyStart;
113 const G4Material* aMaterial = aTrack.GetMaterial();
114 const G4ParticleDefinition* aParticleDef = aTrack.GetDynamicParticle()->GetDefinition();
115 G4ErrorPropagatorData* g4edata = G4ErrorPropagatorData::GetErrorPropagatorData();
116 if (g4edata->GetMode() == G4ErrorMode(G4ErrorMode_PropBackwards)) {
117 kinEnergyLoss = - kinEnergyStart +
118 m_energyLossForExtrapolator->EnergyBeforeStep(kinEnergyStart, currentMinimumStep, aMaterial, aParticleDef);
119 } else if (g4edata->GetMode() == G4ErrorMode(G4ErrorMode_PropForwards)) {
120 kinEnergyLoss = kinEnergyStart -
121 m_energyLossForExtrapolator->EnergyAfterStep(kinEnergyStart, currentMinimumStep, aMaterial, aParticleDef);
122 }
123 B2DEBUG(300, "ExtEnergyLoss::GetContinuousStepLimit() currentMinimumStep " << currentMinimumStep
124 << " kinEnergyLoss " << kinEnergyLoss << " kinEnergyStart " << kinEnergyStart);
125 if (kinEnergyLoss / kinEnergyStart > m_StepLimit) {
126 step = m_StepLimit / (kinEnergyLoss / kinEnergyStart) * currentMinimumStep;
127 B2DEBUG(300, "ExtEnergyLoss::GetContinuousStepLimit() limiting Step " << step
128 << " energy loss fraction " << kinEnergyLoss / kinEnergyStart << " > " << m_StepLimit);
129 }
130 }
131
132 return step;
133
134}
135
Calculate energy loss, fluctuation, and multiple-scattering angle for extrapolator.
G4double EnergyAfterStep(G4double kinEnergy, G4double step, const G4Material *, const G4ParticleDefinition *)
Get kinetic energy after a step in given material by given particle.
G4double EnergyBeforeStep(G4double kinEnergy, G4double step, const G4Material *, const G4ParticleDefinition *)
Get kinetic energy before a step in given material by given particle.
G4double GetContinuousStepLimit(const G4Track &, G4double, G4double, G4double &)
Gets step limit for the particle being swum.
EnergyLossForExtrapolator * m_energyLossForExtrapolator
Pointer to the geant4e-specific energy-loss and mult-scat class.
Definition: ExtEnergyLoss.h:60
G4double m_StepLimit
Step limit for this process (fraction of KE that could be lost in one step)
Definition: ExtEnergyLoss.h:63
G4VParticleChange * AlongStepDoIt(const G4Track &, const G4Step &)
Apply energy loss process along the step.
virtual ~ExtEnergyLoss()
destructor
ExtEnergyLoss(const G4String &processName="ExtEnergyLoss", G4ProcessType aType=fElectromagnetic)
constructor
Abstract base class for different kinds of events.
STL namespace.