removed explicit threshold from energy loss BetheBlochPDG

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;
 

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
   };

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 =

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;

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));

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;