diff --git a/corsika/detail/modules/energy_loss/BetheBlochPDG.inl b/corsika/detail/modules/energy_loss/BetheBlochPDG.inl index 776531540cce808be3fb017c40865d6d1bddc373..1f3d0e9303cb60f40fb588ecad70034c634dd1c8 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 3e64567b974c9dfdb71b6f67c99e7877dbd6d39d..f8e8759376a3833bd3f7df438b2dbe4132c5e3c8 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 eea24768d3aadfd2eef1606684ac2f3c42e4c2da..91bc3f7cc99640a10e067da3360eee6e50911020 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 e20864685a0118f39cc3abc8539fdec7e6447322..24ed9e77d3be401229568a818511c9b2f2735563 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 72e38e49de7c9509ba2ba777b1a1cb922067a56b..83ca7bd57098a6af751f430077dffc7414fa31e1 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 ef8252501703205734ed361e4667df9fabdcbbed..d419efbb86b9cbfb3bc5a2b8db3cc8adb3a6e2bb 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;