improved EnergyLoss, const-correctness & step-length limitation

Processes/EnergyLoss/EnergyLoss.cc

@@ -24,9 +24,8 @@ using namespace std;
 
 using namespace corsika;
 using namespace corsika::units::si;
-using namespace corsika::setup;
-using Particle = Stack::ParticleType;
-using Track = Trajectory;
+using SetupParticle = corsika::setup::Stack::ParticleType;
+using SetupTrack = corsika::setup::Trajectory;
 
 namespace corsika::process::energy_loss {
 
@@ -53,7 +52,7 @@ namespace corsika::process::energy_loss {
    * For shell correction, see Sec 6 of https://www.nap.edu/read/20066/chapter/8#115
    *
    */
-  HEPEnergyType EnergyLoss::BetheBloch(Particle& p, GrammageType const dX) {
+  HEPEnergyType EnergyLoss::BetheBloch(SetupParticle const& p, GrammageType const dX) {
 
     // all these are material constants and have to come through Environment
     // right now: values for nitrogen_D
@@ -70,20 +69,20 @@ namespace corsika::process::energy_loss {
     // end of material constants
 
     // this is the Bethe-Bloch coefficiet 4pi N_A r_e^2 m_e c^2
-    auto const K = 0.307075_MeV / 1_mol * square(1_cm);
+    auto constexpr K = 0.307075_MeV / 1_mol * square(1_cm);
     HEPEnergyType const E = p.GetEnergy();
     HEPMassType const m = p.GetMass();
     double const gamma = E / m;
-    double const Z = p.GetChargeNumber();
-    double const Z2 = pow(Z, 2);
-    HEPMassType const me = particles::Electron::GetMass();
+    int const Z = p.GetChargeNumber();
+    int const Z2 = Z * Z;
+    HEPMassType constexpr me = particles::Electron::GetMass();
     auto const m2 = m * m;
-    auto const me2 = me * me;
-    double const gamma2 = pow(gamma, 2);
+    auto constexpr me2 = me * me;
+    double const gamma2 = gamma * gamma;
 
     double const beta2 = (gamma2 - 1) / gamma2; // 1-1/gamma2    (1-1/gamma)*(1+1/gamma);
                                                 // (gamma_2-1)/gamma_2 = (1-1/gamma2);
-    double const c2 = 1;                        // HEP convention here c=c2=1
+    double constexpr c2 = 1;                    // HEP convention here c=c2=1
     cout << "BetheBloch beta2=" << beta2 << " gamma2=" << gamma2 << endl;
     [[maybe_unused]] double const eta2 = beta2 / (1 - beta2);
     HEPMassType const Wmax =
@@ -148,7 +147,7 @@ namespace corsika::process::energy_loss {
            (0.5 * log(aux) - beta2 - Cadj / Z - delta / 2 + barkas + bloch) * dX;
   }
 
-  process::EProcessReturn EnergyLoss::DoContinuous(Particle& p, Track& t) {
+  process::EProcessReturn EnergyLoss::DoContinuous(SetupParticle& p, SetupTrack& t) {
     if (p.GetChargeNumber() == 0) return process::EProcessReturn::eOk;
     GrammageType const dX =
         p.GetNode()->GetModelProperties().IntegratedGrammage(t, t.GetLength());
@@ -174,8 +173,21 @@ namespace corsika::process::energy_loss {
     return status;
   }
 
-  units::si::LengthType EnergyLoss::MaxStepLength(Particle&, Track&) {
-    return units::si::meter * std::numeric_limits<double>::infinity();
+  units::si::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 Ekin = vParticle.GetEnergy() - vParticle.GetMass();
+    auto const maxLoss = 0.01 * vParticle.GetEnergy();
+    auto const maxGrammage = maxLoss / dE * dX;
+
+    return vParticle.GetNode()->GetModelProperties().ArclengthFromGrammage(vTrack,
+                                                                           maxGrammage) *
+           1.0001; // to make sure particle gets absorbed when DoContinuous() is called
   }
 
   void EnergyLoss::MomentumUpdate(corsika::setup::Stack::ParticleType& vP,

Processes/EnergyLoss/EnergyLoss.h

@@ -36,15 +36,16 @@ namespace corsika::process::energy_loss {
 
     corsika::process::EProcessReturn DoContinuous(corsika::setup::Stack::ParticleType&,
                                                   corsika::setup::Trajectory&);
-    corsika::units::si::LengthType MaxStepLength(corsika::setup::Stack::ParticleType&,
-                                                 corsika::setup::Trajectory&);
+    corsika::units::si::LengthType MaxStepLength(
+        corsika::setup::Stack::ParticleType const&,
+        corsika::setup::Trajectory const&) const;
 
     corsika::units::si::HEPEnergyType GetTotal() const { return fEnergyLossTot; }
     void PrintProfile() const;
 
   private:
-    corsika::units::si::HEPEnergyType BetheBloch(
-        corsika::setup::Stack::ParticleType& p,
+    static corsika::units::si::HEPEnergyType BetheBloch(
+        corsika::setup::Stack::ParticleType const& p,
         const corsika::units::si::GrammageType dX);
 
     int GetXbin(corsika::setup::Stack::ParticleType& p,