From 284a7906dfefa13a711a121e95914a77fe8458d6 Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Mon, 18 Jan 2021 18:15:50 +0000
Subject: [PATCH] removed explicit threshold from energy loss BetheBlochPDG
---
corsika/detail/modules/energy_loss/BetheBlochPDG.inl | 10 +++++-----
corsika/modules/energy_loss/BetheBlochPDG.hpp | 7 ++++---
examples/cascade_example.cpp | 2 +-
examples/cascade_proton_example.cpp | 2 +-
examples/hybrid_MC.cpp | 2 +-
examples/stopping_power.cpp | 2 +-
6 files changed, 13 insertions(+), 12 deletions(-)
diff --git a/corsika/detail/modules/energy_loss/BetheBlochPDG.inl b/corsika/detail/modules/energy_loss/BetheBlochPDG.inl
index 776531540..1f3d0e930 100644
--- a/corsika/detail/modules/energy_loss/BetheBlochPDG.inl
+++ b/corsika/detail/modules/energy_loss/BetheBlochPDG.inl
@@ -26,11 +26,10 @@ namespace corsika {
return sqrt((Elab - m) * (Elab + m));
};
- BetheBlochPDG::BetheBlochPDG(ShowerAxis const& shower_axis, HEPEnergyType emCut)
+ BetheBlochPDG::BetheBlochPDG(ShowerAxis const& shower_axis)
: dX_(10_g / square(1_cm)) // profile binning
, dX_threshold_(0.0001_g / square(1_cm))
, shower_axis_(shower_axis)
- , emCut_(emCut)
, profile_(int(shower_axis.getMaximumX() / dX_) + 1) {}
HEPEnergyType BetheBlochPDG::getBetheBloch(setup::Stack::particle_type const& p,
@@ -170,8 +169,9 @@ namespace corsika {
auto const dEdX = -getTotalEnergyLoss(vParticle, dX) / dX; // dE > 0
//~ auto const Ekin = vParticle.getEnergy() - vParticle.getMass();
- // in any case: never go below 0.99*emCut_ This needs to be
- // slightly smaller than emCut_ since, either this Step is limited
+ auto const emCut = get_energy_threshold(vParticle.getPID());
+ // in any case: never go below 0.99*emCut This needs to be
+ // slightly smaller than emCut since, either this Step is limited
// by energy_lim, then the particle is stopped in a very short
// range (before doing anythin else) and is then removed
// instantly. The exact position where it reaches emCut is not
@@ -179,7 +179,7 @@ namespace corsika {
// afterwards.
//
const auto energy = vParticle.getEnergy();
- auto energy_lim = std::max(0.9 * energy, 0.99 * emCut_);
+ auto energy_lim = std::max(0.9 * energy, 0.99 * emCut);
auto const maxGrammage = (energy - energy_lim) / dEdX;
diff --git a/corsika/modules/energy_loss/BetheBlochPDG.hpp b/corsika/modules/energy_loss/BetheBlochPDG.hpp
index 3e64567b9..f8e875937 100644
--- a/corsika/modules/energy_loss/BetheBlochPDG.hpp
+++ b/corsika/modules/energy_loss/BetheBlochPDG.hpp
@@ -41,11 +41,13 @@ namespace corsika {
using MeVgcm2 = decltype(1e6 * electronvolt / gram * square(1e-2 * meter));
public:
- BetheBlochPDG(ShowerAxis const& showerAxis, HEPEnergyType emCut);
+ BetheBlochPDG(ShowerAxis const& showerAxis);
ProcessReturn doContinuous(setup::Stack::particle_type&, setup::Trajectory const&);
LengthType getMaxStepLength(setup::Stack::particle_type const&,
- setup::Trajectory const&) const;
+ setup::Trajectory const&)
+ const; //! limited by the energy threshold! By default the limit is the particle
+ //! rest mass, i.e. kinetic energy is zero
static HEPEnergyType getBetheBloch(setup::Stack::particle_type const&,
const GrammageType);
static HEPEnergyType getRadiationLosses(setup::Stack::particle_type const&,
@@ -66,7 +68,6 @@ namespace corsika {
GrammageType const dX_ = 10_g / square(1_cm); // profile binning
GrammageType const dX_threshold_ = 0.0001_g / square(1_cm);
ShowerAxis const& shower_axis_;
- HEPEnergyType emCut_;
HEPEnergyType energy_lost_ = HEPEnergyType::zero();
std::vector<HEPEnergyType> profile_; // longitudinal profile
};
diff --git a/examples/cascade_example.cpp b/examples/cascade_example.cpp
index eea24768d..91bc3f7cc 100644
--- a/examples/cascade_example.cpp
+++ b/examples/cascade_example.cpp
@@ -142,7 +142,7 @@ int main() {
ParticleCut cut(80_GeV, true, true);
TrackWriter trackWriter("tracks.dat");
- BetheBlochPDG eLoss{showerAxis, cut.getElectronECut()};
+ BetheBlochPDG eLoss{showerAxis};
// assemble all processes into an ordered process list
auto sequence =
diff --git a/examples/cascade_proton_example.cpp b/examples/cascade_proton_example.cpp
index e20864685..24ed9e77d 100644
--- a/examples/cascade_proton_example.cpp
+++ b/examples/cascade_proton_example.cpp
@@ -132,7 +132,7 @@ int main() {
TrackWriter trackWriter("tracks.dat");
ShowerAxis const showerAxis{injectionPos, Vector{rootCS, 0_m, 0_m, -100_km}, env};
- BetheBlochPDG eLoss{showerAxis, cut.getElectronECut()};
+ BetheBlochPDG eLoss{showerAxis};
// assemble all processes into an ordered process list
// auto sequence = sibyll << decay << hadronicElastic << cut << trackWriter;
diff --git a/examples/hybrid_MC.cpp b/examples/hybrid_MC.cpp
index 72e38e49d..83ca7bd57 100644
--- a/examples/hybrid_MC.cpp
+++ b/examples/hybrid_MC.cpp
@@ -207,7 +207,7 @@ int main(int argc, char** argv) {
decaySibyll.printDecayConfig();
ParticleCut cut{60_GeV, false, true};
- BetheBlochPDG eLoss{showerAxis, cut.getElectronECut()};
+ BetheBlochPDG eLoss{showerAxis};
CONEXhybrid conex_model(center, showerAxis, t, injectionHeight, E0,
get_PDG(Code::Proton));
diff --git a/examples/stopping_power.cpp b/examples/stopping_power.cpp
index ef8252501..d419efbb8 100644
--- a/examples/stopping_power.cpp
+++ b/examples/stopping_power.cpp
@@ -49,7 +49,7 @@ int main() {
112.8_km); // this is the CORSIKA 7 start of atmosphere/universe
ShowerAxis showerAxis{injectionPos, Vector<length_d>{rootCS, 0_m, 0_m, 1_m}, env};
- BetheBlochPDG eLoss{showerAxis, 300_MeV};
+ BetheBlochPDG eLoss{showerAxis};
setup::Stack stack;
--
GitLab