applications/CMakeLists.txt

@@ -12,4 +12,4 @@ endif()
 
 install (
   TARGETS c8_air_shower DESTINATION bin
-)
+)
\ No newline at end of file

applications/c8_air_shower.cpp

@@ -250,7 +250,7 @@ int main(int argc, char** argv) {
       ->group("Misc.");
   app.add_option("-M,--hadronModel", "High-energy hadronic interaction model")
       ->default_val("SIBYLL-2.3d")
-      ->check(CLI::IsMember({"SIBYLL-2.3d", "QGSJet-II.04", "EPOS-LHC"}))
+      ->check(CLI::IsMember({"SIBYLL-2.3d", "QGSJet-II.04", "EPOS-LHC", "Pythia8"}))
       ->group("Misc.");
   app.add_option("-T,--hadronModelTransitionEnergy",
                  "Transition between high-/low-energy hadronic interaction "
@@ -406,6 +406,13 @@ int main(int argc, char** argv) {
   EnergyLossWriter dEdX{showerAxis, dX};
   output.add("energyloss", dEdX);
 
+  // energy threshold for high energy hadronic model. Affects LE/HE switch for
+  // hadron interactions and the hadronic photon model in proposal
+  // valid range varies a bit between models. pythia only has cross sections tabulated
+  // down to 100GeV.
+  HEPEnergyType const heHadronModelThreshold =
+      1_GeV * app["--hadronModelTransitionEnergy"]->as<double>();
+
   DynamicInteractionProcess<StackType> heModel;
 
   auto const all_elements = corsika::get_all_elements_in_universe(env);
@@ -422,6 +429,15 @@ int main(int argc, char** argv) {
   } else if (modelStr == "EPOS-LHC") {
     heModel = DynamicInteractionProcess<StackType>{
         std::make_shared<corsika::epos::Interaction>(corsika::setup::C7trackedParticles)};
+  } else if (modelStr == "Pythia8") {
+    heModel = DynamicInteractionProcess<StackType>{
+        std::make_shared<corsika::pythia8::Interaction>(
+            corsika::setup::C7trackedParticles)};
+    if (heHadronModelThreshold < 100_GeV) {
+      CORSIKA_LOG_CRITICAL(
+          "threshold too low for Pythia! use --hadronModelTransitionEnergy>=100_GeV!");
+      return EXIT_FAILURE;
+    }
   } else {
     CORSIKA_LOG_CRITICAL("invalid choice \"{}\"; also check argument parser", modelStr);
     return EXIT_FAILURE;
@@ -477,11 +493,6 @@ int main(int argc, char** argv) {
   set_energy_production_threshold(Code::TauMinus, prod_threshold);
   set_energy_production_threshold(Code::TauPlus, prod_threshold);
 
-  // energy threshold for high energy hadronic model. Affects LE/HE switch for
-  // hadron interactions and the hadronic photon model in proposal
-  HEPEnergyType const heHadronModelThreshold =
-      1_GeV * app["--hadronModelTransitionEnergy"]->as<double>();
-
   corsika::proposal::Interaction emCascade(
       env, sophia, sibyll->getHadronInteractionModel(), heHadronModelThreshold);
 

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 <tuple>
-#include <stdexcept>
-
-#include <boost/filesystem/path.hpp>
-#include <fmt/core.h>
-
-#include <corsika/modules/pythia8/Interaction.hpp>
-
-#include <corsika/framework/geometry/FourVector.hpp>
-#include <corsika/framework/utility/COMBoost.hpp>
-#include <corsika/framework/core/EnergyMomentumOperations.hpp>
-#include <corsika/media/Environment.hpp>
-#include <corsika/media/NuclearComposition.hpp>
-
-namespace corsika::pythia8 {
-
-  inline Interaction::~Interaction() {
-    CORSIKA_LOG_INFO("Pythia::Interaction n= {}", count_);
-  }
-
-  inline Interaction::Interaction(boost::filesystem::path const& mpiInitFile,
-                                  bool const print_listing)
-      : print_listing_(print_listing)
-      , pythiaMain_{CORSIKA_Pythia8_XML_DIR, false}
-      , pythiaColl_{CORSIKA_Pythia8_XML_DIR, false} {
-    auto rndm = std::make_shared<corsika::pythia8::Random>();
-    pythiaColl_.setRndmEnginePtr(rndm);
-    pythiaMain_.setRndmEnginePtr(rndm);
-
-    CORSIKA_LOG_INFO("Pythia8 MPI init file: {}", mpiInitFile.native());
-    // Main Pythia object for managing the cascade evolution.
-    // Can also do decays, but no hard processes.
-
-    pythiaMain_.readString("ProcessLevel:all = off");
-
-    // Reduce statistics printout to relevant ones.
-    pythiaMain_.readString("Stat:showProcessLevel = off");
-    pythiaMain_.readString("Stat:showPartonLevel = off");
-
-    // Add Argon, since not in Particle data. id:all = name antiName
-    // spinType chargeType colType m0 mWidth mMin mMax tau0.
-    pythiaMain_.readString(
-        "1000180400:all = 40Ar 40Arbar 1 54 0 "
-        "37.22474 0. 0. 0. 0.");
-
-    // we can't test this block, LCOV_EXCL_START
-    if (!pythiaMain_.init())
-      throw std::runtime_error("Pythia::Interaction: Initialization failed!");
-    // LCOV_EXCL_STOP
-
-    // Secondary Pythia object for performing individual collisions.
-    // Variable incoming beam type and energy.
-    pythiaColl_.readString("Beams:allowVariableEnergy = on");
-    pythiaColl_.readString("Beams:allowIDAswitch = on");
-
-    // Set up for fixed-target collisions.
-    pythiaColl_.readString(
-        "Beams:frameType = 3"); // arbitrary frame, need to define full 4-momenta
-    pythiaColl_.settings.parm("Beams:pzA", eMaxLab_ / 1_GeV);
-    pythiaColl_.readString("Beams:pzB = 0.");
-
-    // Must use the soft and low-energy QCD processes.
-    pythiaColl_.readString("SoftQCD:all = on");
-    pythiaColl_.readString("LowEnergyQCD:all = on");
-
-    // Decays to be done by pythiaMain_.
-    pythiaColl_.readString("HadronLevel:Decay = off");
-
-    // Reduce printout and relax energy-momentum conservation.
-    pythiaColl_.readString("Print:quiet = on");
-    pythiaColl_.readString("Check:epTolErr = 0.1");
-    pythiaColl_.readString("Check:epTolWarn = 0.0001");
-    pythiaColl_.readString("Check:mTolErr = 0.01");
-
-    // Redure statistics printout to relevant ones.
-    pythiaColl_.readString("Stat:showProcessLevel = off");
-    pythiaColl_.readString("Stat:showPartonLevel = off");
-
-    bool const reuse = true; // could be made more flexible
-    // Reuse MPI initialization file if it exists; else create a new one.
-    if (reuse)
-      pythiaColl_.readString("MultipartonInteractions:reuseInit = 3");
-    else
-      pythiaColl_.readString("MultipartonInteractions:reuseInit = 1");
-    pythiaColl_.settings.word("MultipartonInteractions:initFile", mpiInitFile.native());
-
-    // initialize
-    // we can't test this block, LCOV_EXCL_START
-    if (!pythiaColl_.init())
-      throw std::runtime_error("Pythia::Interaction: Initialization failed!");
-    // LCOV_EXCL_STOP
-  }
-
-  inline bool Interaction::isValid(Code const projectileId, Code const targetId,
-                                   HEPEnergyType const sqrtS) const {
-    if (is_nucleus(projectileId)) // not yet possible with Pythia
-      return false;
-
-    HEPEnergyType const labE = calculate_lab_energy(
-        static_pow<2>(sqrtS), get_mass(projectileId), get_mass(targetId));
-    if (labE < eKinMinLab_) return false;
-
-    return std::find(validTargets_.begin(), validTargets_.end(), targetId) !=
-           validTargets_.end();
-  }
-
-  inline bool Interaction::canInteract(Code const pCode) const {
-    return is_hadron(pCode) && !is_nucleus(pCode); // should be sufficient
-  }
-
-  inline std::tuple<CrossSectionType, CrossSectionType>
-  Interaction::getCrossSectionInelEla(Code const projectileId, Code const targetId,
-                                      FourMomentum const& projectileP4,
-                                      FourMomentum const& targetP4) const {
-    HEPEnergyType const CoMenergy =
-        is_nucleus(targetId)
-            ? (projectileP4 + targetP4 / get_nucleus_A(targetId)).getNorm()
-            : (projectileP4 + targetP4).getNorm();
-
-    if (is_nucleus(targetId) || is_nucleus(projectileId)) {
-      CORSIKA_LOG_ERROR(
-          "Pythia8::Interaction::getCrossSectionInelEla() called with nuclear projectile "
-          "or target");
-      return std::make_tuple(CrossSectionType::zero(), CrossSectionType::zero());
-    }
-
-    if (!isValid(projectileId, targetId, CoMenergy)) {
-      return std::make_tuple(CrossSectionType::zero(), CrossSectionType::zero());
-    }
-
-    // input particle PDG
-    auto const pdgCodeBeam = static_cast<int>(get_PDG(projectileId));
-    auto const pdgCodeTarget = static_cast<int>(get_PDG(targetId));
-    double const ecm_GeV = CoMenergy * (1 / 1_GeV);
-
-    auto const sigTot =
-        pythiaColl_.getSigmaTotal(pdgCodeBeam, pdgCodeTarget, ecm_GeV) * 1_mb;
-    auto const sigEla =
-        pythiaColl_.getSigmaPartial(pdgCodeBeam, pdgCodeTarget, ecm_GeV, 2) * 1_mb;
-
-    return std::make_tuple(sigTot - sigEla, sigEla);
-  }
-
-  CrossSectionType Interaction::getCrossSection(Code const projectileId,
-                                                Code const targetId,
-                                                FourMomentum const& projectileP4,
-                                                FourMomentum const& targetP4) const {
-
-    if (!is_nucleus(targetId)) {
-      auto const [sigProd, sigEla] =
-          getCrossSectionInelEla(projectileId, targetId, projectileP4, targetP4);
-      return sigProd + sigEla;
-    } else if (get_nucleus_A(targetId) == 1) {
-      auto const [sigProd, sigEla] = getCrossSectionInelEla(
-          projectileId, get_nucleus_Z(targetId) == 1 ? Code::Proton : Code::Neutron,
-          projectileP4, targetP4);
-      return sigProd + sigEla;
-    } else {
-      auto const [sigProd, sigEla] =
-          getCrossSectionInelEla(projectileId, Code::Proton, projectileP4, targetP4);
-      auto const sigTot = sigProd + sigEla;
-      auto const nSubcoll = getAverageSubcollisions(targetId, sigTot);
-      if (nSubcoll == 0) { // no parameterization available -> we can't handle this
-        return CrossSectionType::zero();
-      } else {
-        return sigTot * get_nucleus_A(targetId) / nSubcoll;
-      }
-    }
-  }
-
-  double Interaction::getAverageSubcollisions(Code targetId,
-                                              CrossSectionType sigTot) const {
-    if (targetId == Code::Proton || targetId == Code::Neutron ||
-        targetId == Code::Hydrogen || targetId == Code::AntiProton ||
-        targetId == Code::AntiNeutron)
-      return 1;
-
-    auto const Z = get_nucleus_Z(targetId);
-    auto const A = get_nucleus_A(targetId);
-
-    double nCollAvg;
-
-    if (Z == 7 && A == 14)
-      nCollAvg = (sigTot < 31._mb) // this is from the paper
-                     ? 1. + 0.017 / 1_mb * sigTot
-                     : 1.2 + 0.0105 / 1_mb * sigTot;
-    else if (Z == 8 && A == 16)
-      nCollAvg = (sigTot < 16._mb) // provided by Torbjörn Sjöstrand
-                     ? 1. + 0.0245 / 1_mb * sigTot
-                     : 1.2 + 0.012 / 1_mb * sigTot;
-    else if (Z == 18 && A == 40)
-      nCollAvg =
-          (sigTot < 10._mb) ? 1. + 0.050 / 1_mb * sigTot : 1.28 + 0.022 / 1_mb * sigTot;
-    else {
-      CORSIKA_LOG_ERROR(
-          fmt::format("Pythia8 cross-sections not defined for ({},{}) nucleus", A, Z));
-      nCollAvg = 0;
-    }
-
-    return nCollAvg;
-  }
-
-  template <class TView>
-  inline void Interaction::doInteraction(TView& view, Code const projectileId,
-                                         Code const targetId,
-                                         FourMomentum const& projectileP4,
-                                         FourMomentum const& targetP4) {
-
-    CORSIKA_LOG_DEBUG(
-        "Pythia::Interaction: "
-        "doInteraction: {} interaction? ",
-        projectileId, corsika::pythia8::Interaction::canInteract(projectileId));
-
-    // define system
-    auto const sqrtS2 = (projectileP4 + targetP4).getNormSqr();
-    HEPEnergyType const sqrtS = sqrt(sqrtS2);
-    HEPEnergyType const eProjectileLab = (sqrtS2 - static_pow<2>(get_mass(projectileId)) -
-                                          static_pow<2>(get_mass(targetId))) /
-                                         (2 * get_mass(targetId));
-
-    if (!isValid(projectileId, targetId, sqrtS)) {
-      throw std::runtime_error("invalid target,projectile,energy combination.");
-    }
-
-    CORSIKA_LOG_DEBUG("Interaction: ebeam lab: {} GeV", eProjectileLab / 1_GeV);
-    CORSIKA_LOG_DEBUG(
-        "Interaction: "
-        " doInteraction: E(GeV): {}"
-        " Ecm(GeV): {}",
-        eProjectileLab / 1_GeV, sqrtS / 1_GeV);
-
-    count_++;
-
-    int const idNow = static_cast<int>(get_PDG(projectileId));
-
-    // References to the two event records. Clear main event record.
-    Pythia8::Event& eventMain = pythiaMain_.event;
-    Pythia8::Event& eventColl = pythiaColl_.event;
-
-    COMBoost const labFrameBoost{targetP4.getSpaceLikeComponents(), get_mass(targetId)};
-    auto const proj4MomLab = labFrameBoost.toCoM(projectileP4);
-    auto const& rotCS = labFrameBoost.getRotatedCS();
-    auto const pProjLab = proj4MomLab.getSpaceLikeComponents().getComponents(rotCS);
-
-    auto constexpr invGeV = 1 / 1_GeV;
-
-    Pythia8::Vec4 const pNow{pProjLab.getX() * invGeV, pProjLab.getY() * invGeV,
-                             pProjLab.getZ() * invGeV,
-                             proj4MomLab.getTimeLikeComponent() * invGeV};
-
-    // Insert incoming particle in cleared main event record.
-    int unsuccessful_iterations{0};
-    do { // retry event generation if Pythia gets stuck
-      eventMain.clear();
-      eventMain.append(90, -11, 0, 0, 1, 1, 0, 0, pNow, pNow.mCalc());
-      int const iHad = eventMain.append(idNow, 12, 0, 0, 0, 0, 0, 0, pNow,
-                                        get_mass(projectileId) * (1 / 1_GeV));
-
-      Pythia8::Vec4 const vNow{}; // production vertex
-
-      eventMain[iHad].vProd(vNow);
-
-      auto const [Anow, Znow] = std::invoke([targetId]() {
-        if (targetId == Code::Proton) {
-          return std::make_pair(1, 1);
-        } else if (targetId == Code::Neutron) {
-          return std::make_pair(1, 0);
-        } else if (is_nucleus(targetId)) {
-          return std::make_pair<int, int>(get_nucleus_A(targetId),
-                                          get_nucleus_Z(targetId));
-        } else {
-          // due to the earlier call to isValid(), we shouldn't end up here
-          // LCOV_EXCL_START
-          CORSIKA_LOG_ERROR("invalid target {}; you shouldn't have gotten this far!",
-                            targetId);
-          return std::make_pair(0, 0);
-          // LCOV_EXCL_STOP
-        }
-      });
-
-      // Set up for collisions on a nucleus.
-      int np = Znow;
-      int nn = Anow - Znow;
-      int sizeOld = 0;
-      int sizeNew = 0;
-      Pythia8::Vec4 const dirNow = pNow / pNow.pAbs();
-      Pythia8::Rndm& rndm = pythiaMain_.rndm;
-
-      double constexpr mp = get_mass(Code::Proton) / 1_GeV;
-      double const sqrtSNN_GeV = (pNow + Pythia8::Vec4{0, 0, 0, mp}).mCalc();
-
-      auto const nCollAvg = getAverageSubcollisions(
-          targetId, pythiaColl_.getSigmaTotal(idNow, 2212, sqrtSNN_GeV) * 1_mb);
-      double const probMore = 1. - 1. / nCollAvg;
-
-      // Loop over varying number of hit nucleons in target nucleus.
-      for (int iColl = 1; iColl <= Anow; ++iColl) {
-        if (iColl > 1 && rndm.flat() > probMore) break;
-
-        // Pick incoming projectile: trivial for first subcollision, else ...
-        int iProj = iHad;
-        int procType = 0;
-
-        // ... find highest-pLongitudinal particle from latest subcollision.
-        if (iColl > 1) {
-          iProj = 0;
-          double pMax = 0.;
-          for (int i = sizeOld; i < sizeNew; ++i)
-            if (eventMain[i].isFinal() && eventMain[i].isHadron()) {
-              if (double const pp = Pythia8::dot3(dirNow, eventMain[i].p()); pp > pMax) {
-                iProj = i;
-                pMax = pp;
-              }
-            }
-
-          // No further subcollision if no particle with enough energy.
-          // cannot be reliably provoked in tests
-          // LCOV_EXCL_START
-          if (iProj == 0 ||
-              eventMain[iProj].e() - eventMain[iProj].m() < eKinMinLab_ / 1_GeV)
-            break;
-          // LCOV_EXCL_STOP
-
-          // Choose process; only SD or ND at perturbative energies.
-          double const eCMSub =
-              (eventMain[iProj].p() + Pythia8::Vec4{0, 0, 0, mp}).mCalc();
-          if (eCMSub > 10.) procType = (rndm.flat() < probSD_) ? 4 : 1;
-        }
-
-        // Pick one p or n from target.
-        int const idProj = eventMain[iProj].id();
-        bool const doProton = rndm.flat() < (np / double(np + nn));
-        if (doProton)
-          np--;
-        else
-          nn--;
-        int const idNuc = doProton ? 2212 : 2112;
-
-        // Perform the projectile-nucleon subcollision.
-        pythiaColl_.setBeamIDs(idProj, idNuc);
-        pythiaColl_.setKinematics(eventMain[iProj].p(), Pythia8::Vec4());
-
-        if (!pythiaColl_.next(procType)) {
-          CORSIKA_LOG_WARN("Pythia collision next() failed {} {} {}",
-                           eventMain[iProj].p(), idProj, idNuc);
-          eventMain.clear();
-          ++unsuccessful_iterations;
-          goto event_repeat; // retry event, last remaining good use-case for goto
-        }
-
-        // Insert target nucleon. Mothers are (0,iProj) to mark who it
-        // interacted with. Always use proton mass for simplicity.
-        int const statusNuc = (iColl == 1) ? -181 : -182;
-        int const iNuc =
-            eventMain.append(idNuc, statusNuc, 0, iProj, 0, 0, 0, 0, 0., 0., 0., mp, mp);
-        eventMain[iNuc].vProdAdd(vNow);
-
-        // Update full energy of the event with the proton mass.
-        eventMain[0].e(eventMain[0].e() + mp);
-        eventMain[0].m(eventMain[0].p().mCalc());
-
-        // Insert secondary produced particles (but skip intermediate partons)
-        // into main event record and shift to correct production vertex.
-        sizeOld = eventMain.size();
-        for (int iSub = 3; iSub < eventColl.size(); ++iSub) {
-          if (!eventColl[iSub].isFinal()) continue;
-          int const iNew = eventMain.append(eventColl[iSub]);
-          eventMain[iNew].mothers(iNuc, iProj);
-          eventMain[iNew].vProdAdd(vNow);
-        }
-        sizeNew = eventMain.size();
-
-        // Update daughters of colliding hadrons and other history.
-        eventMain[iProj].daughters(sizeOld, sizeNew - 1);
-        eventMain[iNuc].daughters(sizeOld, sizeNew - 1);
-        eventMain[iProj].statusNeg();
-        double dTau = (iColl == 1) ? (vNow.e() - eventMain[iHad].tProd()) *
-                                         eventMain[iHad].m() / eventMain[iHad].e()
-                                   : 0.;
-        eventMain[iProj].tau(dTau);
-
-        // End of loop over interactions in a nucleus.
-      }
-      break; //
-
-    event_repeat:;
-    } while (unsuccessful_iterations < 100);
-
-    if (unsuccessful_iterations >= 100) {
-      CORSIKA_LOG_CRITICAL(
-          "Pythia event generation failed after 100 trials: projectile: {}, target: {}",
-          projectileId, targetId);
-      throw std::runtime_error{"Pythia event generation failed after 100 trials"};
-    }
-
-    MomentumVector Plab_final{labFrameBoost.getOriginalCS()};
-    auto Elab_final = HEPEnergyType::zero();
-
-    for (Pythia8::Particle const& p8p : eventMain) {
-      // skip particles that have decayed / are initial particles in pythia's event record
-      if (!p8p.isFinal()) continue;
-      try {
-        auto const volatile id = static_cast<PDGCode>(p8p.id());
-        auto const pyId = convert_from_PDG(id);
-
-        MomentumVector const pyPlab(
-            rotCS, {p8p.px() * 1_GeV, p8p.py() * 1_GeV, p8p.pz() * 1_GeV});
-        auto const pyP = labFrameBoost.fromCoM(FourVector{p8p.e() * 1_GeV, pyPlab});
-
-        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, pyPlab.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
-    }
-
-    eventMain.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/InteractionModel.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 <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_.settings.word("MultipartonInteractions:initFile",
+                          dataPath.native() + "/pythia8312.mpi");
+    pythia_.readString("HeavyIon:SasdMpiReuseInit = 2");
+    pythia_.settings.word("HeavyIon:SasdMpiInitFile",
+                          dataPath.native() + "/pythia8312.sasd.mpi");
+    pythia_.readString("HeavyIon:SigFitReuseInit = 2");
+    pythia_.settings.word("HeavyIon:SigFitInitFile",
+                          dataPath.native() + "/pythia8312.sigfit");
+
+    // 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());
+          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());
+    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()};
+    }
+
+    // check if config. is valid
+
+    // NOTE: here we calculate SNN (per nucleon energy) AND S (per particle energy)
+    // because isValid expects sqrtSNN as input but the tables need Elab in the hadron -
+    // nucleus frame. Pending better solution
+    auto const SNN =
+        (projectileP4 / (is_nucleus(projectileId) ? get_nucleus_A(projectileId) : 1.) +
+         targetP4 / (is_nucleus(targetId) ? get_nucleus_A(targetId) : 1.))
+            .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 in the hadron-nucleus
+    // frame) since we cannot assume the input is in the lab. frame we calculate the lab.
+    // energy from S = (projectileP4 + targetP4)**2 and Elab = S/(2. * mTarget)
+    auto const S = (projectileP4 + targetP4).getNormSqr();
+    HEPEnergyType const Elab =
+        calculate_lab_energy(S, get_mass(projectileId), get_mass(targetId));
+
+    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 {
+
+    // check if config is valid
+    auto const SNN =
+        (projectileP4 / (is_nucleus(projectileId) ? get_nucleus_A(projectileId) : 1.) +
+         targetP4 / (is_nucleus(targetId) ? get_nucleus_A(targetId) : 1.))
+            .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 in the hadron-nucleus
+    // frame) since we cannot assume the input is in the lab. frame we calculate the lab.
+    // energy from S = (projectileP4 + targetP4)**2 and Elab = S/(2. * mTarget)
+    auto const S = (projectileP4 + targetP4).getNormSqr();
+    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});
+    // HACK: cross sections were calculated assuming hadron-nucleus but pythia uses
+    // hadron-nucleon, therefore tabulation energies are overestimated by factor A!
+    HEPEnergyType const Elab =
+        calculate_lab_energy(S, get_mass(projectileId), get_mass(targetId));
+    CORSIKA_LOG_DEBUG(
+        "pythia getCrossSection: {}+{} at Elab= {} GeV, sqrtS = {} GeV, sqrtSNN = {} GeV",
+        projectileId, get_name(targetId), Elab / 1_GeV, sqrt(S) / 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/qgsjetII/InteractionModel.inl

@@ -102,6 +102,14 @@ namespace corsika::qgsjetII {
     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,

corsika/detail/modules/sibyll/InteractionModel.inl

@@ -51,6 +51,19 @@ namespace corsika::sibyll {
                                                          target4mom);
   }
 
+  inline std::tuple<CrossSectionType, CrossSectionType>
+  InteractionModel::getCrossSectionInelEla(Code projCode, Code targetCode,
+                                           FourMomentum const& proj4mom,
+                                           FourMomentum const& target4mom) const {
+    if (is_nucleus(projCode))
+      return {getNuclearInteractionModel().getCrossSection(projCode, targetCode, proj4mom,
+                                                           target4mom),
+              CrossSectionType::zero()};
+    else
+      return getHadronInteractionModel().getCrossSectionInelEla(projCode, targetCode,
+                                                                proj4mom, target4mom);
+  }
+
   template <typename TSecondaries>
   inline void InteractionModel::doInteraction(TSecondaries& view, Code projCode,
                                               Code targetCode,

corsika/framework/utility/CrossSectionTable.hpp

+/*
+ * (c) Copyright 2024 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 <algorithm>
+#include <cmath>
+#include <stdexcept>
+#include <utility>
+#include <vector>
+
+#include <boost/functional/identity.hpp>
+#include <boost/iterator/transform_iterator.hpp>
+
+#include <corsika/framework/core/PhysicalUnits.hpp>
+
+namespace corsika {
+  namespace InterpolationTransforms {
+    using Identity = boost::identity;
+
+    struct Log {
+      template <typename T>
+      auto operator()(T val) const {
+        if constexpr (is_quantity_v<T>) {
+          return std::log(val.magnitude());
+        } else {
+          return std::log(val);
+        }
+      }
+    };
+  } // namespace InterpolationTransforms
+
+  template <typename Transform>
+  class CrossSectionTable {
+  public:
+    template <typename InputIt1, typename InputIt2>
+    CrossSectionTable(InputIt1 energiesFirst, InputIt1 energiesLast,
+                      InputIt2 crosssectionsFirst) {
+      static_assert(std::is_same_v<typename InputIt1::value_type, HEPEnergyType>);
+      static_assert(std::is_same_v<typename InputIt2::value_type, CrossSectionType>);
+
+      table_.reserve(std::distance(energiesFirst, energiesLast));
+
+      auto itE = boost::make_transform_iterator(std::move(energiesFirst), transform_);
+      decltype(itE) const itEEnd =
+          boost::make_transform_iterator(std::move(energiesLast), transform_);
+
+      for (auto itXS = std::move(crosssectionsFirst); itE != itEEnd; ++itE, ++itXS) {
+        table_.emplace_back(*itE, *itXS);
+      }
+
+      std::sort(table_.begin(), table_.end(), less_datapoint);
+    }
+
+    CrossSectionType interpolate(HEPEnergyType energy) const {
+      auto const transformed_val = transform_(energy);
+      auto const lb_it =
+          std::lower_bound(table_.cbegin(), table_.cend(), transformed_val, less_x);
+      if (lb_it == table_.cbegin()) {
+        throw std::runtime_error{"CrossSectionTable: value out of bounds (lower limit)"};
+      }
+      if (lb_it == table_.cend()) {
+        throw std::runtime_error{"CrossSectionTable: value out of bounds (upper limit)"};
+      }
+
+      auto const prev_it = std::prev(lb_it);
+      auto const lambda =
+          (transformed_val - prev_it->first) / (lb_it->first - prev_it->first);
+      return lambda * lb_it->second + (1 - lambda) * prev_it->second;
+    }
+
+  private:
+    using key_type = std::invoke_result_t<Transform, HEPEnergyType>;
+    using datapoint_type = std::pair<key_type, CrossSectionType>;
+    std::vector<datapoint_type> table_;
+    Transform transform_{};
+
+    static bool less_datapoint(datapoint_type const& a, datapoint_type const& b) {
+      return a.first < b.first;
+    }
+    static bool less_x(datapoint_type const& a, typename datapoint_type::first_type b) {
+      return a.first < b;
+    }
+  };
+} // namespace corsika

corsika/modules/Pythia8.hpp

@@ -10,3 +10,18 @@
 
 #include <corsika/modules/pythia8/Decay.hpp>
 #include <corsika/modules/pythia8/NeutrinoInteraction.hpp>
+#include <corsika/modules/pythia8/InteractionModel.hpp>
+
+namespace corsika::pythia8 {
+  /**
+   * pythia8::Interaction is the process for ProcessSequence.
+   *
+   * The pythia8::InteractionModel is wrapped as an InteractionProcess here in order
+   * to provide all the functions for ProcessSequence.
+   */
+  class Interaction : public InteractionModel, public InteractionProcess<Interaction> {
+  public:
+    Interaction(std::set<Code> const& stableList = {})
+        : InteractionModel{stableList} {};
+  };
+} // namespace corsika::pythia8

corsika/modules/pythia8/Interaction.hpp → corsika/modules/pythia8/InteractionModel.hpp

@@ -9,41 +9,65 @@
 #pragma once
 
 #include <tuple>
+#include <unordered_map>
 #include <boost/filesystem/path.hpp>
 
 #include <corsika/framework/utility/CorsikaData.hpp>
 #include <corsika/framework/core/ParticleProperties.hpp>
 #include <corsika/framework/core/PhysicalUnits.hpp>
 #include <corsika/framework/random/RNGManager.hpp>
+#include <corsika/framework/utility/CrossSectionTable.hpp>
 #include <corsika/framework/process/InteractionProcess.hpp>
 #include <corsika/modules/pythia8/Pythia8.hpp>
 
 namespace corsika::pythia8 {
-  class Interaction : public InteractionProcess<Interaction> {
+
+  /**
+   * @brief Defines the interface to the PYTHIA8 interaction model. Configured for
+   * hadron-nucleus interactions using Angantyr.
+   *
+   * This is a TModel argument for InteractionProcess<TModel>.
+   */
+
+  class InteractionModel {
 
   public:
-    Interaction(
-        boost::filesystem::path const& mpiInitFile = corsika_data("Pythia/main184.mpi"),
-        bool const print_listing = false);
-    ~Interaction();
+    /**
+     * Constructs the interface for PYTHIA8
+     *
+     * @param stableParticles - set of particle ids to be treated as stable inside
+     * pythia. if empty then all particles are treated as stable (no decays!)
+     * @param dataPath path to the pythia tables
+     * @param printListing switch on/off the pythia printout
+     */
+    InteractionModel(std::set<Code> const& stableParticles = {},
+                     boost::filesystem::path const& dataPath = corsika_data("Pythia"),
+                     bool const printListing = false);
+    ~InteractionModel();
 
     bool canInteract(Code const) const;
 
+    /**
+     * Check if configuration of projectile and target is valid.
+     *
+     * @param projectileId is the Code of the projectile
+     * @param targetId is the Code of the target
+     * @param sqrtS is the squared sum of the 4-momenta of projectile and target
+     */
     bool isValid(Code const projectileId, Code const targetId,
                  HEPEnergyType const sqrtS) const;
+
     /**
      * Returns inelastic AND elastic cross sections.
      *
      * These cross sections must correspond to the process described in doInteraction
-     * AND elastic scattering (sigma_tot = sigma_inel + sigma_el). Allowed targets are:
-     * nuclei or single nucleons (p,n,hydrogen). This "InelEla" method is used since
-     * Sibyll must be useful inside the NuclearInteraction model, which requires that.
+     * (sigma_tot = sigma_inel + sigma_el). Allowed targets are:
+     * nuclei or single nucleons (p,n,hydrogen).
      *
      * @param projectile is the Code of the projectile
      * @param target is the Code of the target
-     * @param sqrtSnn is the center-of-mass energy (per nucleon pair)
-     * @param Aprojectil is the mass number of the projectils, if it is a nucleus
-     * @param Atarget is the mass number of the target, if it is a nucleus
+     * @param projectileP4 is the 4-momentum of the projectile
+     * @param targetP4 is the 4-momentum of the target
      *
      * @return a tuple of: inelastic cross section, elastic cross section
      */
@@ -59,51 +83,72 @@ namespace corsika::pythia8 {
      *
      * @param projectile is the Code of the projectile
      * @param target is the Code of the target
-     * @param sqrtSnn is the center-of-mass energy (per nucleon pair)
-     * @param Aprojectil is the mass number of the projectils, if it is a nucleus
-     * @param Atarget is the mass number of the target, if it is a nucleus
+     * @param projectileP4 is the 4-momentum of the projectile
+     * @param targetP4 is the 4-momentum of the target
      *
      * @return inelastic cross section
      * elastic cross section
      */
-    CrossSectionType getCrossSection(
-        Code const projectile, Code const target, FourMomentum const& projectileP4,
-        FourMomentum const& targetP4) const; // non-const for now
+    CrossSectionType getCrossSection(Code const projectile, Code const target,
+                                     FourMomentum const& projectileP4,
+                                     FourMomentum const& targetP4) const;
 
     /**
-     * In this function PYTHIA is called to produce one event. The
-     * event is copied (and boosted) into the shower lab frame.
+     * In this function PYTHIA8 is called to produce one event. The
+     * event is copied (and boosted) into the reference frame defined by the input
+     * 4-momenta.
      */
     template <typename TView>
     void doInteraction(TView& output, Code const projectileId, Code const targetId,
                        FourMomentum const& projectileP4, FourMomentum const& targetP4);
 
+    using key_type = std::pair<corsika::Code, corsika::Code>;
+
+  private:
+    static auto constexpr GeV_mult = [](float x) { return x * 1_GeV; };
+    static auto constexpr millibarn_mult = [](float x) { return x * 1_mb; };
+
+    static CrossSectionTable<InterpolationTransforms::Log> loadPPTable(
+        boost::filesystem::path const& dataPath, char const* key);
+
+    default_prng_type& RNG_ = RNGManager<>::getInstance().getRandomStream("pythia");
+
     /**
-     * return average number of sub-collisions in a nucleus, using the parameterizations
-     * of Sjöstrand and Utheim, EPJ C 82 (2022) 1, 21, arXiv:2108.03481 [hep-ph]
-     *
-     * @param targetId target (nucleus)
-     * @param sigTot projectile-nucleon (=proton) cross-secion
+     * All particles that we enable as projectiles. Code::Nucleus is not listed as it is
+     * handled separately. The list is unlikely to be exhaustive, but should cover the
+     * mainstream use-cases sufficiently well.
      */
-    double getAverageSubcollisions(Code targetId, CrossSectionType sigTot) const;
 
-    static std::array constexpr validTargets_{Code::Oxygen, Code::Nitrogen,
-                                              Code::Argon,  Code::Hydrogen,
-                                              Code::Proton, Code::Neutron};
+    static std::array constexpr validProjectiles_ = {
+        Code::PiPlus,       Code::PiMinus,   Code::Pi0,           Code::Proton,
+        Code::AntiProton,   Code::Neutron,   Code::AntiNeutron,   Code::KPlus,
+        Code::KMinus,       Code::K0Long,    Code::K0Short,       Code::SigmaMinus,
+        Code::SigmaPlusBar, Code::SigmaPlus, Code::SigmaMinusBar, Code::Xi0,
+        Code::Xi0Bar};
+
+    static std::array constexpr validTargets_ = {
+        Code::Proton, Code::Carbon, Code::Nitrogen, Code::Oxygen, Code::Argon};
+
+    std::unordered_map<corsika::Code, corsika::Code> const xs_map_ = {
+        {Code::SigmaMinus, Code::SigmaPlus},
+        {Code::SigmaMinusBar, Code::SigmaPlus},
+        {Code::SigmaPlusBar, Code::SigmaPlus},
+        {Code::Xi0Bar, Code::Xi0}};
+
+    Pythia8::Pythia pythia_;
+
+    std::unordered_map<key_type, CrossSectionTable<InterpolationTransforms::Log>>
+        crossSectionTables_;
+
+    CrossSectionTable<InterpolationTransforms::Log> crossSectionPPElastic_,
+        crossSectionPPInelastic_;
 
-  private:
-    default_prng_type& RNG_ = RNGManager<>::getInstance().getRandomStream("pythia");
-    //~ bool const internalDecays_ = true;
-    int count_ = 0;
     bool const print_listing_ = false;
-    Pythia8::Pythia pythiaMain_;
-    Pythia8::Pythia mutable pythiaColl_; // Pythia's getSigma...() are not marked const...
-    double const probSD_ =
-        0.3; // Fraction of single diffractive events beyond first collision in nucleus.
-    HEPEnergyType const eMaxLab_ = 1e18_eV;
-    HEPEnergyType const eKinMinLab_ = 0.2_GeV;
+    HEPEnergyType const eMaxLab_ =
+        1e21_eV; // Cross-section tables tabulated up to 10^21 eV
+    HEPEnergyType const eKinMinLab_ = 100_GeV;
   };
 
 } // namespace corsika::pythia8
 
-#include <corsika/detail/modules/pythia8/Interaction.inl>
+#include <corsika/detail/modules/pythia8/InteractionModel.inl>

corsika/modules/qgsjetII/InteractionModel.hpp

@@ -55,6 +55,23 @@ namespace corsika::qgsjetII {
                                      FourMomentum const& projectileP4,
                                      FourMomentum const& targetP4) const;
 
+    /**
+     * Returns inelastic AND elastic cross sections.
+     *
+     * These cross sections must correspond to the process described in doInteraction
+     * AND elastic scattering (sigma_tot = sigma_inel + sigma_el).
+     *
+     * @param projectile is the Code of the projectile
+     * @param target is the Code of the target
+     * @param projectileP4: four-momentum of projectile
+     * @param targetP4: four-momentum of target
+     *
+     * @return a tuple of: inelastic cross section, elastic cross section
+     */
+    std::tuple<CrossSectionType, CrossSectionType> getCrossSectionInelEla(
+        Code const projectile, Code const target, FourMomentum const& projectileP4,
+        FourMomentum const& targetP4) const;
+
     /**
      * In this function QGSJETII is called to produce one event.
      *

corsika/modules/sibyll/InteractionModel.hpp

@@ -30,6 +30,9 @@ namespace corsika::sibyll {
     CrossSectionType getCrossSection(Code, Code, FourMomentum const&,
                                      FourMomentum const&) const;
 
+    std::tuple<CrossSectionType, CrossSectionType> getCrossSectionInelEla(
+        Code, Code, FourMomentum const&, FourMomentum const&) const;
+
     template <typename TSecondaries>
     void doInteraction(TSecondaries&, Code, Code, FourMomentum const&,
                        FourMomentum const&);

examples/CMakeLists.txt

@@ -96,3 +96,11 @@ CORSIKA_REGISTER_EXAMPLE (synchrotron_test_C8tracking)
 add_executable (synchrotron_test_manual_tracking physics_examples/synchrotron_test_manual_tracking.cpp)
 target_link_libraries (synchrotron_test_manual_tracking PRIVATE CORSIKA8::CORSIKA8)
 CORSIKA_REGISTER_EXAMPLE (synchrotron_test_manual_tracking)
+
+add_executable (cross_section physics_examples/cross_section.cpp)
+target_link_libraries (cross_section PRIVATE CORSIKA8::CORSIKA8)
+CORSIKA_REGISTER_EXAMPLE (cross_section RUN_OPTIONS sibyll)
+
+add_executable (had_interactions physics_examples/had_interactions.cpp)
+target_link_libraries (had_interactions PRIVATE CORSIKA8::CORSIKA8)
+CORSIKA_REGISTER_EXAMPLE (had_interactions RUN_OPTIONS sibyll)
\ No newline at end of file

examples/physics_examples/cross_section.cpp

+/*
+ * (c) Copyright 2019 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/Epos.hpp>
+#include <corsika/modules/Sibyll.hpp>
+#include <corsika/modules/QGSJetII.hpp>
+#include <corsika/modules/Pythia8.hpp>
+
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+
+#include <corsika/framework/geometry/RootCoordinateSystem.hpp>
+#include <corsika/framework/geometry/Vector.hpp>
+
+#include <corsika/framework/random/RNGManager.hpp>
+
+#include <corsika/setup/SetupC7trackedParticles.hpp>
+#include <corsika/modules/Random.hpp>
+
+#include <fstream>
+#include <string>
+#include <tuple>
+
+using namespace corsika;
+
+/**
+ * Calculates the inelastic p-p cross section and the production cross section for
+ * p-Oxygen for a given interaction model using the C8 interface getCrossSection(ID1, ID2,
+ * P41, P42) where ID? and P4? are the particle Id and 4-momenta of projectile and target
+ * respectively. The input to getCrossSection can be ANY frame of reference. Internally
+ * the boost to the CM is calculated and the correct cross section is returned.
+ * calculate_cross_sections allows to switch between the CM and lab configurations.
+ *
+ * output is a table in ACSCII file cross-sections-<model-name>.txt
+ *
+ * Energy ranges from 100GeV to 10^12 GeV (10^11 eV to 10^21 eV) in the lab. frame
+ * table format is: "# i, E_lab (GeV), E_cm per nucleon (GeV), cross section p-p (mb),
+ * production "
+ *
+ * @param model the interaction model to use
+ * @param model_name string with model name
+ *
+ * @return a tuple of: inelastic cross section, elastic cross section
+ */
+
+template <typename TModel>
+void calculate_cross_sections(TModel& model, std::string const& model_name) {
+
+  std::stringstream fname;
+  fname << "cross-sections-" << model_name << ".txt";
+  std::ofstream out(fname.str());
+  out << "# cross section table for " << model_name << "\n";
+  out << "# i, P_lab (GeV), E_cm per nucleon (GeV), cross section p-p (mb), cross "
+         "section pi-p (mb), production "
+         "cross section p-O16 (mb), production cross section p-N14 (mb), production "
+         "cross section p-Ar40 (mb), cross section pi-O16 (mb), production cross "
+         "section, pi-N14 (mb), production, cross section pi-Ar40 (mb), cross section "
+         "He-O16 (mb), production cross section, He-N14 (mb), production, cross section "
+         "He-Ar40 (mb)\n";
+
+  CoordinateSystemPtr const& cs = get_root_CoordinateSystem();
+
+  float const amin = 2;
+  float const amax = 11;
+  int const nebins = 15;
+  float const da = (amax - amin) / (float)nebins;
+
+  for (int i = 0; i < nebins; ++i) {
+
+    corsika::units::si::HEPMomentumType const setP = pow(10, da * i + amin) * 1_GeV;
+    CrossSectionType xs_prod_hNuc[3][3] = {};
+    CrossSectionType xs_prod_hp[3] = {};
+    HEPEnergyType setE, comEnn;
+    int l = -1;
+    for (auto projId : {Code::Proton, Code::PiPlus, Code::Helium}) {
+      ++l;
+      setE = calculate_total_energy(setP, get_mass(projId));
+      comEnn = corsika::calculate_com_energy(setE, get_mass(projId),
+                                             corsika::constants::nucleonMass);
+
+      // four momenta in lab frame
+      corsika::FourMomentum p4Proj(setE, MomentumVector(cs, {0_eV, 0_eV, setP}));
+      corsika::FourMomentum p4protonTarg(Proton::mass,
+                                         MomentumVector(cs, {0_eV, 0_eV, 0_eV}));
+
+      // had-p
+      auto const [xs_prod, xs_ela] =
+          model->getCrossSectionInelEla(projId, Code::Proton, p4Proj, p4protonTarg);
+      xs_prod_hp[l] = xs_prod;
+
+      int k = -1;
+      for (auto targNuc : {Code::Oxygen, Code::Nitrogen, Code::Argon}) {
+        ++k;
+        FourMomentum const p4nucTarg = corsika::FourMomentum(
+            get_mass(targNuc), MomentumVector(cs, {0_eV, 0_eV, 0_eV}));
+        xs_prod_hNuc[l][k] = model->getCrossSection(projId, targNuc, p4Proj, p4nucTarg);
+      }
+    }
+    CORSIKA_LOG_INFO(
+        "Elab={:15.2f} GeV,Ecom={:15.2f} GeV, sig-p-p={:6.2f} mb, sig-pi-p={:6.2f} mb, "
+        "p-O={:6.2f} mb, p-N={:6.2f} mb, p-Ar={:6.2f} mb "
+        "pi-O={:6.2f} mb, pi-N={:6.2f} mb, pi-Ar={:6.2f} mb "
+        "He-O={:6.2f} mb, He-N={:6.2f} mb, He-Ar={:6.2f} mb",
+        setP / 1_GeV, comEnn / 1_GeV, xs_prod_hp[0] / 1_mb, xs_prod_hp[1] / 1_mb,
+        xs_prod_hNuc[0][0] / 1_mb, xs_prod_hNuc[0][1] / 1_mb, xs_prod_hNuc[0][2] / 1_mb,
+        xs_prod_hNuc[1][0] / 1_mb, xs_prod_hNuc[1][1] / 1_mb, xs_prod_hNuc[1][2] / 1_mb,
+        xs_prod_hNuc[2][0] / 1_mb, xs_prod_hNuc[2][1] / 1_mb, xs_prod_hNuc[2][2] / 1_mb);
+
+    out << i << " " << setP / 1_GeV << " " << comEnn / 1_GeV << " "
+        << xs_prod_hp[0] / 1_mb << " " << xs_prod_hp[1] / 1_mb << " "
+        << xs_prod_hNuc[0][0] / 1_mb << " " << xs_prod_hNuc[0][1] / 1_mb << " "
+        << xs_prod_hNuc[0][2] / 1_mb << " " << xs_prod_hNuc[1][0] / 1_mb << " "
+        << xs_prod_hNuc[1][1] / 1_mb << " " << xs_prod_hNuc[1][2] / 1_mb << " "
+        << " " << xs_prod_hNuc[2][0] / 1_mb << " " << xs_prod_hNuc[2][1] / 1_mb << " "
+        << xs_prod_hNuc[2][2] / 1_mb << "\n";
+  }
+  out.close();
+  std::cout << "wrote cross section table for " << model_name
+            << " in file: " << fname.str() << std::endl;
+}
+
+int main(int argc, char** argv) {
+
+  logging::set_level(logging::level::info);
+
+  if (argc != 2) {
+    std::cout << "usage: check <interaction model> \n valid models are: sibyll, "
+                 "epos, qgsjet, pythia8"
+              << std::endl;
+    return 1;
+  }
+  std::string int_model_name = std::string(argv[1]);
+
+  if (int_model_name == "sibyll") {
+    RNGManager<>::getInstance().registerRandomStream("sibyll");
+    auto model = std::make_shared<corsika::sibyll::InteractionModel>(
+        std::set<Code>{Code::Proton, Code::Oxygen, Code::Nitrogen},
+        corsika::setup::C7trackedParticles);
+    calculate_cross_sections(model, int_model_name);
+  } else if (int_model_name == "pythia8") {
+    RNGManager<>::getInstance().registerRandomStream("pythia");
+    auto model = std::make_shared<corsika::pythia8::InteractionModel>();
+    calculate_cross_sections(model, int_model_name);
+
+  } else if (int_model_name == "epos") {
+    RNGManager<>::getInstance().registerRandomStream("epos");
+    auto model = std::make_shared<corsika::epos::InteractionModel>(
+        corsika::setup::C7trackedParticles);
+    calculate_cross_sections(model, int_model_name);
+  } else if (int_model_name == "qgsjet") {
+    RNGManager<>::getInstance().registerRandomStream("qgsjet");
+    auto model = std::make_shared<corsika::qgsjetII::InteractionModel>();
+    calculate_cross_sections(model, int_model_name);
+  } else {
+    std::cout << "interaction model should be: sibyll, epos, qgsjet or pythia"
+              << std::endl;
+    return 1;
+  }
+}

examples/physics_examples/had_interactions.cpp

+/*
+ * (c) Copyright 2019 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/Epos.hpp>
+#include <corsika/modules/QGSJetII.hpp>
+#include <corsika/modules/Sibyll.hpp>
+#include <corsika/modules/Pythia8.hpp>
+
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+
+#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/framework/geometry/RootCoordinateSystem.hpp>
+#include <corsika/framework/geometry/Vector.hpp>
+
+#include <corsika/framework/random/RNGManager.hpp>
+
+#include <corsika/framework/core/EnergyMomentumOperations.hpp>
+
+#include <corsika/media/Environment.hpp>
+#include <corsika/media/HomogeneousMedium.hpp>
+#include <corsika/media/MediumPropertyModel.hpp>
+#include <corsika/media/NuclearComposition.hpp>
+#include <corsika/media/ShowerAxis.hpp>
+#include <corsika/media/UniformMagneticField.hpp>
+
+#include <corsika/setup/SetupStack.hpp>
+#include <corsika/setup/SetupC7trackedParticles.hpp>
+
+#include <corsika/modules/Random.hpp>
+
+#include <fstream>
+#include <iostream>
+#include <tuple>
+
+using namespace corsika;
+
+using EnvironmentInterface = IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>;
+using EnvType = Environment<EnvironmentInterface>;
+
+template <typename TModel>
+void create_events(TModel& model, std::string const& model_name, Code const projectileId,
+                   FourMomentum projectileP4, Code const targetId, FourMomentum targetP4,
+                   int const nEvents) {
+
+  logging::set_level(logging::level::info);
+
+  // calculate Ecm per nucleon (sqrtSNN). Used only for output!
+  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);
+  auto const sqrtS = (projectileP4 + targetP4).getNorm();
+  CORSIKA_LOG_INFO("{} + {} interactions at sqrtSNN= {} GeV, sqrtS = {} GeV",
+                   projectileId, targetId, sqrtSNN / 1_GeV, sqrtS / 1_GeV);
+
+  // target material and environment (not used just there to create the C8 stack)
+  EnvType env;
+  auto& universe = *(env.getUniverse());
+  CoordinateSystemPtr const& rootCS = env.getCoordinateSystem();
+  auto world = EnvType::createNode<Sphere>(Point{rootCS, 0_m, 0_m, 0_m}, 150_km);
+  using MyHomogeneousModel =
+      MediumPropertyModel<UniformMagneticField<HomogeneousMedium<EnvironmentInterface>>>;
+  auto const props = world->setModelProperties<MyHomogeneousModel>(
+      Medium::AirDry1Atm, MagneticFieldVector(rootCS, 0_T, 0_T, 0_T),
+      1_kg / (1_m * 1_m * 1_m), NuclearComposition({targetId}, {1.}));
+  world->setModelProperties(props);
+  universe.addChild(std::move(world));
+
+  // projectile (dummy)
+  auto const plab3 = projectileP4.getSpaceLikeComponents();
+  auto p0 = Point(rootCS, 0_m, 0_m, 0_m);
+  setup::Stack<EnvType> stack;
+  stack.addParticle(std::make_tuple(
+      projectileId, calculate_kinetic_energy(plab3.getNorm(), get_mass(projectileId)),
+      plab3.normalized(), p0, 0_ns));
+  auto particle = stack.first();
+  particle.setNode(universe.getContainingNode(p0));
+
+  std::stringstream pname;
+  pname << "particles-" << model_name << "-" << projectileId << "-" << targetId << "-"
+        << sqrtSNN / 1_GeV << "GeV"
+        << ".txt";
+
+  unsigned int count_central_charged = 0;
+  HEPMomentumType total_pt = 0_GeV;
+
+  std::ofstream pout(pname.str());
+  // add header
+  pout << "# particle production in " << model_name << " for " << projectileId << " + "
+       << targetId << " interactions."
+       << "\n";
+  pout << "# lab. momentum of projectile: "
+       << projectileP4.getSpaceLikeComponents().getNorm() / 1_GeV << " (GeV)"
+       << "\n";
+  pout << "# CM energy per nucleon: " << sqrtSNN / 1_GeV << " (GeV)"
+       << "\n";
+  pout << "# i, j, particle name, PDG id, p.rap, rap., energy (GeV), px, py, pz (GeV), "
+          "pt (GeV), charge"
+       << "\n";
+  for (int i = 0; i < nEvents; ++i) {
+
+    // empty stack from previous interaction
+    stack.clear();
+    // add particle to get StackView to write secondaries to C8 stack
+    stack.addParticle(std::make_tuple(
+        projectileId, calculate_kinetic_energy(plab3.getNorm(), get_mass(projectileId)),
+        plab3.normalized(), p0, 0_ns)); // irrelevant
+    auto particle = stack.first();
+    particle.setNode(universe.getContainingNode(p0));
+    setup::Stack<EnvType>::stack_view_type output(particle);
+
+    // generate secondaries and fill output
+    model->doInteraction(output, projectileId, targetId, projectileP4, targetP4);
+
+    // loop through secondaries and calculate higher level variables (rapidity, pt, ...),
+    // write to file
+    int j = -1;
+    HEPMomentumType tpt = 0_GeV;
+    for (auto& secondary : output) {
+      ++j;
+      // pseudorapidity, rapidity
+      auto const p = secondary.getMomentum();
+      auto const pz = p.getZ(rootCS);
+      auto const px = p.getX(rootCS);
+      auto const py = p.getY(rootCS);
+      auto const pabs = p.getNorm();
+      auto const pt = sqrt(p.getSquaredNorm() - square(pz));
+      auto const en = secondary.getEnergy();
+      auto const mt = sqrt(square(pt) + square(get_mass(secondary.getPID())));
+      double prap, rap;
+      if (pz / 1_GeV <= 0) {
+        rap = log(mt / (en - pz));
+        prap = log((pt + 1.0_eV) / (pabs - pz));
+      } else {
+        rap = log((en + pz) / mt);
+        prap = log((pabs + pz) / (pt + 1.0_eV));
+      }
+      pout << i << " " << j << " " << secondary.getPID() << " "
+           << static_cast<int>(get_PDG(secondary.getPID())) << " " << prap << " " << rap
+           << " " << en / 1_GeV << " " << px / 1_GeV << " " << py / 1_GeV << " "
+           << pz / 1_GeV << " " << pt / 1_GeV << " "
+           << get_charge_number(secondary.getPID()) << "\n";
+      // calculate plateau (CMS experiment)
+      if (abs(prap) < 2.5 && get_charge_number(secondary.getPID()) != 0) {
+        count_central_charged++;
+      }
+      tpt += pt;
+    }
+    total_pt += tpt / float(j);
+  }
+  std::cout << "average plateau height: " << count_central_charged / float(nEvents)
+            << std::endl;
+  std::cout << "average pt (GeV): " << total_pt / float(nEvents) / 1_GeV << std::endl;
+  std::cout << "wrote hadronic events for " << model_name << " in file: " << pname.str()
+            << std::endl;
+};
+
+int main(int argc, char** argv) {
+
+  std::string int_model, projectile, target;
+  HEPMomentumType P0;
+  int Nevents = 10;
+  if (argc < 2 || argc > 7) {
+    std::cout << "usage: check <interaction model> \n valid models are: sibyll, "
+                 "epos, qgsjet, pythia8. Other options are (in that order): <projectile> "
+                 "<target> <projectile momentum (lab in GeV)> <no. of events>. Defaults "
+                 "are: proton, proton, 100, 10"
+              << std::endl;
+    return 1;
+  } else if (argc == 2) {
+    int_model = std::string(argv[1]);
+    projectile = "proton";
+    target = "proton";
+    P0 = 100_GeV;
+  } else if (argc == 3) {
+    int_model = std::string(argv[1]);
+    projectile = std::string(argv[2]);
+    target = "proton";
+    P0 = 100_GeV;
+  } else if (argc == 4) {
+    int_model = std::string(argv[1]);
+    projectile = std::string(argv[2]);
+    target = std::string(argv[3]);
+    P0 = 100_GeV;
+  } else if (argc == 5) {
+    int_model = std::string(argv[1]);
+    projectile = std::string(argv[2]);
+    target = std::string(argv[3]);
+    P0 = std::stoi(std::string(argv[4])) * 1_GeV;
+  } else if (argc == 6) {
+    int_model = std::string(argv[1]);
+    projectile = std::string(argv[2]);
+    target = std::string(argv[3]);
+    P0 = std::stoi(std::string(argv[4])) * 1_GeV;
+    Nevents = std::stoi(std::string(argv[5]));
+  }
+
+  Code projectileId;
+  if (projectile == "proton")
+    projectileId = Code::Proton;
+  else if (projectile == "antiproton")
+    projectileId = Code::AntiProton;
+  else if (projectile == "pion")
+    projectileId = Code::PiPlus;
+  else if (projectile == "oxygen")
+    projectileId = Code::Oxygen;
+  else {
+    std::cout << "unknown projectile. pick proton, antiproton, pion or oxygen"
+              << std::endl;
+    return 0;
+  }
+  Code targetId;
+  if (target == "proton")
+    targetId = Code::Proton;
+  else if (target == "nitrogen")
+    targetId = Code::Nitrogen;
+  else {
+    std::cout << "unknown target. pick proton or nitrogen" << std::endl;
+    return 0;
+  }
+
+  auto const rootCS = get_root_CoordinateSystem();
+  FourMomentum const P4proj{
+      calculate_total_energy(get_mass(projectileId), P0),
+      {rootCS, {HEPMomentumType::zero(), HEPMomentumType::zero(), P0}}};
+
+  FourMomentum const P4targ{
+      calculate_total_energy(get_mass(targetId), HEPMomentumType::zero()),
+      {rootCS,
+       {HEPMomentumType::zero(), HEPMomentumType::zero(), HEPMomentumType::zero()}}};
+
+  /*
+   the interface to hadronic models in C8 has two layers. The class InteractionModel
+   (HadronInteractionModel in the case of sibyll) implements the interface to the
+   interaction model. It provides the basic interface functions getCrossSection and
+   doInteraction. The class Interaction extends the interface in InteractionModel to an
+   element C8 process sequence. Here we use only the interface class.
+  */
+  if (int_model == "sibyll") {
+    RNGManager<>::getInstance().registerRandomStream("sibyll");
+    // this is the hadron model in sibyll. can only do hadron-nucleus interactions
+    auto model = std::make_shared<corsika::sibyll::HadronInteractionModel>(
+        corsika::setup::C7trackedParticles);
+    create_events(model, int_model, projectileId, P4proj, targetId, P4targ, Nevents);
+  } else if (int_model == "epos") {
+    RNGManager<>::getInstance().registerRandomStream("epos");
+    auto model = std::make_shared<corsika::epos::InteractionModel>(
+        corsika::setup::C7trackedParticles);
+    create_events(model, int_model, projectileId, P4proj, targetId, P4targ, Nevents);
+  } else if (int_model == "pythia8") {
+    RNGManager<>::getInstance().registerRandomStream("pythia");
+    auto model = std::make_shared<corsika::pythia8::InteractionModel>(
+        corsika::setup::C7trackedParticles);
+
+    create_events(model, int_model, projectileId, P4proj, targetId, P4targ, Nevents);
+  } else {
+    RNGManager<>::getInstance().registerRandomStream("qgsjet");
+    auto model = std::make_shared<corsika::qgsjetII::InteractionModel>();
+    create_events(model, int_model, projectileId, P4proj, targetId, P4targ, Nevents);
+  }
+}

data @ 39446571

-Subproject commit 252a0fe2368cbbd47206f9459a7e5a3657aacbf7
+Subproject commit 39446571780bbe0eccfd38efa32856552a2b64c6

tests/modules/testPythia8.cpp

@@ -199,4 +199,6 @@ TEST_CASE("Pythia8Interface", "modules") {
                                       {Proton::mass, {rootCS, {0_eV, 0_eV, 0_eV}}}) ==
             CrossSectionType::zero());
   }
+
+#include <tests/modules/testPythia8Interaction.inl>
 }

tests/modules/testPythia8Interaction.inl

@@ -13,14 +13,23 @@
   when this is merged the interaction tests can be activated again by including this file
   in testPytha8.cpp eg #include "tests/modules/testPythia8Interaction.inl"
 */
-SECTION("pythia interaction") {
+#include "corsika/framework/core/EnergyMomentumOperations.hpp"
+#include "corsika/framework/core/Logging.hpp"
+#include "corsika/framework/core/PhysicalUnits.hpp"
 
-  // this will be a p-p collision at sqrts=3.5TeV -> no problem for pythia
-  auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
-      Code::Proton, 7_TeV, (DummyEnvironment::BaseNodeType* const)nodePtr, *csPtr);
-  auto& view = *secViewPtr;
+logging::set_level(logging::level::debug);
 
-  corsika::pythia8::Interaction collision;
+// this will be a p-p collision at sqrts=3.5TeV -> no problem for pythia
+auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
+    Code::Proton, 7_TeV, (DummyEnvironment::BaseNodeType* const)nodePtr, *csPtr);
+auto& view = *secViewPtr;
+
+std::set<Code> const tracked_hadrons = {Code::Proton, Code::Neutron, Code::Pi0,
+                                        Code::PiPlus};
+
+corsika::pythia8::InteractionModel collision(tracked_hadrons);
+
+SECTION("pythia interaction configurations") {
 
   REQUIRE(collision.canInteract(Code::Proton));
   REQUIRE(collision.canInteract(Code::AntiProton));
@@ -29,37 +38,46 @@ SECTION("pythia interaction") {
   REQUIRE(collision.canInteract(Code::PiMinus));
   REQUIRE(collision.canInteract(Code::PiPlus));
   REQUIRE_FALSE(collision.canInteract(Code::Electron));
+  REQUIRE_FALSE(collision.canInteract(Code::MuPlus));
+}
 
-  // pi+p
-  REQUIRE(collision.getCrossSection(
-              Code::PiPlus, Code::Proton,
-              {sqrt(static_pow<2>(PiPlus::mass) + static_pow<2>(100_GeV)),
-               {rootCS, {0_eV, 0_eV, 100_GeV}}},
-              {Proton::mass, {rootCS, {0_eV, 0_eV, 0_eV}}}) > 0_mb);
+corsika::units::si::HEPMomentumType P0 = 10_TeV;
 
-  // pi+H
-  REQUIRE(collision.getCrossSection(
-              Code::PiPlus, Code::Hydrogen,
-              {sqrt(static_pow<2>(PiPlus::mass) + static_pow<2>(100_GeV)),
-               {rootCS, {0_eV, 0_eV, 100_GeV}}},
-              {Hydrogen::mass, {rootCS, {0_eV, 0_eV, 0_eV}}}) > 0_mb);
-
-  // K+{p,n,N,O,Ar}
-  Code const target =
-      GENERATE(Code::Proton, Code::Neutron, Code::Nitrogen, Code::Oxygen, Code::Argon);
-  REQUIRE(collision.getCrossSection(
-              Code::KPlus, target,
-              {sqrt(static_pow<2>(KPlus::mass) + static_pow<2>(100_GeV)),
-               {rootCS, {0_eV, 0_eV, 100_GeV}}},
-              {get_mass(target), {rootCS, {0_eV, 0_eV, 0_eV}}}) > 0_mb);
-
-  collision.doInteraction(view, Code::Proton, target,
-                          {sqrt(static_pow<2>(Proton::mass) + static_pow<2>(100_GeV)),
-                           {rootCS, {0_eV, 0_eV, 100_GeV}}},
+SECTION("pythia interaction") {
+  // test some combinations of valid target and projectile particles
+  // so far only hadron-hadron and hadron-Nucleus is allowed
+  Code const target = GENERATE(Code::Nitrogen, Code::Oxygen, Code::Argon);
+  Code const projectile = GENERATE(Code::Proton, Code::PiPlus, Code::KPlus);
+
+  CORSIKA_LOG_INFO("testing: {} - {}", projectile, target);
+  REQUIRE(
+      collision.getCrossSection(
+          projectile, target,
+          {calculate_total_energy(P0, get_mass(projectile)), {rootCS, {0_eV, 0_eV, P0}}},
+          {get_mass(target), {rootCS, {0_eV, 0_eV, 0_eV}}}) > 0_mb);
+
+  collision.doInteraction(view, projectile, target,
+                          {corsika::calculate_total_energy(P0, get_mass(projectile)),
+                           {rootCS, {0_eV, 0_eV, P0}}},
                           {get_mass(target), {rootCS, {0_eV, 0_eV, 0_eV}}});
   REQUIRE(view.getSize() >= 2);
 }
 
+SECTION("pythia interaction: angantyr fail") {
+  // when initializing pythia with angantyr, kaon-proton or pion-proton interactions do
+  // not seem to work
+  Code const target = Code::Proton;
+  Code const projectile = GENERATE(Code::KPlus, Code::PiPlus);
+
+  CORSIKA_LOG_INFO("testing: {}-{}", projectile, target);
+
+  REQUIRE_THROWS(
+      collision.doInteraction(view, projectile, target,
+                              {corsika::calculate_total_energy(P0, get_mass(projectile)),
+                               {rootCS, {0_eV, 0_eV, P0}}},
+                              {get_mass(target), {rootCS, {0_eV, 0_eV, 0_eV}}}));
+}
+
 SECTION("pythia too low energy") {
 
   // this is a projectile neutron with very little energy
@@ -124,26 +142,28 @@ SECTION("pythia wrong projectile") {
       Code::Iron, 1_GeV, (DummyEnvironment::BaseNodeType* const)nodePtr, *csPtr);
   { [[maybe_unused]] auto const& dummy_StackPtr = stackPtr; }
 
+  HEPMomentumType const P0 = 100_GeV;
+
   corsika::pythia8::Interaction collision;
   REQUIRE(collision.getCrossSectionInelEla(
               Code::Electron, Code::Electron,
-              {sqrt(static_pow<2>(Electron::mass) + static_pow<2>(100_GeV)),
-               {rootCS, {0_eV, 0_eV, 100_GeV}}},
+              {sqrt(static_pow<2>(Electron::mass) + static_pow<2>(P0)),
+               {rootCS, {0_eV, 0_eV, P0}}},
               {Proton::mass, {rootCS, {0_eV, 0_eV, 0_eV}}}) == std::tuple{0_mb, 0_mb});
 
-  REQUIRE_THROWS(
-      collision.doInteraction(*secViewPtr, Code::Helium, Code::Nitrogen,
-                              {sqrt(static_pow<2>(Helium::mass) + static_pow<2>(100_GeV)),
-                               {rootCS, {0_eV, 0_eV, 100_GeV}}},
-                              {Nitrogen::mass, {rootCS, {0_eV, 0_eV, 0_eV}}}));
+  REQUIRE_THROWS(collision.doInteraction(
+      *secViewPtr, Code::Helium, Code::Nitrogen,
+      {sqrt(static_pow<2>(Helium::mass) + static_pow<2>(P0)), {rootCS, {0_eV, 0_eV, P0}}},
+      {Nitrogen::mass, {rootCS, {0_eV, 0_eV, 0_eV}}}));
 
   // gamma+p not possible
   REQUIRE(collision.getCrossSection(
-              Code::Photon, Code::Proton, {100_GeV, {rootCS, {0_eV, 0_eV, 100_GeV}}},
+              Code::H0, Code::Proton,
+              {calculate_total_energy(P0, H0::mass), {rootCS, {0_eV, 0_eV, P0}}},
               {Proton::mass, {rootCS, {0_eV, 0_eV, 0_eV}}}) == CrossSectionType::zero());
 
   REQUIRE(collision.getCrossSectionInelEla(
-              Code::Photon, Code::Proton, {100_GeV, {rootCS, {0_eV, 0_eV, 100_GeV}}},
+              Code::Photon, Code::Proton, {P0, {rootCS, {0_eV, 0_eV, P0}}},
               {Proton::mass, {rootCS, {0_eV, 0_eV, 0_eV}}}) ==
           std::make_tuple(CrossSectionType::zero(), CrossSectionType::zero()));
 }
\ No newline at end of file