Merge branch 'radiation_losses' into 'master'

Radiative losses with constant b

See merge request AirShowerPhysics/corsika!142

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

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;

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&,