corsika/detail/modules/conex/CONEXhybrid.inl

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#include <corsika/framework/core/Logging.hpp>
+#include <corsika/media/CORSIKA7Atmospheres.hpp>
+#include <corsika/modules/conex/CONEXhybrid.hpp>
+#include <corsika/modules/conex/CONEXrandom.hpp>
+#include <corsika/modules/Random.hpp>
+#include <corsika/modules/conex/CONEX_f.hpp>
+#include <corsika/framework/random/RNGManager.hpp>
+#include <corsika/framework/core/PhysicalConstants.hpp>
+#include <corsika/framework/core/Logging.hpp>
+
+#include <conexConfig.h>
+
+#include <algorithm>
+#include <fstream>
+#include <numeric>
+#include <utility>
+
+namespace corsika {
+
+  template <typename TOutputE, typename TOutputN>
+  inline CONEXhybrid<TOutputE, TOutputN>::CONEXhybrid(
+      Point const& center, ShowerAxis const& showerAxis, LengthType groundDist,
+      LengthType injectionHeight, HEPEnergyType primaryEnergy, PDGCode primaryPDG,
+      TOutputE& args1, TOutputN& args2)
+      : SubWriter<TOutputE>(args1)
+      , SubWriter<TOutputN>(args2)
+      , center_{center}
+      , showerAxis_{showerAxis}
+      , groundDist_{groundDist}
+      , injectionHeight_{injectionHeight}
+      , primaryEnergy_{primaryEnergy}
+      , primaryPDG_{primaryPDG}
+      , showerCore_{showerAxis_.getStart() + showerAxis_.getDirection() * groundDist_}
+      , conexObservationCS_{std::invoke([&]() {
+        auto const& c8cs = center.getCoordinateSystem();
+        auto const translation = showerCore_ - center;
+        auto const intermediateCS =
+            make_translation(c8cs, translation.getComponents(c8cs));
+        auto const transformCS = make_rotationToZ(intermediateCS, translation);
+
+        CORSIKA_LOG_DEBUG("translation C8/CONEX obs: ", translation.getComponents());
+
+        /*
+        auto const transform = CoordinateSystem::getTransformation(
+            intermediateCS2, c8cs); // either this way or vice versa... TODO: test this!
+        */
+        return transformCS;
+      })}
+      , x_sf_{std::invoke([&]() {
+        Vector<length_d> const a{conexObservationCS_, 0._m, 0._m, 1._m};
+        auto b = a.cross(showerAxis_.getDirection());
+        auto const lengthB = b.getNorm();
+        if (lengthB < 1e-10_m) {
+          b = Vector<length_d>{conexObservationCS_, 1_m, 0_m, 0_m};
+        }
+
+        return b.normalized();
+      })}
+      , y_sf_{showerAxis_.getDirection().cross(x_sf_)} {
+
+    CORSIKA_LOG_DEBUG("x_sf (conexObservationCS): {}",
+                      x_sf_.getComponents(conexObservationCS_));
+    CORSIKA_LOG_DEBUG("x_sf (C8): {}", x_sf_.getComponents(center.getCoordinateSystem()));
+
+    CORSIKA_LOG_DEBUG("y_sf (conexObservationCS): {}",
+                      y_sf_.getComponents(conexObservationCS_));
+
+    CORSIKA_LOG_DEBUG("y_sf (C8): {}", y_sf_.getComponents(center.getCoordinateSystem()));
+
+    CORSIKA_LOG_DEBUG("showerAxisDirection (conexObservationCS): {}",
+                      showerAxis_.getDirection().getComponents(conexObservationCS_));
+    CORSIKA_LOG_DEBUG(
+        "showerAxisDirection (C8): {}",
+        showerAxis_.getDirection().getComponents(center.getCoordinateSystem()));
+
+    CORSIKA_LOG_DEBUG("showerCore (conexObservationCS): {}",
+                      showerCore_.getCoordinates(conexObservationCS_));
+    CORSIKA_LOG_DEBUG("showerCore (C8): {}",
+                      showerCore_.getCoordinates(center.getCoordinateSystem()));
+
+    auto const& components = ::corsika::standardAirComposition.getComponents();
+    auto const& fractions = ::corsika::standardAirComposition.getFractions();
+    if (::corsika::standardAirComposition.getSize() != 3) {
+      throw std::runtime_error{"CONEXhybrid only usable with standard 3-component air"};
+    }
+
+    std::transform(components.cbegin(), components.cend(), ::conex::cxair_.aira.begin(),
+                   get_nucleus_A);
+    std::transform(components.cbegin(), components.cend(), ::conex::cxair_.airz.begin(),
+                   get_nucleus_Z);
+    std::copy(fractions.cbegin(), fractions.cend(), ::conex::cxair_.airw.begin());
+
+    ::conex::cxair_.airava =
+        std::inner_product(::conex::cxair_.airw.cbegin(), ::conex::cxair_.airw.cend(),
+                           ::conex::cxair_.aira.cbegin(), 0.);
+    ::conex::cxair_.airavz =
+        std::inner_product(::conex::cxair_.airw.cbegin(), ::conex::cxair_.airw.cend(),
+                           ::conex::cxair_.airz.cbegin(), 0.);
+
+    // this is the CONEX default but actually unused there
+    ::conex::cxair_.airi = {82.0e-09, 95.0e-09, 188.e-09};
+
+    int randomSeeds[3] = {1234, 0, 0}; // SEEDS ARE NOT USED. All random numbers are
+                                       // obtained from the CORSIKA 8 stream "conex"
+    corsika::connect_random_stream("conex", ::conex::set_rng_function);
+    int heModel = eSibyll23;
+
+    int nShower = 1; // large to avoid final stats.
+    int maxDetail = 0;
+#ifdef CONEX_EXTENSIONS
+    int particleListMode = 0;
+#endif
+
+    std::string configPath = CONEX_CONFIG_PATH;
+    ::conex::initconex_(nShower, randomSeeds, heModel, maxDetail,
+#ifdef CONEX_EXTENSIONS
+                        particleListMode,
+#endif
+                        configPath.c_str(), configPath.size());
+  }
+
+  template <typename TOutputE, typename TOutputN>
+  inline void CONEXhybrid<TOutputE, TOutputN>::initCascadeEquations() {
+
+    // set phi, theta
+    Vector<length_d> ez{conexObservationCS_, {0._m, 0._m, -1_m}};
+    auto const c = showerAxis_.getDirection().dot(ez) / 1_m;
+    double theta = std::acos(c) * 180 / M_PI;
+
+    auto const showerAxisConex =
+        showerAxis_.getDirection().getComponents(conexObservationCS_);
+    double phi = std::atan2(-showerAxisConex.getY().magnitude(),
+                            showerAxisConex.getX().magnitude()) *
+                 180 / M_PI;
+
+    CORSIKA_LOG_DEBUG(
+        "theta (deg) = {}"
+        "; phi (deg) = {}",
+        theta, phi);
+
+    int ipart = static_cast<int>(primaryPDG_);
+
+    double dimpact = 0.; // valid only if shower core is fixed on the observation plane;
+                         // for skimming showers an offset is needed like in CONEX
+
+    // SEEDS ARE NOT USED. All random numbers are obtained from
+    // the CORSIKA 8 stream "conex" and "epos"!
+    std::array<int, 3> ioseed{1, 1, 1};
+
+    double eprima = primaryEnergy_ / 1_GeV;
+    double xminp = injectionHeight_ / 1_m;
+
+    ::conex::conexrun_(ipart, eprima, theta, phi, xminp, dimpact, ioseed.data());
+  }
+
+  template <typename TOutputE, typename TOutputN>
+  template <typename TStackView>
+  inline void CONEXhybrid<TOutputE, TOutputN>::doSecondaries(TStackView& vS) {
+    auto p = vS.begin();
+    while (p != vS.end()) {
+      Code const pid = p.getPID();
+      if (addParticle(pid, p.getEnergy(), p.getMass(), p.getPosition(),
+                      p.getMomentum().normalized(), p.getTime())) {
+        p.erase();
+      }
+      ++p;
+    }
+  }
+
+  template <typename TOutputE, typename TOutputN>
+  inline bool CONEXhybrid<TOutputE, TOutputN>::addParticle(
+      Code pid, HEPEnergyType energy, HEPEnergyType mass, Point const& position,
+      DirectionVector const& direction, TimeType t, double weight) {
+
+    auto const it = std::find_if(egs_em_codes_.cbegin(), egs_em_codes_.cend(),
+                                 [=](auto const& p) { return pid == p.first; });
+    if (it == egs_em_codes_.cend()) { return false; }
+
+    // EM particle
+    auto const egs_pid = it->second;
+    CORSIKA_LOG_DEBUG("position conexObs: {}",
+                      position.getCoordinates(conexObservationCS_));
+
+    auto const coords = position.getCoordinates(conexObservationCS_) / 1_m;
+    double const x = coords[0].magnitude();
+    double const y = coords[1].magnitude();
+
+    double const altitude = ((position - center_).getNorm() - conex::earthRadius) / 1_m;
+    auto const d = position - showerCore_;
+
+    // distance from core to particle projected along shower axis
+    double const slantDistance = -d.dot(showerAxis_.getDirection()) / 1_m;
+
+    // lateral coordinates in CONEX shower frame
+    auto const dShowerPlane = d - d.getParallelProjectionOnto(showerAxis_.getDirection());
+    double const lateralX = dShowerPlane.dot(x_sf_) / 1_m;
+    double const lateralY = dShowerPlane.dot(y_sf_) / 1_m;
+
+    double const slantX = showerAxis_.getProjectedX(position) * (1_cm * 1_cm / 1_g);
+
+    double const time = (t * constants::c - groundDist_) / 1_m;
+
+    // fill u,v,w momentum direction in EGS frame
+    double const u = direction.dot(y_sf_).magnitude();
+    double const v = direction.dot(x_sf_).magnitude();
+    double const w = direction.dot(showerAxis_.getDirection()).magnitude();
+
+    // generation, TO BE CHANGED WHEN WE HAVE THAT INFORMATION AVAILABLE
+    int const latchin = 1;
+
+    double const E = energy / 1_GeV;
+    double const m = mass / 1_GeV;
+
+    CORSIKA_LOG_DEBUG("CONEXhybrid: removing {} {:5e} GeV", egs_pid, energy);
+
+    CORSIKA_LOG_DEBUG("#### parameters to cegs4_() ####");
+    CORSIKA_LOG_DEBUG("egs_pid = {}", egs_pid);
+    CORSIKA_LOG_DEBUG("E = {}", E);
+    CORSIKA_LOG_DEBUG("m = {}", m);
+    CORSIKA_LOG_DEBUG("x = {}", x);
+    CORSIKA_LOG_DEBUG("y = {}", y);
+    CORSIKA_LOG_DEBUG("altitude = {}", altitude);
+    CORSIKA_LOG_DEBUG("slantDistance = {}", slantDistance);
+    CORSIKA_LOG_DEBUG("lateralX = {}", lateralX);
+    CORSIKA_LOG_DEBUG("lateralY = {}", lateralY);
+    CORSIKA_LOG_DEBUG("slantX = {}", slantX);
+    CORSIKA_LOG_DEBUG("time = {}", time);
+    CORSIKA_LOG_DEBUG("u = {}", u);
+    CORSIKA_LOG_DEBUG("v = {}", v);
+    CORSIKA_LOG_DEBUG("w = {}", w);
+
+    ::conex::cxoptl_.dptl[10 - 1] = egs_pid;
+    ::conex::cxoptl_.dptl[4 - 1] = E;
+    ::conex::cxoptl_.dptl[5 - 1] = m;
+    ::conex::cxoptl_.dptl[6 - 1] = x;
+    ::conex::cxoptl_.dptl[7 - 1] = y;
+    ::conex::cxoptl_.dptl[8 - 1] = altitude;
+    ::conex::cxoptl_.dptl[9 - 1] = time;
+    ::conex::cxoptl_.dptl[11 - 1] = weight;
+    ::conex::cxoptl_.dptl[12 - 1] = latchin;
+    ::conex::cxoptl_.dptl[13 - 1] = slantX;
+    ::conex::cxoptl_.dptl[14 - 1] = lateralX;
+    ::conex::cxoptl_.dptl[15 - 1] = lateralY;
+    ::conex::cxoptl_.dptl[16 - 1] = slantDistance;
+    ::conex::cxoptl_.dptl[2 - 1] = u;
+    ::conex::cxoptl_.dptl[1 - 1] = v;
+    ::conex::cxoptl_.dptl[3 - 1] = w;
+
+    int n = 1, i = 1;
+    ::conex::cegs4_(n, i);
+
+    return true;
+  }
+
+  template <typename TOutputE, typename TOutputN>
+  template <typename TStack>
+  inline void CONEXhybrid<TOutputE, TOutputN>::doCascadeEquations(TStack&) {
+
+    ::conex::conexcascade_();
+
+    int nX = ::conex::get_number_of_depth_bins_(); // make sure this works!
+
+    int icut = 1;
+    int icutg = 2;
+    int icute = 3;
+    int icutm = 2;
+    int icuth = 3;
+    int iSec = 0;
+
+    const int maxX = nX;
+
+    auto X = std::make_unique<float[]>(maxX);
+    auto H = std::make_unique<float[]>(maxX);
+    auto D = std::make_unique<float[]>(maxX);
+    auto N = std::make_unique<float[]>(maxX);
+    auto dEdX = std::make_unique<float[]>(maxX);
+    auto Mu = std::make_unique<float[]>(maxX);
+    auto dMu = std::make_unique<float[]>(maxX);
+    auto Photon = std::make_unique<float[]>(maxX);
+    auto Electrons = std::make_unique<float[]>(maxX);
+    auto Hadrons = std::make_unique<float[]>(maxX);
+
+    float EGround[3], fitpars[13];
+
+    ::conex::get_shower_data_(icut, iSec, nX, X[0], N[0], fitpars[0], H[0], D[0]);
+    ::conex::get_shower_edep_(icut, nX, dEdX[0], EGround[0]);
+    ::conex::get_shower_muon_(icutm, nX, Mu[0], dMu[0]);
+    ::conex::get_shower_gamma_(icutg, nX, Photon[0]);
+    ::conex::get_shower_electron_(icute, nX, Electrons[0]);
+    ::conex::get_shower_hadron_(icuth, nX, Hadrons[0]);
+
+    // make sure CONEX binning is same to C8:
+    GrammageType const dX = (X[1] - X[0]) * (1_g / square(1_cm));
+
+    for (int i = 0; i < nX; ++i) {
+      GrammageType const curX = X[i] * (1_g / square(1_cm));
+      SubWriter<TOutputE>::write(curX, curX + dX,
+                                 Code::Unknown, // this is sum of all dEdX
+                                 dEdX[i] * 1_GeV / 1_g * square(1_cm) * dX);
+      SubWriter<TOutputN>::write(curX, curX + dX, Code::Photon, Photon[i]);
+      SubWriter<TOutputN>::write(curX, curX + dX, Code::Proton, Hadrons[i]);
+      SubWriter<TOutputN>::write(curX, curX + dX, Code::Electron, Electrons[i]);
+      SubWriter<TOutputN>::write(curX, curX + dX, Code::MuMinus, Mu[i]);
+    }
+  }
+
+  template <typename TOutputE, typename TOutputN>
+  inline YAML::Node CONEXhybrid<TOutputE, TOutputN>::getConfig() const {
+
+    return YAML::Node();
+  }
+
+} // namespace corsika

corsika/detail/modules/energy_loss/BetheBlochPDG.inl

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/core/ParticleProperties.hpp>
+
+#include <corsika/framework/geometry/Line.hpp>
+#include <corsika/framework/core/Logging.hpp>
+
+#include <cmath>
+#include <fstream>
+#include <limits>
+
+namespace corsika {
+
+  template <typename TOutput>
+  template <typename... TArgs>
+  inline BetheBlochPDG<TOutput>::BetheBlochPDG(TArgs&&... args)
+      : TOutput(std::forward<TArgs>(args)...) {}
+
+  template <typename TOutput>
+  template <typename TParticle>
+  inline HEPEnergyType BetheBlochPDG<TOutput>::getBetheBloch(TParticle const& p,
+                                                             GrammageType const dX) {
+
+    // all these are material constants and have to come through Environment
+    // right now: values for nitrogen_D
+    // 7 nitrogen_gas 82.0 0.49976 D E 0.0011653 0.0 1.7378 4.1323 0.15349 3.2125 10.54
+    auto Ieff = 82.0_eV;
+    [[maybe_unused]] auto Zmat = 7;
+    auto ZoverA = 0.49976_mol / 1_g;
+    const double x0 = 1.7378;
+    const double x1 = 4.1323;
+    const double Cbar = 10.54;
+    const double delta0 = 0.0;
+    const double aa = 0.15349;
+    const double sk = 3.2125;
+    // end of material constants
+
+    // this is the Bethe-Bloch coefficiet 4pi N_A r_e^2 m_e c^2
+    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;
+    int const Z = p.getChargeNumber();
+    int const Z2 = Z * Z;
+    HEPMassType constexpr me = Electron::mass;
+    auto const m2 = m * m;
+    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 constexpr c2 = 1;                    // HEP convention here c=c2=1
+    CORSIKA_LOG_TRACE("BetheBloch beta2={}, gamma2={}", beta2, gamma2);
+    [[maybe_unused]] double const eta2 = beta2 / (1 - beta2);
+    HEPMassType const Wmax =
+        2 * me * c2 * beta2 * gamma2 / (1 + 2 * gamma * me / m + me2 / m2);
+    // approx, but <<1%    HEPMassType const Wmax = 2*me*c2*beta2*gamma2;      for HEAVY
+    // PARTICLES Wmax ~ 2me v2 for non-relativistic particles
+    CORSIKA_LOG_TRACE("BetheBloch Wmax={}", Wmax);
+
+    // Sternheimer parameterization, density corrections towards high energies
+    // NOTE/TODO: when Cbar is 0 it needs to be approximated from parameterization ->
+    // MISSING
+    CORSIKA_LOG_TRACE("BetheBloch p.getMomentum().getNorm()/m={}",
+                      p.getMomentum().getNorm() / m);
+    double const x = log10(p.getMomentum().getNorm() / m);
+    double delta = 0;
+    if (x >= x1) {
+      delta = 2 * (log(10)) * x - Cbar;
+    } else if (x < x1 && x >= x0) {
+      delta = 2 * (log(10)) * x - Cbar + aa * pow((x1 - x), sk);
+    } else if (x < x0) { // and IF conductor (otherwise, this is 0)
+      delta = delta0 * pow(100, 2 * (x - x0));
+    }
+    CORSIKA_LOG_TRACE("BetheBloch delta={}", delta);
+
+    // with further low energies correction, accurary ~1% down to beta~0.05 (1MeV for p)
+
+    // shell correction, <~100MeV
+    // need more clarity about formulas and units
+    const double Cadj = 0;
+    /*
+    // https://www.nap.edu/read/20066/chapter/8#104
+    HEPEnergyType Iadj = 12_eV * Z + 7_eV;  // Iadj<163eV
+    if (Iadj>=163_eV)
+      Iadj = 9.76_eV * Z + 58.8_eV * pow(Z, -0.19);  // Iadj>=163eV
+    double const Cadj = (0.422377/eta2 + 0.0304043/(eta2*eta2) -
+    0.00038106/(eta2*eta2*eta2)) * 1e-6 * Iadj*Iadj + (3.858019/eta2 -
+    0.1667989/(eta2*eta2) + 0.00157955/(eta2*eta2*eta2)) * 1e-9 * Iadj*Iadj*Iadj;
+    */
+
+    // Barkas correction O(Z3) higher-order Born approximation
+    // see Appl. Phys. 85 (1999) 1249
+    // double A = 1;
+    // if (p.getPID() == Code::Nucleus) A = p.getNuclearA();
+    // double const Erel = (p.getEnergy()-p.getMass()) / A / 1_keV;
+    // double const Llow = 0.01 * Erel;
+    // double const Lhigh = 1.5/pow(Erel, 0.4) + 45000./Zmat * pow(Erel, 1.6);
+    // double const barkas = Z * Llow*Lhigh/(Llow+Lhigh); // RU, I think the Z was
+    // missing...
+    double const barkas = 1; // does not work yet
+
+    // Bloch correction for O(Z4) higher-order Born approximation
+    // see Appl. Phys. 85 (1999) 1249
+    const double alpha = 1. / 137.035999173;
+    double const y2 = Z * Z * alpha * alpha / beta2;
+    double const bloch = -y2 * (1.202 - y2 * (1.042 - 0.855 * y2 + 0.343 * y2 * y2));
+
+    double const aux = 2 * me * c2 * beta2 * gamma2 * Wmax / (Ieff * Ieff);
+    return -K * Z2 * ZoverA / beta2 *
+           (0.5 * log(aux) - beta2 - Cadj / Z - delta / 2 + barkas + bloch) * dX;
+  }
+
+  // radiation losses according to PDG 2018, ch. 33 ref. [5]
+  template <typename TOutput>
+  template <typename TParticle>
+  inline HEPEnergyType BetheBlochPDG<TOutput>::getRadiationLosses(
+      TParticle 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;
+  }
+
+  template <typename TOutput>
+  template <typename TParticle>
+  inline HEPEnergyType BetheBlochPDG<TOutput>::getTotalEnergyLoss(
+      TParticle const& vP, GrammageType const vDX) {
+    return getBetheBloch(vP, vDX) + getRadiationLosses(vP, vDX);
+  }
+
+  template <typename TOutput>
+  template <typename TParticle>
+  inline ProcessReturn BetheBlochPDG<TOutput>::doContinuous(Step<TParticle>& step,
+                                                            bool const) {
+
+    // if this step was limiting the CORSIKA stepping, the particle is lost
+    /* see Issue https://gitlab.iap.kit.edu/AirShowerPhysics/corsika/-/issues/389
+    if (limitStep) {
+      fillProfile(t, p.getEnergy());
+      p.setEnergy(p.getMass());
+      return ProcessReturn::ParticleAbsorbed;
+    }
+    */
+
+    if (step.getParticlePre().getChargeNumber() == 0) return ProcessReturn::Ok;
+
+    GrammageType const dX =
+        step.getParticlePre().getNode()->getModelProperties().getIntegratedGrammage(
+            step.getStraightTrack());
+    CORSIKA_LOG_TRACE("EnergyLoss pid={}, z={}, dX={} g/cm2",
+                      step.getParticlePre().getPID(),
+                      step.getParticlePre().getChargeNumber(), dX / 1_g * square(1_cm));
+    HEPEnergyType const dE = getTotalEnergyLoss(step.getParticlePre(), dX);
+    //    if (dE > HEPEnergyType::zero())
+    //      dE = -dE;
+    CORSIKA_LOG_TRACE("EnergyLoss  dE={} MeV, Ekin={} GeV, EkinNew={} GeV", dE / 1_MeV,
+                      step.getEkinPre() / 1_GeV, (step.getEkinPre() + dE) / 1_GeV);
+    step.add_dEkin(dE);
+
+    // also send to output
+    TOutput::write(step.getPositionPre(), step.getPositionPost(),
+                   step.getParticlePre().getPID(),
+                   -step.getParticlePre().getWeight() * dE);
+    return ProcessReturn::Ok;
+  }
+
+  template <typename TOutput>
+  template <typename TParticle, typename TTrajectory>
+  inline LengthType BetheBlochPDG<TOutput>::getMaxStepLength(
+      TParticle const& vParticle, TTrajectory const& vTrack) const {
+    if (vParticle.getChargeNumber() == 0) {
+      return meter * std::numeric_limits<double>::infinity();
+    }
+
+    auto constexpr dX = 1_g / square(1_cm);
+    auto const dEdX = -getTotalEnergyLoss(vParticle, dX) / dX;
+    auto const energy = vParticle.getKineticEnergy();
+    auto const energy_lim =
+        std::max(energy * 0.9, // either 10% relative loss max., or
+                 get_kinetic_energy_propagation_threshold(
+                     vParticle.getPID()) // energy thresholds globally defined
+                                         // for individual particles
+                     * 0.99999 // need to go slightly below global e-cut to assure
+                               // removal in ParticleCut. The 1% does not matter since
+                               // at cut-time the entire energy is removed.
+        );
+    auto const maxGrammage = (energy - energy_lim) / dEdX;
+
+    return vParticle.getNode()->getModelProperties().getArclengthFromGrammage(
+        vTrack, maxGrammage);
+  }
+
+  template <typename TOutput>
+  inline YAML::Node BetheBlochPDG<TOutput>::getConfig() const {
+
+    YAML::Node node;
+    node["type"] = "BetheBlochPDG";
+    return node;
+  }
+
+} // namespace corsika

corsika/detail/modules/epos/InteractionModel.inl

+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/modules/epos/InteractionModel.hpp>
+#include <corsika/modules/epos/EposStack.hpp>
+#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <corsika/framework/utility/COMBoost.hpp>
+#include <corsika/modules/Random.hpp>
+#include <corsika/framework/utility/CorsikaData.hpp>
+
+#include <epos.hpp>
+
+#include <string>
+#include <tuple>
+#include <cmath>
+
+namespace corsika::epos {
+
+  inline InteractionModel::InteractionModel(std::set<Code> vList,
+                                            std::string const& dataPath,
+                                            bool const epos_printout_on)
+      : data_path_(dataPath)
+      , epos_listing_(epos_printout_on) {
+    // initialize Eposlhc
+    corsika::connect_random_stream(RNG_, ::epos::set_rng_function);
+    if (!isInitialized_) {
+      isInitialized_ = true;
+      if (dataPath == "") {
+        data_path_ = (std::string(corsika_data("EPOS").c_str()) + "/").c_str();
+      }
+      initialize();
+      if (vList.empty()) {
+        CORSIKA_LOGGER_DEBUG(logger_,
+                             "set all particles known to CORSIKA stable inside EPOS..");
+        setParticleListStable(get_all_particles());
+      } else {
+        CORSIKA_LOGGER_DEBUG(logger_, "set specific particles stable inside EPOS..");
+        setParticleListStable(vList);
+      }
+    }
+  }
+
+  inline void InteractionModel::setParticleListStable(std::set<Code> vPartList) const {
+    for (auto& p : vPartList) {
+      int const eid = convertToEposRaw(p);
+      if (eid != 0) {
+        // LCOV_EXCL_START
+        // this is only a safeguard against messing up the epos internals by initializing
+        // more than once.
+        unsigned int const n_particles_stable_epos =
+            ::epos::nodcy_.nrnody; // avoid waring -Wsign-compare
+        if (n_particles_stable_epos < ::epos::mxnody) {
+          CORSIKA_LOGGER_TRACE(logger_, "setting {} with EposId={} stable inside EPOS.",
+                               p, eid);
+          ::epos::nodcy_.nrnody = ::epos::nodcy_.nrnody + 1;
+          ::epos::nodcy_.nody[::epos::nodcy_.nrnody - 1] = eid;
+        } else {
+          CORSIKA_LOGGER_ERROR(logger_, "List of stable particles too long for Epos!");
+          throw std::runtime_error("Epos initialization error!");
+        }
+        // LCOV_EXCL_STOP
+      } else {
+        CORSIKA_LOG_TRACE(
+            "particle conversion Corsika-->Epos not known for {}. Using {}. Setting "
+            "unstable in Epos!",
+            p, eid);
+      }
+    }
+    CORSIKA_LOGGER_DEBUG(logger_, "set {} particles stable inside Epos",
+                         ::epos::nodcy_.nrnody);
+  }
+
+  inline bool InteractionModel::isValid(Code const projectileId, Code const targetId,
+                                        HEPEnergyType const sqrtS) const {
+    //! eposlhc only accepts nuclei with X<=A<=Y as targets, or protons aka Hydrogen or
+    //! neutrons (p,n == nucleon)
+    if (!is_nucleus(targetId) && targetId != Code::Neutron && targetId != Code::Proton) {
+      return false;
+    }
+    if (is_nucleus(targetId) && (get_nucleus_A(targetId) >= get_nucleus_A(maxNucleus_))) {
+      return false;
+    }
+    if ((minEnergyCoM_ > sqrtS) || (sqrtS > maxEnergyCoM_)) { return false; }
+    if (!epos::canInteract(projectileId)) { return false; }
+    return true;
+  }
+
+  inline void InteractionModel::initialize() const {
+
+    CORSIKA_LOGGER_DEBUG(logger_, "initializing...");
+
+    // corsika7 ini
+    int iarg = 0;
+    ::epos::aaset_(iarg);
+
+    // debug output settings
+    ::epos::prnt1_.ish = 0;
+    ::epos::prnt3_.iwseed = 0; // 1: printout seeds, 0: off
+    ::epos::files_.ifch = 6;   // output unit, 6: screen
+
+    // dummy set seeds for random number generator in epos. need to fool epos checks...
+    // we will use external generator
+    ::epos::cseed_.seedi = 1;
+    ::epos::cseed_.seedj = 1;
+    ::epos::cseed_.seedc = 1;
+
+    ::epos::enrgy_.egymin = minEnergyCoM_ / 1_GeV; // 6.;
+    ::epos::enrgy_.egymax = maxEnergyCoM_ / 1_GeV; // 2.e6;
+
+    ::epos::lhcparameters_();
+
+    ::epos::hadr6_.isigma = 0; // do not show cross section
+    ::epos::hadr6_.isetcs = 3; /*  !option to obtain pomeron parameters
+      ! 0.....determine parameters but do not use Kfit
+      ! 1.....determine parameters and use Kfit
+      ! else..get from table
+      !         should be sufficiently detailed
+      !          say iclegy1=1,iclegy2=99
+      !         table is always done, more or less detailed!!!
+      !and option to use cross section tables
+      ! 2....tabulation
+      ! 3....simulation
+                               */
+    ::epos::cjinti_.ionudi =
+        1; // !include quasi elastic events but strict calculation of xs
+    ::epos::cjinti_.iorsce = 0; // !color exchange turned on(1) or off(0)
+    ::epos::cjinti_.iorsdf = 3; //  !droplet formation turned on(>0) or off(0)
+    ::epos::cjinti_.iorshh = 0; //    !other hadron-hadron int. turned on(1) or off(0)
+
+    ::epos::othe1_.istore = 0; // do not produce epos output file
+    ::epos::nucl6_.infragm =
+        2; // 0: keep free nucleons in fragmentation,1: one fragment, 2: fragmentation
+
+    ::epos::othe2_.iframe = 11; // cms frame
+
+    // decay settings
+    // activate decays in epos for particles defined by set_stable/set_unstable
+    // ::epos::othe2_.idecay = 0; // no decays in epos
+
+    // set paths to tables in corsika data
+    ::epos::datadir BASE(data_path_);
+    strcpy(::epos::fname_.fnnx, BASE.data);
+    ::epos::nfname_.nfnnx = BASE.length;
+
+    ::epos::datadir TL(data_path_ + "epos.initl");
+    strcpy(::epos::fname_.fnii, TL.data);
+    ::epos::nfname_.nfnii = TL.length;
+
+    ::epos::datadir EV(data_path_ + "epos.iniev");
+    strcpy(::epos::fname_.fnie, EV.data);
+    ::epos::nfname_.nfnie = EV.length;
+
+    ::epos::datadir RJ(data_path_ + "epos.inirj"); // lhcparameters adds ".lhc"
+    strcpy(::epos::fname_.fnrj, RJ.data);
+    ::epos::nfname_.nfnrj = RJ.length;
+
+    ::epos::datadir CS(data_path_ + "epos.inics"); // lhcparameters adds ".lhc"
+    strcpy(::epos::fname_.fncs, CS.data);
+    ::epos::nfname_.nfncs = CS.length;
+
+    // initializes maximum energy and mass
+    initializeEventCoM(
+        maxNucleus_, get_nucleus_A(maxNucleus_), get_nucleus_Z(maxNucleus_), maxNucleus_,
+        get_nucleus_A(maxNucleus_), get_nucleus_Z(maxNucleus_), maxEnergyCoM_);
+  }
+
+  inline void InteractionModel::initializeEventCoM(Code const idBeam, int const iBeamA,
+                                                   int const iBeamZ, Code const idTarget,
+                                                   int const iTargetA, int const iTargetZ,
+                                                   HEPEnergyType const EcmNN) const {
+    CORSIKA_LOGGER_TRACE(logger_,
+                         "initialize event in CoM frame!"
+                         " Ecm={}",
+                         EcmNN);
+    ::epos::lept1_.engy = -1.;
+    ::epos::enrgy_.ecms = -1.;
+    ::epos::enrgy_.elab = -1.;
+    ::epos::enrgy_.ekin = -1.;
+    ::epos::hadr1_.pnll = -1.;
+
+    ::epos::enrgy_.ecms = EcmNN / 1_GeV; // -> c.m.s. frame
+
+    CORSIKA_LOGGER_TRACE(logger_, "inside EPOS: Ecm={}, Elab={}", ::epos::enrgy_.ecms,
+                         ::epos::enrgy_.elab);
+
+    configureParticles(idBeam, iBeamA, iBeamZ, idTarget, iTargetA, iTargetZ);
+    ::epos::ainit_();
+  }
+
+  inline void InteractionModel::configureParticles(Code const idBeam, int const iBeamA,
+                                                   int const iBeamZ, Code const idTarget,
+                                                   int const iTargetA,
+                                                   int const iTargetZ) const {
+    CORSIKA_LOGGER_TRACE(logger_,
+                         "setting "
+                         "Beam={}, "
+                         "BeamA={}, "
+                         "BeamZ={}, "
+                         "Target={}, "
+                         "TargetA={}, "
+                         "TargetZ={} ",
+                         idBeam, iBeamA, iBeamZ, idTarget, iTargetA, iTargetZ);
+
+    if (is_nucleus(idBeam)) {
+      ::epos::hadr25_.idprojin = convertToEposRaw(Code::Proton);
+      ::epos::nucl1_.laproj = iBeamZ;
+      ::epos::nucl1_.maproj = iBeamA;
+    } else {
+      ::epos::hadr25_.idprojin = convertToEposRaw(idBeam);
+      ::epos::nucl1_.laproj = -1;
+      ::epos::nucl1_.maproj = 1;
+    }
+
+    if (is_nucleus(idTarget)) {
+      ::epos::hadr25_.idtargin = convertToEposRaw(Code::Proton);
+      ::epos::nucl1_.matarg = iTargetA;
+      ::epos::nucl1_.latarg = iTargetZ;
+    } else if (idTarget == Code::Proton || idTarget == Code::Hydrogen) {
+      ::epos::hadr25_.idtargin = convertToEposRaw(Code::Proton);
+      ::epos::nucl1_.matarg = 1;
+      ::epos::nucl1_.latarg = -1;
+    } else if (idTarget == Code::Neutron) {
+      ::epos::hadr25_.idtargin = convertToEposRaw(Code::Neutron);
+      ::epos::nucl1_.matarg = 1;
+      ::epos::nucl1_.latarg = -1;
+    }
+
+    CORSIKA_LOGGER_TRACE(logger_,
+                         "inside EPOS: "
+                         "Id beam={}, "
+                         "Z beam={}, "
+                         "A beam={}, "
+                         "XS beam={}, "
+                         "Id target={}, "
+                         "Z target={}, "
+                         "A target={}, "
+                         "XS target={} ",
+                         ::epos::hadr25_.idprojin, ::epos::nucl1_.laproj,
+                         ::epos::nucl1_.maproj, ::epos::had10_.iclpro,
+                         ::epos::hadr25_.idtargin, ::epos::nucl1_.latarg,
+                         ::epos::nucl1_.matarg, ::epos::had10_.icltar);
+  }
+
+  inline InteractionModel::~InteractionModel() {
+    CORSIKA_LOGGER_DEBUG(logger_, "n={} ", count_);
+  }
+
+  inline std::tuple<CrossSectionType, CrossSectionType>
+  InteractionModel::calcCrossSectionCoM(Code const BeamId, int const BeamA,
+                                        int const BeamZ, Code const TargetId,
+                                        int const TargetA, int const TargetZ,
+                                        const HEPEnergyType EnergyCOM) const {
+    CORSIKA_LOGGER_DEBUG(logger_,
+                         "calcCrossSection: input:"
+                         " beamId={}, beamA={}, beamZ={}"
+                         " target={}, targetA={}, targetZ={}"
+                         " Ecm={:4.3f} GeV,",
+                         BeamId, BeamA, BeamZ, TargetId, TargetA, TargetZ,
+                         EnergyCOM / 1_GeV);
+
+    const int iBeam = epos::getEposXSCode(
+        BeamId); // 0 (can not interact, 1: proton-like, 2: pion-like, 3:kaon-like)
+
+    CORSIKA_LOGGER_TRACE(logger_,
+                         "projectile cross section type={} "
+                         "(0: cannot interact, 1:pion, 2:baryon, 3:kaon)",
+                         iBeam);
+    // reset beam particle // (1: pion-like, 2: proton-like, 3:kaon-like)
+    if (iBeam == 1)
+      initializeEventCoM(Code::PiPlus, BeamA, BeamZ, TargetId, TargetA, TargetZ,
+                         EnergyCOM);
+    else if (iBeam == 2)
+      initializeEventCoM(Code::Proton, BeamA, BeamZ, TargetId, TargetA, TargetZ,
+                         EnergyCOM);
+    else if (iBeam == 3)
+      initializeEventCoM(Code::KPlus, BeamA, BeamZ, TargetId, TargetA, TargetZ,
+                         EnergyCOM);
+
+    double sigProd, sigEla = 0;
+    float sigTot1, sigProd1, sigCut1 = 0;
+    if (!is_nucleus(TargetId) && !is_nucleus(BeamId)) {
+      sigProd = ::epos::hadr5_.sigine;
+      sigEla = ::epos::hadr5_.sigela;
+    } else {
+      // calculate from model, SLOW:
+      float sigQEla1 = 0; // target fragmentation/excitation
+      ::epos::crseaaepos_(sigTot1, sigProd1, sigCut1, sigQEla1);
+      sigProd = sigProd1;
+      // sigEla not properly defined here
+    }
+    CORSIKA_LOGGER_DEBUG(logger_,
+                         "calcCrossSectionCoM: output:"
+                         " sigProd={} mb,"
+                         " sigEla={} mb",
+                         sigProd, sigEla);
+
+    return std::make_tuple(sigProd * 1_mb, sigEla * 1_mb);
+  }
+
+  inline std::tuple<CrossSectionType, CrossSectionType>
+  InteractionModel::readCrossSectionTableLab(Code const BeamId, int const BeamA,
+                                             int const BeamZ, Code const TargetId,
+                                             HEPEnergyType const EnergyLab) const {
+    CORSIKA_LOGGER_DEBUG(logger_,
+                         "readCrossSectionTableLab: input: "
+                         "beamId={}, "
+                         "beamA={}, "
+                         "beamZ={} "
+                         "targetId={}, "
+                         "ELab={:12.2f} GeV,",
+                         BeamId, BeamA, BeamZ, TargetId, EnergyLab / 1_GeV);
+
+    // read cross section from epos internal tables
+    int Abeam = 0;
+    float Ekin = -1;
+
+    if (is_nucleus(BeamId)) {
+      Abeam = BeamA;
+      // kinetic energy per nucleon
+      Ekin = (EnergyLab / Abeam - constants::nucleonMass) / 1_GeV;
+    } else {
+      ::epos::hadr2_.idproj = convertToEposRaw(BeamId);
+      int const iBeam = epos::getEposXSCode(
+          BeamId); // 0 (can not interact, 1: pion-like, 2: proton-like, 3:kaon-like)
+      CORSIKA_LOGGER_TRACE(logger_,
+                           "readCrossSectionTableLab: projectile cross section type={} "
+                           "(0: cannot interact, 1:pion, 2:baryon, 3:kaon)",
+                           iBeam);
+
+      ::epos::had10_.iclpro = iBeam;
+      Abeam = 1;
+      Ekin = (EnergyLab - get_mass(BeamId)) / 1_GeV;
+    }
+
+    int Atarget = 1;
+    if (is_nucleus(TargetId)) { Atarget = get_nucleus_A(TargetId); }
+
+    int iMode = 3; // 0: air, >0 not air
+
+    CORSIKA_LOGGER_TRACE(logger_,
+                         "readCrossSectionTableLab: inside Epos "
+                         "beamId={}, beamXS={}",
+                         ::epos::hadr2_.idproj, ::epos::had10_.iclpro);
+
+    CORSIKA_LOGGER_TRACE(logger_,
+                         "readCrossSectionTableLab: calling Epos cross section with"
+                         "Ekin = {}, Abeam = {}, Atarget = {}, iMode = {}",
+                         Ekin, Abeam, Atarget, iMode);
+
+    // cross section from table, FAST
+    float sigProdEpos = ::epos::eposcrse_(Ekin, Abeam, Atarget, iMode);
+    // sig-el from analytic calculation, no fast
+    float sigElaEpos = ::epos::eposelacrse_(Ekin, Abeam, Atarget, iMode);
+
+    CORSIKA_LOGGER_TRACE(logger_,
+                         "readCrossSectionTableLab: result: sigProd = {}, sigEla = {}",
+                         sigProdEpos, sigElaEpos);
+
+    return std::make_tuple(sigProdEpos * 1_mb, sigElaEpos * 1_mb);
+  }
+
+  inline std::tuple<CrossSectionType, CrossSectionType>
+  InteractionModel::getCrossSectionInelEla(Code const projectileId, Code const targetId,
+                                           FourMomentum const& projectileP4,
+                                           FourMomentum const& targetP4) const {
+    auto const sqrtS2 = (projectileP4 + targetP4).getNormSqr();
+    auto const sqrtS = sqrt(sqrtS2);
+
+    if (!isValid(projectileId, targetId, sqrtS)) {
+      return {CrossSectionType::zero(), CrossSectionType::zero()};
+    }
+    HEPEnergyType const Elab = (sqrtS2 - static_pow<2>(get_mass(projectileId)) -
+                                static_pow<2>(get_mass(targetId))) /
+                               (2 * get_mass(targetId));
+    int beamA = 1;
+    int beamZ = 1;
+    if (is_nucleus(projectileId)) {
+      beamA = get_nucleus_A(projectileId);
+      beamZ = get_nucleus_Z(projectileId);
+    }
+
+    CORSIKA_LOGGER_DEBUG(logger_,
+                         "getCrossSectionLab: input:"
+                         " beamId={}, beamA={}, beamZ={}"
+                         " target={}"
+                         " ELab={:4.3f} GeV, sqrtS={}",
+                         projectileId, beamA, beamZ, targetId, Elab / 1_GeV,
+                         sqrtS / 1_GeV);
+    return readCrossSectionTableLab(projectileId, beamA, beamZ, targetId, Elab);
+  }
+
+  template <typename TSecondaryView>
+  inline void InteractionModel::doInteraction(TSecondaryView& view,
+                                              Code const projectileId,
+                                              Code const targetId,
+                                              FourMomentum const& projectileP4,
+                                              FourMomentum const& targetP4) {
+
+    count_ = count_ + 1;
+
+    // define nucleon-nucleon center-of-mass frame
+    auto const projectileP4NN =
+        projectileP4 / (is_nucleus(projectileId) ? get_nucleus_A(projectileId) : 1);
+    auto const targetP4NN =
+        targetP4 / (is_nucleus(targetId) ? get_nucleus_A(targetId) : 1);
+    auto const SNN = (projectileP4NN + targetP4NN).getNormSqr();
+    HEPEnergyType const sqrtSNN = sqrt(SNN);
+    if (!isValid(projectileId, targetId, sqrtSNN)) {
+      throw std::runtime_error("invalid projectile/target/energy combination.");
+    }
+    HEPEnergyType const Elab = (SNN - static_pow<2>(get_mass(projectileId)) -
+                                static_pow<2>(get_mass(targetId))) /
+                               (2 * get_mass(targetId));
+
+    // system of initial-state
+    COMBoost const boost(projectileP4NN, targetP4NN);
+
+    auto const& originalCS = boost.getOriginalCS();
+    auto const& csPrime =
+        boost.getRotatedCS(); // z is along the CM motion (projectile, in Cascade)
+
+    CORSIKA_LOGGER_DEBUG(logger_,
+                         "doInteraction: interaction, projectile id={}, E={}, p3={} ",
+                         projectileId, projectileP4.getTimeLikeComponent(),
+                         projectileP4.getSpaceLikeComponents());
+    CORSIKA_LOGGER_DEBUG(
+        logger_, "doInteraction: projectile per-nucleon ENN={}, p3NN={} ",
+        projectileP4NN.getTimeLikeComponent(), projectileP4NN.getSpaceLikeComponents());
+    CORSIKA_LOGGER_DEBUG(
+        logger_, "doInteraction: interaction, target id={}, E={}, p3={} ", targetId,
+        targetP4.getTimeLikeComponent(), targetP4.getSpaceLikeComponents());
+    CORSIKA_LOGGER_DEBUG(logger_, "doInteraction: target per-nucleon ENN={}, p3NN={} ",
+                         targetP4NN.getTimeLikeComponent(),
+                         targetP4NN.getSpaceLikeComponents());
+    CORSIKA_LOGGER_DEBUG(logger_, "doInteraction: Elab={}, sqrtSNN={} ", Elab, sqrtSNN);
+
+    int beamA = 1;
+    int beamZ = 1;
+    if (is_nucleus(projectileId)) {
+      beamA = get_nucleus_A(projectileId);
+      beamZ = get_nucleus_Z(projectileId);
+      CORSIKA_LOGGER_DEBUG(logger_, "projectile: A={}, Z={} ", beamA, beamZ);
+    }
+
+    // // from corsika7 interface
+    // // NEXLNK-part
+    int targetA = 1;
+    int targetZ = 1;
+    if (is_nucleus(targetId)) {
+      targetA = get_nucleus_A(targetId);
+      targetZ = get_nucleus_Z(targetId);
+      CORSIKA_LOGGER_DEBUG(logger_, "target: A={}, Z={} ", targetA, targetZ);
+    }
+    initializeEventCoM(projectileId, beamA, beamZ, targetId, targetA, targetZ, sqrtSNN);
+
+    // create event
+    int iarg = 1;
+    ::epos::aepos_(iarg);
+    ::epos::afinal_();
+
+    if (epos_listing_) { // LCOV_EXCL_START
+      char nam[9] = "EPOSLHC&";
+      ::epos::alistf_(nam, 9);
+    } // LCOV_EXCL_STOP
+
+    // NSTORE-part
+
+    MomentumVector P_final(originalCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
+    HEPEnergyType E_final = 0_GeV;
+
+    // secondaries
+    EposStack es;
+    CORSIKA_LOGGER_DEBUG(logger_, "number of entries on Epos stack: {}", es.getSize());
+    for (auto& psec : es) {
+      if (!psec.isFinal()) continue;
+
+      auto momentum = psec.getMomentum(csPrime);
+      // transform particle output to frame defined by input 4-momenta
+      auto const P4output = boost.fromCoM(FourVector{psec.getEnergy(), momentum});
+      auto p3output = P4output.getSpaceLikeComponents();
+      p3output.rebase(originalCS); // transform back into standard lab frame
+
+      EposCode const eposId = psec.getPID();
+      Code const pid = epos::convertFromEpos(eposId);
+      HEPEnergyType const mass = get_mass(pid);
+      HEPEnergyType const Ekin = sqrt(p3output.getSquaredNorm() + mass * mass) - mass;
+      CORSIKA_LOGGER_TRACE(logger_,
+                           " id= {}"
+                           " p= {}",
+                           pid, p3output.getComponents() / 1_GeV);
+
+      auto pnew = view.addSecondary(std::make_tuple(pid, Ekin, p3output.normalized()));
+      P_final += pnew.getMomentum();
+      E_final += pnew.getEnergy();
+    }
+    CORSIKA_LOGGER_DEBUG(
+        logger_,
+        "conservation (all GeV): Ecm_final= n/a" /* << Ecm_final / 1_GeV*/
+        ", E_final={} GeV"
+        ", P_final={} GeV"
+        ", no. of particles={}",
+        E_final / 1_GeV, (P_final / 1_GeV).getComponents(), view.getSize());
+  }
+} // namespace corsika::epos

corsika/detail/modules/epos/ParticleConversion.inl

+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <epos.hpp>
+
+namespace corsika::epos {
+
+  inline HEPMassType getEposMass(Code const pCode) {
+    if (is_nucleus(pCode)) throw std::runtime_error("Not suited for Nuclei.");
+    auto sCode = convertToEposRaw(pCode);
+    if (sCode == 0)
+      throw std::runtime_error("getEposMass: unknown particle!");
+    else {
+      float mass2 = 0;
+      ::epos::idmass_(sCode, mass2);
+      return sqrt(mass2) * 1_GeV;
+    }
+  }
+
+  inline PDGCode getEposPDGId(Code const p) {
+    if (!is_nucleus(p)) {
+      int eid = corsika::epos::convertToEposRaw(p);
+      char nxs[4] = "nxs";
+      char pdg[4] = "pdg";
+      return static_cast<PDGCode>(::epos::idtrafo_(nxs, pdg, eid));
+    } else {
+      throw std::runtime_error("Epos id conversion not implemented for nuclei!");
+    }
+  }
+} // namespace corsika::epos

corsika/detail/modules/fluka/InteractionModel.inl

+/*
+ * (c) Copyright 2022 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <algorithm>
+#include <cstdlib>
+#include <vector>
+#include <iterator>
+#include <set>
+#include <utility>
+#include <string_view>
+
+#include <boost/iterator/zip_iterator.hpp>
+#include <Eigen/Dense>
+
+#include <corsika/media/Environment.hpp>
+#include <corsika/media/NuclearComposition.hpp>
+#include <corsika/framework/geometry/FourVector.hpp>
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <corsika/framework/core/EnergyMomentumOperations.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/modules/Random.hpp>
+#include <corsika/modules/fluka/ParticleConversion.hpp>
+
+#include <FLUKA.hpp>
+
+namespace corsika::fluka {
+  inline InteractionModel::InteractionModel(std::set<Code> const& nuccomp)
+      : materials_{genFlukaMaterials(nuccomp)}
+      , cumsgx_{std::make_unique<double[]>(materials_.size() * 3)} {
+    CORSIKA_LOGGER_INFO(logger_, "FLUKA version {}", ::fluka::get_version());
+    for (auto const& [code, matno] : materials_) {
+      CORSIKA_LOGGER_DEBUG(logger_, "FLUKA material initialization: {} -> {}",
+                           get_name(code, full_name{}), matno);
+    }
+
+    if (int const ndmhep = ::fluka::ndmhep_(); ::fluka::nmxhep != ndmhep) {
+      CORSIKA_LOGGER_CRITICAL(logger_, "HEPEVT dimension mismatch. FLUKA reports %d",
+                              ndmhep);
+      throw std::runtime_error{"FLUKA HEPEVT dimension mismatch"};
+    }
+
+    corsika::connect_random_stream(RNG_, ::fluka::set_rng_function);
+  }
+
+  inline bool InteractionModel::isValid(Code projectileID, int material,
+                                        HEPEnergyType /*sqrtS*/) const {
+    if (!fluka::canInteract(projectileID)) {
+      // invalid projectile
+      return false;
+    }
+
+    if (material < 0) {
+      // invalid/uninitialized target
+      return false;
+    }
+
+    // TODO: check validity range of sqrtS
+    return true;
+  }
+
+  inline bool InteractionModel::isValid(Code projectileID, Code targetID,
+                                        HEPEnergyType sqrtS) const {
+    return isValid(projectileID, getMaterialIndex(targetID), sqrtS);
+  }
+
+  inline int InteractionModel::getMaterialIndex(Code targetID) const {
+    auto compare = [=](std::pair<Code, int> const& p) { return p.first == targetID; };
+    if (auto it = std::find_if(materials_.cbegin(), materials_.cend(), compare);
+        it == materials_.cend()) {
+      return -1;
+    } else {
+      return it->second;
+    }
+  }
+
+  inline CrossSectionType InteractionModel::getCrossSection(
+      Code const projectileId, Code const targetId, FourMomentum const& projectileP4,
+      FourMomentum const& targetP4) const {
+    auto const flukaMaterial = getMaterialIndex(targetId);
+
+    HEPEnergyType const sqrtS = (projectileP4 + targetP4).getNorm();
+    if (!isValid(projectileId, flukaMaterial, sqrtS)) { return CrossSectionType::zero(); }
+
+    COMBoost const targetRestBoost{targetP4.getSpaceLikeComponents(), get_mass(targetId)};
+    FourMomentum const projectileLab4mom = targetRestBoost.toCoM(projectileP4);
+    HEPEnergyType const Elab = projectileLab4mom.getTimeLikeComponent();
+    auto constexpr invGeV = 1 / 1_GeV;
+    double const labMomentum =
+        projectileLab4mom.getSpaceLikeComponents().getNorm() * invGeV;
+
+    auto const plab = projectileLab4mom.getSpaceLikeComponents();
+
+    CORSIKA_LOGGER_DEBUG(logger_, fmt::format("Elab = {} GeV", Elab * invGeV));
+    auto const flukaCodeProj =
+        static_cast<FLUKACodeIntType>(convertToFluka(projectileId));
+
+    double const dummyEkin = 0;
+    CrossSectionType const xs = ::fluka::sgmxyz_(&flukaCodeProj, &flukaMaterial,
+                                                 &dummyEkin, &labMomentum, &iflxyz_) *
+                                1_mb;
+    return xs;
+  }
+
+  template <typename TSecondaryView>
+  inline void InteractionModel::doInteraction(TSecondaryView& view,
+                                              Code const projectileId,
+                                              Code const targetId,
+                                              FourMomentum const& projectileP4,
+                                              FourMomentum const& targetP4) {
+
+    auto const flukaCodeProj =
+        static_cast<FLUKACodeIntType>(convertToFluka(projectileId));
+    auto const flukaMaterial = getMaterialIndex(targetId);
+
+    HEPEnergyType const sqrtS = (projectileP4 + targetP4).getNorm();
+    if (!isValid(projectileId, flukaMaterial, sqrtS)) {
+      std::string const errmsg = fmt::format(
+          "Event generation with invalid configuration requested: proj: {}, target: {}",
+          get_name(projectileId, full_name{}), get_name(targetId, full_name{}));
+      CORSIKA_LOGGER_CRITICAL(logger_, errmsg);
+      throw std::runtime_error{errmsg.c_str()};
+    }
+
+    COMBoost const targetRestBoost{targetP4.getSpaceLikeComponents(), get_mass(targetId)};
+    FourMomentum const projectileLab4mom = targetRestBoost.toCoM(projectileP4);
+    auto constexpr invGeV = 1 / 1_GeV;
+
+    auto const plab = projectileLab4mom.getSpaceLikeComponents();
+    auto const& cs = plab.getCoordinateSystem();
+    auto const labMomentum = plab.getNorm();
+    double const labMomentumGeV = labMomentum * invGeV;
+
+    auto const direction = (plab / labMomentum).getComponents().getEigenVector();
+
+    double const dummyEkin = 0;
+    ::fluka::evtxyz_(&flukaCodeProj, &flukaMaterial, &dummyEkin, &labMomentumGeV,
+                     &direction[0], &direction[1], &direction[2], &iflxyz_, cumsgx_.get(),
+                     cumsgx_.get() + materials_.size(),
+                     cumsgx_.get() + materials_.size() * 2);
+
+    for (int i = 0; i < ::fluka::hepevt_.nhep; ++i) {
+      int const status = ::fluka::hepevt_.isthep[i];
+      if (status != 1) // skip non-final-state particles
+        continue;
+
+      auto const pdg = static_cast<corsika::PDGCode>(::fluka::hepevt_.idhep[i]);
+      auto const c8code = corsika::convert_from_PDG(pdg);
+      auto const mom = QuantityVector<hepenergy_d>{
+          Eigen::Map<Eigen::Vector3d>(&::fluka::hepevt_.phep[i][0]) *
+          (1_GeV).magnitude()};
+      auto const pPrime = mom.getNorm();
+      auto const c8mass = corsika::get_mass(c8code);
+
+      auto const fourMomCollisionFrame =
+          FourVector{calculate_total_energy(pPrime, c8mass), MomentumVector{cs, mom}};
+      auto const fourMomOrigFrame = targetRestBoost.fromCoM(fourMomCollisionFrame);
+      auto const momOrigFrame = fourMomOrigFrame.getSpaceLikeComponents();
+      auto const p = momOrigFrame.getNorm();
+
+      view.addSecondary(std::tuple{c8code, corsika::calculate_kinetic_energy(p, c8mass),
+                                   momOrigFrame / p});
+    }
+  }
+
+  inline std::vector<std::pair<Code, int>> InteractionModel::genFlukaMaterials(
+      std::set<Code> const& allElementsInUniverse) {
+
+    /*
+     * We define one material per element/isotope we have in C8. Cross-section averaging
+     * and target isotope selection happen in C8.
+     */
+
+    int const nElements = allElementsInUniverse.size();
+    auto nelmfl = std::make_unique<int[]>(nElements);
+    std::vector<int> izelfl;
+    izelfl.reserve(nElements);
+    auto wfelml = std::make_unique<double[]>(nElements);
+    auto const mxelfl = nElements;
+    double const pptmax = 1e11; // GeV
+    double const ef2dp3 = 0;    // GeV, 0 means default is used
+    double const df2dp3 = -1;   // default
+    bool const lprint = true;
+    auto mtflka = std::make_unique<int[]>(mxelfl);
+    // magic number that FLUKA uses to see if it's the right version
+    char const crvrck[] = "76466879";
+    int const size = 8;
+
+    std::fill(&nelmfl[0], &nelmfl[nElements], 1);
+    std::fill(&wfelml[0], &wfelml[nElements], 1.);
+    std::transform(allElementsInUniverse.cbegin(), allElementsInUniverse.cend(),
+                   std::back_inserter(izelfl), get_nucleus_Z);
+
+    // check if FLUPRO is set
+    if (std::getenv("FLUPRO") == nullptr) {
+      throw std::runtime_error{"FLUPRO not set; required to initialize FLUKA"};
+    }
+
+    // call FLUKA
+    ::fluka::stpxyz_(&nElements, nelmfl.get(), izelfl.data(), wfelml.get(), &mxelfl,
+                     &pptmax, &ef2dp3, &df2dp3, &iflxyz_, &lprint, mtflka.get(), crvrck,
+                     &size);
+
+    // now create & fill vector of (C8 Code, FLUKA mat. no.) pairs
+    std::vector<std::pair<Code, int>> mapping;
+    mapping.reserve(nElements);
+
+    auto it = boost::make_zip_iterator(
+        boost::make_tuple(allElementsInUniverse.begin(), &mtflka[0]));
+    auto const end = boost::make_zip_iterator(
+        boost::make_tuple(allElementsInUniverse.end(), &mtflka[nElements]));
+    for (; it != end; ++it) {
+      boost::tuple<Code const&, int&> tup = *it;
+      mapping.emplace_back(tup.get<0>(), tup.get<1>());
+    }
+
+    return mapping;
+  }
+} // namespace corsika::fluka

corsika/detail/modules/proposal/ContinuousProcess.inl

+
+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#include <PROPOSAL/PROPOSAL.h>
+
+#include <corsika/media/IMediumModel.hpp>
+#include <corsika/modules/proposal/ContinuousProcess.hpp>
+#include <corsika/modules/proposal/InteractionModel.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/utility/COMBoost.hpp>
+#include <corsika/framework/core/Logging.hpp>
+
+namespace corsika::proposal {
+
+  template <typename TOutput>
+  inline void ContinuousProcess<TOutput>::buildCalculator(Code code,
+                                                          size_t const& component_hash) {
+    // search crosssection builder for given particle
+    auto p_cross = cross.find(code);
+    if (p_cross == cross.end())
+      throw std::runtime_error("PROPOSAL could not find corresponding builder");
+    if (code == Code::Photon) return; // no continuous builders needed for photons
+
+    // interpolate the crosssection for given media and energy cut. These may
+    // take some minutes if you have to build the tables and cannot read the tables
+    // from disk
+    auto c = p_cross->second(media.at(component_hash), proposal_energycutsettings[code]);
+
+    // choose multiple scattering model
+    static constexpr auto ms_type = PROPOSAL::MultipleScatteringType::MoliereInterpol;
+
+    // Build displacement integral and scattering object and interpolate them too and
+    // saved in the calc map by a key build out of a hash of composed of the component and
+    // particle code.
+    auto calculator = Calculator{PROPOSAL::make_displacement(c, true),
+                                 PROPOSAL::make_multiple_scattering(
+                                     ms_type, particle[code], media.at(component_hash))};
+    calc[std::make_pair(component_hash, code)] = std::move(calculator);
+  }
+
+  template <typename TOutput>
+  template <typename TEnvironment, typename... TOutputArgs>
+  inline ContinuousProcess<TOutput>::ContinuousProcess(TEnvironment const& _env,
+                                                       TOutputArgs&&... args)
+      : ProposalProcessBase(_env)
+      , TOutput(args...) {
+    //! Initialize PROPOSAL tables for all media and all particles
+    for (auto const& medium : media) {
+      for (auto const particle_code : tracked) {
+        buildCalculator(particle_code, medium.first);
+      }
+    }
+  }
+
+  template <typename TOutput>
+  template <typename TParticle>
+  inline void ContinuousProcess<TOutput>::scatter(Step<TParticle>& step,
+                                                  HEPEnergyType const& loss,
+                                                  GrammageType const& grammage) {
+
+    // get or build corresponding calculators
+    auto c = getCalculator(step.getParticlePre(), calc);
+
+    auto initial_particle_dir = step.getDirectionPre();
+
+    auto E_i_total = step.getEkinPre() + step.getParticlePre().getMass();
+    auto E_f_total = E_i_total - loss;
+
+    // sample scattering_angle, which is the combination sqrt(theta_x^2 + theta_y^2),
+    // where theta_x and theta_y itself are independent angles drawn from the multiple
+    // scattering distribution
+    std::uniform_real_distribution<double> distr(0., 1.);
+    auto scattering_angle = (c->second).scatter->CalculateScatteringAngle2D(
+        grammage / 1_g * square(1_cm), E_i_total / 1_MeV, E_f_total / 1_MeV, distr(RNG_),
+        distr(RNG_));
+
+    auto const& root = initial_particle_dir.getCoordinateSystem();
+
+    // construct vector that is normal to initial direction.
+    DirectionVector normal_vec{root, {0, 0, 0}};
+    if (initial_particle_dir.getX(root) > 0.1) {
+      normal_vec = {
+          root, {-initial_particle_dir.getY(root), initial_particle_dir.getX(root), 0}};
+    } else {
+      // if x is small, use y and z to construct normal vector
+      normal_vec = {
+          root, {0, -initial_particle_dir.getZ(root), initial_particle_dir.getY(root)}};
+    }
+
+    // rotation of zenith by moliere_angle
+    CoordinateSystemPtr rotated1 =
+        make_rotation(root, normal_vec.getComponents(), scattering_angle);
+
+    // rotation of azimuth by random angle between 0 and 2*PI
+    std::uniform_real_distribution<double> distr_azimuth(0., 2 * M_PI);
+    CoordinateSystemPtr rotated2 = make_rotation(
+        rotated1, initial_particle_dir.getComponents(), distr_azimuth(RNG_));
+
+    DirectionVector scattered_particle_dir{root,
+                                           initial_particle_dir.getComponents(rotated2)};
+
+    // update particle direction after continuous loss caused by multiple
+    // scattering
+    DirectionVector diff_dir_ = scattered_particle_dir - initial_particle_dir;
+    step.add_dU(diff_dir_);
+  }
+
+  template <typename TOutput>
+  template <typename TParticle>
+  inline ProcessReturn ContinuousProcess<TOutput>::doContinuous(Step<TParticle>& step,
+                                                                bool const) {
+    if (!canInteract(step.getParticlePre().getPID())) return ProcessReturn::Ok;
+    if (step.getDisplacement().getSquaredNorm() == static_pow<2>(0_m))
+      return ProcessReturn::Ok;
+    if (step.getParticlePre().getPID() == Code::Photon)
+      return ProcessReturn::Ok; // no continuous energy losses, no scattering for photons
+    // calculate passed grammage
+    auto dX = step.getParticlePre().getNode()->getModelProperties().getIntegratedGrammage(
+        step.getStraightTrack());
+
+    // get or build corresponding track integral calculator and solve the
+    // integral
+    auto c = getCalculator(step.getParticlePre(), calc);
+    auto E_i_total = (step.getEkinPre() + step.getParticlePre().getMass());
+    auto E_f_total = (c->second).disp->UpperLimitTrackIntegral(
+                         E_i_total * (1 / 1_MeV), dX * ((1 / 1_g) * 1_cm * 1_cm)) *
+                     1_MeV;
+    auto dE = E_i_total - E_f_total;
+
+    // if the particle has a charge take multiple scattering into account
+    if (step.getParticlePre().getChargeNumber() != 0) scatter(step, dE, dX);
+    step.add_dEkin(-dE); // on the stack, this is just kinetic energy, E-m
+
+    // also send to output
+    TOutput::write(step.getPositionPre(), step.getPositionPost(),
+                   step.getParticlePre().getPID(),
+                   step.getParticlePre().getWeight() * dE);
+
+    return ProcessReturn::Ok;
+  }
+
+  template <typename TOutput>
+  template <typename TParticle, typename TTrajectory>
+  inline LengthType ContinuousProcess<TOutput>::getMaxStepLength(
+      TParticle const& vP, TTrajectory const& track) {
+    auto const code = vP.getPID();
+    if (!canInteract(code)) return meter * std::numeric_limits<double>::infinity();
+    if (code == Code::Photon)
+      return meter *
+             std::numeric_limits<double>::infinity(); // no step limitation for photons
+
+    // Limit the step size of a conitnuous loss. The maximal continuous loss seems to be
+    // a hyper parameter which must be adjusted.
+    //
+    auto const energy = vP.getEnergy();
+    auto const energy_lim = std::max(
+        energy * 0.9, // either 10% relative loss max., or
+        (get_kinetic_energy_propagation_threshold(code) +
+         get_mass(code)) // energy thresholds globally defined for individual particles
+            * 0.9999     // need to go slightly below global e-cut to assure removal
+                         // in ParticleCut. This does not matter since at cut-time
+                         // the entire energy is removed.
+    );
+
+    // solving the track integral for giving energy lim
+    auto c = getCalculator(vP, calc);
+    auto grammage =
+        (c->second).disp->SolveTrackIntegral(energy / 1_MeV, energy_lim / 1_MeV) * 1_g /
+        square(1_cm);
+
+    // return it in distance aequivalent
+    auto dist =
+        vP.getNode()->getModelProperties().getArclengthFromGrammage(track, grammage);
+    CORSIKA_LOG_TRACE("PROPOSAL::getMaxStepLength X={} g/cm2, l={} m ",
+                      grammage / 1_g * square(1_cm), dist / 1_m);
+    if (dist < 0_m) {
+      CORSIKA_LOG_WARN(
+          "PROPOSAL::getMaxStepLength calculated a negative step length of l={} m. "
+          "Return 0_m instead.",
+          dist / 1_m);
+      return 0_m;
+    }
+    return dist;
+  }
+
+  template <typename TOutput>
+  inline YAML::Node ContinuousProcess<TOutput>::getConfig() const {
+    return YAML::Node();
+  }
+
+} // namespace corsika::proposal

corsika/detail/modules/proposal/HadronicPhotonModel.inl

+/*
+ * (c) Copyright 2022 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/core/EnergyMomentumOperations.hpp>
+
+#include <tuple>
+#include <random>
+
+namespace corsika::proposal {
+
+  template <typename THadronicLEModel, typename THadronicHEModel>
+  inline HadronicPhotonModel<THadronicLEModel, THadronicHEModel>::HadronicPhotonModel(
+      THadronicLEModel& _hadintLE, THadronicHEModel& _hadintHE,
+      HEPEnergyType const& _heenthresholdNN)
+      : leHadronicInteraction_(_hadintLE)
+      , heHadronicInteraction_(_hadintHE)
+      , heHadronicModelThresholdLabNN_(_heenthresholdNN) {
+    // check validity of threshold assuming photon-nucleon
+    // sqrtS per target nucleon
+    HEPEnergyType const sqrtS =
+        calculate_com_energy(_heenthresholdNN, Rho0::mass, Proton::mass);
+    if (!heHadronicInteraction_.isValid(Code::Rho0, Code::Proton, sqrtS)) {
+      CORSIKA_LOGGER_CRITICAL(
+          logger_,
+          "Invalid energy threshold for hadron interaction model. threshold_lab= {} GeV, "
+          "threshold_com={} GeV",
+          _heenthresholdNN / 1_GeV, sqrtS / 1_GeV);
+      throw std::runtime_error("Configuration error!");
+    }
+    CORSIKA_LOGGER_DEBUG(
+        logger_, "Threshold for HE hadronic interactions in proposal set to Elab={} GeV",
+        _heenthresholdNN / 1_GeV);
+  }
+
+  template <typename THadronicLEModel, typename THadronicHEModel>
+  template <typename TStackView>
+  inline ProcessReturn
+  HadronicPhotonModel<THadronicLEModel, THadronicHEModel>::doHadronicPhotonInteraction(
+      TStackView& view, CoordinateSystemPtr const& labCS, FourMomentum const& photonP4,
+      Code const& targetId) {
+
+    //  temporarily add to stack, will be removed after interaction in DoInteraction
+    typename TStackView::inner_stack_value_type photonStack;
+    Point const pDummy(labCS, {0_m, 0_m, 0_m});
+    TimeType const tDummy = 0_ns;
+    // target at rest
+    FourMomentum const targetP4(get_mass(targetId),
+                                MomentumVector(labCS, {0_GeV, 0_GeV, 0_GeV}));
+    auto hadronicPhoton = photonStack.addParticle(
+        std::make_tuple(Code::Photon, photonP4.getTimeLikeComponent(),
+                        photonP4.getSpaceLikeComponents().normalized(), pDummy, tDummy));
+    hadronicPhoton.setNode(view.getProjectile().getNode());
+    // create inelastic interaction of the hadronic photon
+    // create new StackView for the photon
+    TStackView photon_secondaries(hadronicPhoton);
+
+    // call inner hadronic event generator
+    CORSIKA_LOGGER_TRACE(logger_, "{} + {} interaction. Ekinlab = {} GeV", Code::Photon,
+                         targetId, photonP4.getTimeLikeComponent() / 1_GeV);
+    // check if had. model can handle configuration
+
+    if (photonP4.getTimeLikeComponent() > heHadronicModelThresholdLabNN_) {
+      CORSIKA_LOGGER_TRACE(logger_, "HE photo-hadronic interaction!");
+      auto const sqrtSNN = (photonP4 + targetP4 / get_nucleus_A(targetId)).getNorm();
+      CORSIKA_LOGGER_DEBUG(logger_, "sqrtS={} GeV", sqrtSNN / 1_GeV);
+      // when Sibyll is used for hadronic interactions Argon cannot be used as target
+      // nucleus. Since PROPOSAL has a non-zero cross section for Argon
+      // targets we have to check here if the model can handle Argon (see Issue #498)
+      if (!heHadronicInteraction_.isValid(Code::Rho0, targetId, sqrtSNN)) {
+        CORSIKA_LOGGER_WARN(
+            logger_,
+            "HE interaction model cannot handle configuration in photo-hadronic "
+            "interaction! projectile={}, target={} (A={}, Z={}), sqrt(S) per "
+            "nuc.={:8.2f} "
+            "GeV. Skipping secondary production!",
+            Code::Rho0, targetId, get_nucleus_A(targetId), get_nucleus_Z(targetId),
+            sqrtSNN / 1_GeV);
+        return ProcessReturn::Ok;
+      }
+      heHadronicInteraction_.doInteraction(photon_secondaries, Code::Rho0, targetId,
+                                           photonP4, targetP4);
+    } else {
+      CORSIKA_LOGGER_TRACE(logger_,
+                           "LE photo-hadronic interaction! implemented via SOPHIA "
+                           "assuming a single nucleon as target");
+      // sample nucleon from nucleus A,Z
+      double const fProtons = get_nucleus_Z(targetId) / double(get_nucleus_A(targetId));
+      double const fNeutrons = 1. - fProtons;
+      std::discrete_distribution<int> nucleonChannelDist{fProtons, fNeutrons};
+      corsika::default_prng_type& rng =
+          corsika::RNGManager<>::getInstance().getRandomStream("proposal");
+      Code const nucleonId = (nucleonChannelDist(rng) ? Code::Neutron : Code::Proton);
+      // target passed to SOPHIA needs to be exactly on-shell!
+      FourMomentum const nucleonP4(get_mass(nucleonId),
+                                   MomentumVector(labCS, {0_GeV, 0_GeV, 0_GeV}));
+      CORSIKA_LOGGER_DEBUG(logger_,
+                           "selected {} as target nucleon (f_proton, f_neutron)={},{}",
+                           nucleonId, fProtons, fNeutrons);
+      auto const sqrtSNN = (photonP4 + nucleonP4).getNorm();
+      CORSIKA_LOGGER_DEBUG(logger_, "sqrtS={} GeV", sqrtSNN / 1_GeV);
+
+      if (!leHadronicInteraction_.isValid(Code::Photon, nucleonId, sqrtSNN)) {
+        CORSIKA_LOGGER_WARN(
+            logger_,
+            "LE interaction model cannot handle configuration in photo-hadronic "
+            "interaction! projectile={}, target={} (A={}, Z={}), sqrt(S) per "
+            "nuc.={:8.2f} "
+            "GeV. Skipping secondary production!",
+            Code::Photon, targetId, get_nucleus_A(targetId), get_nucleus_Z(targetId),
+            sqrtSNN / 1_GeV);
+        return ProcessReturn::Ok;
+      }
+      leHadronicInteraction_.doInteraction(photon_secondaries, Code::Photon, nucleonId,
+                                           photonP4, nucleonP4);
+    }
+    for (const auto& pSec : photon_secondaries) {
+      auto const p3lab = pSec.getMomentum();
+      Code const pid = pSec.getPID();
+      HEPEnergyType const secEkin =
+          calculate_kinetic_energy(p3lab.getNorm(), get_mass(pid));
+      view.addSecondary(std::make_tuple(pid, secEkin, p3lab.normalized()));
+    }
+    return ProcessReturn::Ok;
+  }
+} // namespace corsika::proposal

corsika/detail/modules/proposal/InteractionModel.inl

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#include <corsika/media/IMediumModel.hpp>
+#include <corsika/media/NuclearComposition.hpp>
+#include <corsika/framework/utility/COMBoost.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+
+#include <limits>
+#include <memory>
+#include <random>
+#include <tuple>
+#include <PROPOSAL/particle/Particle.h>
+#include <corsika/modules/proposal/InteractionModel.hpp>
+
+namespace corsika::proposal {
+
+  template <typename THadronicLEModel, typename THadronicHEModel>
+  template <typename TEnvironment>
+  inline InteractionModel<THadronicLEModel, THadronicHEModel>::InteractionModel(
+      TEnvironment const& _env, THadronicLEModel& _hadintLE, THadronicHEModel& _hadintHE,
+      HEPEnergyType const& _enthreshold)
+      : ProposalProcessBase(_env)
+      , HadronicPhotonModel<THadronicLEModel, THadronicHEModel>(_hadintLE, _hadintHE,
+                                                                _enthreshold) {
+    //! Initialize PROPOSAL tables for all media and all particles
+    for (auto const& medium : media) {
+      for (auto const particle_code : tracked) {
+        buildCalculator(particle_code, medium.first);
+      }
+    }
+  }
+
+  template <typename THadronicLEModel, typename THadronicHEModel>
+  inline void InteractionModel<THadronicLEModel, THadronicHEModel>::buildCalculator(
+      Code code, size_t const& component_hash) {
+    // search crosssection builder for given particle
+    auto p_cross = cross.find(code);
+    if (p_cross == cross.end())
+      throw std::runtime_error("PROPOSAL could not find corresponding builder");
+
+    // interpolate the crosssection for given media and energy cut. These may
+    // take some minutes if you have to build the tables and cannot read the tables
+    // from disk
+    auto c = p_cross->second(media.at(component_hash), proposal_energycutsettings[code]);
+
+    // Look which interactions take place and build the corresponding
+    // interaction and secondary builder. The interaction integral will
+    // interpolated too and saved in the calc map by a key build out of a hash
+    // of composed of the component and particle code.
+    auto inter_types = PROPOSAL::CrossSectionVector::GetInteractionTypes(c);
+    calc_[std::make_pair(component_hash, code)] = std::make_tuple(
+        PROPOSAL::make_secondaries(inter_types, particle[code], media.at(component_hash)),
+        PROPOSAL::make_interaction(c, true, true),
+        std::make_unique<LPM_calculator>(media.at(component_hash), code, inter_types));
+  }
+
+  template <typename THadronicLEModel, typename THadronicHEModel>
+  template <typename TStackView>
+  inline ProcessReturn
+  InteractionModel<THadronicLEModel, THadronicHEModel>::doInteraction(
+      TStackView& view, Code const projectileId,
+      [[maybe_unused]] FourMomentum const& projectileP4) {
+
+    auto const projectile = view.getProjectile();
+
+    if (canInteract(projectileId)) {
+
+      // get or build corresponding calculators
+      auto c = getCalculator(projectile, calc_);
+
+      // get the rates of the interaction types for every component.
+      std::uniform_real_distribution<double> distr(0., 1.);
+
+      // sample a interaction-type, loss and component
+      auto rates =
+          std::get<eINTERACTION>(c->second)->Rates(projectile.getEnergy() / 1_MeV);
+      auto [type, target_hash, v] = std::get<eINTERACTION>(c->second)->SampleLoss(
+          projectile.getEnergy() / 1_MeV, rates, distr(RNG_));
+
+      // TODO: This should become obsolete as soon #482 is fixed
+      if (type == PROPOSAL::InteractionType::Undefined) {
+        CORSIKA_LOG_WARN(
+            "PROPOSAL: No particle interaction possible. "
+            "Put initial particle back on stack.");
+        view.addSecondary(std::make_tuple(projectileId,
+                                          projectile.getEnergy() - get_mass(projectileId),
+                                          projectile.getDirection()));
+        return ProcessReturn::Ok;
+      }
+
+      // Read how much random numbers are required to calculate the secondaries.
+      // Calculate the secondaries and deploy them on the corsika stack.
+      auto rnd = std::vector<double>(
+          std::get<eSECONDARIES>(c->second)->RequiredRandomNumbers(type));
+      for (auto& it : rnd) it = distr(RNG_);
+      Point const& place = projectile.getPosition();
+      CoordinateSystemPtr const& labCS = place.getCoordinateSystem();
+
+      auto point = PROPOSAL::Cartesian3D(
+          place.getX(labCS) / 1_cm, place.getY(labCS) / 1_cm, place.getZ(labCS) / 1_cm);
+      auto projectile_dir = projectile.getDirection();
+      auto d = projectile_dir.getComponents(labCS);
+      auto direction = PROPOSAL::Cartesian3D(d.getX().magnitude(), d.getY().magnitude(),
+                                             d.getZ().magnitude());
+
+      auto loss = PROPOSAL::StochasticLoss(
+          static_cast<int>(type), v * projectile.getEnergy() / 1_MeV, point, direction,
+          projectile.getTime() / 1_s, 0., projectile.getEnergy() / 1_MeV);
+      PROPOSAL::Component target;
+      if (type != PROPOSAL::InteractionType::Ioniz)
+        target = PROPOSAL::Component::GetComponentForHash(target_hash);
+
+      auto sec =
+          std::get<eSECONDARIES>(c->second)->CalculateSecondaries(loss, target, rnd);
+
+      // Check for LPM effect suppression!
+
+      if (std::get<eLPM_SUPPRESSION>(c->second)) {
+        auto lpm_suppression = CheckForLPM(
+            *std::get<eLPM_SUPPRESSION>(c->second), projectile.getEnergy(), type, sec,
+            projectile.getNode()->getModelProperties().getMassDensity(place), target, v);
+        if (lpm_suppression) {
+          // Discard interaction - put initial particle back on stack
+          CORSIKA_LOG_DEBUG("LPM suppression detected!");
+          view.addSecondary(std::make_tuple(
+              projectileId, projectile.getEnergy() - get_mass(projectileId),
+              projectile.getDirection()));
+          return ProcessReturn::Ok;
+        }
+      }
+
+      for (auto& s : sec) {
+        auto E = s.energy * 1_MeV;
+        auto vecProposal = s.direction;
+        auto dir = DirectionVector(
+            labCS, {vecProposal.GetX(), vecProposal.GetY(), vecProposal.GetZ()});
+
+        if (static_cast<PROPOSAL::ParticleType>(s.type) ==
+            PROPOSAL::ParticleType::Hadron) {
+          FourMomentum const photonP4(E, E * dir);
+          // for PROPOSAL media target.GetAtomicNum() returns the atomic number in units
+          // of atomic mass, so Nitrogen is 14.0067. When using media in CORSIKA the same
+          // field is filled with the number of nucleons (ie. 14 for Nitrogen) To be sure
+          // we explicitly cast to int
+          auto const A = int(target.GetAtomicNum());
+          auto const Z = int(target.GetNucCharge());
+          Code const targetId = get_nucleus_code(A, Z);
+          CORSIKA_LOGGER_DEBUG(
+              logger_,
+              "photo-hadronic interaction of projectile={} with target={}! Energy={} GeV",
+              projectileId, targetId, E / 1_GeV);
+          this->doHadronicPhotonInteraction(view, labCS, photonP4, targetId);
+        } else {
+          auto sec_code = convert_from_PDG(static_cast<PDGCode>(s.type));
+          // use mass provided by PROPOSAL to ensure correct conversion to kinetic energy
+          auto massProposal =
+              PROPOSAL::ParticleDef::GetParticleDefForType(s.type).mass * 1_MeV;
+          view.addSecondary(std::make_tuple(sec_code, E - massProposal, dir));
+        }
+      }
+    }
+    return ProcessReturn::Ok;
+  }
+
+  template <typename THadronicLEModel, typename THadronicHEModel>
+  template <typename TParticle>
+  inline CrossSectionType
+  InteractionModel<THadronicLEModel, THadronicHEModel>::getCrossSection(
+      TParticle const& projectile, Code const projectileId,
+      FourMomentum const& projectileP4) {
+
+    // ==============================================
+    // this block better diappears. RU 26.10.2021
+    //
+    // determine the volume where the particle is (last) known to be
+    auto const* currentLogicalNode = projectile.getNode();
+    NuclearComposition const& composition =
+        currentLogicalNode->getModelProperties().getNuclearComposition();
+    auto const meanMass = composition.getAverageMassNumber() * constants::u;
+    // ==============================================
+
+    if (canInteract(projectileId)) {
+      auto c = getCalculator(projectile, calc_);
+      return meanMass / (std::get<eINTERACTION>(c->second)->MeanFreePath(
+                             projectileP4.getTimeLikeComponent() / 1_MeV) *
+                         1_g / (1_cm * 1_cm));
+    }
+
+    return CrossSectionType::zero();
+  }
+
+  template <typename THadronicLEModel, typename THadronicHEModel>
+  bool InteractionModel<THadronicLEModel, THadronicHEModel>::CheckForLPM(
+      const LPM_calculator& lpm_calculator, const HEPEnergyType projectile_energy,
+      const PROPOSAL::InteractionType type,
+      const std::vector<PROPOSAL::ParticleState>& sec, const MassDensityType mass_density,
+      const PROPOSAL::Component& target, const double v) {
+
+    double suppression_factor;
+    double current_density = mass_density * ((1_cm * 1_cm * 1_cm) / 1_g);
+
+    if (type == PROPOSAL::InteractionType::Photopair && lpm_calculator.photo_pair_lpm_) {
+      if (sec.size() != 2)
+        throw std::runtime_error(
+            "Invalid number of secondaries for Photopair in CheckForLPM");
+      auto x =
+          sec[0].energy / (sec[0].energy + sec[1].energy); // particle energy asymmetry
+      double density_correction = current_density / lpm_calculator.mass_density_baseline_;
+      suppression_factor = lpm_calculator.photo_pair_lpm_->suppression_factor(
+          projectile_energy / 1_MeV, x, target, density_correction);
+    } else if (type == PROPOSAL::InteractionType::Brems && lpm_calculator.brems_lpm_) {
+      double density_correction = current_density / lpm_calculator.mass_density_baseline_;
+      suppression_factor = lpm_calculator.brems_lpm_->suppression_factor(
+          projectile_energy / 1_MeV, v, target, density_correction);
+    } else if (type == PROPOSAL::InteractionType::Epair && lpm_calculator.epair_lpm_) {
+      if (sec.size() != 3)
+        throw std::runtime_error(
+            "Invalid number of secondaries for EPair in CheckForLPM");
+      double rho = (sec[1].energy - sec[2].energy) /
+                   (sec[1].energy + sec[2].energy); // todo: check
+      // it would be better if these variables are calculated within PROPOSAL
+      double beta = (v * v) / (2 * (1 - v));
+      double xi = (lpm_calculator.particle_mass_ * lpm_calculator.particle_mass_) /
+                  (PROPOSAL::ME * PROPOSAL::ME) * v * v / 4 * (1 - rho * rho) / (1 - v);
+      double density_correction = current_density / lpm_calculator.mass_density_baseline_;
+      suppression_factor = lpm_calculator.epair_lpm_->suppression_factor(
+          projectile_energy / 1_MeV, v, rho * rho, beta, xi, density_correction);
+    } else {
+      return false;
+    }
+
+    // rejection sampling
+    std::uniform_real_distribution<double> distr(0., 1.);
+    auto rnd_num = distr(RNG_);
+
+    CORSIKA_LOGGER_TRACE(logger_,
+                         "Checking for LPM suppression! energy: {}, v: {}, type: {}, "
+                         "target: {}, mass_density: {}, mass_density_baseline: {}, "
+                         "suppression_factor: {}, rnd: {}, result: {}, pmass: {} ",
+                         projectile_energy, v, type, target.GetName(),
+                         mass_density / (1_g / (1_cm * 1_cm * 1_cm)),
+                         lpm_calculator.mass_density_baseline_, suppression_factor,
+                         rnd_num, (rnd_num > suppression_factor),
+                         lpm_calculator.particle_mass_);
+
+    if (rnd_num > suppression_factor)
+      return true;
+    else
+      return false;
+  }
+} // namespace corsika::proposal

corsika/detail/modules/proposal/ProposalProcessBase.inl

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#include <corsika/media/IMediumModel.hpp>
+#include <corsika/media/NuclearComposition.hpp>
+#include <corsika/modules/proposal/ProposalProcessBase.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/utility/COMBoost.hpp>
+#include <corsika/framework/utility/CorsikaData.hpp>
+
+#include <cstdlib>
+#include <iostream>
+#include <limits>
+#include <memory>
+#include <random>
+#include <tuple>
+
+namespace corsika::proposal {
+
+  inline bool ProposalProcessBase::canInteract(Code pcode) const {
+    return std::find(begin(tracked), end(tracked), pcode) != end(tracked);
+  }
+
+  inline HEPEnergyType ProposalProcessBase::getOptimizedEmCut(Code code) const {
+    // get energy above which energy losses need to be considered
+    auto const production_threshold = get_energy_production_threshold(code);
+
+    HEPEnergyType lowest_table_value = 0_GeV;
+
+    // find tables for EnergyCuts closest (but still smaller than) production_threshold
+    for (auto const& table_energy : energycut_table_values) {
+      if (table_energy <= production_threshold && table_energy > lowest_table_value) {
+        lowest_table_value = table_energy;
+      }
+    }
+
+    if (lowest_table_value == 0_GeV) {
+      // no appropriate table available
+      return production_threshold;
+    }
+
+    return lowest_table_value;
+  }
+
+  template <typename TEnvironment>
+  inline ProposalProcessBase::ProposalProcessBase(TEnvironment const& _env) {
+    _env.getUniverse()->walk([&](auto& vtn) {
+      if (vtn.hasModelProperties()) {
+        const auto& prop = vtn.getModelProperties();
+        const auto& medium = mediumData(prop.getMedium());
+
+        auto comp_vec = std::vector<PROPOSAL::Component>();
+        const auto& comp = prop.getNuclearComposition();
+        auto frac_iter = comp.getFractions().cbegin();
+        for (auto const pcode : comp.getComponents()) {
+          comp_vec.emplace_back(std::string(get_name(pcode)), get_nucleus_Z(pcode),
+                                get_nucleus_A(pcode), *frac_iter);
+          ++frac_iter;
+        }
+
+        media[comp.getHash()] = PROPOSAL::Medium(
+            medium.getName(), medium.getIeff(), -medium.getCbar(), medium.getAA(),
+            medium.getSK(), medium.getX0(), medium.getX1(), medium.getDlt0(),
+            medium.getCorrectedDensity(), comp_vec);
+      }
+    });
+
+    //! If corsika data exist store interpolation tables to the corresponding
+    //! path, otherwise interpolation tables would only stored in main memory if
+    //! no explicit intrpolation def is specified.
+    PROPOSAL::InterpolationSettings::TABLES_PATH = corsika_data("PROPOSAL").c_str();
+
+    //! Initialize EnergyCutSettings
+    for (auto const particle_code : tracked) {
+      if (particle_code == Code::Photon) {
+        // no EnergyCut for photon, only-stochastic propagation
+        continue;
+      }
+      // use optimized emcut for which tables should be available
+      auto optimized_ecut = getOptimizedEmCut(particle_code);
+      CORSIKA_LOG_INFO("PROPOSAL: Use tables with EnergyCut of {} for particle {}.",
+                       optimized_ecut, particle_code);
+      proposal_energycutsettings[particle_code] = optimized_ecut;
+    }
+
+    //! Initialize PROPOSAL tables for all media and all particles
+    for (auto const& medium : media) {
+      for (auto const particle_code : tracked) {
+        buildTables(medium.second, particle_code,
+                    proposal_energycutsettings[particle_code]);
+      }
+    }
+  }
+
+  void ProposalProcessBase::buildTables(PROPOSAL::Medium medium, Code particle_code,
+                                        HEPEnergyType emcut) {
+    CORSIKA_LOG_DEBUG("PROPOSAL: Initialize tables for particle {} in {}", particle_code,
+                      medium.GetName());
+    cross[particle_code](medium, emcut);
+  }
+
+  inline size_t ProposalProcessBase::hash::operator()(const calc_key_t& p) const
+      noexcept {
+    return p.first ^ std::hash<Code>{}(p.second);
+  }
+
+} // namespace corsika::proposal

corsika/detail/modules/pythia8/Decay.inl

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#include <corsika/modules/pythia8/Pythia8.hpp>
+#include <corsika/modules/pythia8/Decay.hpp>
+#include <corsika/modules/pythia8/Random.hpp>
+
+#include <corsika/framework/utility/COMBoost.hpp>
+
+namespace corsika::pythia8 {
+
+  inline Decay::Decay(bool const print_listing)
+      : Pythia8::Pythia(CORSIKA_Pythia8_XML_DIR)
+      , print_listing_(print_listing) {
+    init();
+  }
+
+  inline Decay::Decay(std::set<Code> const& those)
+      : handleAllDecays_(false)
+      , handledDecays_(those) {
+    init();
+  }
+
+  inline Decay::~Decay() { CORSIKA_LOG_INFO("Pythia::Decay n={}", count_); }
+
+  inline void Decay::init() {
+
+    // run this only once during construction
+
+    // link random number generator in pythia to CORSIKA8
+    Pythia8::Pythia::setRndmEnginePtr(std::make_shared<corsika::pythia8::Random>());
+
+    Pythia8::Pythia::readString("Next:numberShowInfo = 0");
+    Pythia8::Pythia::readString("Next:numberShowProcess = 0");
+    Pythia8::Pythia::readString("Next:numberShowEvent = 0");
+
+    Pythia8::Pythia::readString("Print:quiet = on");
+    Pythia8::Pythia::readString("Check:particleData = off");
+
+    /*
+       switching off event check in pythia is needed to allow decays that are off-shell
+       according to the mass definition in pythia.
+       the consistency of particle masses between event generators is an unsolved issues
+    */
+    CORSIKA_LOG_INFO("Pythia::Init: switching off event checking in pythia..");
+    Pythia8::Pythia::readString("Check:event = 1");
+
+    Pythia8::Pythia::readString("ProcessLevel:all = off");
+    Pythia8::Pythia::readString("ProcessLevel:resonanceDecays = on");
+
+    // making sure
+    setStable(Code::Pi0);
+
+    //    Pythia8::Pythia::particleData.readString("59:m0 = 101.00");
+
+    // LCOV_EXCL_START, we don't validate pythia8 internals
+    if (!Pythia8::Pythia::init())
+      throw std::runtime_error("Pythia::Decay: Initialization failed!");
+    // LCOV_EXCL_STOP
+  }
+
+  inline bool Decay::canHandleDecay(Code const vParticleCode) {
+    // if known to pythia and not proton, electron or neutrino it can decay
+    if (vParticleCode == Code::Proton || vParticleCode == Code::AntiProton ||
+        vParticleCode == Code::NuE || vParticleCode == Code::NuMu ||
+        vParticleCode == Code::NuTau || vParticleCode == Code::NuEBar ||
+        vParticleCode == Code::NuMuBar || vParticleCode == Code::NuTauBar ||
+        vParticleCode == Code::Electron || vParticleCode == Code::Positron)
+      return false;
+    else if (canDecay(vParticleCode)) // check pythia8 internal
+      return true;
+    else
+      return false;
+  }
+
+  inline void Decay::setHandleDecay(Code const vParticleCode) {
+    handleAllDecays_ = false;
+    CORSIKA_LOG_DEBUG("Pythia::Decay: set to handle decay of {} ", vParticleCode);
+    if (Decay::canHandleDecay(vParticleCode))
+      handledDecays_.insert(vParticleCode);
+    else
+      throw std::runtime_error("this decay can not be handled by pythia!");
+  }
+
+  inline void Decay::setHandleDecay(std::vector<Code> const& vParticleList) {
+    handleAllDecays_ = false;
+    for (auto const p : vParticleList) setHandleDecay(p);
+  }
+
+  inline bool Decay::isDecayHandled(Code const vParticleCode) {
+    if (handleAllDecays_ && canHandleDecay(vParticleCode))
+      return true;
+    else
+      return handledDecays_.find(vParticleCode) != Decay::handledDecays_.end();
+  }
+
+  inline void Decay::setStable(std::vector<Code> const& particleList) {
+    for (auto const p : particleList) Decay::setStable(p);
+  }
+
+  inline void Decay::setUnstable(Code const pCode) {
+    CORSIKA_LOG_DEBUG("Pythia::Decay: setting {} unstable..", pCode);
+    Pythia8::Pythia::particleData.mayDecay(static_cast<int>(get_PDG(pCode)), true);
+  }
+
+  inline void Decay::setStable(Code const pCode) {
+    CORSIKA_LOG_DEBUG("Pythia::Decay: setting {} stable..", pCode);
+    Pythia8::Pythia::particleData.mayDecay(static_cast<int>(get_PDG(pCode)), false);
+  }
+
+  inline bool Decay::isStable(Code const vCode) {
+    return Pythia8::Pythia::particleData.canDecay(static_cast<int>(get_PDG(vCode)));
+  }
+
+  inline bool Decay::canDecay(Code const pCode) {
+    bool const ans =
+        Pythia8::Pythia::particleData.canDecay(static_cast<int>(get_PDG(pCode)));
+    CORSIKA_LOG_DEBUG("Pythia::Decay: checking if particle: {} can decay in PYTHIA? {} ",
+                      pCode, ans);
+    return ans;
+  }
+
+  inline void Decay::printDecayConfig(Code const vCode) {
+    CORSIKA_LOG_INFO("Decay: Pythia decay configuration:");
+    CORSIKA_LOG_INFO(" {} is {} ", vCode, (isStable(vCode) ? "stable" : "unstable"));
+  }
+
+  inline void Decay::printDecayConfig() {
+    CORSIKA_LOG_INFO("Pythia::Decay: decay configuration:");
+    if (handleAllDecays_)
+      CORSIKA_LOG_INFO(" all particles known to Pythia are handled by Pythia::Decay!");
+    else
+      for (auto const pCode : handledDecays_)
+        CORSIKA_LOG_INFO("Decay of {} is handled by Pythia!", pCode);
+  }
+
+  template <typename TParticle>
+  inline TimeType Decay::getLifetime(TParticle const& particle) {
+
+    auto const pid = particle.getPID();
+    if (canDecay(pid)) {
+      HEPEnergyType E = particle.getEnergy();
+      HEPMassType m = particle.getMass();
+
+      double const gamma = E / m;
+
+      TimeType const t0 = get_lifetime(pid);
+      auto const lifetime = gamma * t0;
+      CORSIKA_LOG_TRACE("Pythia::Decay: code: {}", particle.getPID());
+      CORSIKA_LOG_TRACE("Pythia::Decay: MinStep: t0: {}", t0);
+      CORSIKA_LOG_TRACE("Pythia::Decay: MinStep: energy: {} GeV", E / 1_GeV);
+      CORSIKA_LOG_TRACE("Pythia::Decay: momentum: {} GeV",
+                        particle.getMomentum().getComponents() / 1_GeV);
+      CORSIKA_LOG_TRACE("Pythia::Decay: MinStep: gamma: {}", gamma);
+      CORSIKA_LOG_TRACE("Pythia::Decay: MinStep: tau: {} ", lifetime);
+
+      return lifetime;
+    } else
+      return std::numeric_limits<double>::infinity() * 1_s;
+  }
+
+  template <typename TView>
+  inline void Decay::doDecay(TView& view) {
+
+    auto projectile = view.getProjectile();
+
+    auto const& labMomentum = projectile.getMomentum();
+    [[maybe_unused]] CoordinateSystemPtr const& labCS = labMomentum.getCoordinateSystem();
+
+    // define target kinematics in lab frame
+    // define boost to and from CoM frame
+    // CoM frame definition in Pythia projectile: +z
+    COMBoost const boost(labMomentum, projectile.getMass());
+    auto const& rotatedCS = boost.getRotatedCS();
+
+    count_++;
+
+    // pythia stack
+    Pythia8::Event& event = Pythia8::Pythia::event;
+    event.reset();
+
+    auto const particleId = projectile.getPID();
+
+    // set particle unstable
+    Decay::setUnstable(particleId);
+
+    // input particle PDG
+    auto const pdgCode = static_cast<int>(get_PDG(particleId));
+
+    double constexpr px = 0;
+    double constexpr py = 0;
+    double constexpr pz = 0;
+    double const en = projectile.getMass() / 1_GeV;
+    double const m = en;
+
+    // add particle to pythia stack
+    event.append(pdgCode, 1, 0, 0, // PID, status, col, acol
+                 px, py, pz, en, m);
+
+    // LCOV_EXCL_START, we don't validate pythia8 internals
+    if (!Pythia8::Pythia::next())
+      throw std::runtime_error("Pythia::Decay: decay failed!");
+    // LCOV_EXCL_STOP
+
+    CORSIKA_LOG_DEBUG("Pythia::Decay: particles after decay: {} ", event.size());
+
+    // LCOV_EXCL_START, we don't validate pythia8 internals
+    if (print_listing_) {
+      // list final state
+      event.list();
+    }
+    // LCOV_EXCL_STOP
+
+    // loop over final state
+    for (int i = 0; i < event.size(); ++i) {
+      if (event[i].isFinal()) {
+        try {
+          auto const id = static_cast<PDGCode>(event[i].id());
+          auto const pyId = convert_from_PDG(id);
+          HEPEnergyType const Erest = event[i].e() * 1_GeV;
+          MomentumVector const pRest(
+              rotatedCS,
+              {event[i].px() * 1_GeV, event[i].py() * 1_GeV, event[i].pz() * 1_GeV});
+          FourVector const fourMomRest{Erest, pRest};
+          auto const fourMomLab = boost.fromCoM(fourMomRest);
+          auto const p3 = fourMomLab.getSpaceLikeComponents();
+
+          HEPEnergyType const mass = get_mass(pyId);
+          HEPEnergyType const Ekin = sqrt(p3.getSquaredNorm() + mass * mass) - mass;
+
+          CORSIKA_LOG_TRACE(
+              "particle: id={} momentum={} energy={} ", pyId,
+              fourMomLab.getSpaceLikeComponents().getComponents(labCS) / 1_GeV,
+              fourMomLab.getTimeLikeComponent());
+
+          view.addSecondary(std::make_tuple(pyId, Ekin, p3.normalized()));
+        } catch (std::out_of_range const& ex) {
+          CORSIKA_LOG_CRITICAL("Pythia ID {} unknown in C8", event[i].id());
+          throw ex;
+        }
+      }
+    }
+
+    // set particle stable
+    Decay::setStable(particleId);
+  }
+} // namespace corsika::pythia8

corsika/detail/modules/pythia8/InteractionModel.inl

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <utility>
+#include <stdexcept>
+
+#include <boost/filesystem/path.hpp>
+#include <boost/iterator/transform_iterator.hpp>
+#include <fmt/core.h>
+#include <cnpy.hpp>
+
+#include <corsika/modules/pythia8/InteractionModel.hpp>
+
+#include <corsika/framework/geometry/FourVector.hpp>
+#include <corsika/framework/utility/COMBoost.hpp>
+#include <corsika/framework/utility/CrossSectionTable.hpp>
+#include <corsika/framework/core/EnergyMomentumOperations.hpp>
+#include <corsika/media/Environment.hpp>
+#include <corsika/media/NuclearComposition.hpp>
+
+namespace corsika::pythia8 {
+
+  inline InteractionModel::~InteractionModel() {}
+
+  inline InteractionModel::InteractionModel(std::set<Code> const& stableParticles,
+                                            boost::filesystem::path const& dataPath,
+                                            bool const printListing)
+      : crossSectionPPElastic_{loadPPTable(dataPath, "sig_el")}
+      , crossSectionPPInelastic_{loadPPTable(dataPath, "sig_inel")}
+      , print_listing_(printListing) {
+    auto rndm = std::make_shared<corsika::pythia8::Random>();
+    pythia_.setRndmEnginePtr(rndm);
+
+    CORSIKA_LOG_INFO("Pythia8 reuse files: {}", dataPath.native());
+
+    // projectile and target init to p Nitrogen to initialize pythia with angantyr
+    pythia_.readString("Beams:allowIDASwitch = on");
+    pythia_.settings.mode("Beams:idA",
+                          static_cast<PDGCodeIntType>(get_PDG(Code::Proton)));
+    pythia_.settings.mode("Beams:idB",
+                          static_cast<PDGCodeIntType>(get_PDG(Code::Oxygen)));
+
+    pythia_.readString("Beams:allowVariableEnergy = on");
+    pythia_.readString("Beams:frameType = 1"); // CoM
+    // initialize a maximum energy event
+    pythia_.settings.parm("Beams:eCM", 400_TeV / 1_GeV);
+
+    // needed when regular Pythia (not Angantyr) takes over for pp
+    pythia_.readString("SoftQCD:all = on");
+
+    /* reuse file settings
+    reuseInit = 2: initialization data is reuse file, but if the file is not found,
+    initialization fails reuseInit = 3: same, but if the file is not found, it will be
+    generated and saved after normal initialization */
+    pythia_.readString("MultipartonInteractions:reuseInit = 2");
+    pythia_.readFile(dataPath.native() + "/pythia8313.cmnd");
+    // switch on decays for all hadrons except the ones defined as tracked by C8
+    if (!stableParticles.empty()) {
+      pythia_.readString("HadronLevel:Decay = on");
+      for (auto pCode : stableParticles)
+        pythia_.particleData.mayDecay(static_cast<int>(get_PDG(pCode)), false);
+    } else {
+      // all hadrons stable
+      pythia_.readString("HadronLevel:Decay = on");
+    }
+    // Reduce printout and relax energy-momentum conservation.
+    pythia_.readString("Print:quiet = on");
+    pythia_.readString("Check:epTolErr = 0.1");
+    pythia_.readString("Check:epTolWarn = 0.0001");
+    pythia_.readString("Check:mTolErr = 0.01");
+
+    // Reduce statistics printout to relevant ones.
+    pythia_.readString("Stat:showProcessLevel = off");
+    pythia_.readString("Stat:showPartonLevel = off");
+
+    // initialize
+    // we can't test this block, LCOV_EXCL_START
+    if (!pythia_.init())
+      throw std::runtime_error("Pythia::InteractionModel: Initialization failed!");
+    // LCOV_EXCL_STOP
+
+    // create/load tables of cross sections
+    for (Code const projectile : validProjectiles_) {
+      for (Code const target : validTargets_) {
+        if (xs_map_.find(projectile) == xs_map_.end()) {
+          // projectile not mapped, load table directly
+          auto const tablePath =
+              dataPath / "pythia8_xs_npz" /
+              fmt::format("xs_{}_{}.npz",
+                          static_cast<PDGCodeIntType>(get_PDG(projectile)),
+                          static_cast<PDGCodeIntType>(get_PDG(target)));
+          auto const tables = cnpy::npz_load(tablePath.native());
+          // NOTE: tables are calculated for plab. In C8 we we assume elab. starting at
+          // plab=100GeV the difference is at most (1+m**2/p**2)=1.000088035
+          auto const energies = tables.at("plab").as_vec<float>();
+          auto const total_xs = tables.at("sig_tot").as_vec<float>();
+
+          if (auto const e_size = energies.size(), xs_size = total_xs.size();
+              xs_size != e_size) {
+            throw std::runtime_error{
+                fmt::format("Pythia8 table corrupt (plab size = {}; sig_tot size = {})",
+                            e_size, xs_size)};
+          }
+
+          auto xs_it = boost::make_transform_iterator(total_xs.cbegin(), millibarn_mult);
+          auto e_it_beg = boost::make_transform_iterator(energies.cbegin(), GeV_mult);
+          decltype(e_it_beg) e_it_end =
+              boost::make_transform_iterator(energies.cend(), GeV_mult);
+          crossSectionTables_.insert(
+              {std::pair{projectile, target},
+               CrossSectionTable<InterpolationTransforms::Log>{
+                   std::move(e_it_beg), std::move(e_it_end), std::move(xs_it)}});
+        }
+      }
+    }
+  }
+
+  inline CrossSectionTable<InterpolationTransforms::Log> InteractionModel::loadPPTable(
+      boost::filesystem::path const& dataPath, char const* key) {
+    auto const ppTablePath =
+        dataPath / "pythia8_xs_npz" /
+        fmt::format("xs_{}_{}.npz", static_cast<PDGCodeIntType>(get_PDG(Code::Proton)),
+                    static_cast<PDGCodeIntType>(get_PDG(Code::Proton)));
+    auto const tables = cnpy::npz_load(ppTablePath.native());
+    // NOTE: tables are calculated for plab. In C8 we we assume elab. starting at
+    // plab=100GeV the difference is at most (1+m**2/p**2)=1.000088035
+    auto const energies = tables.at("plab").as_vec<float>();
+    auto const xs = tables.at(key).as_vec<float>();
+
+    if (auto e_size = energies.size(), xs_size = xs.size(); xs_size != e_size) {
+      throw std::runtime_error{fmt::format(
+          "Pythia8 table corrupt (plab size = {}; xs size = {})", e_size, xs_size)};
+    }
+
+    auto xs_it = boost::make_transform_iterator(xs.cbegin(), millibarn_mult);
+    auto e_it_beg = boost::make_transform_iterator(energies.cbegin(), GeV_mult);
+    decltype(e_it_beg) e_it_end =
+        boost::make_transform_iterator(energies.cend(), GeV_mult);
+    return CrossSectionTable<InterpolationTransforms::Log>{
+        std::move(e_it_beg), std::move(e_it_end), std::move(xs_it)};
+  }
+
+  inline bool InteractionModel::isValid(Code const projectileId, Code const targetId,
+                                        HEPEnergyType const sqrtSNN) const {
+
+    CORSIKA_LOG_DEBUG("pythia isValid: {} + {} at sqrtSNN = {} GeV", projectileId,
+                      targetId, sqrtSNN / 1_GeV);
+    if (is_nucleus(targetId))
+      CORSIKA_LOG_DEBUG("nucleus = {}-{}", get_nucleus_A(targetId),
+                        get_nucleus_Z(targetId));
+    if (is_nucleus(projectileId)) // not yet possible with Pythia
+      return false;
+
+    if (!canInteract(projectileId)) return false;
+
+    auto const mass_target =
+        get_mass(targetId) / (is_nucleus(targetId) ? get_nucleus_A(targetId) : 1.);
+    HEPEnergyType const labE =
+        calculate_lab_energy(static_pow<2>(sqrtSNN), get_mass(projectileId), mass_target);
+    if (labE < eKinMinLab_) return false;
+
+    bool const validProjectile =
+        std::find(validProjectiles_.begin(), validProjectiles_.end(), projectileId) !=
+        validProjectiles_.end();
+    bool const validTarget = std::find(validTargets_.begin(), validTargets_.end(),
+                                       targetId) != validTargets_.end();
+    return validProjectile && validTarget;
+  }
+
+  inline bool InteractionModel::canInteract(Code const pCode) const {
+    return is_hadron(pCode) && !is_nucleus(pCode);
+  }
+
+  inline std::tuple<CrossSectionType, CrossSectionType>
+  InteractionModel::getCrossSectionInelEla(Code const projectileId, Code const targetId,
+                                           FourMomentum const& projectileP4,
+                                           FourMomentum const& targetP4) const {
+    // elastic cross sections only available for proton
+    if (projectileId != Code::Proton || targetId != Code::Proton) {
+      return std::tuple{CrossSectionType::zero(), CrossSectionType::zero()};
+    }
+
+    // NOTE: here we calculate SNN (per nucleon energy) AND S (per particle energy)
+    // because isValid expects sqrtSNN as input but the tables need Elab
+    auto const Aprojectile =
+        (is_nucleus(projectileId) ? get_nucleus_A(projectileId) : 1.);
+    auto const Atarget = (is_nucleus(targetId) ? get_nucleus_A(targetId) : 1.);
+    auto const SNN = (projectileP4 / Aprojectile + targetP4 / Atarget).getNormSqr();
+    HEPEnergyType const sqrtSNN = sqrt(SNN);
+    if (!isValid(projectileId, targetId, sqrtSNN))
+      return std::tuple{CrossSectionType::zero(), CrossSectionType::zero()};
+
+    // read cross section from table
+    // define system (the tables were created for lab. energies) since we cannot assume
+    // the input is in the lab. frame we calculate the lab. energy from SNN =
+    // (projectileP4 / Aprojectile + targetP4 / Atarget)**2 and Elab = S/(2. * mTarget /
+    // Atarget) where A* is the number of nucleons in a nucleus. A* is 1 if projectile or
+    // target are simple hadrons.
+    HEPEnergyType const Elab = calculate_lab_energy(
+        SNN, get_mass(projectileId) / Aprojectile, get_mass(targetId) / Atarget);
+
+    CrossSectionType const xsEl = crossSectionPPElastic_.interpolate(Elab);
+    CrossSectionType const xsInel = crossSectionPPInelastic_.interpolate(Elab);
+    return std::pair{xsInel, xsEl};
+  }
+
+  inline CrossSectionType InteractionModel::getCrossSection(
+      Code const projectileId, Code const targetId, FourMomentum const& projectileP4,
+      FourMomentum const& targetP4) const {
+
+    auto const Aprojectile =
+        (is_nucleus(projectileId) ? get_nucleus_A(projectileId) : 1.);
+    auto const Atarget = (is_nucleus(targetId) ? get_nucleus_A(targetId) : 1.);
+
+    auto const SNN = (projectileP4 / Aprojectile + targetP4 / Atarget).getNormSqr();
+    HEPEnergyType const sqrtSNN = sqrt(SNN);
+    if (!isValid(projectileId, targetId, sqrtSNN)) return CrossSectionType::zero();
+
+    // read cross section from table
+    // define system (the tables were created for lab. energies) since we cannot assume
+    // the input is in the lab. frame we calculate the lab. energy from SNN =
+    // (projectileP4 / Aprojectile + targetP4 / Atarget)**2 and Elab = S/(2. * mTarget /
+    // Atarget) where A* is the number of nucleons in a nucleus. A* is 1 if projectile or
+    // target are simple hadrons.
+    auto const it = xs_map_.find(projectileId);
+    auto const mapped_projectile = (it == xs_map_.end()) ? projectileId : it->second;
+    auto const& table = crossSectionTables_.at(std::pair{mapped_projectile, targetId});
+    HEPEnergyType const Elab = calculate_lab_energy(
+        SNN, get_mass(projectileId) / Aprojectile, get_mass(targetId) / Atarget);
+    CORSIKA_LOG_DEBUG("pythia getCrossSection: {}+{} at Elab= {} GeV, sqrtSNN = {} GeV",
+                      projectileId, get_name(targetId), Elab / 1_GeV, sqrtSNN / 1_GeV);
+    return table.interpolate(Elab);
+  }
+
+  template <class TView>
+  inline void InteractionModel::doInteraction(TView& view, Code const projectileId,
+                                              Code const targetId,
+                                              FourMomentum const& projectileP4,
+                                              FourMomentum const& targetP4) {
+
+    CORSIKA_LOG_DEBUG(
+        "Pythia::InteractionModel: "
+        "doInteraction: {} interaction? {} ",
+        projectileId, corsika::pythia8::InteractionModel::canInteract(projectileId));
+
+    // define system nucleon-nucleon
+    auto const projectileP4NN =
+        projectileP4 / (is_nucleus(projectileId) ? get_nucleus_A(projectileId) : 1);
+    auto const targetP4NN =
+        targetP4 / (is_nucleus(targetId) ? get_nucleus_A(targetId) : 1);
+    auto const SNN = (projectileP4NN + targetP4NN).getNormSqr();
+    auto const sqrtSNN = sqrt(SNN);
+    HEPEnergyType const eProjectileLab = calculate_lab_energy(
+        SNN, get_mass(projectileId), get_mass(targetId) / get_nucleus_A(targetId));
+
+    CORSIKA_LOG_DEBUG("InteractionModel: ebeam lab: {} GeV", eProjectileLab / 1_GeV);
+    CORSIKA_LOG_DEBUG(
+        "InteractionModel: "
+        " doInteraction: E(GeV): {}"
+        " Ecm_NN(GeV): {}",
+        eProjectileLab / 1_GeV, sqrtSNN / 1_GeV);
+
+    if (!isValid(projectileId, targetId, sqrtSNN)) {
+      throw std::runtime_error("invalid target,projectile,energy combination.");
+    }
+
+    COMBoost const boost{projectileP4NN, targetP4NN};
+    auto const& rotCS = boost.getRotatedCS();
+
+    // variables to talk to Pythia
+    double const eCM_pythia = sqrtSNN / 1_GeV; // CoM energy hadron-Nucleon
+    double const idA_pythia = static_cast<double>(get_PDG(projectileId));
+    double const idB_pythia = static_cast<double>(get_PDG(targetId));
+    CORSIKA_LOG_DEBUG("re-configuring pythia idA={}, idB={}, ecm={} GeV", idA_pythia,
+                      idB_pythia, eCM_pythia);
+    // configure event
+    pythia_.setBeamIDs(idA_pythia, idB_pythia);
+    pythia_.setKinematics(eCM_pythia);
+
+    // generate event (one chance only!)
+    if (!pythia_.next()) {
+      throw std::runtime_error("Pythia::InteractionModel: interaction failed!");
+    }
+
+    // LCOV_EXCL_START, we don't validate pythia8 internals
+    if (print_listing_) {
+      // print final state
+      pythia_.event.list();
+    }
+    // LCOV_EXCL_STOP
+
+    MomentumVector Plab_final{boost.getOriginalCS()};
+    auto Elab_final = HEPEnergyType::zero();
+
+    for (Pythia8::Particle const& p8p : pythia_.event) {
+      // skip particles that have decayed and initial particles
+      if (!p8p.isFinal()) continue;
+      try {
+        auto const volatile id = static_cast<PDGCode>(p8p.id());
+        auto const pyId = convert_from_PDG(id);
+
+        // skip nuclear remnants
+        if (is_nucleus(pyId)) continue;
+
+        MomentumVector const pyPcom(
+            rotCS, {p8p.px() * 1_GeV, p8p.py() * 1_GeV, p8p.pz() * 1_GeV});
+        auto const pyP = boost.fromCoM(FourVector{p8p.e() * 1_GeV, pyPcom});
+
+        HEPEnergyType const mass = get_mass(pyId);
+        HEPEnergyType const Ekin =
+            sqrt(pyP.getSpaceLikeComponents().getSquaredNorm() + mass * mass) - mass;
+
+        // add to corsika stack
+        auto pnew = view.addSecondary(
+            std::make_tuple(pyId, Ekin, pyP.getSpaceLikeComponents().normalized()));
+
+        Plab_final += pnew.getMomentum();
+        Elab_final += pnew.getEnergy();
+      }
+      // irreproducible in tests, LCOV_EXCL_START
+      catch (std::out_of_range const& ex) {
+        CORSIKA_LOG_CRITICAL("Pythia ID {} unknown in C8", p8p.id());
+        throw ex;
+      }
+      // LCOV_EXCL_STOP
+    }
+
+    pythia_.event.clear();
+
+    CORSIKA_LOG_DEBUG(
+        "conservation (all GeV): "
+        "Elab_final= {}"
+        ", Plab_final= {}",
+        Elab_final / 1_GeV, (Plab_final / 1_GeV).getComponents());
+  }
+
+} // namespace corsika::pythia8

corsika/detail/modules/pythia8/NeutrinoInteraction.inl

+/*
+ * (c) Copyright 2024 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <boost/filesystem/path.hpp>
+#include "corsika/framework/core/EnergyMomentumOperations.hpp"
+
+namespace corsika::pythia8 {
+
+  inline NeutrinoInteraction::NeutrinoInteraction(std::set<Code> const& list,
+                                                  bool const& handleNC,
+                                                  bool const& handleCC)
+      : stable_particles_(list)
+      , handle_nc_(handleNC)
+      , handle_cc_(handleCC)
+      , pythiaMain_{CORSIKA_Pythia8_XML_DIR, false} {
+
+    pythiaMain_.setRndmEnginePtr(std::make_shared<corsika::pythia8::Random>());
+  }
+
+  inline NeutrinoInteraction::~NeutrinoInteraction() {
+    CORSIKA_LOGGER_DEBUG(logger_, "Pythia::NeutrinoInteraction n= {}", count_);
+  }
+
+  template <class TView>
+  void NeutrinoInteraction::doInteraction(TView& view, Code const projectileId,
+                                          Code const targetId,
+                                          FourMomentum const& projectileP4,
+                                          FourMomentum const& targetP4) {
+
+    CORSIKA_LOGGER_DEBUG(logger_, "Primary {} - {} interaction with E_nu = {} GeV",
+                         projectileId, targetId,
+                         projectileP4.getTimeLikeComponent() / 1_GeV);
+    CORSIKA_LOGGER_DEBUG(
+        logger_, "configure Pythia for primary neutrino interactions. NC={}, CC={}",
+        handle_nc_, handle_cc_);
+    if (!handle_nc_ && !handle_cc_) {
+      CORSIKA_LOGGER_ERROR(
+          logger_,
+          "no neutrino interaction channel configured! Select either NC, CC or both!");
+      throw std::runtime_error("Configuration error!");
+    }
+    CORSIKA_LOGGER_DEBUG(logger_, "minimal Q2 in DIS: {} GeV2", minQ2_ / 1_GeV / 1_GeV);
+
+    if (!isValid(projectileId, targetId, projectileP4, targetP4)) {
+      CORSIKA_LOGGER_ERROR(logger_, "wrong projectile, target or energy configuration!");
+      throw std::runtime_error("Configuration error!");
+    }
+    // sample nucleon from nucleus A,Z
+    double const fProtons = get_nucleus_Z(targetId) / double(get_nucleus_A(targetId));
+    double const fNeutrons = 1. - fProtons;
+    std::discrete_distribution<int> nucleonChannelDist{fProtons, fNeutrons};
+    corsika::default_prng_type& rng =
+        corsika::RNGManager<>::getInstance().getRandomStream("pythia");
+    Code const targetNucleonId = (nucleonChannelDist(rng) ? Code::Neutron : Code::Proton);
+    int const idTarget_pythia = static_cast<int>(get_PDG(targetNucleonId));
+    auto const targetP4NN =
+        targetP4 / (is_nucleus(targetId) ? get_nucleus_A(targetId) : 1);
+    CORSIKA_LOGGER_DEBUG(logger_, "selected {} target. P4NN={}", targetNucleonId,
+                         targetP4NN.getSpaceLikeComponents());
+
+    // set projectile
+    double const Q2min = minQ2_ / 1_GeV / 1_GeV;
+    int const idProjectile_pythia = static_cast<int>(get_PDG(projectileId));
+
+    // calculate CoM energy of the neutrino nucleon interaction
+    double const ecm_pythia =
+        calculate_com_energy(projectileP4.getTimeLikeComponent(), get_mass(projectileId),
+                             get_mass(targetNucleonId)) /
+        1_GeV;
+
+    CORSIKA_LOGGER_DEBUG(logger_, "center-of-mass energy: {} GeV", ecm_pythia);
+    // Set up CoM interaction
+    pythiaMain_.readString("Beams:frameType = 1");
+    // center-of-mass energy
+    pythiaMain_.settings.parm("Beams:eCM", ecm_pythia);
+    // BeamA (+z in CoM)
+    pythiaMain_.settings.mode("Beams:idA", idProjectile_pythia);
+    // BeamB (-z direction in CoM)
+    pythiaMain_.settings.mode("Beams:idB", idTarget_pythia);
+
+    // Set up DIS process within some phase space.
+    // Neutral current (with gamma/Z interference).
+    if (handle_nc_) pythiaMain_.readString("WeakBosonExchange:ff2ff(t:gmZ) = on");
+    // Charged current.
+    if (handle_cc_) pythiaMain_.readString("WeakBosonExchange:ff2ff(t:W) = on");
+    // Phase-space cut: minimal Q2 of process.
+    pythiaMain_.settings.parm("PhaseSpace:Q2Min", Q2min);
+
+    // Set dipole recoil on. Necessary for DIS + shower.
+    pythiaMain_.readString("SpaceShower:dipoleRecoil = on");
+
+    // Allow emissions up to the kinematical limit,
+    // since rate known to match well to matrix elements everywhere.
+    pythiaMain_.readString("SpaceShower:pTmaxMatch = 2");
+
+    // QED radiation off lepton not handled yet by the new procedure.
+    pythiaMain_.readString("PDF:lepton = off");
+    pythiaMain_.readString("TimeShower:QEDshowerByL = off");
+
+    // switch on decays for all hadrons except the ones defined as tracked by C8
+    if (!stable_particles_.empty()) {
+      pythiaMain_.readString("HadronLevel:Decay = on");
+      for (auto pCode : stable_particles_)
+        pythiaMain_.particleData.mayDecay(static_cast<int>(get_PDG(pCode)), false);
+    } else {
+      // all hadrons stable
+      pythiaMain_.readString("HadronLevel:Decay = off");
+    }
+
+    pythiaMain_.readString("Stat:showProcessLevel = off");
+    pythiaMain_.readString("Stat:showPartonLevel = off");
+
+    pythiaMain_.readString("Print:quiet = on");
+    pythiaMain_.readString("Check:epTolErr = 0.1");
+    pythiaMain_.readString("Check:epTolWarn = 0.0001");
+    pythiaMain_.readString("Check:mTolErr = 0.01");
+
+    pythiaMain_.init();
+
+    // References to the event record
+    Pythia8::Event& eventMain = pythiaMain_.event;
+    // define boost assuming target nucleon is at rest
+    COMBoost const boost{projectileP4, targetP4NN};
+    // the boost is along the total momentum axis and does not include a rotation
+    // get rotated frame where momentum of projectile in the CoM is along +z
+    auto const& rotCS = boost.getRotatedCS();
+
+    if (!pythiaMain_.next()) {
+      throw std::runtime_error("Pythia neutrino collision failed ");
+    } else {
+      CORSIKA_LOGGER_DEBUG(logger_, "pythia neutrino interaction done!");
+    }
+
+    MomentumVector Plab_final{boost.getOriginalCS()};
+    auto Elab_final = HEPEnergyType::zero();
+    CORSIKA_LOGGER_DEBUG(logger_, "particles generated in neutrino interaction:");
+    for (int i = 0; i < eventMain.size(); ++i) {
+      auto const& p8p = eventMain[i];
+      if (p8p.isFinal()) {
+        try {
+          // get particle ids from pythia stack and convert to corsika id
+          auto const volatile id = static_cast<PDGCode>(p8p.id());
+          auto const pyId = convert_from_PDG(id);
+          // get pythia momentum in CoM (rotated cs!)
+          MomentumVector const pyPcom(
+              rotCS, {p8p.px() * 1_GeV, p8p.py() * 1_GeV, p8p.pz() * 1_GeV});
+          // boost pythia 4Momentum in CoM to lab. frame. This includes the rotation.
+          // calculate 3-momentum and kinetic energy
+          CORSIKA_LOGGER_DEBUG(logger_, "momentum in CoM: {} GeV",
+                               pyPcom.getComponents() / 1_GeV);
+          auto const pyP4lab = boost.fromCoM(FourVector{p8p.e() * 1_GeV, pyPcom});
+          auto const pyPlab = pyP4lab.getSpaceLikeComponents();
+          HEPEnergyType const pyEkinlab =
+              calculate_kinetic_energy(pyPlab.getNorm(), get_mass(pyId));
+
+          // add to corsika stack
+          auto pnew =
+              view.addSecondary(std::make_tuple(pyId, pyEkinlab, pyPlab.normalized()));
+
+          CORSIKA_LOGGER_DEBUG(logger_, "id = {}, E = {} GeV, p = {} GeV", pyId,
+                               pyEkinlab / 1_GeV, pyPlab.getComponents() / 1_GeV);
+          Plab_final += pnew.getMomentum();
+          Elab_final += pnew.getEnergy();
+        }
+        // irreproducible in tests, LCOV_EXCL_START
+        catch (std::out_of_range const& ex) {
+          CORSIKA_LOGGER_CRITICAL(logger_, "Pythia ID {} unknown in C8", p8p.id());
+          throw ex;
+        }
+        // LCOV_EXCL_STOP
+      }
+    }
+
+    CORSIKA_LOGGER_DEBUG(logger_,
+                         "conservation (all GeV): "
+                         "Elab_final= {}"
+                         ", Plab_final= {}",
+                         Elab_final / 1_GeV, (Plab_final / 1_GeV).getComponents());
+
+    count_++;
+  }
+
+} // namespace corsika::pythia8
\ No newline at end of file

corsika/detail/modules/pythia8/Random.inl

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/modules/pythia8/Random.hpp>
+
+namespace corsika::pythia8 {
+
+  inline double Random::flat() { return Dist_(RNG_); }
+
+} // namespace corsika::pythia8

corsika/detail/modules/qgsjetII/InteractionModel.inl

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#include <corsika/modules/Random.hpp>
+#include <corsika/modules/qgsjetII/InteractionModel.hpp>
+#include <corsika/modules/qgsjetII/ParticleConversion.hpp>
+#include <corsika/modules/qgsjetII/QGSJetIIFragmentsStack.hpp>
+#include <corsika/modules/qgsjetII/QGSJetIIStack.hpp>
+
+#include <corsika/framework/geometry/FourVector.hpp>
+#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/framework/core/Logging.hpp>
+#include <corsika/framework/core/EnergyMomentumOperations.hpp>
+#include <corsika/framework/utility/COMBoost.hpp>
+
+#include <sstream>
+#include <tuple>
+
+#include <qgsjet-II-04.hpp>
+
+namespace corsika::qgsjetII {
+
+  inline InteractionModel::InteractionModel(boost::filesystem::path const dataPath) {
+    // initialize QgsjetII
+    corsika::connect_random_stream(rng_, ::qgsjetII::set_rng_function);
+
+    static bool initialized = false;
+    if (!initialized) {
+      CORSIKA_LOG_DEBUG("Reading QGSJetII data tables from {}", dataPath);
+      qgset_();
+      datadir DIR(dataPath.string() + "/");
+      qgaini_(DIR.data);
+      initialized = true;
+    }
+  }
+
+  inline InteractionModel::~InteractionModel() {
+    CORSIKA_LOG_DEBUG("QgsjetII::InteractionModel n= {}", count_);
+  }
+
+  inline bool InteractionModel::isValid(Code const projectileId, Code const targetId,
+                                        HEPEnergyType const sqrtS) const {
+
+    if (sqrtS < sqrtSmin_) { return false; }
+    if (is_nucleus(targetId)) {
+      size_t iTarget = get_nucleus_A(targetId);
+      if (iTarget > int(maxMassNumber_) || iTarget <= 0) { return false; }
+    } else if (targetId != Proton::code) {
+      return false;
+    }
+
+    if (is_nucleus(projectileId)) {
+      size_t iProjectile = get_nucleus_A(projectileId);
+      if (iProjectile > int(maxMassNumber_) || iProjectile <= 0) { return false; }
+    } else if (!is_hadron(projectileId)) {
+      return false;
+    }
+    return true;
+  }
+
+  inline CrossSectionType InteractionModel::getCrossSection(
+      Code const projectileId, Code const targetId, FourMomentum const& projectileP4,
+      FourMomentum const& targetP4) const {
+
+    if (!corsika::qgsjetII::canInteract(projectileId)) {
+      return CrossSectionType::zero();
+    }
+
+    auto const AfactorProjectile =
+        is_nucleus(projectileId) ? get_nucleus_A(projectileId) : 1;
+    auto const AfactorTarget = is_nucleus(targetId) ? get_nucleus_A(targetId) : 1;
+
+    // define projectile, in lab frame
+    auto const S = (projectileP4 + targetP4).getNormSqr();
+    auto const SNN =
+        (projectileP4 / AfactorProjectile + targetP4 / AfactorTarget).getNormSqr();
+    auto const sqrtSNN = sqrt(SNN);
+    if (!isValid(projectileId, targetId, sqrtSNN)) { return CrossSectionType::zero(); }
+
+    auto const projMass = get_mass(projectileId);
+    auto const targetMass = get_mass(targetId);
+
+    // lab-frame energy per projectile nucleon as required by qgsect()
+    HEPEnergyType const ElabN =
+        calculate_lab_energy(S, projMass, targetMass) / AfactorProjectile;
+
+    int const iBeam = static_cast<QgsjetIIXSClassIntType>(
+        corsika::qgsjetII::getQgsjetIIXSCode(projectileId));
+
+    CORSIKA_LOG_DEBUG(
+        "QgsjetII::getCrossSection Elab= {} GeV iBeam= {}"
+        " iProjectile= {} iTarget= {}",
+        ElabN / 1_GeV, iBeam, AfactorProjectile, AfactorTarget);
+    double const ElabNGeV{ElabN * (1 / 1_GeV)};
+    double const sigProd = qgsect_(ElabNGeV, iBeam, AfactorProjectile, AfactorTarget);
+    CORSIKA_LOG_DEBUG("QgsjetII::getCrossSection sigProd= {} mb", sigProd);
+    return sigProd * 1_mb;
+  }
+
+  inline std::tuple<CrossSectionType, CrossSectionType>
+  InteractionModel::getCrossSectionInelEla(Code projCode, Code targetCode,
+                                           FourMomentum const& proj4mom,
+                                           FourMomentum const& target4mom) const {
+    return {getCrossSection(projCode, targetCode, proj4mom, target4mom),
+            CrossSectionType::zero()};
+  }
+
+  template <typename TSecondaries>
+  inline void InteractionModel::doInteraction(TSecondaries& view, Code const projectileId,
+                                              Code const targetId,
+                                              FourMomentum const& projectileP4,
+                                              FourMomentum const& targetP4) {
+
+    CORSIKA_LOG_DEBUG(
+        "ProcessQgsjetII: "
+        "doInteraction: {} interaction possible? {}",
+        projectileId, corsika::qgsjetII::canInteract(projectileId));
+
+    // define projectile, in lab frame
+    auto const AfactorProjectile =
+        is_nucleus(projectileId) ? get_nucleus_A(projectileId) : 1;
+    auto const AfactorTarget = is_nucleus(targetId) ? get_nucleus_A(targetId) : 1;
+
+    auto const S = (projectileP4 + targetP4).getNormSqr();
+    auto const SNN =
+        (projectileP4 / AfactorProjectile + targetP4 / AfactorTarget).getNormSqr();
+    auto const sqrtSNN = sqrt(SNN);
+
+    if (!corsika::qgsjetII::canInteract(projectileId) ||
+        !isValid(projectileId, targetId, sqrtSNN)) {
+      throw std::runtime_error(fmt::format(
+          "invalid target [{}]/ projectile [{}] /energy [{} GeV] combination.",
+          get_name(targetId, full_name{}), get_name(projectileId, full_name{}),
+          sqrtSNN / 1_GeV));
+    }
+
+    auto const projMass = get_mass(projectileId);
+    auto const targetMass = get_mass(targetId);
+
+    // lab-frame energy per projectile nucleon
+    HEPEnergyType const Elab = calculate_lab_energy(S, projMass, targetMass);
+    auto const ElabN = Elab / AfactorProjectile;
+
+    CORSIKA_LOG_DEBUG("ebeam lab: {} GeV per projectile nucleon", ElabN / 1_GeV);
+
+    int const targetMassNumber = AfactorTarget;
+    CORSIKA_LOG_DEBUG("target: {}, qgsjetII code/A: {}", targetId, targetMassNumber);
+
+    // select QGSJetII internal projectile type
+    QgsjetIIHadronType qgsjet_hadron_type = qgsjetII::getQgsjetIIHadronType(projectileId);
+    if (qgsjet_hadron_type == QgsjetIIHadronType::NucleusType) {
+      qgsjet_hadron_type = bernoulli_(rng_) ? QgsjetIIHadronType::ProtonType
+                                            : QgsjetIIHadronType::NeutronType;
+    } else if (qgsjet_hadron_type == QgsjetIIHadronType::NeutralLightMesonType) {
+      // from conex: replace pi0 or rho0 with pi+/pi- in alternating sequence
+      qgsjet_hadron_type = alternate_;
+      alternate_ =
+          (alternate_ == QgsjetIIHadronType::PiPlusType ? QgsjetIIHadronType::PiMinusType
+                                                        : QgsjetIIHadronType::PiPlusType);
+    }
+
+    count_++;
+    int const qgsjet_hadron_type_int =
+        static_cast<QgsjetIICodeIntType>(qgsjet_hadron_type);
+    CORSIKA_LOG_DEBUG(
+        "qgsjet_hadron_type_int={} projectileMassNumber={} targetMassNumber={}",
+        qgsjet_hadron_type_int, AfactorProjectile, AfactorTarget);
+    qgini_(ElabN / 1_GeV, qgsjet_hadron_type_int, AfactorProjectile, AfactorTarget);
+    qgconf_();
+
+    CoordinateSystemPtr const& rootCS = get_root_CoordinateSystem();
+
+    // bookkeeping
+    MomentumVector Plab_final(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
+    HEPEnergyType Elab_final = 0_GeV;
+
+    // to read the secondaries
+    // define rotation to and from CoM frame
+    // CoM frame definition in QgsjetII projectile: +z
+
+    // QGSJetII, both, in input and output only considers the lab frame with a target at
+    // rest.
+
+    // system of initial-state
+    COMBoost boost(projectileP4 / AfactorProjectile, targetP4 / AfactorTarget);
+
+    auto const& originalCS = boost.getOriginalCS();
+    auto const& csPrime =
+        boost.getRotatedCS(); // z is along the CM motion (projectile, in Cascade)
+
+    HEPMomentumType const pLabMag = calculate_momentum(Elab, get_mass(projectileId));
+    MomentumVector const pLab{csPrime, {0_eV, 0_eV, pLabMag}};
+
+    // internal QGSJetII system: hadron-nucleon lab. frame!
+    COMBoost const boostInternal({Elab / AfactorProjectile, pLab / AfactorProjectile},
+                                 targetMass / AfactorTarget); // felix
+
+    // fragments
+    QGSJetIIFragmentsStack qfs;
+    for (auto& fragm : qfs) {
+      int const A = fragm.getFragmentSize();
+      if (A == 1) { // nucleon
+        Code const idFragm = bernoulli_(rng_) ? Code::Proton : Code::Neutron;
+
+        HEPMassType const nucleonMass = get_mass(idFragm);
+        // no pT, fragments just go forward
+        MomentumVector const momentum{
+            csPrime, {0_eV, 0_eV, calculate_momentum(ElabN, nucleonMass)}};
+
+        // this is not "CoM" here, but rather the system defined by projectile+target,
+        // which in Cascade-mode is already lab
+        auto const P4com = boostInternal.toCoM(FourVector{ElabN, momentum});
+        auto const P4output = boost.fromCoM(P4com);
+        auto p3output = P4output.getSpaceLikeComponents();
+        p3output.rebase(originalCS); // transform back into standard lab frame
+
+        HEPEnergyType const Ekin =
+            calculate_kinetic_energy(p3output.getNorm(), nucleonMass);
+
+        CORSIKA_LOG_DEBUG(
+            "secondary fragment> id= {}"
+            " p={}",
+            idFragm, p3output.getComponents());
+
+        auto pnew =
+            view.addSecondary(std::make_tuple(idFragm, Ekin, p3output.normalized()));
+        Plab_final += pnew.getMomentum();
+        Elab_final += pnew.getEnergy();
+
+      } else { // nucleus, A>1
+
+        int Z = 0;
+        switch (A) {
+          case 2: // deuterium
+            Z = 1;
+            break;
+          case 3: // tritium
+            Z = 1;
+            break;
+          case 4: // helium
+            Z = 2;
+            break;
+          default: // nucleus
+          {
+            Z = int(A / 2.15 + 0.7);
+          }
+        }
+
+        Code const idFragm = get_nucleus_code(A, Z);
+        HEPEnergyType const mass = get_mass(idFragm);
+        // no pT, frgments just go forward
+        MomentumVector momentum{csPrime,
+                                {0.0_GeV, 0.0_GeV, calculate_momentum(ElabN * A, mass)}};
+
+        // this is not "CoM" here, but rather the system defined by projectile+target,
+        // which in Cascade-mode is already lab
+        auto const P4com = boostInternal.toCoM(FourVector{ElabN * A, momentum});
+        auto const P4output = boost.fromCoM(P4com);
+        auto p3output = P4output.getSpaceLikeComponents();
+        p3output.rebase(originalCS); // transform back into standard lab frame
+
+        HEPEnergyType const Ekin = calculate_kinetic_energy(p3output.getNorm(), mass);
+
+        CORSIKA_LOG_DEBUG(
+            "secondary fragment> id={}"
+            " p={}"
+            " A={}"
+            " Z={}",
+            idFragm, p3output.getComponents(), A, Z);
+
+        auto pnew =
+            view.addSecondary(std::make_tuple(idFragm, Ekin, p3output.normalized()));
+        Plab_final += pnew.getMomentum();
+        Elab_final += pnew.getEnergy();
+      }
+    }
+
+    // secondaries
+    QGSJetIIStack qs;
+    for (auto& psec : qs) {
+      auto momentum = psec.getMomentum(csPrime);
+      // this is not "CoM" here, but rather the system defined by projectile+target,
+      // which in Cascade-mode is already lab
+      auto const P4com = boostInternal.toCoM(FourVector{psec.getEnergy(), momentum});
+      auto const P4output = boost.fromCoM(P4com);
+      auto p3output = P4output.getSpaceLikeComponents();
+      p3output.rebase(originalCS); // transform back into standard lab frame
+
+      Code const pid = corsika::qgsjetII::convertFromQgsjetII(psec.getPID());
+      HEPEnergyType const mass = get_mass(pid);
+      HEPEnergyType const Ekin = calculate_kinetic_energy(p3output.getNorm(), mass);
+
+      CORSIKA_LOG_DEBUG("secondary> id= {}, p= {}", pid, p3output.getComponents());
+      auto pnew = view.addSecondary(std::make_tuple(pid, Ekin, p3output.normalized()));
+      Plab_final += pnew.getMomentum();
+      Elab_final += pnew.getEnergy();
+    }
+    CORSIKA_LOG_DEBUG(
+        "conservation (all GeV): Ecm_final= n/a " /* << Ecm_final / 1_GeV*/
+        ", Elab_final={} "
+        ", Plab_final={}"
+        ", N_wounded,targ={}"
+        ", N_wounded,proj={}"
+        ", N_fragm,proj={}",
+        Elab_final / 1_GeV, (Plab_final / 1_GeV).getComponents(),
+        QGSJetIIFragmentsStackData::getWoundedNucleonsTarget(),
+        QGSJetIIFragmentsStackData::getWoundedNucleonsProjectile(), qfs.getSize());
+  }
+} // namespace corsika::qgsjetII

corsika/detail/modules/qgsjetII/ParticleConversion.inl

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <corsika/modules/qgsjetII/ParticleConversion.hpp>
+
+using namespace corsika::qgsjetII;

corsika/detail/modules/qgsjetII/QGSJetIIStack.inl

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+namespace corsika::qgsjetII {
+
+  inline void QGSJetIIStackData::clear() {
+    qgarr12_.nsp = 0;
+    qgarr13_.nsf = 0;
+    qgarr55_.nwt = 0;
+  }
+  inline unsigned int QGSJetIIStackData::getSize() const { return qgarr12_.nsp; }
+  inline unsigned int QGSJetIIStackData::getCapacity() const { return nptmax; }
+
+  inline void QGSJetIIStackData::setId(const unsigned int i, const int v) {
+    qgarr14_.ich[i] = v;
+  }
+  inline void QGSJetIIStackData::setEnergy(const unsigned int i, const HEPEnergyType v) {
+    qgarr14_.esp[i][0] = v / 1_GeV;
+  }
+
+  inline void QGSJetIIStackData::setMomentum(const unsigned int i,
+                                             const MomentumVector& v) {
+    auto tmp = v.getComponents();
+    qgarr14_.esp[i][2] = tmp[0] / 1_GeV;
+    qgarr14_.esp[i][3] = tmp[1] / 1_GeV;
+    qgarr14_.esp[i][1] = tmp[2] / 1_GeV;
+  }
+
+  inline int QGSJetIIStackData::getId(const unsigned int i) const {
+    return qgarr14_.ich[i];
+  }
+  inline HEPEnergyType QGSJetIIStackData::getEnergy(const int i) const {
+    return qgarr14_.esp[i][0] * 1_GeV;
+  }
+  inline MomentumVector QGSJetIIStackData::getMomentum(
+      const unsigned int i, const CoordinateSystemPtr& CS) const {
+    QuantityVector<hepmomentum_d> components = {qgarr14_.esp[i][2] * 1_GeV,
+                                                qgarr14_.esp[i][3] * 1_GeV,
+                                                qgarr14_.esp[i][1] * 1_GeV};
+    return MomentumVector(CS, components);
+  }
+
+  inline void QGSJetIIStackData::copy(const unsigned int i1, const unsigned int i2) {
+    qgarr14_.ich[i2] = qgarr14_.ich[i1];
+    for (unsigned int i = 0; i < 4; ++i) qgarr14_.esp[i2][i] = qgarr14_.esp[i1][i];
+  }
+
+  inline void QGSJetIIStackData::swap(const unsigned int i1, const unsigned int i2) {
+    std::swap(qgarr14_.ich[i1], qgarr14_.ich[i2]);
+    for (unsigned int i = 0; i < 4; ++i)
+      std::swap(qgarr14_.esp[i1][i], qgarr14_.esp[i2][i]);
+  }
+
+  inline void QGSJetIIStackData::incrementSize() { qgarr12_.nsp++; }
+  inline void QGSJetIIStackData::decrementSize() {
+    if (qgarr12_.nsp > 0) { qgarr12_.nsp--; }
+  }
+
+  template <typename StackIteratorInterface>
+  inline void ParticleInterface<StackIteratorInterface>::setParticleData(
+      const int vID, const HEPEnergyType vE, const MomentumVector& vP,
+      const HEPMassType) {
+    setPID(vID);
+    setEnergy(vE);
+    setMomentum(vP);
+  }
+
+  template <typename StackIteratorInterface>
+  inline void ParticleInterface<StackIteratorInterface>::setParticleData(
+      ParticleInterface<StackIteratorInterface>& /*parent*/, const int vID,
+      const HEPEnergyType vE, const MomentumVector& vP, const HEPMassType) {
+    setPID(vID);
+    setEnergy(vE);
+    setMomentum(vP);
+  }
+
+  template <typename StackIteratorInterface>
+  inline void ParticleInterface<StackIteratorInterface>::setEnergy(
+      const HEPEnergyType v) {
+    getStackData().setEnergy(getIndex(), v);
+  }
+
+  template <typename StackIteratorInterface>
+  inline HEPEnergyType ParticleInterface<StackIteratorInterface>::getEnergy() const {
+    return getStackData().getEnergy(getIndex());
+  }
+
+  template <typename StackIteratorInterface>
+  inline void ParticleInterface<StackIteratorInterface>::setPID(const int v) {
+    getStackData().setId(getIndex(), v);
+  }
+
+  template <typename StackIteratorInterface>
+  inline corsika::qgsjetII::QgsjetIICode
+  ParticleInterface<StackIteratorInterface>::getPID() const {
+    return static_cast<corsika::qgsjetII::QgsjetIICode>(getStackData().getId(getIndex()));
+  }
+
+  template <typename StackIteratorInterface>
+  inline MomentumVector ParticleInterface<StackIteratorInterface>::getMomentum(
+      const CoordinateSystemPtr& CS) const {
+    return getStackData().getMomentum(getIndex(), CS);
+  }
+
+  template <typename StackIteratorInterface>
+  inline void ParticleInterface<StackIteratorInterface>::setMomentum(
+      const MomentumVector& v) {
+    getStackData().setMomentum(getIndex(), v);
+  }
+
+} // namespace corsika::qgsjetII

corsika/detail/modules/qgsjetII/qgsjet-II-04.inl

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/modules/qgsjetII/qgsjet-II-04.hpp>
+
+#include <corsika/framework/random/RNGManager.hpp>
+#include <corsika/framework/core/Logging.hpp>
+
+#include <iostream>
+#include <random>
+
+inline datadir::datadir(const std::string& dir) {
+  if (dir.length() > 130) {
+    CORSIKA_LOG_ERROR("QGSJetII error, will cut datadir \"{}\" to 130 characters: ", {});
+  }
+  int i = 0;
+  for (i = 0; i < std::min(130, int(dir.length())); ++i) data[i] = dir[i];
+  data[i + 0] = ' ';
+  data[i + 1] = '\0';
+}
+
+inline void lzmaopenfile_(const char*, int) {}
+inline void lzmaclosefile_() {}
+inline void lzmafillarray_(const double&, const int&) {}

corsika/detail/modules/radio/CoREAS.inl

+/*
+ * (c) Copyright 2022 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/modules/radio/CoREAS.hpp>
+
+namespace corsika {
+
+  template <typename TRadioDetector, typename TPropagator>
+  template <typename Particle>
+  inline ProcessReturn CoREAS<TRadioDetector, TPropagator>::simulate(
+      Step<Particle> const& step) {
+
+    // get the global simulation time for that track.
+    auto const startTime_{
+        step.getTimePre()}; // time at the start point of the track. I
+                            // should use something similar to fCoreHitTime (?)
+    auto const endTime_{step.getTimePost()}; // time at end point of track.
+
+    if (startTime_ == endTime_) {
+      return ProcessReturn::Ok;
+    } else {
+
+      // get start and end position of the track
+      Point const startPoint_{step.getPositionPre()};
+      Point const endPoint_{step.getPositionPost()};
+      // get the coordinate system of the startpoint and hence the track
+      auto const cs_{startPoint_.getCoordinateSystem()};
+
+      auto const currDirection{(endPoint_ - startPoint_).normalized()};
+      // calculate the track length
+      auto const tracklength_{(endPoint_ - startPoint_).getNorm()};
+
+      // beta is velocity / speed of light. Start & end should be the same in endpoints!
+      auto const corrBetaValue{(endPoint_ - startPoint_).getNorm() /
+                               (constants::c * (endTime_ - startTime_))};
+      auto const beta_{currDirection * corrBetaValue};
+
+      // get particle charge
+      auto const charge_{get_charge(step.getParticlePre().getPID())};
+
+      // get thinning weight
+      auto const thinningWeight{step.getParticlePre().getWeight()};
+
+      // constants for electric field vector calculation
+      auto const constants_{charge_ * emConstant_ * thinningWeight};
+
+      // set threshold for application of ZHS-like approximation.
+      const double approxThreshold_{1.0e-3};
+
+      // loop over each observer in the observer collection (detector)
+      for (auto& observer : observers_.getObservers()) {
+
+        // get the SignalPathCollection (path1) from the start "endpoint" to the observer.
+        auto paths1{this->propagator_.propagate(step.getParticlePre(), startPoint_,
+                                                observer.getLocation())};
+
+        // get the SignalPathCollection (path2) from the end "endpoint" to the observer.
+        auto paths2{this->propagator_.propagate(step.getParticlePre(), endPoint_,
+                                                observer.getLocation())};
+
+        // LCOV_EXCL_START
+        // This should never happen unless someone implements a bad propagator
+        // good to catch this, hard to test without writing a custom and broken propagator
+        if (paths1.size() != paths2.size()) {
+          CORSIKA_LOG_CRITICAL(
+              "Signal Paths do not have the same size in CoREAS. Path starts: {} and "
+              "path ends: {}",
+              paths1.size(), paths2.size());
+        }
+        // LCOV_EXCL_STOP
+
+        // loop over both paths at once and directly compare 'start' and 'end'
+        // attributes
+        for (size_t i = 0; (i < paths1.size() && i < paths2.size()); i++) {
+
+          // calculate preDoppler factor
+          double preDoppler_{1.0 - paths1[i].refractive_index_source_ *
+                                       beta_.dot(paths1[i].emit_)};
+
+          // check if preDoppler has become zero in case of refractive index of unity
+          // because of numerical limitations here you might need std::fabs(preDoppler)
+          // in the if statement - same with post & mid
+          // LCOV_EXCL_START
+          if (preDoppler_ == 0) {
+            CORSIKA_LOG_WARN("preDoppler factor numerically zero in COREAS");
+            // redo calculation with higher precision
+            auto const& beta_components_{beta_.getComponents(cs_)};
+            auto const& emit_components_{paths1[i].emit_.getComponents(cs_)};
+            long double const indexL_{paths1[i].refractive_index_source_};
+            long double const betaX_{static_cast<double>(beta_components_.getX())};
+            long double const betaY_{static_cast<double>(beta_components_.getY())};
+            long double const betaZ_{static_cast<double>(beta_components_.getZ())};
+            long double const startX_{static_cast<double>(emit_components_.getX())};
+            long double const startY_{static_cast<double>(emit_components_.getY())};
+            long double const startZ_{static_cast<double>(emit_components_.getZ())};
+            long double const doppler =
+                1.0l - indexL_ * (betaX_ * startX_ + betaY_ * startY_ + betaZ_ * startZ_);
+            preDoppler_ = doppler;
+          }
+          // LCOV_EXCL_STOP
+
+          // calculate postDoppler factor
+          double postDoppler_{1.0 - paths2[i].refractive_index_source_ *
+                                        beta_.dot(paths2[i].emit_)};
+
+          // check if postDoppler has become zero in case of refractive index of unity
+          // because of numerical limitations
+          // LCOV_EXCL_START
+          if (postDoppler_ == 0) {
+            CORSIKA_LOG_WARN("postDoppler factor numerically zero in CoREAS");
+            // redo calculation with higher precision
+            auto const& beta_components_{beta_.getComponents(cs_)};
+            auto const& emit_components_{paths2[i].emit_.getComponents(cs_)};
+            long double const indexL_{paths2[i].refractive_index_source_};
+            long double const betaX_{static_cast<double>(beta_components_.getX())};
+            long double const betaY_{static_cast<double>(beta_components_.getY())};
+            long double const betaZ_{static_cast<double>(beta_components_.getZ())};
+            long double const endX_{static_cast<double>(emit_components_.getX())};
+            long double const endY_{static_cast<double>(emit_components_.getY())};
+            long double const endZ_{static_cast<double>(emit_components_.getZ())};
+            long double const doppler =
+                1.0l - indexL_ * (betaX_ * endX_ + betaY_ * endY_ + betaZ_ * endZ_);
+            postDoppler_ = doppler;
+          }
+          // LCOV_EXCL_STOP
+
+          // calculate receive time for startpoint (aka time delay)
+          auto startPointReceiveTime_{
+              startTime_ +
+              paths1[i].propagation_time_}; // TODO: time 0 is when the imaginary
+                                            // primary hits the ground
+
+          // calculate receive time for endpoint
+          auto endPointReceiveTime_{endTime_ + paths2[i].propagation_time_};
+
+          // get unit vector for startpoint at observer location
+          auto ReceiveVectorStart_{paths1[i].receive_};
+
+          // get unit vector for endpoint at observer location
+          auto ReceiveVectorEnd_{paths2[i].receive_};
+
+          // perform ZHS-like calculation close to Cherenkov angle and for refractive
+          // index at observer location greater than 1
+          if ((paths1[i].refractive_index_destination_ > 1) &&
+              ((std::fabs(preDoppler_) < approxThreshold_) ||
+               (std::fabs(postDoppler_) < approxThreshold_))) {
+            CORSIKA_LOG_DEBUG("Used ZHS-like approximation in CoREAS - radio");
+
+            // clear the existing paths for this particle and track, since we don't need
+            // them anymore
+            paths1.clear();
+            paths2.clear();
+
+            auto const halfVector_{(startPoint_ - endPoint_) * 0.5};
+            auto const midPoint_{endPoint_ + halfVector_};
+
+            // get global simulation time for the middle point of that track.
+            TimeType const midTime_{(startTime_ + endTime_) * 0.5};
+
+            // get the SignalPathCollection (path3) from the middle "endpoint" to the
+            // observer.
+            auto paths3{this->propagator_.propagate(step.getParticlePre(), midPoint_,
+                                                    observer.getLocation())};
+
+            // now loop over the paths for endpoint that we got above
+            for (auto const& path : paths3) {
+
+              auto const midPointReceiveTime_{midTime_ + path.propagation_time_};
+              double midDoppler_{1.0 -
+                                 path.refractive_index_source_ * beta_.dot(path.emit_)};
+
+              // check if midDoppler has become zero because of numerical limitations
+              // LCOV_EXCL_START
+              if (midDoppler_ == 0) {
+                CORSIKA_LOG_WARN("midDoppler factor numerically zero in COREAS");
+                // redo calculation with higher precision
+                auto const& beta_components_{beta_.getComponents(cs_)};
+                auto const& emit_components_{path.emit_.getComponents(cs_)};
+                long double const indexL_{path.refractive_index_source_};
+                long double const betaX_{static_cast<double>(beta_components_.getX())};
+                long double const betaY_{static_cast<double>(beta_components_.getY())};
+                long double const betaZ_{static_cast<double>(beta_components_.getZ())};
+                long double const midX_{static_cast<double>(emit_components_.getX())};
+                long double const midY_{static_cast<double>(emit_components_.getY())};
+                long double const midZ_{static_cast<double>(emit_components_.getZ())};
+                long double const doppler =
+                    1.0l - indexL_ * (betaX_ * midX_ + betaY_ * midY_ + betaZ_ * midZ_);
+                midDoppler_ = doppler;
+              }
+              // LCOV_EXCL_STOP
+
+              // change the values of the receive unit vectors of start and end
+              ReceiveVectorStart_ = path.receive_;
+              ReceiveVectorEnd_ = path.receive_;
+
+              // CoREAS calculation -> get ElectricFieldVector for "midPoint"
+              ElectricFieldVector EVmid_ = (path.emit_.cross(path.emit_.cross(beta_))) /
+                                           midDoppler_ / path.R_distance_ * constants_ *
+                                           observer.getSampleRate();
+
+              ElectricFieldVector EV1_{EVmid_};
+              ElectricFieldVector EV2_{EVmid_ * (-1.0)};
+
+              TimeType deltaT_{tracklength_ / (constants::c * corrBetaValue) *
+                               std::fabs(midDoppler_)}; // TODO: Caution with this!
+
+              if (startPointReceiveTime_ <
+                  endPointReceiveTime_) // EVstart_ arrives earlier
+              {
+                startPointReceiveTime_ = midPointReceiveTime_ - 0.5 * deltaT_;
+                endPointReceiveTime_ = midPointReceiveTime_ + 0.5 * deltaT_;
+              } else // EVend_ arrives earlier
+              {
+                startPointReceiveTime_ = midPointReceiveTime_ + 0.5 * deltaT_;
+                endPointReceiveTime_ = midPointReceiveTime_ - 0.5 * deltaT_;
+              }
+
+              TimeType const gridResolution_{1 / observer.getSampleRate()};
+              deltaT_ = endPointReceiveTime_ - startPointReceiveTime_;
+
+              // redistribute contributions over time scale defined by the observation
+              // time resolution
+              if (abs(deltaT_) < (gridResolution_)) {
+
+                EV1_ *= std::fabs((deltaT_ / gridResolution_));
+                EV2_ *= std::fabs((deltaT_ / gridResolution_));
+
+                // ToDO: be careful with times in C8!!! where is the zero (time). Is it
+                // close-by?
+                long const startBin = static_cast<long>(
+                    std::floor(startPointReceiveTime_ / gridResolution_ + 0.5l));
+                long const endBin = static_cast<long>(
+                    std::floor(endPointReceiveTime_ / gridResolution_ + 0.5l));
+                double const startBinFraction =
+                    (startPointReceiveTime_ / gridResolution_) -
+                    std::floor(startPointReceiveTime_ / gridResolution_);
+                double const endBinFraction =
+                    (endPointReceiveTime_ / gridResolution_) -
+                    std::floor(endPointReceiveTime_ / gridResolution_);
+
+                // only do timing modification if contributions would land in same bin
+                if (startBin == endBin) {
+
+                  // if startE arrives before endE
+                  if ((deltaT_) >= 0_s) {
+                    if ((startBinFraction >= 0.5) &&
+                        (endBinFraction >= 0.5)) // both points left of bin center
+                    {
+                      startPointReceiveTime_ -=
+                          gridResolution_; // shift EV1_ to previous gridpoint
+                    } else if ((startBinFraction < 0.5) &&
+                               (endBinFraction < 0.5)) // both points right of bin center
+                    {
+                      endPointReceiveTime_ +=
+                          gridResolution_; // shift EV2_ to next gridpoint
+                    } else                 // points on both sides of bin center
+                    {
+                      double const leftDist = 1.0 - startBinFraction;
+                      double const rightDist = endBinFraction;
+                      // check if asymmetry to right or left
+                      if (rightDist >= leftDist) {
+                        endPointReceiveTime_ +=
+                            gridResolution_; // shift EV2_ to next gridpoint
+                      } else {
+                        startPointReceiveTime_ -=
+                            gridResolution_; // shift EV1_ to previous gridpoint
+                      }
+                    }
+                  } else // if endE arrives before startE
+                  {
+                    if ((startBinFraction >= 0.5) &&
+                        (endBinFraction >= 0.5)) // both points left of bin center
+                    {
+                      endPointReceiveTime_ -=
+                          gridResolution_; // shift EV2_ to previous gridpoint
+                    } else if ((startBinFraction < 0.5) &&
+                               (endBinFraction < 0.5)) // both points right of bin center
+                    {
+                      startPointReceiveTime_ +=
+                          gridResolution_; // shift EV1_ to next gridpoint
+                    } else                 // points on both sides of bin center
+                    {
+                      double const leftDist = 1.0 - endBinFraction;
+                      double const rightDist = startBinFraction;
+                      // check if asymmetry to right or left
+                      if (rightDist >= leftDist) {
+                        startPointReceiveTime_ +=
+                            gridResolution_; // shift EV1_ to next gridpoint
+                      } else {
+                        endPointReceiveTime_ -=
+                            gridResolution_; // shift EV2_ to previous gridpoint
+                      }
+                    }
+                  } // End of else statement
+                }   // End of if for startbin == endbin
+              }     // End of if deltaT < gridresolution
+
+              // TODO: Be very careful with this. Maybe the EVs should be fed after the
+              // for loop of paths3
+              observer.receive(startPointReceiveTime_, ReceiveVectorStart_, EV1_);
+              observer.receive(endPointReceiveTime_, ReceiveVectorEnd_, EV2_);
+            } // End of looping over paths3
+
+          } // end of ZHS-like approximation
+          else {
+
+            // calculate electric field vector for startpoint
+            ElectricFieldVector EV1_ =
+                (paths1[i].emit_.cross(paths1[i].emit_.cross(beta_))) / preDoppler_ /
+                paths1[i].R_distance_ * constants_ * observer.getSampleRate();
+
+            // calculate electric field vector for endpoint
+            ElectricFieldVector EV2_ =
+                (paths2[i].emit_.cross(paths2[i].emit_.cross(beta_))) / postDoppler_ /
+                paths2[i].R_distance_ * constants_ * (-1.0) * observer.getSampleRate();
+
+            if ((preDoppler_ < 1.e-9) || (postDoppler_ < 1.e-9)) {
+
+              TimeType const gridResolution_{1 / observer.getSampleRate()};
+              TimeType deltaT_{endPointReceiveTime_ - startPointReceiveTime_};
+
+              if (abs(deltaT_) < (gridResolution_)) {
+
+                EV1_ *= std::fabs(deltaT_ / gridResolution_); // Todo: rename EV1 and 2
+                EV2_ *= std::fabs(deltaT_ / gridResolution_);
+
+                long const startBin = static_cast<long>(
+                    std::floor(startPointReceiveTime_ / gridResolution_ + 0.5l));
+                long const endBin = static_cast<long>(
+                    std::floor(endPointReceiveTime_ / gridResolution_ + 0.5l));
+                double const startBinFraction =
+                    (startPointReceiveTime_ / gridResolution_) -
+                    std::floor(startPointReceiveTime_ / gridResolution_);
+                double const endBinFraction =
+                    (endPointReceiveTime_ / gridResolution_) -
+                    std::floor(endPointReceiveTime_ / gridResolution_);
+
+                // only do timing modification if contributions would land in same bin
+                if (startBin == endBin) {
+
+                  if ((startBinFraction >= 0.5) &&
+                      (endBinFraction >= 0.5)) // both points left of bin center
+                  {
+                    startPointReceiveTime_ -=
+                        gridResolution_; // shift EV1_ to previous gridpoint
+                  } else if ((startBinFraction < 0.5) &&
+                             (endBinFraction < 0.5)) // both points right of bin center
+                  {
+                    endPointReceiveTime_ +=
+                        gridResolution_; // shift EV2_ to next gridpoint
+                  } else                 // points on both sides of bin center
+                  {
+                    double const leftDist = 1.0 - startBinFraction;
+                    double const rightDist = endBinFraction;
+                    // check if asymmetry to right or left
+                    if (rightDist >= leftDist) {
+                      endPointReceiveTime_ +=
+                          gridResolution_; // shift EV2_ to next gridpoint
+                    } else {
+                      startPointReceiveTime_ -=
+                          gridResolution_; // shift EV1_ to previous gridpoint
+                    }
+                  }
+
+                } // End of if for startbin == endbin
+              }   // End of if deltaT < gridresolution
+            }     // End of if that checks small doppler factors
+            observer.receive(startPointReceiveTime_, ReceiveVectorStart_, EV1_);
+            observer.receive(endPointReceiveTime_, ReceiveVectorEnd_, EV2_);
+          } // End of else that does not perform ZHS-like approximation
+
+        } // End of loop over both paths to get signal info
+      }   // End of looping over observer
+
+      return ProcessReturn::Ok;
+    }
+  } // End of simulate method
+
+} // namespace corsika

corsika/detail/modules/radio/LimitedRadioProcess.inl

+/*
+ * (c) Copyright 2024 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+#pragma once
+
+namespace corsika {
+
+  template <class TUnderlyingProcess>
+  inline LimitedRadioProcess<TUnderlyingProcess>::LimitedRadioProcess(
+      TUnderlyingProcess&& process,
+      std::function<LimitedRadioProcess::functor_signature> runThis)
+      : process_{std::forward<TUnderlyingProcess>(process)}
+      , runThis_{std::move(runThis)} {}
+
+  template <class TUnderlyingProcess>
+  template <typename Particle>
+  inline ProcessReturn LimitedRadioProcess<TUnderlyingProcess>::doContinuous(
+      const Step<Particle>& step, const bool arg) {
+
+    if ((step.getParticlePre().getPID() != Code::Electron) &&
+        (step.getParticlePre().getPID() != Code::Positron)) {
+      CORSIKA_LOG_DEBUG("Not e+/-, will not run RadioProcess");
+      return ProcessReturn::Ok;
+    }
+
+    // Check for running the wrapped RadioProcess function
+    if (runThis_(step.getPositionPre())) {
+      CORSIKA_LOG_TRACE("Will run RadioProcess for particle at {}",
+                        step.getPositionPre());
+      return process_.doContinuous(step, arg);
+    }
+    CORSIKA_LOG_DEBUG("Particle at {} is not selected by LimitedRadioProcess",
+                      step.getPositionPre());
+    return ProcessReturn::Ok;
+  }
+
+  template <class TUnderlyingProcess>
+  template <typename Particle, typename Track>
+  inline LengthType LimitedRadioProcess<TUnderlyingProcess>::getMaxStepLength(
+      Particle const& vParticle, Track const& vTrack) const {
+    return process_.getMaxStepLength(vParticle, vTrack);
+  }
+
+} // namespace corsika
\ No newline at end of file

corsika/detail/modules/radio/RadioProcess.inl

+/*
+ * (c) Copyright 2022 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+#pragma once
+
+#include <corsika/modules/radio/RadioProcess.hpp>
+
+namespace corsika {
+
+  template <typename TObserverCollection, typename TRadioImpl, typename TPropagator>
+  inline TRadioImpl&
+  RadioProcess<TObserverCollection, TRadioImpl, TPropagator>::implementation() {
+    return static_cast<TRadioImpl&>(*this);
+  }
+
+  template <typename TObserverCollection, typename TRadioImpl, typename TPropagator>
+  inline TRadioImpl const&
+  RadioProcess<TObserverCollection, TRadioImpl, TPropagator>::implementation() const {
+    return static_cast<TRadioImpl const&>(*this);
+  }
+
+  template <typename TObserverCollection, typename TRadioImpl, typename TPropagator>
+  inline RadioProcess<TObserverCollection, TRadioImpl, TPropagator>::RadioProcess(
+      TObserverCollection& observers, TPropagator& propagator)
+      : observers_(observers)
+      , propagator_(propagator) {}
+
+  template <typename TObserverCollection, typename TRadioImpl, typename TPropagator>
+  template <typename Particle>
+  inline ProcessReturn RadioProcess<TObserverCollection, TRadioImpl,
+                                    TPropagator>::doContinuous(const Step<Particle>& step,
+                                                               const bool) {
+
+    // return immediately if radio process does not have any observers
+    if (observers_.size() == 0) return ProcessReturn::Ok;
+
+    // we want the following particles:
+    // Code::Electron & Code::Positron
+
+    // we wrap Simulate() in doContinuous as the plan is to add particle level
+    // filtering or thinning for calculation of the radio emission. This is
+    // important for controlling the runtime of radio (by ignoring particles
+    // that aren't going to contribute i.e. heavy hadrons)
+    // if (valid(step)) {
+    auto const particleID_{step.getParticlePre().getPID()};
+    if ((particleID_ == Code::Electron) || (particleID_ == Code::Positron)) {
+      CORSIKA_LOG_DEBUG("Particle for radio calculation: {} ", particleID_);
+      return this->implementation().simulate(step);
+    } else {
+      CORSIKA_LOG_DEBUG("Particle {} is irrelevant for radio", particleID_);
+      return ProcessReturn::Ok;
+    }
+    //}
+  }
+
+  template <typename TObserverCollection, typename TRadioImpl, typename TPropagator>
+  template <typename Particle, typename Track>
+  inline LengthType
+  RadioProcess<TObserverCollection, TRadioImpl, TPropagator>::getMaxStepLength(
+      [[maybe_unused]] const Particle& vParticle,
+      [[maybe_unused]] const Track& vTrack) const {
+    return meter * std::numeric_limits<double>::infinity();
+  }
+
+  template <typename TObserverCollection, typename TRadioImpl, typename TPropagator>
+  inline void RadioProcess<TObserverCollection, TRadioImpl, TPropagator>::startOfLibrary(
+      const boost::filesystem::path& directory) {
+
+    // setup the streamer
+    output_.initStreamer((directory / ("observers.parquet")).string());
+
+    // enable compression with the default level
+    output_.enableCompression();
+
+    // LCOV_EXCL_START
+    // build the schema
+    output_.addField("Time", parquet::Repetition::REQUIRED, parquet::Type::DOUBLE,
+                     parquet::ConvertedType::NONE);
+
+    output_.addField("Ex", parquet::Repetition::REQUIRED, parquet::Type::DOUBLE,
+                     parquet::ConvertedType::NONE);
+
+    output_.addField("Ey", parquet::Repetition::REQUIRED, parquet::Type::DOUBLE,
+                     parquet::ConvertedType::NONE);
+
+    output_.addField("Ez", parquet::Repetition::REQUIRED, parquet::Type::DOUBLE,
+                     parquet::ConvertedType::NONE);
+    // LCOV_EXCL_STOP
+    // and build the streamer
+    output_.buildStreamer();
+  }
+
+  template <typename TObserverCollection, typename TRadioImpl, typename TPropagator>
+  inline void RadioProcess<TObserverCollection, TRadioImpl, TPropagator>::endOfShower(
+      const unsigned int) {
+
+    // loop over every observer and instruct them to
+    // flush data to disk, and then reset the observer
+    // before the next event
+    for (auto& observer : observers_.getObservers()) {
+
+      auto const sampleRate = observer.getSampleRate() * 1_s;
+      auto const radioImplementation =
+          static_cast<std::string>(this->implementation().algorithm);
+
+      // get the axis labels for this observer and write the first row.
+      axistype axis = observer.implementation().getAxis();
+
+      // get the copy of the waveform data for this event
+      std::vector<double> const& dataX = observer.implementation().getWaveformX();
+      std::vector<double> const& dataY = observer.implementation().getWaveformY();
+      std::vector<double> const& dataZ = observer.implementation().getWaveformZ();
+
+      // check for the axis name
+      std::string label = "Unknown";
+      if (observer.getDomainLabel() == "Time") {
+        label = "Time";
+      }
+      // LCOV_EXCL_START
+      else if (observer.getDomainLabel() == "Frequency") {
+        label = "Frequency";
+      }
+      // LCOV_EXCL_STOP
+      if (radioImplementation == "ZHS" && label == "Time") {
+        for (size_t i = 0; i < axis.size() - 1; i++) {
+          auto time = (axis.at(i + 1) + axis.at(i)) / 2.;
+          auto Ex = -(dataX.at(i + 1) - dataX.at(i)) * sampleRate;
+          auto Ey = -(dataY.at(i + 1) - dataY.at(i)) * sampleRate;
+          auto Ez = -(dataZ.at(i + 1) - dataZ.at(i)) * sampleRate;
+
+          *(output_.getWriter())
+              << showerId_ << static_cast<double>(time) << static_cast<double>(Ex)
+              << static_cast<double>(Ey) << static_cast<double>(Ez) << parquet::EndRow;
+        }
+      } else if (radioImplementation == "CoREAS" && label == "Time") {
+        for (size_t i = 0; i < axis.size() - 1; i++) {
+          *(output_.getWriter())
+              << showerId_ << static_cast<double>(axis[i])
+              << static_cast<double>(dataX[i]) << static_cast<double>(dataY[i])
+              << static_cast<double>(dataZ[i]) << parquet::EndRow;
+        }
+      }
+
+      observer.reset();
+    }
+    output_.closeStreamer();
+
+    // increment our event counter
+    showerId_++;
+  }
+
+  template <typename TObserverCollection, typename TRadioImpl, typename TPropagator>
+  inline YAML::Node
+  RadioProcess<TObserverCollection, TRadioImpl, TPropagator>::getConfig() const {
+
+    // top-level YAML node
+    YAML::Node config;
+
+    // fill in some basics
+    config["type"] = "RadioProcess";
+    config["algorithm"] = this->implementation().algorithm;
+    config["units"]["time"] = "ns";
+    config["units"]["frequency"] = "GHz";
+    config["units"]["electric field"] = "V/m";
+    config["units"]["distance"] = "m";
+
+    for (auto& observer : observers_.getObservers()) {
+      // get the name/location of this observer
+      auto name = observer.getName();
+      auto location = observer.getLocation().getCoordinates();
+
+      // get the observers config
+      config["observers"][name] = observer.getConfig();
+
+      // write the location of this observer
+      config["observers"][name]["location"].push_back(location.getX() / 1_m);
+      config["observers"][name]["location"].push_back(location.getY() / 1_m);
+      config["observers"][name]["location"].push_back(location.getZ() / 1_m);
+    }
+
+    return config;
+  }
+
+} // namespace corsika