corsika/detail/modules/StackInspector.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
@@ -27,7 +26,9 @@ namespace corsika {
       : StackProcess<StackInspector<TStack>>(vNStep)
       , ReportStack_(vReportStack)
       , E0_(vE0)
-      , StartTime_(std::chrono::system_clock::now()) {}
+      , StartTime_(std::chrono::system_clock::now())
+      , energyPostInit_(HEPEnergyType::zero())
+      , timePostInit_(std::chrono::system_clock::now()) {}
 
   template <typename TStack>
   inline StackInspector<TStack>::~StackInspector() {}
@@ -54,50 +55,62 @@ namespace corsika {
       }
     }
 
-    auto const now = std::chrono::system_clock::now();
-    std::chrono::duration<double> const elapsed_seconds = now - StartTime_;
-    std::time_t const now_time = std::chrono::system_clock::to_time_t(now);
-    auto const dE = E0_ - Etot;
-    if (dE < dE_threshold_) return;
-    double const progress = dE / E0_;
+    if ((E0_ - Etot) < dE_threshold_) return;
 
-    std::time_t const start_time = std::chrono::system_clock::to_time_t(StartTime_);
+    // limit number of printouts
+    if (PrintoutCounter_ < MaxNumberOfPrintouts_) {
+      std::chrono::system_clock::time_point const now = std::chrono::system_clock::now();
+      std::chrono::duration<double> const elapsed_seconds = now - StartTime_; // seconds
 
-    if (progress > 0) {
+      // Select reference times and energies using either the true start
+      // or the delayed start to avoid counting overhead in ETA
+      bool const usePostVals = (energyPostInit_ != HEPEnergyType::zero());
+      auto const dE = (usePostVals ? energyPostInit_ : E0_) - Etot;
+      std::chrono::duration<double> const usedSeconds = now - timePostInit_;
+      double const dT = usedSeconds.count();
 
-      double const eta_seconds = elapsed_seconds.count() / progress;
-      std::time_t const eta_time = std::chrono::system_clock::to_time_t(
-          StartTime_ + std::chrono::seconds((int)eta_seconds));
+      double const progress = (E0_ - Etot) / E0_;
+      // for printout
+      std::time_t const now_time = std::chrono::system_clock::to_time_t(now);
 
-      int const yday0 = std::localtime(&start_time)->tm_yday;
-      int const yday1 = std::localtime(&eta_time)->tm_yday;
+      double const ETA_seconds = (1.0 - progress) * dT * (E0_ / dE);
+      std::chrono::system_clock::time_point const ETA =
+          now + std::chrono::seconds((int)ETA_seconds);
+
+      std::time_t const ETA_time = std::chrono::system_clock::to_time_t(ETA);
+
+      int const yday0 = std::localtime(&now_time)->tm_yday;
+      int const yday1 = std::localtime(&ETA_time)->tm_yday;
       int const dyday = yday1 - yday0;
 
-      CORSIKA_LOG_INFO(
+      std::stringstream ETA_string;
+      ETA_string << std::put_time(std::localtime(&ETA_time), "%T %Z");
+
+      CORSIKA_LOG_DEBUG(
           "StackInspector: "
           " time={}"
-          ", running={} seconds"
+          ", running={:.1f} seconds"
           " ( {:.1f}%)"
           ", nStep={}"
           ", stackSize={}"
-          ", Estack={} GeV"
+          ", Estack={:.1f} GeV"
           ", ETA={}{}",
-          std::put_time(std::localtime(&now_time), "%T"), elapsed_seconds.count(),
+          std::put_time(std::localtime(&now_time), "%T %Z"), elapsed_seconds.count(),
           (progress * 100), getStep(), vS.getSize(), Etot / 1_GeV,
-          (dyday == 0 ? "" : fmt::format("+{}d ", dyday)),
-          std::put_time(std::localtime(&eta_time), "%T"));
-    } else {
-      CORSIKA_LOG_INFO(
-          "StackInspector: "
-          " time={}"
-          ", running={} seconds"
-          " ( {:.1f}%)"
-          ", nStep={}"
-          ", stackSize={}"
-          ", Estack={} GeV"
-          ", ETA={}{}",
-          std::put_time(std::localtime(&now_time), "%T"), elapsed_seconds.count(),
-          (progress * 100), getStep(), vS.getSize(), Etot / 1_GeV, "---", "---");
+          (dyday == 0 ? "" : fmt::format("+{}d ", dyday)), ETA_string.str());
+
+      PrintoutCounter_++;
+
+      if (PrintoutCounter_ == MaxNumberOfPrintouts_) {
+        CORSIKA_LOG_DEBUG("StackInspector reached allowed maximum of {} lines printout",
+                          MaxNumberOfPrintouts_);
+      }
+
+      // Change reference time once the shower has begin (avoid counting overhead time)
+      if (progress > 0.02 && energyPostInit_ == HEPEnergyType::zero()) {
+        energyPostInit_ = Etot;
+        timePostInit_ = std::chrono::system_clock::now();
+      }
     }
   }
 

corsika/detail/modules/TimeCut.inl

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/modules/TimeCut.hpp>
+
+namespace corsika {
+
+  inline TimeCut::TimeCut(const TimeType time)
+      : time_(time) {}
+
+  template <typename Particle>
+  inline ProcessReturn TimeCut::doContinuous(Step<Particle> const& step, bool const) {
+    CORSIKA_LOG_TRACE("TimeCut::doContinuous");
+    if (step.getTimePost() >= time_) {
+      CORSIKA_LOG_TRACE("stopping continuous process");
+      return ProcessReturn::ParticleAbsorbed;
+    }
+    return ProcessReturn::Ok;
+  }
+
+} // namespace corsika

corsika/detail/modules/TrackWriter.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
@@ -15,41 +14,41 @@
 
 namespace corsika {
 
-  template <typename TOutputWriter>
-  inline TrackWriter<TOutputWriter>::TrackWriter() {}
+  template <typename TOutput>
+  inline TrackWriter<TOutput>::TrackWriter() {}
 
-  template <typename TOutputWriter>
-  template <typename TParticle, typename TTrack>
-  inline ProcessReturn TrackWriter<TOutputWriter>::doContinuous(TParticle const& vP,
-                                                                TTrack const& vT,
-                                                                bool const) {
+  template <typename TOutput>
+  template <typename TParticle>
+  inline ProcessReturn TrackWriter<TOutput>::doContinuous(Step<TParticle> const& step,
+                                                          bool const) {
 
-    auto const start = vT.getPosition(0).getCoordinates();
-    auto const end = vT.getPosition(1).getCoordinates();
+    auto const start = step.getPositionPre().getCoordinates();
+    auto const end = step.getPositionPost().getCoordinates();
 
     // write the track to the file
-    this->write(vP.getPID(), vP.getEnergy(), start, vP.getTime() - vT.getDuration(), end,
-                vP.getTime());
+    TOutput::write(step.getParticlePre().getPID(), step.getEkinPre(),
+                   step.getParticlePre().getWeight(), start, step.getTimePre(), end,
+                   step.getEkinPost(), step.getTimePost());
 
     return ProcessReturn::Ok;
   }
 
-  template <typename TOutputWriter>
+  template <typename TOutput>
   template <typename TParticle, typename TTrack>
-  inline LengthType TrackWriter<TOutputWriter>::getMaxStepLength(TParticle const&,
-                                                                 TTrack const&) {
+  inline LengthType TrackWriter<TOutput>::getMaxStepLength(TParticle const&,
+                                                           TTrack const&) {
     return meter * std::numeric_limits<double>::infinity();
   }
 
-  template <typename TOutputWriter>
-  YAML::Node TrackWriter<TOutputWriter>::getConfig() const {
+  template <typename TOutput>
+  YAML::Node TrackWriter<TOutput>::getConfig() const {
     using namespace units::si;
 
     YAML::Node node;
 
     // add default units for values
     node["type"] = "TrackWriter";
-    node["units"] = "GeV | m | s";
+    node["units"] = "GeV | m | ns";
 
     return node;
   }

corsika/detail/modules/conex/CONEXhybrid.inl

 /*
  * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
  *
- * This software is distributed under the terms of the 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/Logging.hpp>
+#include <corsika/media/CORSIKA7Atmospheres.hpp>
 #include <corsika/modules/conex/CONEXhybrid.hpp>
+#include <corsika/modules/conex/CONEXrandom.hpp>
+#include <corsika/modules/Random.hpp>
 #include <corsika/modules/conex/CONEX_f.hpp>
 #include <corsika/framework/random/RNGManager.hpp>
 #include <corsika/framework/core/PhysicalConstants.hpp>
@@ -17,15 +19,19 @@
 
 #include <algorithm>
 #include <fstream>
-#include <iomanip>
+#include <numeric>
 #include <utility>
 
 namespace corsika {
 
-  inline CONEXhybrid::CONEXhybrid(Point center, ShowerAxis const& showerAxis,
-                                  LengthType groundDist, LengthType injectionHeight,
-                                  HEPEnergyType primaryEnergy, PDGCode primaryPDG)
-      : center_{center}
+  template <typename TOutputE, typename TOutputN>
+  inline CONEXhybrid<TOutputE, TOutputN>::CONEXhybrid(
+      Point const& center, ShowerAxis const& showerAxis, LengthType groundDist,
+      LengthType injectionHeight, HEPEnergyType primaryEnergy, PDGCode primaryPDG,
+      TOutputE& args1, TOutputN& args2)
+      : SubWriter<TOutputE>(args1)
+      , SubWriter<TOutputN>(args2)
+      , center_{center}
       , showerAxis_{showerAxis}
       , groundDist_{groundDist}
       , injectionHeight_{injectionHeight}
@@ -57,8 +63,7 @@ namespace corsika {
 
         return b.normalized();
       })}
-      , y_sf_{showerAxis_.getDirection().cross(x_sf_)}
-      , energy_em_(0_GeV) {
+      , y_sf_{showerAxis_.getDirection().cross(x_sf_)} {
 
     CORSIKA_LOG_DEBUG("x_sf (conexObservationCS): {}",
                       x_sf_.getComponents(conexObservationCS_));
@@ -80,9 +85,31 @@ namespace corsika {
     CORSIKA_LOG_DEBUG("showerCore (C8): {}",
                       showerCore_.getCoordinates(center.getCoordinateSystem()));
 
-    int randomSeeds[3] = {1234, 0,
-                          0}; // SEEDS ARE NOT USED. All random numbers are obtained from
-                              // the CORSIKA 8 stream "conex" and "epos"!
+    auto const& components = ::corsika::standardAirComposition.getComponents();
+    auto const& fractions = ::corsika::standardAirComposition.getFractions();
+    if (::corsika::standardAirComposition.getSize() != 3) {
+      throw std::runtime_error{"CONEXhybrid only usable with standard 3-component air"};
+    }
+
+    std::transform(components.cbegin(), components.cend(), ::conex::cxair_.aira.begin(),
+                   get_nucleus_A);
+    std::transform(components.cbegin(), components.cend(), ::conex::cxair_.airz.begin(),
+                   get_nucleus_Z);
+    std::copy(fractions.cbegin(), fractions.cend(), ::conex::cxair_.airw.begin());
+
+    ::conex::cxair_.airava =
+        std::inner_product(::conex::cxair_.airw.cbegin(), ::conex::cxair_.airw.cend(),
+                           ::conex::cxair_.aira.cbegin(), 0.);
+    ::conex::cxair_.airavz =
+        std::inner_product(::conex::cxair_.airw.cbegin(), ::conex::cxair_.airw.cend(),
+                           ::conex::cxair_.airz.cbegin(), 0.);
+
+    // this is the CONEX default but actually unused there
+    ::conex::cxair_.airi = {82.0e-09, 95.0e-09, 188.e-09};
+
+    int randomSeeds[3] = {1234, 0, 0}; // SEEDS ARE NOT USED. All random numbers are
+                                       // obtained from the CORSIKA 8 stream "conex"
+    corsika::connect_random_stream("conex", ::conex::set_rng_function);
     int heModel = eSibyll23;
 
     int nShower = 1; // large to avoid final stats.
@@ -99,7 +126,8 @@ namespace corsika {
                         configPath.c_str(), configPath.size());
   }
 
-  inline void CONEXhybrid::initCascadeEquations() {
+  template <typename TOutputE, typename TOutputN>
+  inline void CONEXhybrid<TOutputE, TOutputN>::initCascadeEquations() {
 
     // set phi, theta
     Vector<length_d> ez{conexObservationCS_, {0._m, 0._m, -1_m}};
@@ -132,8 +160,9 @@ namespace corsika {
     ::conex::conexrun_(ipart, eprima, theta, phi, xminp, dimpact, ioseed.data());
   }
 
+  template <typename TOutputE, typename TOutputN>
   template <typename TStackView>
-  inline void CONEXhybrid::doSecondaries(TStackView& vS) {
+  inline void CONEXhybrid<TOutputE, TOutputN>::doSecondaries(TStackView& vS) {
     auto p = vS.begin();
     while (p != vS.end()) {
       Code const pid = p.getPID();
@@ -145,9 +174,10 @@ namespace corsika {
     }
   }
 
-  inline bool CONEXhybrid::addParticle(Code pid, HEPEnergyType energy, HEPEnergyType mass,
-                                       Point const& position,
-                                       DirectionVector const& direction, TimeType t) {
+  template <typename TOutputE, typename TOutputN>
+  inline bool CONEXhybrid<TOutputE, TOutputN>::addParticle(
+      Code pid, HEPEnergyType energy, HEPEnergyType mass, Point const& position,
+      DirectionVector const& direction, TimeType t, double weight) {
 
     auto const it = std::find_if(egs_em_codes_.cbegin(), egs_em_codes_.cend(),
                                  [=](auto const& p) { return pid == p.first; });
@@ -182,14 +212,11 @@ namespace corsika {
     double const v = direction.dot(x_sf_).magnitude();
     double const w = direction.dot(showerAxis_.getDirection()).magnitude();
 
-    double const weight = 1; // NEEDS TO BE CHANGED WHEN WE HAVE WEIGHTS!
-
     // generation, TO BE CHANGED WHEN WE HAVE THAT INFORMATION AVAILABLE
     int const latchin = 1;
 
     double const E = energy / 1_GeV;
     double const m = mass / 1_GeV;
-    energy_em_ += energy;
 
     CORSIKA_LOG_DEBUG("CONEXhybrid: removing {} {:5e} GeV", egs_pid, energy);
 
@@ -232,8 +259,9 @@ namespace corsika {
     return true;
   }
 
+  template <typename TOutputE, typename TOutputN>
   template <typename TStack>
-  inline void CONEXhybrid::doCascadeEquations(TStack&) {
+  inline void CONEXhybrid<TOutputE, TOutputN>::doCascadeEquations(TStack&) {
 
     ::conex::conexcascade_();
 
@@ -268,33 +296,25 @@ namespace corsika {
     ::conex::get_shower_electron_(icute, nX, Electrons[0]);
     ::conex::get_shower_hadron_(icuth, nX, Hadrons[0]);
 
-    std::ofstream file{"conex_output.txt"};
-    file << fmt::format("#{:>10} {:>13} {:>13} {:>13} {:>13} {:>13} {:>13} {:>13}\n", "X",
-                        "N", "dEdX", "Mu", "dMu", "Photon", "El", "Had");
+    // make sure CONEX binning is same to C8:
+    GrammageType const dX = (X[1] - X[0]) * (1_g / square(1_cm));
+
     for (int i = 0; i < nX; ++i) {
-      file << fmt::format(" {:>10.2f} {:.5e} {:.5e} {:.5e} {:.5e} {:.5e} {:.5e} {:.5e}\n",
-                          X[i], N[i], dEdX[i], Mu[i], dMu[i], Photon[i], Electrons[i],
-                          Hadrons[i]);
+      GrammageType const curX = X[i] * (1_g / square(1_cm));
+      SubWriter<TOutputE>::write(curX, curX + dX,
+                                 Code::Unknown, // this is sum of all dEdX
+                                 dEdX[i] * 1_GeV / 1_g * square(1_cm) * dX);
+      SubWriter<TOutputN>::write(curX, curX + dX, Code::Photon, Photon[i]);
+      SubWriter<TOutputN>::write(curX, curX + dX, Code::Proton, Hadrons[i]);
+      SubWriter<TOutputN>::write(curX, curX + dX, Code::Electron, Electrons[i]);
+      SubWriter<TOutputN>::write(curX, curX + dX, Code::MuMinus, Mu[i]);
     }
-
-    std::ofstream fitout{"conex_fit.txt"};
-    fitout << fitpars[1 - 1] << " # log10(eprima/eV)" << std::endl;
-    fitout << fitpars[2 - 1] << " # theta" << std::endl;
-    fitout << fitpars[3 - 1] << " # X1 (first interaction)" << std::endl;
-    fitout << fitpars[4 - 1] << " # Nmax" << std::endl;
-    fitout << fitpars[5 - 1] << " # X0" << std::endl;
-    fitout << fitpars[6 - 1] << " # P1" << std::endl;
-    fitout << fitpars[7 - 1] << " # P2" << std::endl;
-    fitout << fitpars[8 - 1] << " # P3" << std::endl;
-    fitout << fitpars[9 - 1] << " # chi^2 / sqrt(Nmax)" << std::endl;
-    fitout << fitpars[10 - 1] << " # Xmax" << std::endl;
-    fitout << fitpars[11 - 1] << " # phi" << std::endl;
-    fitout << fitpars[12 - 1] << " # inelasticity 1st int." << std::endl;
-    fitout << fitpars[13 - 1] << " # ???" << std::endl;
   }
 
-  inline HEPEnergyType CONEXhybrid::getEnergyEM() const { return energy_em_; }
+  template <typename TOutputE, typename TOutputN>
+  inline YAML::Node CONEXhybrid<TOutputE, TOutputN>::getConfig() const {
 
-  inline void CONEXhybrid::reset() { energy_em_ = 0_GeV; }
+    return YAML::Node();
+  }
 
 } // namespace corsika

corsika/detail/modules/energy_loss/BetheBlochPDG.inl

 /*
  * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
  *
- * This software is distributed under the terms of the 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
@@ -19,19 +18,15 @@
 
 namespace corsika {
 
-  auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
-    return sqrt((Elab - m) * (Elab + m));
-  };
-
-  inline BetheBlochPDG::BetheBlochPDG(ShowerAxis const& shower_axis)
-      : dX_(10_g / square(1_cm)) // profile binning
-      , dX_threshold_(0.0001_g / square(1_cm))
-      , shower_axis_(shower_axis)
-      , profile_(int(shower_axis.getMaximumX() / dX_) + 1) {}
+  template <typename TOutput>
+  template <typename... TArgs>
+  inline BetheBlochPDG<TOutput>::BetheBlochPDG(TArgs&&... args)
+      : TOutput(std::forward<TArgs>(args)...) {}
 
+  template <typename TOutput>
   template <typename TParticle>
-  inline HEPEnergyType BetheBlochPDG::getBetheBloch(TParticle const& p,
-                                                    GrammageType const dX) {
+  inline HEPEnergyType BetheBlochPDG<TOutput>::getBetheBloch(TParticle const& p,
+                                                             GrammageType const dX) {
 
     // all these are material constants and have to come through Environment
     // right now: values for nitrogen_D
@@ -124,27 +119,30 @@ namespace corsika {
   }
 
   // radiation losses according to PDG 2018, ch. 33 ref. [5]
+  template <typename TOutput>
   template <typename TParticle>
-  inline HEPEnergyType BetheBlochPDG::getRadiationLosses(TParticle const& vP,
-                                                         GrammageType const vDX) {
+  inline HEPEnergyType BetheBlochPDG<TOutput>::getRadiationLosses(
+      TParticle const& vP, GrammageType const vDX) {
     // simple-minded hard-coded value for b(E) inspired by data from
     // http://pdg.lbl.gov/2018/AtomicNuclearProperties/ for N and O.
     auto constexpr b = 3.0 * 1e-6 * square(1_cm) / 1_g;
     return -vP.getEnergy() * b * vDX;
   }
 
+  template <typename TOutput>
   template <typename TParticle>
-  inline HEPEnergyType BetheBlochPDG::getTotalEnergyLoss(TParticle const& vP,
-                                                         GrammageType const vDX) {
+  inline HEPEnergyType BetheBlochPDG<TOutput>::getTotalEnergyLoss(
+      TParticle const& vP, GrammageType const vDX) {
     return getBetheBloch(vP, vDX) + getRadiationLosses(vP, vDX);
   }
 
-  template <typename TParticle, typename TTrajectory>
-  inline ProcessReturn BetheBlochPDG::doContinuous(TParticle& p, TTrajectory const& t,
-                                                   bool const) {
+  template <typename TOutput>
+  template <typename TParticle>
+  inline ProcessReturn BetheBlochPDG<TOutput>::doContinuous(Step<TParticle>& step,
+                                                            bool const) {
 
     // if this step was limiting the CORSIKA stepping, the particle is lost
-    /* see Issue https://gitlab.ikp.kit.edu/AirShowerPhysics/corsika/-/issues/389
+    /* see Issue https://gitlab.iap.kit.edu/AirShowerPhysics/corsika/-/issues/389
     if (limitStep) {
       fillProfile(t, p.getEnergy());
       p.setEnergy(p.getMass());
@@ -152,25 +150,32 @@ namespace corsika {
     }
     */
 
-    if (p.getChargeNumber() == 0) return ProcessReturn::Ok;
-
-    GrammageType const dX = p.getNode()->getModelProperties().getIntegratedGrammage(t);
-    CORSIKA_LOG_TRACE("EnergyLoss pid={}, z={}, dX={} g/cm2", p.getPID(),
-                      p.getChargeNumber(), dX / 1_g * square(1_cm));
-    HEPEnergyType dE = getTotalEnergyLoss(p, dX);
-    auto E = p.getEnergy();
-    [[maybe_unused]] const auto Ekin = E - p.getMass();
-    auto Enew = E + dE;
-    CORSIKA_LOG_TRACE("EnergyLoss  dE={} MeV, E={} GeV, Ekin={} GeV, Enew={} GeV",
-                      dE / 1_MeV, E / 1_GeV, Ekin / 1_GeV, Enew / 1_GeV);
-    p.setEnergy(Enew); // kinetic energy on stack
-    fillProfile(t, dE);
+    if (step.getParticlePre().getChargeNumber() == 0) return ProcessReturn::Ok;
+
+    GrammageType const dX =
+        step.getParticlePre().getNode()->getModelProperties().getIntegratedGrammage(
+            step.getStraightTrack());
+    CORSIKA_LOG_TRACE("EnergyLoss pid={}, z={}, dX={} g/cm2",
+                      step.getParticlePre().getPID(),
+                      step.getParticlePre().getChargeNumber(), dX / 1_g * square(1_cm));
+    HEPEnergyType const dE = getTotalEnergyLoss(step.getParticlePre(), dX);
+    //    if (dE > HEPEnergyType::zero())
+    //      dE = -dE;
+    CORSIKA_LOG_TRACE("EnergyLoss  dE={} MeV, Ekin={} GeV, EkinNew={} GeV", dE / 1_MeV,
+                      step.getEkinPre() / 1_GeV, (step.getEkinPre() + dE) / 1_GeV);
+    step.add_dEkin(dE);
+
+    // also send to output
+    TOutput::write(step.getPositionPre(), step.getPositionPost(),
+                   step.getParticlePre().getPID(),
+                   -step.getParticlePre().getWeight() * dE);
     return ProcessReturn::Ok;
   }
 
+  template <typename TOutput>
   template <typename TParticle, typename TTrajectory>
-  inline LengthType BetheBlochPDG::getMaxStepLength(TParticle const& vParticle,
-                                                    TTrajectory const& vTrack) const {
+  inline LengthType BetheBlochPDG<TOutput>::getMaxStepLength(
+      TParticle const& vParticle, TTrajectory const& vTrack) const {
     if (vParticle.getChargeNumber() == 0) {
       return meter * std::numeric_limits<double>::infinity();
     }
@@ -180,7 +185,7 @@ namespace corsika {
     auto const energy = vParticle.getKineticEnergy();
     auto const energy_lim =
         std::max(energy * 0.9, // either 10% relative loss max., or
-                 calculate_kinetic_energy_threshold(
+                 get_kinetic_energy_propagation_threshold(
                      vParticle.getPID()) // energy thresholds globally defined
                                          // for individual particles
                      * 0.99999 // need to go slightly below global e-cut to assure
@@ -193,80 +198,12 @@ namespace corsika {
         vTrack, maxGrammage);
   }
 
-  template <typename TTrajectory>
-  inline void BetheBlochPDG::fillProfile(TTrajectory const& vTrack,
-                                         const HEPEnergyType dE) {
-
-    GrammageType grammageStart = shower_axis_.getProjectedX(vTrack.getPosition(0));
-    GrammageType grammageEnd = shower_axis_.getProjectedX(vTrack.getPosition(1));
-
-    if (grammageStart > grammageEnd) { // particle going upstream
-      std::swap(grammageStart, grammageEnd);
-    }
-
-    GrammageType const deltaX = grammageEnd - grammageStart;
-
-    if (deltaX < dX_threshold_) return;
-
-    // only register the range that is covered by the profile
-    int const maxBin = int(profile_.size() - 1);
-    int binStart = grammageStart / dX_;
-    if (binStart < 0) binStart = 0;
-    if (binStart > maxBin) binStart = maxBin;
-    int binEnd = grammageEnd / dX_;
-    if (binEnd < 0) binEnd = 0;
-    if (binEnd > maxBin) binEnd = maxBin;
-
-    CORSIKA_LOG_TRACE("energy deposit of -dE={} between {} and {}", -dE, grammageStart,
-                      grammageEnd);
-
-    auto energyCount = HEPEnergyType::zero();
-
-    auto const factor = -dE / deltaX;
-    auto fill = [&](int const bin, GrammageType const weight) {
-      auto const increment = factor * weight;
-      profile_[bin] += increment;
-      energyCount += increment;
-
-      CORSIKA_LOG_TRACE("filling bin {} with weight {} : {} ", bin, weight, increment);
-    };
-
-    // fill longitudinal profile
-    if (binStart == binEnd) {
-      fill(binStart, deltaX);
-    } else {
-      fill(binStart, ((1 + binStart) * dX_ - grammageStart));
-      fill(binEnd, (grammageEnd - binEnd * dX_));
-      for (int bin = binStart + 1; bin < binEnd; ++bin) { fill(bin, dX_); }
-    }
-
-    CORSIKA_LOG_TRACE("total energy added to histogram: {} ", energyCount);
-  }
-
-  inline void BetheBlochPDG::printProfile() const {
-    std::ofstream file("EnergyLossProfile.dat");
-    file << "# EnergyLoss profile" << std::endl
-         << "# lower X bin edge [g/cm2]  dE/dX [GeV/g/cm2]\n";
-    double const deltaX = dX_ / 1_g * square(1_cm);
-    for (size_t i = 0; i < profile_.size(); ++i) {
-      file << std::scientific << std::setw(15) << i * deltaX << std::setw(15)
-           << profile_.at(i) / (deltaX * 1_GeV) << '\n';
-    }
-    file.close();
-  }
-
-  inline HEPEnergyType BetheBlochPDG::getTotal() const {
-    return std::accumulate(profile_.cbegin(), profile_.cend(), HEPEnergyType::zero());
-  }
+  template <typename TOutput>
+  inline YAML::Node BetheBlochPDG<TOutput>::getConfig() const {
 
-  inline void BetheBlochPDG::showResults() const {
-    CORSIKA_LOG_INFO(
-        " ******************************\n"
-        " PROCESS::ContinuousProcess: \n"
-        " energy lost dE (GeV)      : {}\n  ",
-        energy_lost_ / 1_GeV);
+    YAML::Node node;
+    node["type"] = "BetheBlochPDG";
+    return node;
   }
 
-  inline void BetheBlochPDG::reset() { energy_lost_ = 0_GeV; }
-
 } // namespace corsika

corsika/detail/modules/epos/InteractionModel.inl

 /*
  * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
  *
- * This software is distributed under the terms of the 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
 
 #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>
@@ -24,32 +23,58 @@
 
 namespace corsika::epos {
 
-  inline InteractionModel::InteractionModel(std::string const& dataPath,
+  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);
+      }
     }
-    setParticlesStable();
   }
 
-  inline void InteractionModel::setParticlesStable() const {
-    CORSIKA_LOGGER_DEBUG(logger_,
-                         "set all particles known to CORSIKA stable inside EPOS..");
-    for (auto& p : get_all_particles()) {
-      if (!is_hadron(p)) continue;
+  inline void InteractionModel::setParticleListStable(std::set<Code> vPartList) const {
+    for (auto& p : vPartList) {
       int const eid = convertToEposRaw(p);
       if (eid != 0) {
-        ::epos::nodcy_.nrnody = ::epos::nodcy_.nrnody + 1;
-        ::epos::nodcy_.nody[::epos::nodcy_.nrnody - 1] = eid;
+        // 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,
@@ -59,7 +84,7 @@ namespace corsika::epos {
     if (!is_nucleus(targetId) && targetId != Code::Neutron && targetId != Code::Proton) {
       return false;
     }
-    if (is_nucleus(targetId) && (get_nucleus_A(targetId) >= maxTargetMassNumber_)) {
+    if (is_nucleus(targetId) && (get_nucleus_A(targetId) >= get_nucleus_A(maxNucleus_))) {
       return false;
     }
     if ((minEnergyCoM_ > sqrtS) || (sqrtS > maxEnergyCoM_)) { return false; }
@@ -115,6 +140,10 @@ namespace corsika::epos {
 
     ::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);
@@ -136,9 +165,10 @@ namespace corsika::epos {
     strcpy(::epos::fname_.fncs, CS.data);
     ::epos::nfname_.nfncs = CS.length;
 
-    // initialiazes maximum energy and mass
-    initializeEventCoM(Code::Lead, Lead::nucleus_A, Lead::nucleus_Z, Code::Lead,
-                       Lead::nucleus_A, Lead::nucleus_Z, 1_PeV);
+    // 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,
@@ -284,7 +314,7 @@ namespace corsika::epos {
                          "beamA={}, "
                          "beamZ={} "
                          "targetId={}, "
-                         "ELab={:4.3f} GeV,",
+                         "ELab={:12.2f} GeV,",
                          BeamId, BeamA, BeamZ, TargetId, EnergyLab / 1_GeV);
 
     // read cross section from epos internal tables
@@ -300,7 +330,7 @@ namespace corsika::epos {
       int const iBeam = epos::getEposXSCode(
           BeamId); // 0 (can not interact, 1: pion-like, 2: proton-like, 3:kaon-like)
       CORSIKA_LOGGER_TRACE(logger_,
-                           "projectile cross section type={} "
+                           "readCrossSectionTableLab: projectile cross section type={} "
                            "(0: cannot interact, 1:pion, 2:baryon, 3:kaon)",
                            iBeam);
 
@@ -314,16 +344,25 @@ namespace corsika::epos {
 
     int iMode = 3; // 0: air, >0 not air
 
-    CORSIKA_LOGGER_DEBUG(logger_,
-                         "inside Epos "
+    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);
   }
 
@@ -417,7 +456,7 @@ namespace corsika::epos {
     if (is_nucleus(targetId)) {
       targetA = get_nucleus_A(targetId);
       targetZ = get_nucleus_Z(targetId);
-      CORSIKA_LOGGER_DEBUG(logger_, "target: A={}, Z={} ", beamA, beamZ);
+      CORSIKA_LOGGER_DEBUG(logger_, "target: A={}, Z={} ", targetA, targetZ);
     }
     initializeEventCoM(projectileId, beamA, beamZ, targetId, targetA, targetZ, sqrtSNN);
 

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

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

@@ -2,133 +2,170 @@
 /*
  * (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>
 
 namespace corsika::proposal {
 
-  inline 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 const emCut =
-        calculate_kinetic_energy_threshold(code) +
-        get_mass(code); //! energy thresholds globally defined for individual particles
-    auto c = p_cross->second(media.at(comp.getHash()), emCut);
+    auto c = p_cross->second(media.at(component_hash), proposal_energycutsettings[code]);
 
-    // Use higland multiple scattering and deactivate stochastic deflection by
-    // passing an empty vector
-    static constexpr auto ms_type = PROPOSAL::MultipleScatteringType::Highland;
-    auto s_type = std::vector<PROPOSAL::InteractionType>();
+    // choose multiple scattering model
+    static constexpr auto ms_type = PROPOSAL::MultipleScatteringType::MoliereInterpol;
 
     // Build displacement integral and scattering object and interpolate them too and
     // saved in the calc map by a key build out of a hash of composed of the component and
     // particle code.
-    auto calculator =
-        Calculator{PROPOSAL::make_displacement(c, true),
-                   PROPOSAL::make_scattering(ms_type, s_type, particle[code],
-                                             media.at(comp.getHash()))};
-    calc[std::make_pair(comp.getHash(), code)] = std::move(calculator);
+    auto calculator = Calculator{PROPOSAL::make_displacement(c, true),
+                                 PROPOSAL::make_multiple_scattering(
+                                     ms_type, particle[code], media.at(component_hash))};
+    calc[std::make_pair(component_hash, code)] = std::move(calculator);
   }
 
-  template <typename TEnvironment>
-  inline ContinuousProcess::ContinuousProcess(TEnvironment const& _env)
-      : ProposalProcessBase(_env) {}
+  template <typename TOutput>
+  template <typename TEnvironment, typename... TOutputArgs>
+  inline ContinuousProcess<TOutput>::ContinuousProcess(TEnvironment const& _env,
+                                                       TOutputArgs&&... args)
+      : ProposalProcessBase(_env)
+      , TOutput(args...) {
+    //! Initialize PROPOSAL tables for all media and all particles
+    for (auto const& medium : media) {
+      for (auto const particle_code : tracked) {
+        buildCalculator(particle_code, medium.first);
+      }
+    }
+  }
 
+  template <typename TOutput>
   template <typename TParticle>
-  inline void ContinuousProcess::scatter(TParticle& vP, HEPEnergyType const& loss,
-                                         GrammageType const& grammage) {
+  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::Cartesian3D(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 = std::array<double, 4>();
-    for (auto& it : rnd) it = distr(RNG_);
-
-    // calculate deflection based on particle energy, loss
-    auto deflection = (c->second).scatter->CalculateMultipleScattering(
-        grammage / 1_g * square(1_cm), vP.getEnergy() / 1_MeV, E_f / 1_MeV, rnd);
-
-    [[maybe_unused]] auto [unused1, final_direction] =
-        PROPOSAL::multiple_scattering::ScatterInitialDirection(direction, deflection);
+    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
-    vP.setDirection(
-        {vP_dir.getCoordinateSystem(),
-         {final_direction.GetX(), final_direction.GetY(), final_direction.GetZ()}});
+    DirectionVector diff_dir_ = scattered_particle_dir - initial_particle_dir;
+    step.add_dU(diff_dir_);
   }
 
-  template <typename TParticle, typename TTrajectory>
-  inline ProcessReturn ContinuousProcess::doContinuous(TParticle& vP,
-                                                       TTrajectory const& vT,
-                                                       bool const) {
-    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);
+    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 = (c->second).disp->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); // on the stack, this is just kinetic energy, E-m
+    if (step.getParticlePre().getChargeNumber() != 0) scatter(step, dE, dX);
+    step.add_dEkin(-dE); // on the stack, this is just kinetic energy, E-m
+
+    // also send to output
+    TOutput::write(step.getPositionPre(), step.getPositionPost(),
+                   step.getParticlePre().getPID(),
+                   step.getParticlePre().getWeight() * dE);
+
     return ProcessReturn::Ok;
   }
 
+  template <typename TOutput>
   template <typename TParticle, typename TTrajectory>
-  inline LengthType ContinuousProcess::getMaxStepLength(TParticle const& vP,
-                                                        TTrajectory const& vT) {
+  inline LengthType ContinuousProcess<TOutput>::getMaxStepLength(
+      TParticle const& vP, TTrajectory const& track) {
     auto const code = vP.getPID();
     if (!canInteract(code)) return meter * std::numeric_limits<double>::infinity();
+    if (code == Code::Photon)
+      return meter *
+             std::numeric_limits<double>::infinity(); // no step limitation for photons
 
     // Limit the step size of a conitnuous loss. The maximal continuous loss seems to be
     // a hyper parameter which must be adjusted.
     //
     auto const energy = vP.getEnergy();
-    auto const energy_lim =
-        std::max(energy * 0.9, // either 10% relative loss max., or
-                 calculate_kinetic_energy_threshold(
-                     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.
-        );
+    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);
@@ -137,20 +174,23 @@ namespace corsika::proposal {
         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;
   }
 
-  inline void ContinuousProcess::showResults() const {
-    CORSIKA_LOG_DEBUG(
-        " ******************************\n"
-        " PROCESS::ContinuousProcess: \n"
-        " energy lost dE (GeV)      :  {}",
-        energy_lost_ / 1_GeV);
+  template <typename TOutput>
+  inline YAML::Node ContinuousProcess<TOutput>::getConfig() const {
+    return YAML::Node();
   }
 
-  inline 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 → 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 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>
 #include <corsika/media/NuclearComposition.hpp>
-#include <corsika/modules/proposal/Interaction.hpp>
 #include <corsika/framework/utility/COMBoost.hpp>
 #include <corsika/framework/core/PhysicalUnits.hpp>
 
@@ -16,42 +14,57 @@
 #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 Interaction::Interaction(TEnvironment const& _env)
-      : ProposalProcessBase(_env) {}
+  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);
+      }
+    }
+  }
 
-  inline void Interaction::buildCalculator(Code code, NuclearComposition const& comp) {
+  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
+    // take some minutes if you have to build the tables and cannot read the tables
     // from disk
-    auto const emCut =
-        get_kinetic_energy_threshold(code) +
-        get_mass(code); //! energy thresholds globally defined for individual particles
-
-    auto c = p_cross->second(media.at(comp.getHash()), emCut);
+    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(comp.getHash(), code)] = std::make_tuple(
-        PROPOSAL::make_secondaries(inter_types, particle[code], media.at(comp.getHash())),
-        PROPOSAL::make_interaction(c, true));
+    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 Interaction::doInteraction(TStackView& view,
-                                                  Code const projectileId,
-                                                  FourMomentum const& projectileP4) {
+  inline ProcessReturn
+  InteractionModel<THadronicLEModel, THadronicHEModel>::doInteraction(
+      TStackView& view, Code const projectileId,
+      [[maybe_unused]] FourMomentum const& projectileP4) {
 
     auto const projectile = view.getProjectile();
 
@@ -69,6 +82,17 @@ namespace corsika::proposal {
       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>(
@@ -87,25 +111,68 @@ namespace corsika::proposal {
       auto loss = PROPOSAL::StochasticLoss(
           static_cast<int>(type), v * projectile.getEnergy() / 1_MeV, point, direction,
           projectile.getTime() / 1_s, 0., projectile.getEnergy() / 1_MeV);
-      auto target = PROPOSAL::Component::component_map->operator[](target_hash);
+      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()});
-        auto sec_code = convert_from_PDG(static_cast<PDGCode>(s.type));
-        view.addSecondary(std::make_tuple(sec_code, E - get_mass(sec_code), dir));
+
+        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 Interaction::getCrossSection(TParticle const& projectile,
-                                                       Code const projectileId,
-                                                       FourMomentum const& projectileP4) {
+  inline CrossSectionType
+  InteractionModel<THadronicLEModel, THadronicHEModel>::getCrossSection(
+      TParticle const& projectile, Code const projectileId,
+      FourMomentum const& projectileP4) {
 
     // ==============================================
     // this block better diappears. RU 26.10.2021
@@ -126,4 +193,64 @@ namespace corsika::proposal {
 
     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>
@@ -23,8 +22,28 @@
 namespace corsika::proposal {
 
   inline bool ProposalProcessBase::canInteract(Code pcode) const {
-    if (std::find(begin(tracked), end(tracked), pcode) != end(tracked)) return true;
-    return false;
+    return std::find(begin(tracked), end(tracked), pcode) != end(tracked);
+  }
+
+  inline HEPEnergyType ProposalProcessBase::getOptimizedEmCut(Code code) const {
+    // get energy above which energy losses need to be considered
+    auto const production_threshold = get_energy_production_threshold(code);
+
+    HEPEnergyType lowest_table_value = 0_GeV;
+
+    // find tables for EnergyCuts closest (but still smaller than) production_threshold
+    for (auto const& table_energy : energycut_table_values) {
+      if (table_energy <= production_threshold && table_energy > lowest_table_value) {
+        lowest_table_value = table_energy;
+      }
+    }
+
+    if (lowest_table_value == 0_GeV) {
+      // no appropriate table available
+      return production_threshold;
+    }
+
+    return lowest_table_value;
   }
 
   template <typename TEnvironment>
@@ -37,7 +56,7 @@ namespace corsika::proposal {
         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;
@@ -54,6 +73,34 @@ namespace corsika::proposal {
     //! path, otherwise interpolation tables would only stored in main memory if
     //! no explicit intrpolation def is specified.
     PROPOSAL::InterpolationSettings::TABLES_PATH = corsika_data("PROPOSAL").c_str();
+
+    //! Initialize EnergyCutSettings
+    for (auto const particle_code : tracked) {
+      if (particle_code == Code::Photon) {
+        // no EnergyCut for photon, only-stochastic propagation
+        continue;
+      }
+      // use optimized emcut for which tables should be available
+      auto optimized_ecut = getOptimizedEmCut(particle_code);
+      CORSIKA_LOG_INFO("PROPOSAL: Use tables with EnergyCut of {} for particle {}.",
+                       optimized_ecut, particle_code);
+      proposal_energycutsettings[particle_code] = optimized_ecut;
+    }
+
+    //! Initialize PROPOSAL tables for all media and all particles
+    for (auto const& medium : media) {
+      for (auto const particle_code : tracked) {
+        buildTables(medium.second, particle_code,
+                    proposal_energycutsettings[particle_code]);
+      }
+    }
+  }
+
+  void ProposalProcessBase::buildTables(PROPOSAL::Medium medium, Code particle_code,
+                                        HEPEnergyType emcut) {
+    CORSIKA_LOG_DEBUG("PROPOSAL: Initialize tables for particle {} in {}", particle_code,
+                      medium.GetName());
+    cross[particle_code](medium, emcut);
   }
 
   inline size_t ProposalProcessBase::hash::operator()(const calc_key_t& p) const

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>
@@ -33,8 +32,7 @@ namespace corsika::pythia8 {
     // run this only once during construction
 
     // link random number generator in pythia to CORSIKA8
-    Pythia8::RndmEngine* rndm = new corsika::pythia8::Random();
-    Pythia8::Pythia::setRndmEnginePtr(rndm);
+    Pythia8::Pythia::setRndmEnginePtr(std::make_shared<corsika::pythia8::Random>());
 
     Pythia8::Pythia::readString("Next:numberShowInfo = 0");
     Pythia8::Pythia::readString("Next:numberShowProcess = 0");
@@ -90,7 +88,7 @@ namespace corsika::pythia8 {
 
   inline void Decay::setHandleDecay(std::vector<Code> const& vParticleList) {
     handleAllDecays_ = false;
-    for (auto p : vParticleList) setHandleDecay(p);
+    for (auto const p : vParticleList) setHandleDecay(p);
   }
 
   inline bool Decay::isDecayHandled(Code const vParticleCode) {
@@ -101,7 +99,7 @@ namespace corsika::pythia8 {
   }
 
   inline void Decay::setStable(std::vector<Code> const& particleList) {
-    for (auto p : particleList) Decay::setStable(p);
+    for (auto const p : particleList) Decay::setStable(p);
   }
 
   inline void Decay::setUnstable(Code const pCode) {
@@ -136,7 +134,7 @@ namespace corsika::pythia8 {
     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);
   }
 
@@ -218,27 +216,34 @@ namespace corsika::pythia8 {
     // 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);
-        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()));
+        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 <tuple>
-
-namespace corsika::pythia8 {
-
-  inline Interaction::~Interaction() {
-    CORSIKA_LOG_INFO("Pythia::Interaction n= {}", count_);
-  }
-
-  inline Interaction::Interaction(bool const print_listing)
-      : Pythia8::Pythia(CORSIKA_Pythia8_XML_DIR)
-      , print_listing_(print_listing) {
-
-    CORSIKA_LOG_INFO("Configuring Pythia8 from: {}", CORSIKA_Pythia8_XML_DIR);
-
-    // initialize Pythia
-
-    // reduce output from pythia if set to "on"
-    Pythia8::Pythia::readString("Print:quiet = on");
-    // check if data in particle data file is minimally consistent. Very verbose! set to
-    // "off"! we do not change the basic file provided by pythia.
-    Pythia8::Pythia::readString("Check:particleData = off");
-    Pythia8::Pythia::readString("Check:event = on");             // default: on
-    Pythia8::Pythia::readString("Check:levelParticleData = 12"); // 1 is default
-    /** \TODO: proper process initialization for MinBias needed, see
-        also Issue https://gitlab.ikp.kit.edu/AirShowerPhysics/corsika/-/issues/369 **/
-    Pythia8::Pythia::readString("HardQCD:all = on");
-    Pythia8::Pythia::readString("ProcessLevel:resonanceDecays = off");
-
-    // we can't test this block, LCOV_EXCL_START
-    if (!Pythia8::Pythia::init())
-      throw std::runtime_error("Pythia::Interaction: Initialization failed!");
-    // LCOV_EXCL_STOP
-
-    // 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
-      std::vector<Code> const 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(&(Pythia8::Pythia::info), Pythia8::Pythia::settings,
-                &(Pythia8::Pythia::particleData), &(Pythia8::Pythia::rndm));
-  }
-
-  inline void Interaction::setStable(std::vector<Code> const& particleList) {
-    for (auto p : particleList) Interaction::setStable(p);
-  }
-
-  inline void Interaction::setUnstable(Code const pCode) {
-    CORSIKA_LOG_DEBUG("Pythia::Interaction: setting {} unstable..", pCode);
-    Pythia8::Pythia::particleData.mayDecay(static_cast<int>(get_PDG(pCode)), true);
-  }
-
-  inline void Interaction::setStable(Code const pCode) {
-    CORSIKA_LOG_DEBUG("Pythia::Interaction: setting {} stable..", pCode);
-    Pythia8::Pythia::particleData.mayDecay(static_cast<int>(get_PDG(pCode)), false);
-  }
-
-  inline bool Interaction::isValid(Code const projectileId, Code const targetId,
-                                   HEPEnergyType const sqrtS) const {
-
-    if ((10_GeV > sqrtS) || (sqrtS > 1_PeV)) { return false; }
-
-    if (targetId != Code::Hydrogen && targetId != Code::Neutron &&
-        targetId != Code::Proton) {
-      return false;
-    }
-
-    if (is_nucleus(projectileId)) { return false; }
-
-    if (!canInteract(projectileId)) { return false; }
-    return true;
-  }
-
-  inline void Interaction::configureLabFrameCollision(Code const projectileId,
-                                                      Code const targetId,
-                                                      HEPEnergyType const 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(projectileId));
-    std::stringstream stBeam;
-    stBeam << "Beams:idA = " << pdgBeam;
-    Pythia8::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;
-    Pythia8::Pythia::readString(stTarget.str());
-    // set frame to lab. frame
-    Pythia8::Pythia::readString("Beams:frameType = 2");
-    // set beam energy
-    double const Elab = BeamEnergy / 1_GeV;
-    std::stringstream stEnergy;
-    stEnergy << "Beams:eA = " << Elab;
-    Pythia8::Pythia::readString(stEnergy.str());
-    // target at rest
-    Pythia8::Pythia::readString("Beams:eB = 0.");
-    // initialize this config
-
-    // we can't test this block, LCOV_EXCL_START
-    if (!Pythia8::Pythia::init())
-      throw std::runtime_error("Pythia::Interaction: Initialization failed!");
-    // LCOV_EXCL_STOP
-  }
-
-  inline bool Interaction::canInteract(Code const pCode) const {
-    return pCode == Code::Proton || pCode == Code::Neutron || pCode == Code::AntiProton ||
-           pCode == Code::AntiNeutron || pCode == Code::PiMinus || pCode == Code::PiPlus;
-  }
-
-  inline std::tuple<CrossSectionType, CrossSectionType>
-  Interaction::getCrossSectionInelEla(Code const projectileId, Code const targetId,
-                                      FourMomentum const& projectileP4,
-                                      FourMomentum const& targetP4) const {
-
-    HEPEnergyType const CoMenergy = (projectileP4 + targetP4).getNorm();
-
-    if (!isValid(projectileId, targetId, CoMenergy)) {
-      return {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 = CoMenergy / 1_GeV;
-
-    //! @todo: remove this const_cast, when Pythia8 becomes const-correct! CHECK!
-    Pythia8::SigmaTotal& sigma = *const_cast<Pythia8::SigmaTotal*>(&sigma_);
-
-    // calculate cross section
-    sigma.calc(pdgCodeBeam, pdgCodeTarget, ecm);
-    if (sigma.hasSigmaTot()) {
-      double const sigEla = sigma.sigmaEl();
-      double const sigProd = sigma.sigmaTot() - sigEla;
-
-      return std::make_tuple(sigProd * (1_fm * 1_fm), sigEla * (1_fm * 1_fm));
-
-    } else {
-      // we can't test pythia8 internals, LCOV_EXCL_START
-      throw std::runtime_error("pythia cross section init failed");
-      // we can't test pythia8 internals, LCOV_EXCL_STOP
-    }
-  }
-
-  template <class TView>
-  inline void Interaction::doInteraction(TView& view, Code const projectileId,
-                                         Code const targetId,
-                                         FourMomentum const& projectileP4,
-                                         FourMomentum const& targetP4) {
-
-    auto projectile = view.getProjectile();
-
-    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.");
-    }
-
-    // position and time of interaction
-    Point const& pOrig = projectile.getPosition();
-    TimeType const tOrig = projectile.getTime();
-
-    CORSIKA_LOG_DEBUG("Interaction: ebeam lab: {} GeV", eProjectileLab / 1_GeV);
-
-    // define target kinematics in lab frame
-    // define boost to and from CoM frame
-    // CoM frame definition in Pythia projectile: +z
-    COMBoost const boost(projectileP4, constants::nucleonMass);
-    auto const& labCS = boost.getOriginalCS();
-
-    CORSIKA_LOG_DEBUG("Interaction: position of interaction: ", pOrig.getCoordinates());
-    CORSIKA_LOG_DEBUG("Interaction: time: {}", tOrig);
-
-    CORSIKA_LOG_DEBUG(
-        "Interaction: "
-        " doInteraction: E(GeV): {}"
-        " Ecm(GeV): {}",
-        eProjectileLab / 1_GeV, sqrtS / 1_GeV);
-
-    count_++;
-
-    configureLabFrameCollision(projectileId, targetId, eProjectileLab);
-
-    // create event in pytia. LCOV_EXCL_START: we don't validate pythia8 internals
-    if (!Pythia8::Pythia::next())
-      throw std::runtime_error("Pythia::DoInteraction: failed!");
-    // LCOV_EXCL_STOP
-
-    // link to pythia stack
-    Pythia8::Event& event = Pythia8::Pythia::event;
-
-    // LCOV_EXCL_START, we don't validate pythia8 internals
-    if (print_listing_) {
-      // print final state
-      event.list();
-    }
-    // LCOV_EXCL_STOP
-
-    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 const& 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()));
-
-      MomentumVector const pyPlab(labCS,
-                                  {p8p.px() * 1_GeV, p8p.py() * 1_GeV, p8p.pz() * 1_GeV});
-
-      HEPEnergyType const mass = get_mass(pyId);
-      HEPEnergyType const Ekin = sqrt(pyPlab.getSquaredNorm() + mass * mass) - mass;
-
-      // add to corsika stack
-      auto pnew =
-          projectile.addSecondary(std::make_tuple(pyId, Ekin, pyPlab.normalized()));
-
-      Plab_final += pnew.getMomentum();
-      Elab_final += pnew.getEnergy();
-    }
-
-    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 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

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 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/Random.hpp>
 #include <corsika/modules/qgsjetII/InteractionModel.hpp>
 #include <corsika/modules/qgsjetII/ParticleConversion.hpp>
 #include <corsika/modules/qgsjetII/QGSJetIIFragmentsStack.hpp>
@@ -25,11 +25,12 @@
 namespace corsika::qgsjetII {
 
   inline InteractionModel::InteractionModel(boost::filesystem::path const dataPath) {
-    CORSIKA_LOG_DEBUG("Reading QGSJetII data tables from {}", 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);
@@ -69,32 +70,45 @@ namespace corsika::qgsjetII {
       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 / (is_nucleus(projectileId) ? get_nucleus_A(projectileId) : 1) +
-         targetP4 / (is_nucleus(targetId) ? get_nucleus_A(targetId) : 1))
-            .getNormSqr();
+        (projectileP4 / AfactorProjectile + targetP4 / AfactorTarget).getNormSqr();
     auto const sqrtSNN = sqrt(SNN);
     if (!isValid(projectileId, targetId, sqrtSNN)) { return CrossSectionType::zero(); }
-    HEPEnergyType const Elab =
-        calculate_lab_energy(SNN, get_mass(projectileId), get_mass(targetId));
+
+    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));
-    int iTarget = 1;
-    if (is_nucleus(targetId)) { iTarget = get_nucleus_A(targetId); }
-    int iProjectile = 1;
-    if (is_nucleus(projectileId)) { iProjectile = get_nucleus_A(projectileId); }
 
     CORSIKA_LOG_DEBUG(
         "QgsjetII::getCrossSection Elab= {} GeV iBeam= {}"
         " iProjectile= {} iTarget= {}",
-        Elab / 1_GeV, iBeam, iProjectile, iTarget);
-    double sigProd = qgsect_(Elab / 1_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,
@@ -107,46 +121,40 @@ namespace corsika::qgsjetII {
         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 / (is_nucleus(projectileId) ? get_nucleus_A(projectileId) : 1) +
-         targetP4 / (is_nucleus(targetId) ? get_nucleus_A(targetId) : 1))
-            .getNormSqr();
-    auto const sqrtS = sqrt(SNN);
+        (projectileP4 / AfactorProjectile + targetP4 / AfactorTarget).getNormSqr();
+    auto const sqrtSNN = sqrt(SNN);
+
     if (!corsika::qgsjetII::canInteract(projectileId) ||
-        !isValid(projectileId, targetId, sqrtS)) {
-      throw std::runtime_error("invalid target/projectile/energy combination.");
+        !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 projectileMass =
-        (is_nucleus(projectileId) ? constants::nucleonMass : get_mass(projectileId));
-    auto const targetMass =
-        (projectileId == Code::Proton
-             ? get_mass(Code::Proton)
-             : constants::nucleonMass); // qgsjet target is always proton or nucleon.
-                                        // always nucleon??
 
-    // lab energy/hadron
-    HEPEnergyType const Elab = calculate_lab_energy(SNN, projectileMass, targetMass);
+    auto const projMass = get_mass(projectileId);
+    auto const targetMass = get_mass(targetId);
 
-    int beamA = 0;
-    if (is_nucleus(projectileId)) { beamA = get_nucleus_A(projectileId); }
+    // 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 ", Elab / 1_GeV);
+    CORSIKA_LOG_DEBUG("ebeam lab: {} GeV per projectile nucleon", ElabN / 1_GeV);
 
-    int targetMassNumber = 1;   // proton
-    if (is_nucleus(targetId)) { // nucleus
-      targetMassNumber = get_nucleus_A(targetId);
-    }
+    int const targetMassNumber = AfactorTarget;
     CORSIKA_LOG_DEBUG("target: {}, qgsjetII code/A: {}", targetId, targetMassNumber);
 
     // select QGSJetII internal projectile type
-    int projectileMassNumber = 1; // "1" means "hadron"
     QgsjetIIHadronType qgsjet_hadron_type = qgsjetII::getQgsjetIIHadronType(projectileId);
     if (qgsjet_hadron_type == QgsjetIIHadronType::NucleusType) {
-      projectileMassNumber = get_nucleus_A(projectileId);
-      std::array<QgsjetIIHadronType, 2> constexpr nucleons = {
-          QgsjetIIHadronType::ProtonType, QgsjetIIHadronType::NeutronType};
-      std::uniform_int_distribution select(0, 1);
-      qgsjet_hadron_type = nucleons[select(rng_)];
+      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_;
@@ -156,11 +164,12 @@ namespace corsika::qgsjetII {
     }
 
     count_++;
-    int qgsjet_hadron_type_int = static_cast<QgsjetIICodeIntType>(qgsjet_hadron_type);
+    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, projectileMassNumber, targetMassNumber);
-    qgini_(Elab / 1_GeV, 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();
@@ -177,46 +186,40 @@ namespace corsika::qgsjetII {
     // rest.
 
     // system of initial-state
-    COMBoost boost(projectileP4 / projectileMassNumber, targetP4 / targetMassNumber);
+    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 =
-        sqrt((Elab - projectileMass) * (Elab + projectileMass));
-    MomentumVector pLab(csPrime, {0_eV, 0_eV, pLabMag});
+    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 boostInternal({Elab, pLab}, targetMass);
+    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
-        std::uniform_real_distribution<double> select;
-        Code idFragm = Code::Proton;
-        if (select(rng_) > 0.5) { idFragm = Code::Neutron; }
+        Code const idFragm = bernoulli_(rng_) ? Code::Proton : Code::Neutron;
 
-        const HEPMassType nucleonMass = get_mass(idFragm);
+        HEPMassType const nucleonMass = get_mass(idFragm);
         // no pT, fragments just go forward
-        HEPEnergyType const projectileEnergyLabPerNucleon = Elab / beamA;
-        MomentumVector momentum{csPrime,
-                                {0.0_GeV, 0.0_GeV,
-                                 sqrt((projectileEnergyLabPerNucleon + nucleonMass) *
-                                      (projectileEnergyLabPerNucleon - nucleonMass))}};
+        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{projectileEnergyLabPerNucleon, momentum});
+        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 mass = get_mass(idFragm);
-        HEPEnergyType const Ekin = sqrt(p3output.getSquaredNorm() + mass * mass) - mass;
+        HEPEnergyType const Ekin =
+            calculate_kinetic_energy(p3output.getNorm(), nucleonMass);
 
         CORSIKA_LOG_DEBUG(
             "secondary fragment> id= {}"
@@ -247,24 +250,19 @@ namespace corsika::qgsjetII {
           }
         }
 
-        HEPMassType const nucleusMass = Proton::mass * Z + Neutron::mass * (A - Z);
+        Code const idFragm = get_nucleus_code(A, Z);
+        HEPEnergyType const mass = get_mass(idFragm);
         // no pT, frgments just go forward
-        HEPEnergyType const projectileEnergyLabPerNucleon = Elab / beamA;
-        MomentumVector momentum{
-            csPrime,
-            {0.0_GeV, 0.0_GeV,
-             calculate_momentum(projectileEnergyLabPerNucleon * A, nucleusMass)}};
+        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{projectileEnergyLabPerNucleon * A, momentum});
+        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
 
-        Code const idFragm = get_nucleus_code(A, Z);
-        HEPEnergyType const mass = get_mass(idFragm);
         HEPEnergyType const Ekin = calculate_kinetic_energy(p3output.getNorm(), mass);
 
         CORSIKA_LOG_DEBUG(

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