From 2c8b209a787fc4534462960ddc8723a0629140ed Mon Sep 17 00:00:00 2001
From: Maximilian Reininghaus <maximilian.reininghaus@kit.edu>
Date: Thu, 6 Jun 2019 14:58:22 -0300
Subject: [PATCH] added radiative losses with constant b
---
Documentation/Examples/stopping_power.cc | 2 +-
Processes/EnergyLoss/EnergyLoss.cc | 22 ++++++++++++++++++----
Processes/EnergyLoss/EnergyLoss.h | 8 ++++++--
3 files changed, 25 insertions(+), 7 deletions(-)
diff --git a/Documentation/Examples/stopping_power.cc b/Documentation/Examples/stopping_power.cc
index 0ca9f8be4..edbd2bc32 100644
--- a/Documentation/Examples/stopping_power.cc
+++ b/Documentation/Examples/stopping_power.cc
@@ -76,7 +76,7 @@ int main() {
beamCode, E0, plab, pos, 0_ns});
auto const p = stack.GetNextParticle();
- HEPEnergyType dE = eLoss.BetheBloch(p, 1_g / square(1_cm));
+ HEPEnergyType dE = eLoss.TotalEnergyLoss(p, 1_g / square(1_cm));
file << P0 / mass << "\t" << -dE / 1_eV << std::endl;
}
}
diff --git a/Processes/EnergyLoss/EnergyLoss.cc b/Processes/EnergyLoss/EnergyLoss.cc
index b49030a57..febabc776 100644
--- a/Processes/EnergyLoss/EnergyLoss.cc
+++ b/Processes/EnergyLoss/EnergyLoss.cc
@@ -146,6 +146,20 @@ namespace corsika::process::energy_loss {
(0.5 * log(aux) - beta2 - Cadj / Z - delta / 2 + barkas + bloch) * dX;
}
+ // radiation losses according to PDG 2018, ch. 33 ref. [5]
+ HEPEnergyType EnergyLoss::RadiationLosses(SetupParticle const& vP,
+ GrammageType const vDX) {
+ // simple-minded hard-coded value for b(E) inspired by data from
+ // http://pdg.lbl.gov/2018/AtomicNuclearProperties/ for N and O.
+ auto constexpr b = 3.0 * 1e-6 * square(1_cm) / 1_g;
+ return -vP.GetEnergy() * b * vDX;
+ }
+
+ HEPEnergyType EnergyLoss::TotalEnergyLoss(SetupParticle const& vP,
+ GrammageType const vDX) {
+ return BetheBloch(vP, vDX) + RadiationLosses(vP, vDX);
+ }
+
process::EProcessReturn EnergyLoss::DoContinuous(SetupParticle& p,
SetupTrack const& t) {
if (p.GetChargeNumber() == 0) return process::EProcessReturn::eOk;
@@ -153,7 +167,7 @@ namespace corsika::process::energy_loss {
p.GetNode()->GetModelProperties().IntegratedGrammage(t, t.GetLength());
cout << "EnergyLoss " << p.GetPID() << ", z=" << p.GetChargeNumber()
<< ", dX=" << dX / 1_g * square(1_cm) << "g/cm2" << endl;
- HEPEnergyType dE = BetheBloch(p, dX);
+ HEPEnergyType dE = TotalEnergyLoss(p, dX);
auto E = p.GetEnergy();
const auto Ekin = E - p.GetMass();
auto Enew = E + dE;
@@ -173,14 +187,14 @@ namespace corsika::process::energy_loss {
return status;
}
- units::si::LengthType EnergyLoss::MaxStepLength(SetupParticle const& vParticle,
- SetupTrack const& vTrack) const {
+ LengthType EnergyLoss::MaxStepLength(SetupParticle const& vParticle,
+ SetupTrack const& vTrack) const {
if (vParticle.GetChargeNumber() == 0) {
return units::si::meter * std::numeric_limits<double>::infinity();
}
auto constexpr dX = 1_g / square(1_cm);
- auto const dE = -BetheBloch(vParticle, dX); // dE > 0
+ auto const dE = -TotalEnergyLoss(vParticle, dX); // dE > 0
//~ auto const Ekin = vParticle.GetEnergy() - vParticle.GetMass();
auto const maxLoss = 0.01 * vParticle.GetEnergy();
auto const maxGrammage = maxLoss / dE * dX;
diff --git a/Processes/EnergyLoss/EnergyLoss.h b/Processes/EnergyLoss/EnergyLoss.h
index 36d1ee6b4..4380a8ddc 100644
--- a/Processes/EnergyLoss/EnergyLoss.h
+++ b/Processes/EnergyLoss/EnergyLoss.h
@@ -38,8 +38,12 @@ namespace corsika::process::energy_loss {
units::si::HEPEnergyType GetTotal() const { return fEnergyLossTot; }
void PrintProfile() const;
- static units::si::HEPEnergyType BetheBloch(setup::Stack::ParticleType const& p,
- const units::si::GrammageType dX);
+ static units::si::HEPEnergyType BetheBloch(setup::Stack::ParticleType const&,
+ const units::si::GrammageType);
+ static units::si::HEPEnergyType RadiationLosses(setup::Stack::ParticleType const&,
+ const units::si::GrammageType);
+ static units::si::HEPEnergyType TotalEnergyLoss(setup::Stack::ParticleType const&,
+ const units::si::GrammageType);
private:
int GetXbin(setup::Stack::ParticleType const&, setup::Trajectory const&,
--
GitLab