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 GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
+ * 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/Interaction.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>
 
-#include <corsika/setup/SetupEnvironment.hpp>
-#include <corsika/setup/SetupStack.hpp>
-#include <corsika/setup/SetupTrajectory.hpp>
-
-#include <iostream>
-
 namespace corsika::proposal {
 
-  void ContinuousProcess::buildCalculator(Code code, NuclearComposition const& comp) {
+  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
+    // take some minutes if you have to build the tables and cannot read the tables
     // from disk
-    auto c = p_cross->second(media.at(comp.getHash()), emCut_);
+    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 disp = PROPOSAL::make_displacement(c, true);
-    auto scatter =
-        PROPOSAL::make_scattering("highland", particle[code], media.at(comp.getHash()));
-    calc[std::make_pair(comp.getHash(), code)] =
-        std::make_tuple(std::move(disp), std::move(scatter));
+    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 <>
-  ContinuousProcess::ContinuousProcess(setup::Environment const& _env,
-                                       HEPEnergyType _emCut)
-      : ProposalProcessBase(_env, _emCut) {}
+  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 <>
-  void ContinuousProcess::scatter(setup::Stack::particle_type& vP,
-                                  HEPEnergyType const& loss,
-                                  GrammageType const& grammage) {
+  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(vP, calc);
+    auto c = getCalculator(step.getParticlePre(), calc);
 
-    // Cast corsika vector to proposal vector
-    auto vP_dir = vP.getDirection();
-    auto d = vP_dir.getComponents();
-    auto direction = PROPOSAL::Vector3D(d.getX().magnitude(), d.getY().magnitude(),
-                                        d.getZ().magnitude());
+    auto initial_particle_dir = step.getDirectionPre();
 
-    auto E_f = vP.getEnergy() - loss;
+    auto E_i_total = step.getEkinPre() + step.getParticlePre().getMass();
+    auto E_f_total = E_i_total - loss;
 
-    // draw random numbers required for scattering process
+    // 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 rnd = array<double, 4>();
-    for (auto& it : rnd) it = distr(RNG_);
-
-    // calculate deflection based on particle energy, loss
-    auto [mean_dir, final_dir] = get<eSCATTERING>(c->second)->Scatter(
-        grammage / 1_g * square(1_cm), vP.getEnergy() / 1_MeV, E_f / 1_MeV, direction,
-        rnd);
-
-    // TODO: neglect mean direction deflection because Trajectory is a const ref
-    (void)mean_dir;
+    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
-    auto vec = QuantityVector(final_dir.GetX() * E_f, final_dir.GetY() * E_f,
-                              final_dir.GetZ() * E_f);
-    vP.setMomentum(MomentumVector(vP_dir.getCoordinateSystem(), vec));
+    DirectionVector diff_dir_ = scattered_particle_dir - initial_particle_dir;
+    step.add_dU(diff_dir_);
   }
 
-  template <>
-  ProcessReturn ContinuousProcess::doContinuous(setup::Stack::particle_type& vP,
-                                                setup::Trajectory const& vT) {
-
-    if (!canInteract(vP.getPID())) return ProcessReturn::Ok;
-    if (vT.getLength() == 0_m) return ProcessReturn::Ok;
-
+  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 =
-        vP.getNode()->getModelProperties().getIntegratedGrammage(vT, vT.getLength());
+    auto dX = step.getParticlePre().getNode()->getModelProperties().getIntegratedGrammage(
+        step.getStraightTrack());
 
     // get or build corresponding track integral calculator and solve the
     // integral
-    auto c = getCalculator(vP, calc);
-    auto final_energy = get<eDISPLACEMENT>(c->second)->UpperLimitTrackIntegral(
-                            vP.getEnergy() / 1_MeV, dX / 1_g * 1_cm * 1_cm) *
-                        1_MeV;
-    auto dE = vP.getEnergy() - final_energy;
-    energy_lost_ += dE;
+    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 (vP.getChargeNumber() != 0) scatter(vP, dE, dX);
-    vP.setEnergy(final_energy);
-    vP.setMomentum(vP.getMomentum() * vP.getEnergy() / vP.getMomentum().getNorm());
-    return ProcessReturn::Ok;
-  }
+    if (step.getParticlePre().getChargeNumber() != 0) scatter(step, dE, dX);
+    step.add_dEkin(-dE); // on the stack, this is just kinetic energy, E-m
 
-  template <>
-  LengthType ContinuousProcess::getMaxStepLength(setup::Stack::particle_type const& vP,
-                                                 setup::Trajectory const& vT) {
+    // also send to output
+    TOutput::write(step.getPositionPre(), step.getPositionPost(),
+                   step.getParticlePre().getPID(),
+                   step.getParticlePre().getWeight() * dE);
 
-    if (!canInteract(vP.getPID())) return meter * std::numeric_limits<double>::infinity();
+    return ProcessReturn::Ok;
+  }
 
-    // Limit the step size of a conitnuous loss. The maximal continuous loss seems to be a
-    // hyper parameter which must be adjusted.
+  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.
     //
-    // in any case: never go below 0.99*emCut_ This needs to be
-    // slightly smaller than emCut_ since, either this Step is limited
-    // by energy_lim, then the particle is stopped in a very short
-    // range (before doing anythin else) and is then removed
-    // instantly. The exact position where it reaches emCut is not
-    // important, the important fact is that its E_kin is zero
-    // afterwards.
-    //
-    auto energy_lim = std::max(0.9 * vP.getEnergy(), 0.99 * emCut_);
+    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 = get<eDISPLACEMENT>(c->second)->SolveTrackIntegral(
-                        vP.getEnergy() / 1_MeV, energy_lim / 1_MeV) *
-                    1_g / square(1_cm);
+    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(vT, grammage);
+    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;
   }
 
-  void ContinuousProcess::showResults() const {
-    std::cout << " ******************************" << std::endl
-              << " PROCESS::ContinuousProcess: " << std::endl;
-    std::cout << " energy lost dE (GeV)      :  " << energy_lost_ / 1_GeV << std::endl;
+  template <typename TOutput>
+  inline YAML::Node ContinuousProcess<TOutput>::getConfig() const {
+    return YAML::Node();
   }
 
-  void ContinuousProcess::reset() { energy_lost_ = 0_GeV; }
-
 } // 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/Interaction.inl

-/*
- * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#include <corsika/media/IMediumModel.hpp>
-#include <corsika/media/NuclearComposition.hpp>
-#include <corsika/modules/proposal/Interaction.hpp>
-#include <corsika/framework/utility/COMBoost.hpp>
-#include <corsika/framework/core/PhysicalUnits.hpp>
-
-#include <corsika/setup/SetupEnvironment.hpp>
-#include <corsika/setup/SetupStack.hpp>
-#include <corsika/setup/SetupTrajectory.hpp>
-
-#include <limits>
-#include <memory>
-#include <random>
-#include <tuple>
-
-namespace corsika::proposal {
-
-  template <>
-  Interaction::Interaction(setup::Environment const& _env, HEPEnergyType _emCut)
-      : ProposalProcessBase(_env, _emCut) {}
-
-  void Interaction::buildCalculator(Code code, NuclearComposition const& comp) {
-    // 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
-    // from disk
-    auto c = p_cross->second(media.at(comp.getHash()), emCut_);
-
-    // Look which interactions take place and build the corresponding
-    // interaction and secondarie 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(comp.getHash(), code)] = std::make_tuple(
-        PROPOSAL::make_secondaries(inter_types, particle[code], media.at(comp.getHash())),
-        PROPOSAL::make_interaction(c, true));
-  }
-
-  template <>
-  ProcessReturn Interaction::doInteraction(setup::StackView& view) {
-
-    auto const projectile = view.getProjectile();
-
-    if (canInteract(projectile.getPID())) {
-
-      // 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 = get<eINTERACTION>(c->second)->Rates(projectile.getEnergy() / 1_MeV);
-      auto [type, comp_ptr, v] = get<eINTERACTION>(c->second)->SampleLoss(
-          projectile.getEnergy() / 1_MeV, rates, distr(RNG_));
-
-      // Read how much random numbers are required to calculate the secondaries.
-      // Calculate the secondaries and deploy them on the corsika stack.
-      auto rnd =
-          vector<double>(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::Vector3D(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::Vector3D(d.getX().magnitude(), d.getY().magnitude(),
-                                          d.getZ().magnitude());
-      auto loss = make_tuple(static_cast<int>(type), point, direction,
-                             v * projectile.getEnergy() / 1_MeV, 0.);
-      auto sec = get<eSECONDARIES>(c->second)->CalculateSecondaries(
-          projectile.getEnergy() / 1_MeV, loss, *comp_ptr, rnd);
-      for (auto& s : sec) {
-        auto E = get<PROPOSAL::Loss::ENERGY>(s) * 1_MeV;
-        auto vecProposal = get<PROPOSAL::Loss::DIRECTION>(s);
-        auto vec = QuantityVector(vecProposal.GetX() * E, vecProposal.GetY() * E,
-                                  vecProposal.GetZ() * E);
-        auto p = MomentumVector(labCS, vec);
-        auto sec_code =
-            convert_from_PDG(static_cast<PDGCode>(get<PROPOSAL::Loss::TYPE>(s)));
-        view.addSecondary(
-            make_tuple(sec_code, E, p, projectile.getPosition(), projectile.getTime()));
-      }
-    }
-    return ProcessReturn::Ok;
-  }
-
-  template <>
-  GrammageType Interaction::getInteractionLength(
-      setup::Stack::particle_type const& projectile) {
-
-    if (canInteract(projectile.getPID())) {
-      auto c = getCalculator(projectile, calc);
-      return get<eINTERACTION>(c->second)->MeanFreePath(projectile.getEnergy() / 1_MeV) *
-             1_g / (1_cm * 1_cm);
-    }
-    return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
-  }
-} // 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 GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
+ * 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>
@@ -11,10 +10,7 @@
 #include <corsika/modules/proposal/ProposalProcessBase.hpp>
 #include <corsika/framework/core/PhysicalUnits.hpp>
 #include <corsika/framework/utility/COMBoost.hpp>
-
-#include <corsika/setup/SetupEnvironment.hpp>
-#include <corsika/setup/SetupStack.hpp>
-#include <corsika/setup/SetupTrajectory.hpp>
+#include <corsika/framework/utility/CorsikaData.hpp>
 
 #include <cstdlib>
 #include <iostream>
@@ -25,25 +21,42 @@
 
 namespace corsika::proposal {
 
-  bool ProposalProcessBase::canInteract(Code pcode) const {
-    if (std::find(begin(tracked), end(tracked), pcode) != end(tracked)) return true;
-    return false;
+  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;
   }
 
-  ProposalProcessBase::ProposalProcessBase(setup::Environment const& _env,
-                                           HEPEnergyType _emCut)
-      : emCut_(_emCut)
-      , RNG_(RNGManager::getInstance().getRandomStream("proposal")) {
+  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(Point(_env.getCoordinateSystem(), 0_cm, 0_cm, 0_cm)));
+        const auto& medium = mediumData(prop.getMedium());
 
-        auto comp_vec = std::vector<PROPOSAL::Components::Component>();
+        auto comp_vec = std::vector<PROPOSAL::Component>();
         const auto& comp = prop.getNuclearComposition();
         auto frac_iter = comp.getFractions().cbegin();
-        for (auto& pcode : comp.getComponents()) {
+        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;
@@ -56,24 +69,42 @@ namespace corsika::proposal {
       }
     });
 
-    PROPOSAL::InterpolationDef::order_of_interpolation = 2;
-    PROPOSAL::InterpolationDef::nodes_cross_section = 100;
-    PROPOSAL::InterpolationDef::nodes_propagate = 1000;
-
     //! 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.
-    if (auto data_path = std::getenv("CORSIKA_DATA")) {
-      PROPOSAL::InterpolationDef::path_to_tables = std::string(data_path) + "/PROPOSAL";
-    } else {
-      throw std::runtime_error(
-          "It is not recommended to run PROPOSAL without its tables in "
-          "$CORSIKA_DATA/PROPOSAL. This would be extremely slow. Please provide the "
-          "table directory. ");
+    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);
   }
 
-  size_t ProposalProcessBase::hash::operator()(const calc_key_t& p) const noexcept {
+  inline size_t ProposalProcessBase::hash::operator()(const calc_key_t& p) const
+      noexcept {
     return p.first ^ std::hash<Code>{}(p.second);
   }
 

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 GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
+ * 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>
@@ -12,47 +11,35 @@
 
 #include <corsika/framework/utility/COMBoost.hpp>
 
-#include <corsika/setup/SetupStack.hpp>
-#include <corsika/setup/SetupTrajectory.hpp>
-
 namespace corsika::pythia8 {
 
-  Decay::Decay(const bool print_listing)
-      : print_listing_(print_listing) {
+  inline Decay::Decay(bool const print_listing)
+      : Pythia8::Pythia(CORSIKA_Pythia8_XML_DIR)
+      , print_listing_(print_listing) {
     init();
   }
 
-  Decay::Decay(std::set<Code> const& those)
+  inline Decay::Decay(std::set<Code> const& those)
       : handleAllDecays_(false)
       , handledDecays_(those) {
     init();
   }
 
-  Decay::~Decay() { CORSIKA_LOG_INFO("Pythia::Decay n={}", count_); }
+  inline Decay::~Decay() { CORSIKA_LOG_INFO("Pythia::Decay n={}", count_); }
 
-  void Decay::init() {
+  inline void Decay::init() {
 
     // run this only once during construction
 
-    // set random number generator in pythia
-    Pythia8::RndmEngine* rndm = new corsika::pythia8::Random();
-    pythia_.setRndmEnginePtr(rndm);
-
-    /*
-       issue xyz: definition of particles and decay channels use the same mechanism in
-       corsika and pythia we should force pythia to use the file in corsika.
-     */
-    // bool ParticleData::reInit(string startFile, bool xmlFormat = true)
-    // read in particle data from Corsika 8
-    // pythia_.particleData.reInit("/home/felix/ngcorsika/corsika-build/include/corsika/particles/ParticleData.xml");
-    // pythia_.particleData.checkTable();
+    // link random number generator in pythia to CORSIKA8
+    Pythia8::Pythia::setRndmEnginePtr(std::make_shared<corsika::pythia8::Random>());
 
-    pythia_.readString("Next:numberShowInfo = 0");
-    pythia_.readString("Next:numberShowProcess = 0");
-    pythia_.readString("Next:numberShowEvent = 0");
+    Pythia8::Pythia::readString("Next:numberShowInfo = 0");
+    Pythia8::Pythia::readString("Next:numberShowProcess = 0");
+    Pythia8::Pythia::readString("Next:numberShowEvent = 0");
 
-    pythia_.readString("Print:quiet = on");
-    pythia_.readString("Check:particleData = 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
@@ -60,21 +47,23 @@ namespace corsika::pythia8 {
        the consistency of particle masses between event generators is an unsolved issues
     */
     CORSIKA_LOG_INFO("Pythia::Init: switching off event checking in pythia..");
-    pythia_.readString("Check:event = 1");
+    Pythia8::Pythia::readString("Check:event = 1");
 
-    pythia_.readString("ProcessLevel:all = off");
-    pythia_.readString("ProcessLevel:resonanceDecays = off");
+    Pythia8::Pythia::readString("ProcessLevel:all = off");
+    Pythia8::Pythia::readString("ProcessLevel:resonanceDecays = on");
 
     // making sure
     setStable(Code::Pi0);
 
-    //    pythia_.particleData.readString("59:m0 = 101.00");
+    //    Pythia8::Pythia::particleData.readString("59:m0 = 101.00");
 
-    if (!pythia_.init())
+    // LCOV_EXCL_START, we don't validate pythia8 internals
+    if (!Pythia8::Pythia::init())
       throw std::runtime_error("Pythia::Decay: Initialization failed!");
+    // LCOV_EXCL_STOP
   }
 
-  bool Decay::canHandleDecay(Code const vParticleCode) {
+  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 ||
@@ -82,91 +71,92 @@ namespace corsika::pythia8 {
         vParticleCode == Code::NuMuBar || vParticleCode == Code::NuTauBar ||
         vParticleCode == Code::Electron || vParticleCode == Code::Positron)
       return false;
-    else if (canDecay(vParticleCode)) // non-zero for particles known to sibyll
+    else if (canDecay(vParticleCode)) // check pythia8 internal
       return true;
     else
       return false;
   }
 
-  void Decay::setHandleDecay(Code const vParticleCode) {
+  inline void Decay::setHandleDecay(Code const vParticleCode) {
     handleAllDecays_ = false;
-    CORSIKA_LOG_INFO("Pythia::Decay: set to handle decay of {} ", vParticleCode);
+    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!");
   }
 
-  void Decay::setHandleDecay(std::vector<Code> const& vParticleList) {
+  inline void Decay::setHandleDecay(std::vector<Code> const& vParticleList) {
     handleAllDecays_ = false;
-    for (auto p : vParticleList) setHandleDecay(p);
+    for (auto const p : vParticleList) setHandleDecay(p);
   }
 
-  bool Decay::isDecayHandled(Code const vParticleCode) {
+  inline bool Decay::isDecayHandled(Code const vParticleCode) {
     if (handleAllDecays_ && canHandleDecay(vParticleCode))
       return true;
     else
       return handledDecays_.find(vParticleCode) != Decay::handledDecays_.end();
   }
 
-  void Decay::setStable(std::vector<Code> const& particleList) {
-    for (auto p : particleList) Decay::setStable(p);
+  inline void Decay::setStable(std::vector<Code> const& particleList) {
+    for (auto const p : particleList) Decay::setStable(p);
   }
 
-  void Decay::setUnstable(Code const pCode) {
-    CORSIKA_LOG_INFO("Pythia::Decay: setting {} unstable..", pCode);
-    pythia_.particleData.mayDecay(static_cast<int>(get_PDG(pCode)), true);
+  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);
   }
 
-  void Decay::setStable(Code const pCode) {
-    CORSIKA_LOG_INFO("Pythia::Decay: setting {} stable..", pCode);
-    pythia_.particleData.mayDecay(static_cast<int>(get_PDG(pCode)), false);
+  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);
   }
 
-  bool Decay::isStable(Code const vCode) {
-    return pythia_.particleData.canDecay(static_cast<int>(get_PDG(vCode)));
+  inline bool Decay::isStable(Code const vCode) {
+    return Pythia8::Pythia::particleData.canDecay(static_cast<int>(get_PDG(vCode)));
   }
 
-  bool Decay::canDecay(Code const pCode) {
-    const bool ans = pythia_.particleData.canDecay(static_cast<int>(get_PDG(pCode)));
-    CORSIKA_LOG_INFO("Pythia::Decay: checking if particle: {} can decay in PYTHIA? {} ",
-                     pCode, ans);
+  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;
   }
 
-  void Decay::printDecayConfig(const Code vCode) {
+  inline void Decay::printDecayConfig(Code const vCode) {
     CORSIKA_LOG_INFO("Decay: Pythia decay configuration:");
     CORSIKA_LOG_INFO(" {} is {} ", vCode, (isStable(vCode) ? "stable" : "unstable"));
   }
 
-  void Decay::printDecayConfig() {
+  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& pCode : handledDecays_)
+      for (auto const pCode : handledDecays_)
         CORSIKA_LOG_INFO("Decay of {} is handled by Pythia!", pCode);
   }
 
   template <typename TParticle>
-  TimeType Decay::getLifetime(TParticle const& particle) {
+  inline TimeType Decay::getLifetime(TParticle const& particle) {
 
-    const auto pid = particle.getPID();
+    auto const pid = particle.getPID();
     if (canDecay(pid)) {
       HEPEnergyType E = particle.getEnergy();
       HEPMassType m = particle.getMass();
 
-      const double gamma = E / m;
+      double const gamma = E / m;
 
-      const TimeType t0 = get_lifetime(pid);
+      TimeType const t0 = get_lifetime(pid);
       auto const lifetime = gamma * t0;
-      CORSIKA_LOG_INFO("Pythia::Decay: code: {}", particle.getPID());
-      CORSIKA_LOG_INFO("Pythia::Decay: MinStep: t0: {}", t0);
-      CORSIKA_LOG_INFO("Pythia::Decay: MinStep: energy: {} GeV", E / 1_GeV);
-      CORSIKA_LOG_INFO("Pythia::Decay: momentum: {} GeV",
-                       particle.getMomentum().getComponents() / 1_GeV);
-      CORSIKA_LOG_INFO("Pythia::Decay: MinStep: gamma: {}", gamma);
-      CORSIKA_LOG_INFO("Pythia::Decay: MinStep: tau: {} ", lifetime);
+      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
@@ -174,15 +164,12 @@ namespace corsika::pythia8 {
   }
 
   template <typename TView>
-  void Decay::doDecay(TView& view) {
+  inline void Decay::doDecay(TView& view) {
 
     auto projectile = view.getProjectile();
 
-    auto const& decayPoint = projectile.getPosition();
-    auto const t0 = projectile.getTime();
-
     auto const& labMomentum = projectile.getMomentum();
-    CoordinateSystemPtr const& labCS = labMomentum.getCoordinateSystem();
+    [[maybe_unused]] CoordinateSystemPtr const& labCS = labMomentum.getCoordinateSystem();
 
     // define target kinematics in lab frame
     // define boost to and from CoM frame
@@ -193,7 +180,7 @@ namespace corsika::pythia8 {
     count_++;
 
     // pythia stack
-    Pythia8::Event& event = pythia_.event;
+    Pythia8::Event& event = Pythia8::Pythia::event;
     event.reset();
 
     auto const particleId = projectile.getPID();
@@ -211,37 +198,52 @@ namespace corsika::pythia8 {
     double const m = en;
 
     // add particle to pythia stack
-    event.append(pdgCode, 1, 0, 0, px, py, pz, en, m);
+    event.append(pdgCode, 1, 0, 0, // PID, status, col, acol
+                 px, py, pz, en, m);
 
-    if (!pythia_.next())
+    // LCOV_EXCL_START, we don't validate pythia8 internals
+    if (!Pythia8::Pythia::next())
       throw std::runtime_error("Pythia::Decay: decay failed!");
-    else
-      CORSIKA_LOG_INFO("Pythia::Decay: particles after decay: {} ", event.size());
+    // 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)
+    for (int i = 0; i < event.size(); ++i) {
       if (event[i].isFinal()) {
-        auto const pyId = convert_from_PDG(static_cast<PDGCode>(event[i].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);
-
-        CORSIKA_LOG_INFO("particle: id={} momentum={} energy={} ", pyId,
-                         fourMomLab.getSpaceLikeComponents().getComponents(labCS) / 1_GeV,
-                         fourMomLab.getTimeLikeComponent());
-
-        view.addSecondary(std::make_tuple(pyId, fourMomLab.getTimeLikeComponent(),
-                                          fourMomLab.getSpaceLikeComponents(), decayPoint,
-                                          t0));
+        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);

corsika/detail/modules/pythia8/Interaction.inl

-/*
- * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#pragma once
-
-#include <corsika/modules/pythia8/Interaction.hpp>
-
-#include <corsika/framework/geometry/FourVector.hpp>
-#include <corsika/framework/utility/COMBoost.hpp>
-#include <corsika/media/Environment.hpp>
-#include <corsika/media/NuclearComposition.hpp>
-
-#include <corsika/setup/SetupStack.hpp>
-
-#include <tuple>
-
-namespace corsika::pythia8 {
-
-  Interaction::~Interaction() {
-    std::cout << "Pythia::Interaction n=" << count_ << std::endl;
-  }
-
-  Interaction::Interaction(const bool print_listing)
-      : print_listing_(print_listing) {
-
-    std::cout << "Pythia::Interaction n=" << count_ << std::endl;
-
-    // initialize Pythia
-    if (!initialized_) {
-
-      pythia_.readString("Print:quiet = off");
-      pythia_.readString("Check:particleData = on");      // during init
-      pythia_.readString("Check:event = on");             // default: on
-      pythia_.readString("Check:levelParticleData = 12"); // 1 is default
-      // TODO: proper process initialization for MinBias needed
-      pythia_.readString("HardQCD:all = on");
-      pythia_.readString("ProcessLevel:resonanceDecays = off");
-
-      pythia_.init();
-
-      // any decays in pythia? if yes need to define which particles
-      if (internalDecays_) {
-        // define which particles are passed to corsika, i.e. which particles make it into
-        // history even very shortlived particles like charm or pi0 are of interest here
-        const std::vector<Code> HadronsWeWantTrackedByCorsika = {
-            Code::PiPlus,     Code::PiMinus, Code::Pi0,     Code::KMinus,
-            Code::KPlus,      Code::K0Long,  Code::K0Short, Code::SigmaPlus,
-            Code::SigmaMinus, Code::Lambda0, Code::Xi0,     Code::XiMinus,
-            Code::OmegaMinus, Code::DPlus,   Code::DMinus,  Code::D0,
-            Code::D0Bar};
-
-        Interaction::setStable(HadronsWeWantTrackedByCorsika);
-      }
-
-      // basic initialization of cross section routines
-      sigma_.init(&pythia_.info, pythia_.settings, &pythia_.particleData, &pythia_.rndm);
-
-      initialized_ = true;
-    }
-  }
-
-  void Interaction::setStable(std::vector<Code> const& particleList) {
-    for (auto p : particleList) Interaction::setStable(p);
-  }
-
-  void Interaction::setUnstable(Code const pCode) {
-    std::cout << "Pythia::Interaction: setting " << pCode << " unstable.." << std::endl;
-    pythia_.particleData.mayDecay(static_cast<int>(get_PDG(pCode)), true);
-  }
-
-  void Interaction::setStable(Code const pCode) {
-    std::cout << "Pythia::Interaction: setting " << pCode << " stable.." << std::endl;
-    pythia_.particleData.mayDecay(static_cast<int>(get_PDG(pCode)), false);
-  }
-
-  void Interaction::configureLabFrameCollision(const Code BeamId, const Code TargetId,
-                                               const HEPEnergyType BeamEnergy) {
-    // Pythia configuration of the current event
-    // very clumsy. I am sure this can be done better..
-
-    // set beam
-    // beam id for pythia
-    auto const pdgBeam = static_cast<int>(get_PDG(BeamId));
-    std::stringstream stBeam;
-    stBeam << "Beams:idA = " << pdgBeam;
-    pythia_.readString(stBeam.str());
-    // set target
-    auto pdgTarget = static_cast<int>(get_PDG(TargetId));
-    // replace hydrogen with proton, otherwise pythia goes into heavy ion mode!
-    if (TargetId == Code::Hydrogen) pdgTarget = static_cast<int>(get_PDG(Code::Proton));
-    std::stringstream stTarget;
-    stTarget << "Beams:idB = " << pdgTarget;
-    pythia_.readString(stTarget.str());
-    // set frame to lab. frame
-    pythia_.readString("Beams:frameType = 2");
-    // set beam energy
-    const double Elab = BeamEnergy / 1_GeV;
-    std::stringstream stEnergy;
-    stEnergy << "Beams:eA = " << Elab;
-    pythia_.readString(stEnergy.str());
-    // target at rest
-    pythia_.readString("Beams:eB = 0.");
-    // initialize this config
-    pythia_.init();
-  }
-
-  bool Interaction::canInteract(Code const pCode) {
-    return pCode == Code::Proton || pCode == Code::Neutron || pCode == Code::AntiProton ||
-           pCode == Code::AntiNeutron || pCode == Code::PiMinus || pCode == Code::PiPlus;
-  }
-
-  std::tuple<CrossSectionType, CrossSectionType> Interaction::getCrossSection(
-      const Code BeamId, const Code TargetId, const HEPEnergyType CoMenergy) {
-    // interaction possible in pythia?
-    if (TargetId == Code::Proton || TargetId == Code::Hydrogen) {
-      if (canInteract(BeamId) && isValidCoMEnergy(CoMenergy)) {
-        // input particle PDG
-        auto const pdgCodeBeam = static_cast<int>(get_PDG(BeamId));
-        auto const pdgCodeTarget = static_cast<int>(get_PDG(TargetId));
-        const double ecm = CoMenergy / 1_GeV;
-
-        // calculate cross section
-        sigma_.calc(pdgCodeBeam, pdgCodeTarget, ecm);
-        if (sigma_.hasSigmaTot()) {
-          const double sigEla = sigma_.sigmaEl();
-          const double sigProd = sigma_.sigmaTot() - sigEla;
-
-          return std::make_tuple(sigProd * (1_fm * 1_fm), sigEla * (1_fm * 1_fm));
-
-        } else
-          throw std::runtime_error("pythia cross section init failed");
-
-      } else {
-        return std::make_tuple(std::numeric_limits<double>::infinity() * 1_mb,
-                               std::numeric_limits<double>::infinity() * 1_mb);
-      }
-    } else {
-      throw std::runtime_error("invalid target for pythia");
-    }
-  }
-
-  //  template <>
-  GrammageType Interaction::getInteractionLength(
-      corsika::setup::Stack::particle_type const& particle) {
-
-    // coordinate system, get global frame of reference
-    MomentumVector const& pMomentum = particle.getMomentum();
-    CoordinateSystemPtr const& labCS = pMomentum.getCoordinateSystem();
-
-    Code const corsikaBeamId = particle.getPID();
-
-    // beam particles for pythia : 1, 2, 3 for p, pi, k
-    // read from cross section code table
-    bool const kInteraction = canInteract(corsikaBeamId);
-
-    // FOR NOW: assume target is at rest
-    MomentumVector pTarget(labCS, {0_GeV, 0_GeV, 0_GeV});
-
-    // total momentum and energy
-    HEPEnergyType Elab = particle.getEnergy() + constants::nucleonMass;
-    MomentumVector pTotLab(labCS, {0_GeV, 0_GeV, 0_GeV});
-    pTotLab += pMomentum;
-    pTotLab += pTarget;
-    auto const pTotLabNorm = pTotLab.getNorm();
-    // calculate cm. energy
-    const HEPEnergyType ECoM = sqrt(
-        (Elab + pTotLabNorm) * (Elab - pTotLabNorm)); // binomial for numerical accuracy
-
-    std::cout << "Interaction: LambdaInt: \n"
-              << " input energy: " << particle.getEnergy() / 1_GeV << std::endl
-              << " beam can interact:" << kInteraction << std::endl
-              << " beam pid:" << particle.getPID() << std::endl;
-
-    // TODO: move limits into variables
-    if (kInteraction && Elab >= 8.5_GeV && isValidCoMEnergy(ECoM)) {
-
-      // get target from environment
-      /*
-        the target should be defined by the Environment,
-        ideally as full particle object so that the four momenta
-        and the boosts can be defined..
-      */
-      const auto* currentNode = particle.getNode();
-      const auto mediumComposition =
-          currentNode->getModelProperties().getNuclearComposition();
-      // determine average interaction length
-
-      auto const weightedProdCrossSection =
-          mediumComposition.getWeightedSum([=](auto vTargetID) {
-            return std::get<0>(this->getCrossSection(corsikaBeamId, vTargetID, ECoM));
-          });
-
-      std::cout << "Interaction: IntLength: weighted CrossSection (mb): "
-                << weightedProdCrossSection / 1_mb << std::endl
-                << "Interaction: IntLength: average mass number: "
-                << mediumComposition.getAverageMassNumber() << std::endl;
-
-      // calculate interaction length in medium
-      GrammageType const int_length = mediumComposition.getAverageMassNumber() *
-                                      constants::u / weightedProdCrossSection;
-      std::cout << "Interaction: "
-                << "interaction length (g/cm2): " << int_length / (0.001_kg) * 1_cm * 1_cm
-                << std::endl;
-
-      return int_length;
-    }
-
-    return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
-  }
-
-  /**
-     In this function PYTHIA is called to produce one event. The
-     event is copied (and boosted) into the shower lab frame.
-   */
-
-  template <class TView>
-  void Interaction::doInteraction(TView& view) {
-
-    auto projectile = view.getProjectile();
-
-    const auto corsikaBeamId = projectile.getPID();
-    std::cout << "Pythia::Interaction: "
-              << "DoInteraction: " << corsikaBeamId << " interaction? "
-              << corsika::pythia8::Interaction::canInteract(corsikaBeamId) << std::endl;
-
-    if (is_nucleus(corsikaBeamId)) {
-      // nuclei handled by different process, this should not happen
-      throw std::runtime_error("Nuclear projectile are not handled by PYTHIA!");
-    }
-
-    if (corsika::pythia8::Interaction::canInteract(corsikaBeamId)) {
-
-      // define projectile
-      HEPEnergyType const eProjectileLab = projectile.getEnergy();
-      auto const pProjectileLab = projectile.getMomentum();
-      CoordinateSystemPtr const& labCS = pProjectileLab.getCoordinateSystem();
-
-      // position and time of interaction, not used in Sibyll
-      Point pOrig = projectile.getPosition();
-      TimeType tOrig = projectile.getTime();
-
-      // define target
-      // FOR NOW: target is always at rest
-      const auto eTargetLab = 0_GeV + constants::nucleonMass;
-      const auto pTargetLab = MomentumVector(labCS, 0_GeV, 0_GeV, 0_GeV);
-      const FourVector PtargLab(eTargetLab, pTargetLab);
-
-      std::cout << "Interaction: ebeam lab: " << eProjectileLab / 1_GeV << std::endl
-                << "Interaction: pbeam lab: " << pProjectileLab.getComponents() / 1_GeV
-                << std::endl;
-      std::cout << "Interaction: etarget lab: " << eTargetLab / 1_GeV << std::endl
-                << "Interaction: ptarget lab: " << pTargetLab.getComponents() / 1_GeV
-                << std::endl;
-
-      const FourVector PprojLab(eProjectileLab, pProjectileLab);
-
-      // define target kinematics in lab frame
-      // define boost to and from CoM frame
-      // CoM frame definition in Pythia projectile: +z
-      COMBoost const boost(PprojLab, constants::nucleonMass);
-
-      // just for show:
-      // boost projecticle
-      auto const PprojCoM = boost.toCoM(PprojLab);
-
-      // boost target
-      auto const PtargCoM = boost.toCoM(PtargLab);
-
-      std::cout << "Interaction: ebeam CoM: " << PprojCoM.getTimeLikeComponent() / 1_GeV
-                << std::endl
-                << "Interaction: pbeam CoM: "
-                << PprojCoM.getSpaceLikeComponents().getComponents() / 1_GeV << std::endl;
-      std::cout << "Interaction: etarget CoM: " << PtargCoM.getTimeLikeComponent() / 1_GeV
-                << std::endl
-                << "Interaction: ptarget CoM: "
-                << PtargCoM.getSpaceLikeComponents().getComponents() / 1_GeV << std::endl;
-
-      std::cout << "Interaction: position of interaction: " << pOrig.getCoordinates()
-                << std::endl;
-      std::cout << "Interaction: time: " << tOrig << std::endl;
-
-      HEPEnergyType Etot = eProjectileLab + eTargetLab;
-      MomentumVector Ptot = projectile.getMomentum();
-      // invariant mass, i.e. cm. energy
-      HEPEnergyType Ecm = sqrt(Etot * Etot - Ptot.getSquaredNorm());
-
-      // sample target mass number
-      const auto* currentNode = projectile.getNode();
-      const auto& mediumComposition =
-          currentNode->getModelProperties().getNuclearComposition();
-      // get cross sections for target materials
-      /*
-        Here we read the cross section from the interaction model again,
-        should be passed from getInteractionLength if possible
-       */
-      //#warning reading interaction cross section again, should not be necessary
-      auto const& compVec = mediumComposition.getComponents();
-      std::vector<si::CrossSectionType> cross_section_of_components(compVec.size());
-
-      for (size_t i = 0; i < compVec.size(); ++i) {
-        auto const targetId = compVec[i];
-        const auto [sigProd, sigEla] = getCrossSection(corsikaBeamId, targetId, Ecm);
-        [[maybe_unused]] const auto& dummy_sigEla = sigEla;
-        cross_section_of_components[i] = sigProd;
-      }
-
-      const auto corsikaTargetId =
-          mediumComposition.sampleTarget(cross_section_of_components, RNG_);
-      std::cout << "Interaction: target selected: " << corsikaTargetId << std::endl;
-
-      if (corsikaTargetId != Code::Hydrogen && corsikaTargetId != Code::Neutron &&
-          corsikaTargetId != Code::Proton)
-        throw std::runtime_error("DoInteraction: wrong target for PYTHIA");
-
-      std::cout << "Interaction: "
-                << " DoInteraction: E(GeV):" << eProjectileLab / 1_GeV
-                << " Ecm(GeV): " << Ecm / 1_GeV << std::endl;
-
-      if (eProjectileLab < 8.5_GeV || !isValidCoMEnergy(Ecm)) {
-        std::cout << "Interaction: "
-                  << " DoInteraction: should have dropped particle.. "
-                  << "THIS IS AN ERROR" << std::endl;
-        throw std::runtime_error("energy too low for PYTHIA");
-
-      } else {
-        count_++;
-
-        configureLabFrameCollision(corsikaBeamId, corsikaTargetId, eProjectileLab);
-
-        // create event in pytia
-        if (!pythia_.next()) throw std::runtime_error("Pythia::DoInteraction: failed!");
-
-        // link to pythia stack
-        Pythia8::Event& event = pythia_.event;
-
-	if (print_listing_) {
-	  // print final state
-	  event.list();
-	}
-
-        MomentumVector Plab_final(labCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
-        HEPEnergyType Elab_final = 0_GeV;
-        for (int i = 0; i < event.size(); ++i) {
-          Pythia8::Particle& p8p = event[i];
-          // skip particles that have decayed in pythia
-          if (!p8p.isFinal()) continue;
-
-          auto const pyId = convert_from_PDG(static_cast<PDGCode>(p8p.id()));
-
-          const MomentumVector pyPlab(
-              labCS, {p8p.px() * 1_GeV, p8p.py() * 1_GeV, p8p.pz() * 1_GeV});
-          HEPEnergyType const pyEn = p8p.e() * 1_GeV;
-
-          // add to corsika stack
-          auto pnew =
-              projectile.addSecondary(std::make_tuple(pyId, pyEn, pyPlab, pOrig, tOrig));
-
-          Plab_final += pnew.getMomentum();
-          Elab_final += pnew.getEnergy();
-        }
-        std::cout << "conservation (all GeV): "
-                  << "Elab_final=" << Elab_final / 1_GeV
-                  << ", Plab_final=" << (Plab_final / 1_GeV).getComponents() << std::endl;
-      }
-    }
-  }
-
-} // 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 GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
+ * 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
@@ -12,6 +11,6 @@
 
 namespace corsika::pythia8 {
 
-  double Random::flat() { return Dist_(RNG_); }
+  inline double Random::flat() { return Dist_(RNG_); }
 
 } // namespace corsika::pythia8

corsika/detail/modules/qgsjetII/Interaction.inl

-/*
- * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#include <corsika/modules/qgsjetII/Interaction.hpp>
-
-#include <corsika/media/Environment.hpp>
-#include <corsika/media/NuclearComposition.hpp>
-#include <corsika/framework/geometry/QuantityVector.hpp>
-#include <corsika/framework/geometry/FourVector.hpp>
-#include <corsika/modules/qgsjetII/ParticleConversion.hpp>
-#include <corsika/modules/qgsjetII/QGSJetIIFragmentsStack.hpp>
-#include <corsika/modules/qgsjetII/QGSJetIIStack.hpp>
-#include <corsika/setup/SetupStack.hpp>
-#include <corsika/setup/SetupTrajectory.hpp>
-#include <corsika/framework/utility/COMBoost.hpp>
-
-#include <sstream>
-#include <string>
-#include <tuple>
-
-#include <qgsjet-II-04.hpp>
-
-namespace corsika::qgsjetII {
-
-  Interaction::Interaction(const std::string& dataPath)
-      : data_path_(dataPath) {
-    if (dataPath == "") {
-      if (std::getenv("CORSIKA_DATA")) {
-        data_path_ = std::string(std::getenv("CORSIKA_DATA")) + "/QGSJetII/";
-        std::cout << "Searching for QGSJetII data tables in " << data_path_ << std::endl;
-      }
-    }
-
-    // initialize QgsjetII
-    static bool initialized = false;
-    if (!initialized) {
-      qgset_();
-      datadir DIR(data_path_);
-      qgaini_(DIR.data);
-      initialized = true;
-    }
-  }
-
-  Interaction::~Interaction() {
-    std::cout << "QgsjetII::Interaction n=" << count_ << std::endl;
-  }
-
-  CrossSectionType Interaction::getCrossSection(const Code beamId, const Code targetId,
-                                                const HEPEnergyType Elab,
-                                                const unsigned int Abeam,
-                                                const unsigned int targetA) const {
-    double sigProd = std::numeric_limits<double>::infinity();
-
-    if (corsika::qgsjetII::canInteract(beamId)) {
-
-      int const iBeam = static_cast<QgsjetIIXSClassIntType>(
-          corsika::qgsjetII::getQgsjetIIXSCode(beamId));
-      int iTarget = 1;
-      if (is_nucleus(targetId)) {
-        iTarget = targetA;
-        if (iTarget > maxMassNumber_ || iTarget <= 0) {
-          std::ostringstream txt;
-          txt << "QgsjetII target outside range. iTarget=" << iTarget;
-          throw std::runtime_error(txt.str().c_str());
-        }
-      }
-      int iProjectile = 1;
-      if (is_nucleus(beamId)) {
-        iProjectile = Abeam;
-        if (iProjectile > maxMassNumber_ || iProjectile <= 0)
-          throw std::runtime_error("QgsjetII target outside range. ");
-      }
-
-      std::cout << "QgsjetII::getCrossSection Elab=" << Elab << " iBeam=" << iBeam
-                << " iProjectile=" << iProjectile << " iTarget=" << iTarget << std::endl;
-      sigProd = qgsect_(Elab / 1_GeV, iBeam, iProjectile, iTarget);
-      std::cout << "QgsjetII::getCrossSection sigProd=" << sigProd << std::endl;
-    }
-
-    return sigProd * 1_mb;
-  }
-
-  template <typename TParticle>
-  GrammageType Interaction::getInteractionLength(const TParticle& vP) const {
-
-    // coordinate system, get global frame of reference
-    CoordinateSystemPtr const& rootCS = get_root_CoordinateSystem();
-
-    const Code corsikaBeamId = vP.getPID();
-
-    // beam particles for qgsjetII : 1, 2, 3 for p, pi, k
-    // read from cross section code table
-    const bool kInteraction = corsika::qgsjetII::canInteract(corsikaBeamId);
-
-    // FOR NOW: assume target is at rest
-    MomentumVector pTarget(rootCS, {0_GeV, 0_GeV, 0_GeV});
-
-    // total momentum and energy
-    HEPEnergyType Elab = vP.getEnergy();
-
-    std::cout << "Interaction: LambdaInt: \n"
-              << " input energy: " << vP.getEnergy() / 1_GeV << std::endl
-              << " beam can interact:" << kInteraction << std::endl
-              << " beam pid:" << vP.getPID() << std::endl;
-
-    if (kInteraction) {
-
-      int Abeam = 0;
-      if (is_nucleus(vP.getPID())) Abeam = vP.getNuclearA();
-
-      // get target from environment
-      /*
-        the target should be defined by the Environment,
-        ideally as full particle object so that the four momenta
-        and the boosts can be defined..
-      */
-
-      auto const* currentNode = vP.getNode();
-      const auto& mediumComposition =
-          currentNode->getModelProperties().getNuclearComposition();
-
-      CrossSectionType weightedProdCrossSection =
-          mediumComposition.getWeightedSum([=](Code targetID) -> CrossSectionType {
-            int targetA = 0;
-            if (is_nucleus(targetID)) targetA = get_nucleus_A(targetID);
-            return getCrossSection(corsikaBeamId, targetID, Elab, Abeam, targetA);
-          });
-
-      std::cout << "Interaction: "
-                << "IntLength: weighted CrossSection (mb): "
-                << weightedProdCrossSection / 1_mb << std::endl;
-
-      // calculate interaction length in medium
-      GrammageType const int_length = mediumComposition.getAverageMassNumber() *
-                                      constants::u / weightedProdCrossSection;
-      std::cout << "Interaction: "
-                << "interaction length (g/cm2): " << int_length / (0.001_kg) * 1_cm * 1_cm
-                << std::endl;
-
-      return int_length;
-    }
-
-    return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
-  }
-
-  /**
-     In this function QGSJETII is called to produce one event. The
-     event is copied (and boosted) into the shower lab frame.
-   */
-
-  template <typename TView>
-  void Interaction::doInteraction(TView& view) {
-
-    auto const projectile = view.getProjectile();
-
-    const auto corsikaBeamId = projectile.getPID();
-    std::cout << "ProcessQgsjetII: "
-              << "DoInteraction: " << corsikaBeamId << " interaction? "
-              << corsika::qgsjetII::canInteract(corsikaBeamId) << std::endl;
-
-    if (corsika::qgsjetII::canInteract(corsikaBeamId)) {
-
-      CoordinateSystemPtr const& rootCS = get_root_CoordinateSystem();
-
-      // position and time of interaction, not used in QgsjetII
-      Point pOrig = projectile.getPosition();
-      TimeType tOrig = projectile.getTime();
-
-      // define target
-      // for QgsjetII is always a single nucleon
-      // FOR NOW: target is always at rest
-      auto const targetEnergyLab = 0_GeV + constants::nucleonMass;
-      auto const targetMomentumLab = MomentumVector(rootCS, 0_GeV, 0_GeV, 0_GeV);
-      FourVector const PtargLab(targetEnergyLab, targetMomentumLab);
-
-      // define projectile
-      HEPEnergyType const projectileEnergyLab = projectile.getEnergy();
-      auto const projectileMomentumLab = projectile.getMomentum();
-
-      int beamA = 0;
-      if (is_nucleus(corsikaBeamId)) beamA = projectile.getNuclearA();
-
-      HEPEnergyType const projectileEnergyLabPerNucleon = projectileEnergyLab / beamA;
-
-      std::cout << "Interaction: ebeam lab: " << projectileEnergyLab / 1_GeV << std::endl
-                << "Interaction: pbeam lab: "
-                << projectileMomentumLab.getComponents() / 1_GeV << std::endl;
-      std::cout << "Interaction: etarget lab: " << targetEnergyLab / 1_GeV << std::endl
-                << "Interaction: ptarget lab: "
-                << targetMomentumLab.getComponents() / 1_GeV << std::endl;
-
-      std::cout << "Interaction: position of interaction: " << pOrig.getCoordinates()
-                << std::endl;
-      std::cout << "Interaction: time: " << tOrig << std::endl;
-
-      // sample target mass number
-      auto const* currentNode = projectile.getNode();
-      auto const& mediumComposition =
-          currentNode->getModelProperties().getNuclearComposition();
-      // get cross sections for target materials
-      /*
-        Here we read the cross section from the interaction model again,
-        should be passed from getInteractionLength if possible
-       */
-      auto const& compVec = mediumComposition.getComponents();
-      std::vector<CrossSectionType> cross_section_of_components(compVec.size());
-
-      for (size_t i = 0; i < compVec.size(); ++i) {
-        auto const targetId = compVec[i];
-        int targetA = 0;
-        if (is_nucleus(targetId)) targetA = get_nucleus_A(targetId);
-        const auto sigProd =
-            getCrossSection(corsikaBeamId, targetId, projectileEnergyLab, beamA, targetA);
-        cross_section_of_components[i] = sigProd;
-      }
-
-      const auto targetCode =
-          mediumComposition.sampleTarget(cross_section_of_components, rng_);
-      std::cout << "Interaction: target selected: " << targetCode << std::endl;
-
-      int targetMassNumber = 1;               // proton
-      if (is_nucleus(targetCode)) { // nucleus
-        targetMassNumber = get_nucleus_A(targetCode);
-        if (targetMassNumber > maxMassNumber_)
-          throw std::runtime_error("QgsjetII target mass outside range.");
-      } else {
-        if (targetCode != Proton::code)
-          throw std::runtime_error("QgsjetII Taget not possible.");
-      }
-      std::cout << "Interaction: target qgsjetII code/A: " << targetMassNumber << std::endl;
-
-      int projectileMassNumber = 1; // "1" means "hadron"
-      QgsjetIIHadronType qgsjet_hadron_type =
-          qgsjetII::getQgsjetIIHadronType(corsikaBeamId);
-      if (qgsjet_hadron_type == QgsjetIIHadronType::NucleusType) {
-        projectileMassNumber = projectile.getNuclearA();
-        if (projectileMassNumber > maxMassNumber_)
-          throw std::runtime_error("QgsjetII projectile mass outside range.");
-        std::array<QgsjetIIHadronType, 2> constexpr nucleons = {
-            QgsjetIIHadronType::ProtonType, QgsjetIIHadronType::NeutronType};
-        std::uniform_int_distribution select(0, 1);
-        qgsjet_hadron_type = nucleons[select(rng_)];
-      } else {
-        // from conex: replace pi0 or rho0 with pi+/pi- in alternating sequence
-        if (qgsjet_hadron_type == QgsjetIIHadronType::NeutralLightMesonType) {
-          qgsjet_hadron_type = alternate_;
-          alternate_ = (alternate_ == QgsjetIIHadronType::PiPlusType
-                            ? QgsjetIIHadronType::PiMinusType
-                            : QgsjetIIHadronType::PiPlusType);
-        }
-      }
-
-      // beam id for qgsjetII
-      int kBeam = 2; // default: proton Shouldn't we randomize neutron/proton for nuclei?
-      if (corsikaBeamId != Code::Nucleus) {
-        kBeam = corsika::qgsjetII::convertToQgsjetIIRaw(corsikaBeamId);
-        // from conex
-        if (kBeam == 0) { // replace pi0 or rho0 with pi+/pi-
-          static int select = 1;
-          kBeam = select;
-          select *= -1;
-        }
-        // replace lambda by neutron
-        if (kBeam == 6)
-          kBeam = 3;
-        else if (kBeam == -6)
-          kBeam = -3;
-        // else if (abs(kBeam)>6) -> throw
-      }
-
-      std::cout << "Interaction: "
-                << " DoInteraction: E(GeV):" << projectileEnergyLab / 1_GeV << std::endl;
-      count_++;
-      int qgsjet_hadron_type_int = static_cast<QgsjetIICodeIntType>(qgsjet_hadron_type);
-      qgini_(projectileEnergyLab / 1_GeV, qgsjet_hadron_type_int, projectileMassNumber,
-             targetMassNumber);
-      qgconf_();
-
-      // 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
-      auto const& originalCS = projectileMomentumLab.getCoordinateSystem();
-      CoordinateSystemPtr const zAxisFrame =
-          make_rotationToZ(originalCS, projectileMomentumLab);
-
-      // fragments
-      QGSJetIIFragmentsStack qfs;
-      for (auto& fragm : qfs) {
-        Code idFragm = Code::Nucleus;
-        int A = fragm.getFragmentSize();
-        int Z = 0;
-        switch (A) {
-          case 1: { // proton/neutron
-            std::uniform_real_distribution<double> select;
-            if (select(rng_) > 0.5) {
-              idFragm = Code::Proton;
-              Z = 1;
-            } else {
-              idFragm = Code::Neutron;
-              Z = 0;
-            }
-
-            const HEPMassType nucleonMass = get_mass(idFragm);
-
-            auto momentum = Vector(
-                zAxisFrame, corsika::QuantityVector<hepmomentum_d>{
-                                0.0_GeV, 0.0_GeV,
-                                sqrt((projectileEnergyLabPerNucleon + nucleonMass) *
-                                     (projectileEnergyLabPerNucleon - nucleonMass))});
-
-            auto const energy = sqrt(momentum.getSquaredNorm() + square(nucleonMass));
-            momentum.rebase(originalCS); // transform back into standard lab frame
-            std::cout << "secondary fragment> id=" << idFragm
-                      << " p=" << momentum.getComponents() << std::endl;
-            auto pnew = view.addSecondary(
-                std::make_tuple(idFragm, energy, momentum,
-                                                            pOrig, tOrig));
-            Plab_final += pnew.getMomentum();
-            Elab_final += pnew.getEnergy();
-          } break;
-          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);
-          }
-        }
-
-        if (idFragm == Code::Nucleus) { // thus: not p or n
-          const HEPMassType nucleusMass =
-              Proton::mass * Z + Neutron::mass * (A - Z);
-          auto momentum = Vector(
-              zAxisFrame, QuantityVector<hepmomentum_d>{
-                              0.0_GeV, 0.0_GeV,
-                              sqrt((projectileEnergyLabPerNucleon * A + nucleusMass) *
-                                   (projectileEnergyLabPerNucleon * A - nucleusMass))});
-
-          auto const energy = sqrt(momentum.getSquaredNorm() + square(nucleusMass));
-          momentum.rebase(originalCS); // transform back into standard lab frame
-          std::cout << "secondary fragment> id=" << idFragm
-                    << " p=" << momentum.getComponents() << " A=" << A << " Z=" << Z
-                    << std::endl;
-          auto pnew = view.addSecondary(
-              std::make_tuple(
-                  idFragm, energy, momentum, pOrig, tOrig, A, Z));
-          Plab_final += pnew.getMomentum();
-          Elab_final += pnew.getEnergy();
-        }
-      }
-
-      // secondaries
-      QGSJetIIStack qs;
-      for (auto& psec : qs) {
-
-        auto momentum = psec.getMomentum(zAxisFrame);
-        auto const energy = psec.getEnergy();
-
-        momentum.rebase(originalCS); // transform back into standard lab frame
-        std::cout << "secondary fragment> id="
-                  << corsika::qgsjetII::convertFromQgsjetII(psec.getPID())
-                  << " p=" << momentum.getComponents() << std::endl;
-        auto pnew = view.addSecondary(
-            std::make_tuple(corsika::qgsjetII::convertFromQgsjetII(psec.getPID()), energy,
-                            momentum, pOrig, tOrig));
-        Plab_final += pnew.getMomentum();
-        Elab_final += pnew.getEnergy();
-      }
-      std::cout << "conservation (all GeV): Ecm_final= n/a" /* << Ecm_final / 1_GeV*/
-                << std::endl
-                << "Elab_final=" << Elab_final / 1_GeV
-                << ", Plab_final=" << (Plab_final / 1_GeV).getComponents()
-                << ", N_wounded,targ="
-                << QGSJetIIFragmentsStackData::getWoundedNucleonsTarget()
-                << ", N_wounded,proj="
-                << QGSJetIIFragmentsStackData::getWoundedNucleonsProjectile()
-                << ", N_fragm,proj=" << qfs.getSize() << std::endl;
-    }
-  }
-
-} // namespace corsika::qgsjetII

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 GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
+ * 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>

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 GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
+ * 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 {
 
-  void QGSJetIIStackData::clear() {
+  inline void QGSJetIIStackData::clear() {
     qgarr12_.nsp = 0;
     qgarr13_.nsf = 0;
     qgarr55_.nwt = 0;
   }
-  unsigned int QGSJetIIStackData::getSize() const { return qgarr12_.nsp; }
-  unsigned int QGSJetIIStackData::getCapacity() const { return nptmax; }
+  inline unsigned int QGSJetIIStackData::getSize() const { return qgarr12_.nsp; }
+  inline unsigned int QGSJetIIStackData::getCapacity() const { return nptmax; }
 
-  void QGSJetIIStackData::setId(const unsigned int i, const int v) {
+  inline void QGSJetIIStackData::setId(const unsigned int i, const int v) {
     qgarr14_.ich[i] = v;
   }
-  void QGSJetIIStackData::setEnergy(const unsigned int i, const HEPEnergyType v) {
+  inline void QGSJetIIStackData::setEnergy(const unsigned int i, const HEPEnergyType v) {
     qgarr14_.esp[i][0] = v / 1_GeV;
   }
 
-  void QGSJetIIStackData::setMomentum(const unsigned int i, const MomentumVector& v) {
+  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;
   }
 
-  int QGSJetIIStackData::getId(const unsigned int i) const { return qgarr14_.ich[i]; }
-  HEPEnergyType QGSJetIIStackData::getEnergy(const int i) const {
+  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;
   }
-  MomentumVector QGSJetIIStackData::getMomentum(const unsigned int i,
-                                                const CoordinateSystemPtr& CS) const {
+  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);
   }
 
-  void QGSJetIIStackData::copy(const unsigned int i1, const unsigned int i2) {
+  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];
   }
 
-  void QGSJetIIStackData::swap(const unsigned int i1, const unsigned int i2) {
+  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]);
   }
 
-  void QGSJetIIStackData::incrementSize() { qgarr12_.nsp++; }
-  void QGSJetIIStackData::decrementSize() {
+  inline void QGSJetIIStackData::incrementSize() { qgarr12_.nsp++; }
+  inline void QGSJetIIStackData::decrementSize() {
     if (qgarr12_.nsp > 0) { qgarr12_.nsp--; }
   }
 
   template <typename StackIteratorInterface>
-  void ParticleInterface<StackIteratorInterface>::setParticleData(
+  inline void ParticleInterface<StackIteratorInterface>::setParticleData(
       const int vID, const HEPEnergyType vE, const MomentumVector& vP,
       const HEPMassType) {
     setPID(vID);
@@ -70,7 +72,7 @@ namespace corsika::qgsjetII {
   }
 
   template <typename StackIteratorInterface>
-  void ParticleInterface<StackIteratorInterface>::setParticleData(
+  inline void ParticleInterface<StackIteratorInterface>::setParticleData(
       ParticleInterface<StackIteratorInterface>& /*parent*/, const int vID,
       const HEPEnergyType vE, const MomentumVector& vP, const HEPMassType) {
     setPID(vID);
@@ -79,34 +81,36 @@ namespace corsika::qgsjetII {
   }
 
   template <typename StackIteratorInterface>
-  void ParticleInterface<StackIteratorInterface>::setEnergy(const HEPEnergyType v) {
+  inline void ParticleInterface<StackIteratorInterface>::setEnergy(
+      const HEPEnergyType v) {
     getStackData().setEnergy(getIndex(), v);
   }
 
   template <typename StackIteratorInterface>
-  HEPEnergyType ParticleInterface<StackIteratorInterface>::getEnergy() const {
+  inline HEPEnergyType ParticleInterface<StackIteratorInterface>::getEnergy() const {
     return getStackData().getEnergy(getIndex());
   }
 
   template <typename StackIteratorInterface>
-  void ParticleInterface<StackIteratorInterface>::setPID(const int v) {
+  inline void ParticleInterface<StackIteratorInterface>::setPID(const int v) {
     getStackData().setId(getIndex(), v);
   }
 
   template <typename StackIteratorInterface>
-  corsika::qgsjetII::QgsjetIICode ParticleInterface<StackIteratorInterface>::getPID()
-      const {
+  inline corsika::qgsjetII::QgsjetIICode
+  ParticleInterface<StackIteratorInterface>::getPID() const {
     return static_cast<corsika::qgsjetII::QgsjetIICode>(getStackData().getId(getIndex()));
   }
 
   template <typename StackIteratorInterface>
-  MomentumVector ParticleInterface<StackIteratorInterface>::getMomentum(
+  inline MomentumVector ParticleInterface<StackIteratorInterface>::getMomentum(
       const CoordinateSystemPtr& CS) const {
     return getStackData().getMomentum(getIndex(), CS);
   }
 
   template <typename StackIteratorInterface>
-  void ParticleInterface<StackIteratorInterface>::setMomentum(const MomentumVector& v) {
+  inline void ParticleInterface<StackIteratorInterface>::setMomentum(
+      const MomentumVector& v) {
     getStackData().setMomentum(getIndex(), v);
   }
 

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 GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
+ * 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
@@ -11,14 +10,14 @@
 #include <corsika/modules/qgsjetII/qgsjet-II-04.hpp>
 
 #include <corsika/framework/random/RNGManager.hpp>
+#include <corsika/framework/core/Logging.hpp>
 
 #include <iostream>
 #include <random>
 
-datadir::datadir(const std::string& dir) {
+inline datadir::datadir(const std::string& dir) {
   if (dir.length() > 130) {
-    std::cerr << "QGSJetII error, will cut datadir \"" << dir
-              << "\" to 130 characters: " << std::endl;
+    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];
@@ -26,14 +25,6 @@ datadir::datadir(const std::string& dir) {
   data[i + 1] = '\0';
 }
 
-double qgran_(int&) {
-  static corsika::default_prng_type& rng =
-      corsika::RNGManager::getInstance().GetRandomStream("qgran");
-
-  std::uniform_real_distribution<double> dist;
-  return dist(rng);
-}
-
-void lzmaopenfile_(const char*, int) {}
-void lzmaclosefile_() {}
-void lzmafillarray_(const double&, const int&) {}
+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