corsika/modules/fluka/InteractionModel.hpp

+/*
+ * (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 <vector>
+#include <utility>
+#include <memory>
+
+#include <corsika/media/Environment.hpp>
+#include <corsika/media/NuclearComposition.hpp>
+#include <corsika/framework/geometry/FourVector.hpp>
+#include <corsika/framework/utility/COMBoost.hpp>
+#include <corsika/framework/core/Logging.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/random/RNGManager.hpp>
+
+namespace corsika::fluka {
+  /**
+   * This class exposes the (hadronic) interactions of FLUKA. FLUKA needs to be
+   * initialized with a predefined set of target materials and a flag describing the type
+   * of interactions (elastic, inelastic, electromagnetic dissociation). Currently, only
+   * inelastic events are supported.
+   */
+  class InteractionModel {
+  public:
+    /**
+     * Create a new InteractionModel. The FLUKA materials are collected from the elements
+     * present in the environment. Each element is its own FLUKA material, no FLUKA
+     * compounds are used.
+     */
+    InteractionModel(std::set<Code> const&);
+
+    //! Return the cross-section of a given combination of projectile/target.
+    CrossSectionType getCrossSection(Code projectileId, Code targetId,
+                                     FourMomentum const& projectileP4,
+                                     FourMomentum const& targetP4) const;
+
+    bool isValid(Code projectileID, Code targetID, HEPEnergyType sqrtS) const;
+    bool isValid(Code projectileID, int material, HEPEnergyType sqrtS) const;
+
+    //! convert target Code to FLUKA material number
+    int getMaterialIndex(Code targetID) const;
+
+    /**
+     * Perform an interaction. Since FLUKA expects a fixed-target configuration, we
+     * perform a Lorentz transform into the rest frame of the target.
+     */
+    template <typename TSecondaryView>
+    void doInteraction(TSecondaryView& view, Code const projectileId, Code const targetId,
+                       FourMomentum const& projectileP4, FourMomentum const& targetP4);
+
+  private:
+    default_prng_type& RNG_ = RNGManager<>::getInstance().getRandomStream("fluka");
+    std::vector<std::pair<Code, int>> const
+        materials_; //!< map target Code to FLUKA material no.
+    std::shared_ptr<spdlog::logger> logger_ = get_logger("corsika_FLUKA_Interaction");
+    std::unique_ptr<double[]> cumsgx_; //!< dump for evtxyz cumsg*, never read again
+
+    static std::vector<std::pair<Code, int>> genFlukaMaterials(std::set<Code> const&);
+  };
+
+  inline static int const iflxyz_ = 1;
+} // namespace corsika::fluka
+
+#include <corsika/detail/modules/fluka/InteractionModel.inl>

corsika/modules/fluka/ParticleConversion.hpp

+/*
+ * (c) Copyright 2023 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 <array>
+#include <cstdint>
+#include <type_traits>
+#include <vector>
+
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <fmt/core.h>
+
+#include <string>
+
+namespace corsika::fluka {
+
+  enum class FLUKACode : int;
+  using FLUKACodeIntType = std::underlying_type<FLUKACode>::type;
+
+#include <corsika/modules/fluka/Generated.inc>
+
+  FLUKACode constexpr convertToFluka(Code const c8id) {
+    if (is_nucleus(c8id)) {
+      throw std::runtime_error{"nucleus conversion to FLUKA not implemented"};
+    }
+
+    FLUKACode const flukaID = corsika2fluka[static_cast<corsika::CodeIntType>(c8id)];
+    if (flukaID == FLUKACode::Unknown) {
+      throw std::runtime_error{
+          fmt::format("no correspondig FLUKA id for {}", get_name(c8id)).c_str()};
+    }
+
+    return flukaID;
+  }
+
+  int constexpr convertToFlukaRaw(Code const code) {
+    return static_cast<FLUKACodeIntType>(convertToFluka(code));
+  }
+
+  bool const canInteract(Code const code) {
+    if (is_nucleus(code)) return false; // nuclei support not yet implemented
+    return flukaCanInteract[static_cast<CodeIntType>(code)];
+  }
+} // namespace corsika::fluka

corsika/modules/proposal/ContinuousProcess.hpp

 /*
  * (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
@@ -14,8 +13,11 @@
 
 #include <corsika/framework/core/ParticleProperties.hpp>
 #include <corsika/framework/process/ContinuousProcess.hpp>
+#include <corsika/framework/process/ProcessReturn.hpp>
 #include <corsika/framework/random/RNGManager.hpp>
 #include <corsika/framework/random/UniformRealDistribution.hpp>
+#include <corsika/modules/writers/WriterOff.hpp>
+#include <corsika/framework/core/Step.hpp>
 
 #include <corsika/modules/proposal/ProposalProcessBase.hpp>
 
@@ -28,31 +30,32 @@ namespace corsika::proposal {
   //! use of interpolation tables which are runtime intensive calculation, but can be
   //! reused by setting the \param PROPOSAL::InterpolationDef::path_to_tables variable.
   //!
+  template <typename TOutput = WriterOff>
   class ContinuousProcess
-      : public corsika::ContinuousProcess<proposal::ContinuousProcess>,
-        ProposalProcessBase {
+      : public corsika::ContinuousProcess<proposal::ContinuousProcess<TOutput>>,
+        public ProposalProcessBase,
+        public TOutput {
 
-    enum { eDISPLACEMENT, eSCATTERING };
-    using calc_t = std::tuple<std::unique_ptr<PROPOSAL::Displacement>,
-                              std::unique_ptr<PROPOSAL::Scattering>>;
+    struct Calculator {
+      std::unique_ptr<PROPOSAL::Displacement> disp;
+      std::unique_ptr<PROPOSAL::multiple_scattering::Parametrization> scatter;
+    };
 
-    std::unordered_map<calc_key_t, calc_t, hash>
+    std::unordered_map<calc_key_t, Calculator, hash>
         calc; //!< Stores the displacement and scattering calculators.
 
-    HEPEnergyType energy_lost_ = 0 * electronvolt;
-
     //!
     //! Build the displacement and scattering calculators and add it to calc.
     //!
-    void buildCalculator(Code, NuclearComposition const&) final;
+    void buildCalculator(Code, size_t const&) final;
 
   public:
     //!
     //! Produces the continuous loss calculator for leptons based on nuclear
     //! compositions and stochastic description limited by the particle cut.
     //!
-    template <typename TEnvironment>
-    ContinuousProcess(TEnvironment const&);
+    template <typename TEnvironment, typename... TOutputArgs>
+    ContinuousProcess(TEnvironment const&, TOutputArgs&&...);
 
     //!
     //! Multiple Scattering of the lepton. Stochastic deflection is not yet taken into
@@ -60,7 +63,7 @@ namespace corsika::proposal {
     //! because of the constant referernce. The final direction will be updated anyway.
     //!
     template <typename TParticle>
-    void scatter(TParticle&, HEPEnergyType const&, GrammageType const&);
+    void scatter(Step<TParticle>&, HEPEnergyType const&, GrammageType const&);
 
     //!
     //! Produces the loss and deflection after given distance for the particle.
@@ -70,8 +73,8 @@ namespace corsika::proposal {
     //! \param limitFlag is true, if the track was actually limited by
     //! proposal::ContinuousProcess::getMaxStepLength
     //!
-    template <typename TParticle, typename TTrack>
-    ProcessReturn doContinuous(TParticle&, TTrack const& track, bool const limitFlag);
+    template <typename TParticle>
+    ProcessReturn doContinuous(Step<TParticle>&, bool const limitFlag);
 
     //!
     //! Calculates maximal step length of process.
@@ -79,9 +82,10 @@ namespace corsika::proposal {
     template <typename TParticle, typename TTrack>
     LengthType getMaxStepLength(TParticle const&, TTrack const&);
 
-    void showResults() const;
-    void reset();
-    HEPEnergyType getEnergyLost() const { return energy_lost_; }
+    /**
+     * Provide the config as YAML object to be stored on disk as output.
+     */
+    YAML::Node getConfig() const;
   };
 } // namespace corsika::proposal
 

corsika/modules/proposal/HadronicPhotonModel.hpp

+/*
+ * (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/ParticleProperties.hpp>
+#include <corsika/framework/geometry/FourVector.hpp>
+#include <corsika/framework/process/ProcessReturn.hpp>
+
+namespace corsika::proposal {
+
+  //! Implements the production of secondary hadrons for the hadronic interaction of real
+  //! and virtual photons for PROPOSAL. The model distinguishes between resonance
+  //! production (LE) and continuum (HE). HE production replaces the photon with a rho0.
+  //! External models are needed that implement the hadronic particle production via the
+  //! doInteraction(TSecondaries& view, Code const projectile, Code const target,
+  //! FourMomentum const& projectileP4, FourMomentum const& targetP4) routine. The
+  //! threshold between LE and HE interactions is defined in lab energy.
+  //! @tparam THadronicModel
+
+  template <class THadronicLEModel, class THadronicHEModel>
+  class HadronicPhotonModel {
+  public:
+    HadronicPhotonModel(THadronicLEModel&, THadronicHEModel&, HEPEnergyType const&);
+    //!
+    //! Calculate produce the hadronic secondaries in a hadronic photon interaction and
+    //! store them on the particle stack.
+    //!
+    template <typename TSecondaryView>
+    ProcessReturn doHadronicPhotonInteraction(TSecondaryView&, CoordinateSystemPtr const&,
+                                              FourMomentum const&, Code const&);
+
+  private:
+    inline static auto logger_{get_logger("corsika_proposal_HadronicPhotonModel")};
+    THadronicLEModel& leHadronicInteraction_;
+    THadronicHEModel& heHadronicInteraction_;
+    //! threshold for high energy hadronic interaction model. Lab. energy per nucleon
+    HEPEnergyType const heHadronicModelThresholdLabNN_;
+  };
+} // namespace corsika::proposal
+
+#include <corsika/detail/modules/proposal/HadronicPhotonModel.inl>
\ No newline at end of file

corsika/modules/proposal/Interaction.hpp

-/*
- * (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 <PROPOSAL/PROPOSAL.h>
-
-#include <corsika/framework/core/ParticleProperties.hpp>
-#include <corsika/framework/process/InteractionProcess.hpp>
-#include <corsika/framework/random/RNGManager.hpp>
-#include <corsika/framework/random/UniformRealDistribution.hpp>
-
-#include <corsika/media/Environment.hpp>
-
-#include <corsika/modules/proposal/ProposalProcessBase.hpp>
-
-#include <array>
-
-namespace corsika::proposal {
-
-  //!
-  //! Electro-magnetic and photon stochastic losses produced by proposal. It makes
-  //! use of interpolation tables which are runtime intensive calculation, but can be
-  //! reused by setting the \param PROPOSAL::InterpolationDef::path_to_tables variable.
-  //!
-  class Interaction : public InteractionProcess<Interaction>, ProposalProcessBase {
-
-    enum { eSECONDARIES, eINTERACTION };
-    using calculator_t = tuple<unique_ptr<PROPOSAL::SecondariesCalculator>,
-                               unique_ptr<PROPOSAL::Interaction>>;
-
-    std::unordered_map<calc_key_t, calculator_t, hash>
-        calc; //!< Stores the secondaries and interaction calculators.
-
-    //!
-    //! Build the secondaries and interaction calculators and add it to calc.
-    //!
-    void buildCalculator(Code, NuclearComposition const&) final;
-
-  public:
-    //!
-    //! Produces the stoachastic loss calculator for leptons based on nuclear
-    //! compositions and stochastic description limited by the particle cut.
-    //!
-    template <typename TEnvironment>
-    Interaction(TEnvironment const& env);
-
-    //!
-    //! Calculate the rates for the different targets and interactions. Sample a
-    //! pair of interaction-type, component and rate, followed by sampling a loss and
-    //! produce the corresponding secondaries and store them on the particle stack.
-    //!
-    template <typename TSecondaryView>
-    ProcessReturn doInteraction(TSecondaryView&);
-
-    //!
-    //! Calculates the  mean free path length
-    //!
-    template <typename TParticle>
-    GrammageType getInteractionLength(TParticle const& p);
-  };
-} // namespace corsika::proposal
-
-#include <corsika/detail/modules/proposal/Interaction.inl>

corsika/modules/proposal/InteractionModel.hpp

+/*
+ * (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 <PROPOSAL/PROPOSAL.h>
+
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <corsika/framework/process/InteractionProcess.hpp>
+#include <corsika/framework/process/ProcessReturn.hpp>
+#include <corsika/framework/geometry/FourVector.hpp>
+#include <corsika/framework/random/RNGManager.hpp>
+#include <corsika/framework/random/UniformRealDistribution.hpp>
+#include <corsika/modules/proposal/ProposalProcessBase.hpp>
+#include <corsika/modules/proposal/HadronicPhotonModel.hpp>
+
+namespace corsika::proposal {
+
+  //!
+  //! Electro-magnetic and photon stochastic losses produced by proposal. It makes
+  //! use of interpolation tables which are runtime intensive calculation, but can be
+  //! reused by setting the \param PROPOSAL::InterpolationDef::path_to_tables variable.
+  //! Hadroninc interactions of photons with nuclei are included. The cross section is
+  //! calculated by PROPOSAL. For the production of hadronic secondaries an external model
+  //! is needed that implements the
+  //! doInteraction(TSecondaries& view, Code const projectile, Code const
+  //! target,FourMomentum const& projectileP4, FourMomentum const& targetP4) routine.
+  //! @tparam THadronicModel
+  //!
+
+  template <class THadronicLEModel, class THadronicHEModel>
+  class InteractionModel
+      : public ProposalProcessBase,
+        public HadronicPhotonModel<THadronicLEModel, THadronicHEModel> {
+
+    enum { eSECONDARIES, eINTERACTION, eLPM_SUPPRESSION };
+    struct LPM_calculator;
+    using calculator_t = std::tuple<std::unique_ptr<PROPOSAL::SecondariesCalculator>,
+                                    std::unique_ptr<PROPOSAL::Interaction>,
+                                    std::unique_ptr<LPM_calculator>>;
+
+    std::unordered_map<calc_key_t, calculator_t, hash>
+        calc_; //!< Stores the secondaries and interaction calculators.
+
+    //!
+    //! Build the secondaries and interaction calculators and add it to calc.
+    //!
+    void buildCalculator(Code, size_t const&) final;
+
+    inline static auto logger_{get_logger("corsika_proposal_InteractionModel")};
+
+    // Calculators for the LPM effect
+    struct LPM_calculator {
+      std::unique_ptr<PROPOSAL::crosssection::PhotoPairLPM> photo_pair_lpm_ = nullptr;
+      std::unique_ptr<PROPOSAL::crosssection::BremsLPM> brems_lpm_ = nullptr;
+      std::unique_ptr<PROPOSAL::crosssection::EpairLPM> epair_lpm_ = nullptr;
+
+      const double mass_density_baseline_; // base mass density, note PROPOSAL units
+      const double particle_mass_;         // particle mass, note PROPOSAL units
+
+      LPM_calculator(const PROPOSAL::Medium& medium, const Code code,
+                     std::vector<PROPOSAL::InteractionType> inter_types)
+          : mass_density_baseline_(medium.GetMassDensity())
+          , particle_mass_(particle[code].mass) {
+        // at the moment, we always calculate the LPM effect based on the default cross
+        // sections. one could check for the parametrizations used in the other
+        // calculators.
+        if (std::find(inter_types.begin(), inter_types.end(),
+                      PROPOSAL::InteractionType::Photopair) != inter_types.end())
+          photo_pair_lpm_ = std::make_unique<PROPOSAL::crosssection::PhotoPairLPM>(
+              particle[code], medium, PROPOSAL::crosssection::PhotoPairKochMotz());
+        if (std::find(inter_types.begin(), inter_types.end(),
+                      PROPOSAL::InteractionType::Brems) != inter_types.end())
+          brems_lpm_ = std::make_unique<PROPOSAL::crosssection::BremsLPM>(
+              particle[code], medium, PROPOSAL::crosssection::BremsElectronScreening());
+        if (std::find(inter_types.begin(), inter_types.end(),
+                      PROPOSAL::InteractionType::Epair) != inter_types.end())
+          epair_lpm_ =
+              std::make_unique<PROPOSAL::crosssection::EpairLPM>(particle[code], medium);
+      }
+    };
+
+    //!
+    //! Checks whether the interaction is suppressed by the LPM effect
+    //!
+    bool CheckForLPM(const LPM_calculator&, const HEPEnergyType,
+                     const PROPOSAL::InteractionType,
+                     const std::vector<PROPOSAL::ParticleState>&, const MassDensityType,
+                     const PROPOSAL::Component&, const double v);
+
+  public:
+    //!
+    //! Produces the stoachastic loss calculator for leptons based on nuclear
+    //! compositions and stochastic description limited by the particle cut.
+    //!
+    template <typename TEnvironment>
+    InteractionModel(TEnvironment const& env, THadronicLEModel&, THadronicHEModel&,
+                     HEPEnergyType const&);
+
+    //!
+    //! Calculate the rates for the different targets and interactions. Sample a
+    //! pair of interaction-type, component and rate, followed by sampling a loss and
+    //! produce the corresponding secondaries and store them on the particle stack.
+    //! interactions in PROPOSAL are:
+    //!
+    //! InteractionType::Particle
+    //! InteractionType::Brems
+    //! InteractionType::Ioniz
+    //! InteractionType::Epair
+    //! InteractionType::Photonuclear
+    //! InteractionType::MuPair
+    //! InteractionType::Hadrons
+    //! InteractionType::ContinuousEnergyLoss
+    //! InteractionType::WeakInt
+    //! InteractionType::Compton
+    //! InteractionType::Decay
+    //! InteractionType::Annihilation
+    //! InteractionType::Photopair
+    //! InteractionType::Photoproduction
+    //! InteractionType::Photoeffect
+    //!
+    //! more information can be found at:
+    //! https://github.com/tudo-astroparticlephysics/PROPOSAL
+    template <typename TSecondaryView>
+    ProcessReturn doInteraction(TSecondaryView&, Code const projectileId,
+                                FourMomentum const& projectileP4);
+
+    //!
+    //! Calculates and returns the cross section.
+    //!
+    template <typename TParticle>
+    CrossSectionType getCrossSection(TParticle const& p, Code const projectileId,
+                                     FourMomentum const& projectileP4);
+  };
+
+} // namespace corsika::proposal
+
+#include <corsika/detail/modules/proposal/InteractionModel.inl>

corsika/modules/proposal/ProposalProcessBase.hpp

 /*
  * (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
@@ -13,8 +12,6 @@
 #include <corsika/framework/core/ParticleProperties.hpp>
 #include <corsika/framework/random/RNGManager.hpp>
 
-#include <corsika/setup/SetupEnvironment.hpp>
-
 #include <array>
 
 namespace corsika::proposal {
@@ -28,6 +25,20 @@ namespace corsika::proposal {
       Code::MuPlus, Code::TauPlus,  Code::TauMinus,
   };
 
+  //!
+  //! Relative EnergyCut value for PROPOSAL. Is defines that PROPOSAL will treat energy
+  //! losses with a relative energy loss larger than v_cut stochastically.
+  //!
+  static constexpr double v_cut = 0.01;
+
+  //!
+  //! List of EnergyCut values for which CORSIKA will, by default, create and provide
+  //! PROPOSAL tables.
+  //!
+  static constexpr std::array<HEPEnergyType, 10> energycut_table_values{
+      1000_MeV, 100_MeV, 20_MeV,   10_MeV,   3_MeV,
+      1_MeV,    0.4_MeV, 0.25_MeV, 0.15_MeV, 0.05_MeV};
+
   //!
   //! Internal map from particle codes to particle properties required for
   //! crosssections, decay and scattering algorithms. In the future the
@@ -42,25 +53,162 @@ namespace corsika::proposal {
   //!
   //! Crosssection factories for different particle types.
   //!
+
+  // Cross sections for muons and taus
   template <typename T>
-  static auto cross_builder =
-      [](PROPOSAL::Medium& m,
+  auto cross_builder = [](PROPOSAL::Medium& medium,
+                          corsika::units::si::HEPEnergyType
+                              emCut) { //!< Stochastic losses smaller than the given cut
+                                       //!< will be handeled continuously.
+    auto particle_def = T();
+    CORSIKA_LOG_DEBUG("PROPOSAL: Generate general cross sections for particle type {} ",
+                      particle_def.name);
+    auto interpolate = true;
+
+    auto p_cut =
+        std::make_shared<const PROPOSAL::EnergyCutSettings>(emCut / 1_MeV, v_cut, false);
+
+    // do not use v_cut for ionization due to numerical problems
+    auto p_cut_no_vcut =
+        std::make_shared<const PROPOSAL::EnergyCutSettings>(emCut / 1_MeV, 1, false);
+
+    auto cross_vec = std::vector<std::shared_ptr<PROPOSAL::CrossSectionBase>>();
+    cross_vec.push_back(PROPOSAL::make_crosssection(
+        PROPOSAL::crosssection::BremsKelnerKokoulinPetrukhin{false}, particle_def, medium,
+        p_cut, interpolate));
+    cross_vec.push_back(PROPOSAL::make_crosssection(
+        PROPOSAL::crosssection::EpairKelnerKokoulinPetrukhin{false}, particle_def, medium,
+        p_cut, interpolate));
+    cross_vec.push_back(PROPOSAL::make_crosssection(
+        PROPOSAL::crosssection::IonizBetheBlochRossi{*p_cut_no_vcut}, particle_def,
+        medium, p_cut_no_vcut, interpolate));
+    cross_vec.push_back(PROPOSAL::make_crosssection(
+        PROPOSAL::crosssection::PhotoAbramowiczLevinLevyMaor97{
+            std::make_unique<PROPOSAL::crosssection::ShadowButkevichMikheyev>()},
+        particle_def, medium, p_cut, interpolate));
+
+    return cross_vec;
+  };
+
+  // cross sections for photons
+  template <>
+  auto cross_builder<PROPOSAL::GammaDef> =
+      [](PROPOSAL::Medium& medium,
+         corsika::units::si::HEPEnergyType) { //!< Only-stochastic propagation
+        auto particle_def = PROPOSAL::GammaDef();
+        CORSIKA_LOG_DEBUG(
+            "PROPOSAL: Generate photon cross sections for particle type {} ",
+            particle_def.name);
+
+        auto interpolate = true;
+        auto cross_vec = std::vector<std::shared_ptr<PROPOSAL::CrossSectionBase>>();
+        auto photopair =
+            PROPOSAL::make_crosssection(PROPOSAL::crosssection::PhotoPairKochMotz{false},
+                                        particle_def, medium, nullptr, interpolate);
+        cross_vec.push_back(std::move(photopair));
+        auto compton =
+            PROPOSAL::make_crosssection(PROPOSAL::crosssection::ComptonKleinNishina{},
+                                        particle_def, medium, nullptr, interpolate);
+        cross_vec.push_back(std::move(compton));
+        auto photoproduction = PROPOSAL::make_crosssection(
+            PROPOSAL::crosssection::PhotoproductionHeckC7Shadowing{}, particle_def,
+            medium, nullptr, interpolate);
+        cross_vec.push_back(std::move(photoproduction));
+        auto photoeffect =
+            PROPOSAL::make_crosssection(PROPOSAL::crosssection::PhotoeffectSauter{},
+                                        particle_def, medium, nullptr, interpolate);
+        cross_vec.push_back(std::move(photoeffect));
+        auto muonpair = PROPOSAL::make_crosssection(
+            PROPOSAL::crosssection::PhotoMuPairBurkhardtKelnerKokoulin{}, particle_def,
+            medium, nullptr, interpolate);
+        cross_vec.push_back(std::move(muonpair));
+        return cross_vec;
+      };
+
+  // cross sections for electrons
+  template <>
+  auto cross_builder<PROPOSAL::EMinusDef> =
+      [](PROPOSAL::Medium& medium,
          corsika::units::si::HEPEnergyType
              emCut) { //!< Stochastic losses smaller than the given cut
                       //!< will be handeled continuously.
-        using namespace corsika::units::si;
-        auto p_cut =
-            std::make_shared<const PROPOSAL::EnergyCutSettings>(emCut / 1_MeV, 1, true);
-        return PROPOSAL::DefaultCrossSections<T>::template Get<std::false_type>(
-            T(), m, p_cut, true);
+        auto particle_def = PROPOSAL::EMinusDef();
+        CORSIKA_LOG_DEBUG(
+            "PROPOSAL: Generate electron cross sections for particle type {} ",
+            particle_def.name);
+        auto interpolate = true;
+
+        auto p_cut = std::make_shared<const PROPOSAL::EnergyCutSettings>(emCut / 1_MeV,
+                                                                         v_cut, false);
+
+        // do not use v_cut for ionization due to numerical problems
+        auto p_cut_no_vcut =
+            std::make_shared<const PROPOSAL::EnergyCutSettings>(emCut / 1_MeV, 1, false);
+
+        auto cross_vec = std::vector<std::shared_ptr<PROPOSAL::CrossSectionBase>>();
+        cross_vec.push_back(PROPOSAL::make_crosssection(
+            PROPOSAL::crosssection::BremsElectronScreening{false}, particle_def, medium,
+            p_cut, interpolate));
+        cross_vec.push_back(PROPOSAL::make_crosssection(
+            PROPOSAL::crosssection::EpairForElectronPositron{false}, particle_def, medium,
+            p_cut, interpolate));
+        cross_vec.push_back(PROPOSAL::make_crosssection(
+            PROPOSAL::crosssection::IonizBergerSeltzerMoller{*p_cut_no_vcut},
+            particle_def, medium, p_cut_no_vcut, interpolate));
+        cross_vec.push_back(PROPOSAL::make_crosssection(
+            PROPOSAL::crosssection::PhotoAbramowiczLevinLevyMaor97{
+                std::make_unique<PROPOSAL::crosssection::ShadowButkevichMikheyev>()},
+            particle_def, medium, p_cut, interpolate));
+
+        return cross_vec;
+      };
+
+  // cross sections for positrons
+  template <>
+  auto cross_builder<PROPOSAL::EPlusDef> =
+      [](PROPOSAL::Medium& medium,
+         corsika::units::si::HEPEnergyType
+             emCut) { //!< Stochastic losses smaller than the given cut
+                      //!< will be handeled continuously.
+        auto particle_def = PROPOSAL::EPlusDef();
+        CORSIKA_LOG_DEBUG(
+            "PROPOSAL: Generate positron cross sections for particle type {} ",
+            particle_def.name);
+        auto interpolate = true;
+
+        auto p_cut = std::make_shared<const PROPOSAL::EnergyCutSettings>(emCut / 1_MeV,
+                                                                         v_cut, false);
+
+        // do not use v_cut for ionization due to numerical problems
+        auto p_cut_no_vcut =
+            std::make_shared<const PROPOSAL::EnergyCutSettings>(emCut / 1_MeV, 1, false);
+
+        auto cross_vec = std::vector<std::shared_ptr<PROPOSAL::CrossSectionBase>>();
+        cross_vec.push_back(PROPOSAL::make_crosssection(
+            PROPOSAL::crosssection::BremsElectronScreening{false}, particle_def, medium,
+            p_cut, interpolate));
+        cross_vec.push_back(PROPOSAL::make_crosssection(
+            PROPOSAL::crosssection::EpairForElectronPositron{false}, particle_def, medium,
+            p_cut, interpolate));
+        cross_vec.push_back(PROPOSAL::make_crosssection(
+            PROPOSAL::crosssection::IonizBergerSeltzerBhabha{*p_cut_no_vcut},
+            particle_def, medium, p_cut_no_vcut, interpolate));
+        cross_vec.push_back(PROPOSAL::make_crosssection(
+            PROPOSAL::crosssection::PhotoAbramowiczLevinLevyMaor97{
+                std::make_unique<PROPOSAL::crosssection::ShadowButkevichMikheyev>()},
+            particle_def, medium, p_cut, interpolate));
+        cross_vec.push_back(
+            PROPOSAL::make_crosssection(PROPOSAL::crosssection::AnnihilationHeitler{},
+                                        particle_def, medium, nullptr, interpolate));
+
+        return cross_vec;
       };
 
   //!
   //! PROPOSAL default crosssections are maped to corresponding corsika particle
   //! code.
   //!
-  static std::map<Code, std::function<PROPOSAL::crosssection_list_t<PROPOSAL::ParticleDef,
-                                                                    PROPOSAL::Medium>(
+  static std::map<Code, std::function<PROPOSAL::crosssection_list_t(
                             PROPOSAL::Medium&, corsika::units::si::HEPEnergyType)>>
       cross = {{Code::Photon, cross_builder<PROPOSAL::GammaDef>},
                {Code::Electron, cross_builder<PROPOSAL::EMinusDef>},
@@ -76,23 +224,34 @@ namespace corsika::proposal {
   //!
   class ProposalProcessBase {
   protected:
-    RNGManager::prng_type RNG_; //!< random number generator used by proposal
+    default_prng_type& RNG_ = RNGManager<>::getInstance().getRandomStream("proposal");
 
     std::unordered_map<std::size_t, PROPOSAL::Medium>
         media; //!< maps nuclear composition from univers to media to produce
                //!< crosssections, which requires further ionization constants.
 
+    //!< save emcut for tracked particles
+    std::unordered_map<Code, corsika::units::si::HEPEnergyType>
+        proposal_energycutsettings;
+
     //!
     //! Store cut and  nuclear composition of the whole universe in media which are
     //! required for creating crosssections by proposal.
     //!
-    ProposalProcessBase(corsika::setup::Environment const& _env);
+    template <typename TEnvironment>
+    ProposalProcessBase(TEnvironment const& _env);
 
     //!
     //! Checks if a particle can be processed by proposal
     //!
     bool canInteract(Code pcode) const;
 
+    //!
+    //! Finds the optimal EnergyCut for which PROPOSAL tables should (by default) be
+    //! available.
+    //!
+    HEPEnergyType getOptimizedEmCut(Code code) const;
+
     using calc_key_t = std::pair<std::size_t, Code>;
 
     //!
@@ -106,7 +265,12 @@ namespace corsika::proposal {
     //!
     //! Builds the calculator to the corresponding class
     //!
-    virtual void buildCalculator(Code, NuclearComposition const&) = 0;
+    virtual void buildCalculator(Code, size_t const&) = 0;
+
+    //!
+    //! Initialize PROPOSAL tables for given medium, code, and energy cut
+    //!
+    void buildTables(PROPOSAL::Medium, Code, HEPEnergyType);
 
     //!
     //! Searches the particle dependet calculator dependent of actuall medium composition
@@ -118,7 +282,7 @@ namespace corsika::proposal {
       const auto& comp = vP.getNode()->getModelProperties().getNuclearComposition();
       auto calc_it = calc.find(std::make_pair(comp.getHash(), vP.getPID()));
       if (calc_it != calc.end()) return calc_it;
-      buildCalculator(vP.getPID(), comp);
+      buildCalculator(vP.getPID(), comp.getHash());
       return getCalculator(vP, calc);
     }
   };

corsika/modules/pythia8/Decay.hpp

 /*
  * (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/modules/pythia8/Interaction.hpp

-/*
- * (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/framework/core/ParticleProperties.hpp>
-#include <corsika/framework/core/PhysicalUnits.hpp>
-#include <corsika/framework/random/RNGManager.hpp>
-#include <corsika/framework/process/InteractionProcess.hpp>
-#include <corsika/modules/pythia8/Pythia8.hpp>
-
-#include <corsika/setup/SetupStack.hpp>
-
-#include <tuple>
-
-namespace corsika::pythia8 {
-
-  class Interaction : public InteractionProcess<Interaction>, public Pythia8::Pythia {
-
-  public:
-    Interaction(const bool print_listing = false);
-    ~Interaction();
-
-    void setStable(std::vector<Code> const&);
-    void setUnstable(const Code);
-    void setStable(const Code);
-
-    bool isValidCoMEnergy(HEPEnergyType ecm) { return (10_GeV < ecm) && (ecm < 1_PeV); }
-
-    bool canInteract(const Code);
-    void configureLabFrameCollision(const Code, const Code, const HEPEnergyType);
-
-    std::tuple<CrossSectionType, CrossSectionType> getCrossSection(
-        const Code BeamId, const Code TargetId, const HEPEnergyType CoMenergy);
-
-    GrammageType getInteractionLength(corsika::setup::Stack::particle_type const&);
-
-    /**
-       In this function PYTHIA is called to produce one event. The
-       event is copied (and boosted) into the shower lab frame.
-     */
-    template <typename TView>
-    void doInteraction(TView&);
-
-  private:
-    default_prng_type& RNG_ = RNGManager::getInstance().getRandomStream("pythia");
-    Pythia8::SigmaTotal sigma_;
-    const bool internalDecays_ = true;
-    int count_ = 0;
-    bool print_listing_ = false;
-  };
-
-} // namespace corsika::pythia8
-
-#include <corsika/detail/modules/pythia8/Interaction.inl>

corsika/modules/pythia8/InteractionModel.hpp

+/*
+ * (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 <tuple>
+#include <unordered_map>
+#include <boost/filesystem/path.hpp>
+
+#include <corsika/framework/utility/CorsikaData.hpp>
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/random/RNGManager.hpp>
+#include <corsika/framework/utility/CrossSectionTable.hpp>
+#include <corsika/framework/process/InteractionProcess.hpp>
+#include <corsika/modules/pythia8/Pythia8.hpp>
+
+namespace corsika::pythia8 {
+
+  /**
+   * @brief Defines the interface to the PYTHIA8 interaction model. Configured for
+   * hadron-nucleus interactions using Angantyr.
+   *
+   * This is a TModel argument for InteractionProcess<TModel>.
+   */
+
+  class InteractionModel {
+
+  public:
+    /**
+     * Constructs the interface for PYTHIA8
+     *
+     * @param stableParticles - set of particle ids to be treated as stable inside
+     * pythia. if empty then all particles are treated as stable (no decays!)
+     * @param dataPath path to the pythia tables
+     * @param printListing switch on/off the pythia printout
+     */
+    InteractionModel(std::set<Code> const& stableParticles = {},
+                     boost::filesystem::path const& dataPath = corsika_data("Pythia"),
+                     bool const printListing = false);
+    ~InteractionModel();
+
+    bool canInteract(Code const) const;
+
+    /**
+     * Check if configuration of projectile and target is valid.
+     *
+     * @param projectileId is the Code of the projectile
+     * @param targetId is the Code of the target
+     * @param sqrtS is the squared sum of the 4-momenta of projectile and target
+     */
+    bool isValid(Code const projectileId, Code const targetId,
+                 HEPEnergyType const sqrtS) const;
+
+    /**
+     * Returns inelastic AND elastic cross sections.
+     *
+     * These cross sections must correspond to the process described in doInteraction
+     * (sigma_tot = sigma_inel + sigma_el). Allowed targets are:
+     * nuclei or single nucleons (p,n,hydrogen).
+     *
+     * @param projectile is the Code of the projectile
+     * @param target is the Code of the target
+     * @param projectileP4 is the 4-momentum of the projectile
+     * @param targetP4 is the 4-momentum of the target
+     *
+     * @return a tuple of: inelastic cross section, elastic cross section
+     */
+    std::tuple<CrossSectionType, CrossSectionType> getCrossSectionInelEla(
+        Code const projectile, Code const target, FourMomentum const& projectileP4,
+        FourMomentum const& targetP4) const;
+
+    /**
+     * Returns inelastic (production) cross section.
+     *
+     * This cross section must correspond to the process described in doInteraction.
+     * Allowed targets are: nuclei or single nucleons (p,n,hydrogen).
+     *
+     * @param projectile is the Code of the projectile
+     * @param target is the Code of the target
+     * @param projectileP4 is the 4-momentum of the projectile
+     * @param targetP4 is the 4-momentum of the target
+     *
+     * @return inelastic cross section
+     * elastic cross section
+     */
+    CrossSectionType getCrossSection(Code const projectile, Code const target,
+                                     FourMomentum const& projectileP4,
+                                     FourMomentum const& targetP4) const;
+
+    /**
+     * In this function PYTHIA8 is called to produce one event. The
+     * event is copied (and boosted) into the reference frame defined by the input
+     * 4-momenta.
+     */
+    template <typename TView>
+    void doInteraction(TView& output, Code const projectileId, Code const targetId,
+                       FourMomentum const& projectileP4, FourMomentum const& targetP4);
+
+    using key_type = std::pair<corsika::Code, corsika::Code>;
+
+  private:
+    static auto constexpr GeV_mult = [](float x) { return x * 1_GeV; };
+    static auto constexpr millibarn_mult = [](float x) { return x * 1_mb; };
+
+    static CrossSectionTable<InterpolationTransforms::Log> loadPPTable(
+        boost::filesystem::path const& dataPath, char const* key);
+
+    default_prng_type& RNG_ = RNGManager<>::getInstance().getRandomStream("pythia");
+
+    /**
+     * All particles that we enable as projectiles. Code::Nucleus is not listed as it is
+     * handled separately. The list is unlikely to be exhaustive, but should cover the
+     * mainstream use-cases sufficiently well.
+     */
+
+    static std::array constexpr validProjectiles_ = {
+        Code::PiPlus,       Code::PiMinus,   Code::Pi0,           Code::Proton,
+        Code::AntiProton,   Code::Neutron,   Code::AntiNeutron,   Code::KPlus,
+        Code::KMinus,       Code::K0Long,    Code::K0Short,       Code::SigmaMinus,
+        Code::SigmaPlusBar, Code::SigmaPlus, Code::SigmaMinusBar, Code::Xi0,
+        Code::Xi0Bar};
+
+    static std::array constexpr validTargets_ = {
+        Code::Proton, Code::Carbon, Code::Nitrogen, Code::Oxygen, Code::Argon};
+
+    std::unordered_map<corsika::Code, corsika::Code> const xs_map_ = {
+        {Code::SigmaMinus, Code::SigmaPlus},
+        {Code::SigmaMinusBar, Code::SigmaPlus},
+        {Code::SigmaPlusBar, Code::SigmaPlus},
+        {Code::Xi0Bar, Code::Xi0}};
+
+    Pythia8::Pythia pythia_;
+
+    std::unordered_map<key_type, CrossSectionTable<InterpolationTransforms::Log>>
+        crossSectionTables_;
+
+    CrossSectionTable<InterpolationTransforms::Log> crossSectionPPElastic_,
+        crossSectionPPInelastic_;
+
+    bool const print_listing_ = false;
+    HEPEnergyType const eMaxLab_ =
+        1e21_eV; // Cross-section tables tabulated up to 10^21 eV
+    HEPEnergyType const eKinMinLab_ = 100_GeV;
+  };
+
+} // namespace corsika::pythia8
+
+#include <corsika/detail/modules/pythia8/InteractionModel.inl>

corsika/modules/pythia8/NeutrinoInteraction.hpp

+/*
+ * (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 <tuple>
+#include <boost/filesystem/path.hpp>
+
+#include <corsika/framework/utility/CorsikaData.hpp>
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/random/RNGManager.hpp>
+#include <corsika/framework/process/InteractionProcess.hpp>
+#include <corsika/modules/pythia8/Pythia8.hpp>
+
+namespace corsika::pythia8 {
+
+  using HEPEnergyTypeSqr = decltype(1_GeV * 1_GeV);
+
+  /**
+   * @brief Defines the interface to PYTHIA8. Configured for
+   * DIS neutrino - nucleon interactions.
+   *
+   */
+
+  class NeutrinoInteraction : public InteractionProcess<NeutrinoInteraction> {
+
+  public:
+    /**
+     * Constructs the interface for PYTHIA8 neutrino interactions
+     *
+     * @param listOfStableParticles - Switch on Pythia internal decay for all hadrons
+     * except those in the list
+     * @param handleNC - Switch on/off neutral current interactions
+     * @param handleCC - Switch on/off charged current interactions
+     */
+    NeutrinoInteraction(std::set<Code> const& listOfStableParticles,
+                        bool const& handleNC = true, bool const& handleCC = true);
+    ~NeutrinoInteraction();
+    /**
+     * Returns inelastic (production) cross section.
+     *
+     * This cross section must correspond to the process described in doInteraction.
+     * Allowed targets are: nuclei or single nucleons (p,n,hydrogen).
+     *
+     * @param projectile is the Code of the projectile
+     * @param target is the Code of the target
+     * @param projectileP4 is the four momentum of the projectile
+     * @param targetP4 is the four momentum of the target
+     *
+     * @return inelastic cross section
+     */
+    /* NOTE: the cross section for neutrino interactions is fixed to an arbitrary small
+     * value. This is needed just so the interaction process is not skipped in the process
+     * loop when the interaction is forced. If needed the exact value could be made
+     * accessible. */
+    CrossSectionType getCrossSection(Code const projectile, Code const target,
+                                     FourMomentum const& projectileP4,
+                                     FourMomentum const& targetP4) const {
+      if (isValid(projectile, target, projectileP4, targetP4))
+        return cross_section_;
+      else
+        return CrossSectionType::zero();
+    };
+
+    bool isValid(Code const projectileId, [[maybe_unused]] Code const targetId,
+                 FourMomentum const& projectileP4, FourMomentum const& targetP4) const {
+      auto const S = (projectileP4 + targetP4).getNormSqr();
+      return is_neutrino(projectileId) &&
+             (is_nucleus(targetId) || targetId == Code::Proton ||
+              targetId == Code::Neutron) &&
+             S >= minQ2_;
+    };
+
+    /**
+     * In this function PYTHIA is called to produce one event. The
+     * event is copied (and boosted) into the shower lab frame.
+     */
+    template <typename TView>
+    void doInteraction(TView& output, Code const projectileId, Code const targetId,
+                       FourMomentum const& projectileP4, FourMomentum const& targetP4);
+
+  private:
+    std::shared_ptr<spdlog::logger> logger_ =
+        get_logger("corsika_pythia8_NeutrinoInteraction");
+    CrossSectionType const cross_section_ = 4_nb;
+    int count_ = 0;
+    std::set<Code> const stable_particles_ = {};
+    bool const handle_nc_ = true;
+    bool const handle_cc_ = true;
+    HEPEnergyTypeSqr minQ2_ = 25_GeV * 1_GeV;
+    Pythia8::Pythia pythiaMain_;
+  };
+} // namespace corsika::pythia8
+
+#include <corsika/detail/modules/pythia8/NeutrinoInteraction.inl>
\ No newline at end of file

corsika/modules/pythia8/Pythia8.hpp

 /*
  * (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/modules/pythia8/Random.hpp

 /*
  * (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,7 +18,7 @@ namespace corsika::pythia8 {
 
   private:
     std::uniform_real_distribution<double> Dist_;
-    default_prng_type& RNG_ = RNGManager::getInstance().getRandomStream("pythia");
+    default_prng_type& RNG_ = RNGManager<>::getInstance().getRandomStream("pythia");
   };
 
 } // namespace corsika::pythia8

corsika/modules/qgsjetII/Interaction.hpp

-/*
- * (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/framework/core/ParticleProperties.hpp>
-#include <corsika/framework/core/PhysicalUnits.hpp>
-#include <corsika/framework/random/RNGManager.hpp>
-#include <corsika/framework/process/InteractionProcess.hpp>
-#include <corsika/modules/qgsjetII/ParticleConversion.hpp>
-#include <corsika/framework/utility/CorsikaData.hpp>
-
-#include <boost/filesystem/path.hpp>
-
-#include <qgsjet-II-04.hpp>
-
-#include <string>
-
-namespace corsika::qgsjetII {
-
-  class Interaction : public corsika::InteractionProcess<Interaction> {
-
-  public:
-    Interaction(boost::filesystem::path dataPath = corsika_data("QGSJetII"));
-    ~Interaction();
-
-    bool wasInitialized() { return initialized_; }
-    int getMaxTargetMassNumber() const { return maxMassNumber_; }
-    bool isValidTarget(corsika::Code TargetId) const {
-      return is_nucleus(TargetId) && (get_nucleus_A(TargetId) < maxMassNumber_);
-    }
-
-    CrossSectionType getCrossSection(const Code, const Code, const HEPEnergyType,
-                                     const unsigned int Abeam = 0,
-                                     const unsigned int Atarget = 0) const;
-
-    template <typename TParticle>
-    GrammageType getInteractionLength(TParticle const&) const;
-
-    /**
-       In this function QGSJETII is called to produce one event. The
-       event is copied (and boosted) into the shower lab frame.
-     */
-
-    template <typename TSecondaryView>
-    void doInteraction(TSecondaryView&);
-
-  private:
-    int count_ = 0;
-    bool initialized_ = false;
-    QgsjetIIHadronType alternate_ =
-        QgsjetIIHadronType::PiPlusType; // for pi0, rho0 projectiles
-
-    corsika::default_prng_type& rng_ =
-        corsika::RNGManager::getInstance().getRandomStream("qgsjet");
-    const int maxMassNumber_ = 208;
-  };
-
-} // namespace corsika::qgsjetII
-
-#include <corsika/detail/modules/qgsjetII/Interaction.inl>

corsika/modules/qgsjetII/InteractionModel.hpp

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/modules/qgsjetII/ParticleConversion.hpp>
+#include <qgsjet-II-04.hpp>
+
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/random/RNGManager.hpp>
+#include <corsika/framework/utility/COMBoost.hpp>
+#include <corsika/framework/utility/CorsikaData.hpp>
+
+#include <boost/filesystem/path.hpp>
+
+#include <random>
+
+namespace corsika::qgsjetII {
+
+  class InteractionModel {
+
+  public:
+    InteractionModel(boost::filesystem::path dataPath = corsika_data("QGSJetII"));
+    ~InteractionModel();
+
+    /**
+     * Throws exception if invalid system is passed.
+     *
+     * @param beamId
+     * @param targetId
+     */
+    bool isValid(Code const beamId, Code const targetId, HEPEnergyType const sqrtS) const;
+
+    /**
+     * Return the QGSJETII inelastic/production cross section.
+     *
+     * This cross section must correspond to the process described in doInteraction.
+     * Allowed targets are: nuclei or single nucleons (p,n,hydrogen).
+     *
+     * @param projectile is the Code of the projectile
+     * @param target is the Code of the target
+     * @param sqrtSnn is the center-of-mass energy (per nucleon pair)
+     * @param Aprojectile is the mass number of the projectils, if it is a nucleus
+     * @param Atarget is the mass number of the target, if it is a nucleus
+     *
+     * @return inelastic cross section.
+     */
+    CrossSectionType getCrossSection(Code const projectile, Code const target,
+                                     FourMomentum const& projectileP4,
+                                     FourMomentum const& targetP4) const;
+
+    /**
+     * Returns inelastic AND elastic cross sections.
+     *
+     * These cross sections must correspond to the process described in doInteraction
+     * AND elastic scattering (sigma_tot = sigma_inel + sigma_el).
+     *
+     * @param projectile is the Code of the projectile
+     * @param target is the Code of the target
+     * @param projectileP4: four-momentum of projectile
+     * @param targetP4: four-momentum of target
+     *
+     * @return a tuple of: inelastic cross section, elastic cross section
+     */
+    std::tuple<CrossSectionType, CrossSectionType> getCrossSectionInelEla(
+        Code const projectile, Code const target, FourMomentum const& projectileP4,
+        FourMomentum const& targetP4) const;
+
+    /**
+     * In this function QGSJETII is called to produce one event.
+     *
+     * The event is copied (and boosted) into the shower lab frame.
+     */
+    template <typename TSecondaries>
+    void doInteraction(TSecondaries&, Code const projectile, Code const target,
+                       FourMomentum const& projectileP4, FourMomentum const& targetP4);
+
+  private:
+    int count_ = 0;
+    QgsjetIIHadronType alternate_ =
+        QgsjetIIHadronType::PiPlusType; // for pi0, rho0 projectiles
+
+    corsika::default_prng_type& rng_ =
+        corsika::RNGManager<>::getInstance().getRandomStream("qgsjet");
+    std::bernoulli_distribution bernoulli_;
+    static size_t constexpr maxMassNumber_ = 208;
+    static HEPEnergyType constexpr sqrtSmin_ = 10_GeV;
+  };
+
+} // namespace corsika::qgsjetII
+
+#include <corsika/detail/modules/qgsjetII/InteractionModel.inl>

corsika/modules/qgsjetII/ParticleConversion.hpp

 /*
  * (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,13 +14,13 @@
 namespace corsika::qgsjetII {
 
   /**
-     These are the possible secondaries produced by QGSJetII
+   * These are the possible secondaries produced by QGSJetII.
    */
   enum class QgsjetIICode : int8_t;
   using QgsjetIICodeIntType = std::underlying_type<QgsjetIICode>::type;
 
   /**
-     These are the possible projectile for which QGSJetII knwos cross section
+   * These are the possible projectile for which QGSJetII knwos cross section.
    */
   enum class QgsjetIIXSClass : int8_t {
     CannotInteract = 0,
@@ -32,7 +31,7 @@ namespace corsika::qgsjetII {
   using QgsjetIIXSClassIntType = std::underlying_type<QgsjetIIXSClass>::type;
 
   /**
-     These are the only possible projectile types in QGSJetII
+   *  These are the only possible projectile types in QGSJetII.
    */
   enum class QgsjetIIHadronType : int8_t {
     UndefinedType = 0,
@@ -58,41 +57,44 @@ namespace corsika::qgsjetII {
 
 namespace corsika::qgsjetII {
 
-  QgsjetIICode constexpr convertToQgsjetII(Code pCode) {
+  QgsjetIICode constexpr convertToQgsjetII(Code const pCode) {
     return corsika2qgsjetII[static_cast<CodeIntType>(pCode)];
   }
 
-  Code constexpr convertFromQgsjetII(QgsjetIICode pCode) {
-    auto const pCodeInt = static_cast<QgsjetIICodeIntType>(pCode);
-    auto const corsikaCode = qgsjetII2corsika[pCodeInt - minQgsjetII];
+  Code constexpr convertFromQgsjetII(QgsjetIICode const code) {
+    auto const codeInt = static_cast<QgsjetIICodeIntType>(code);
+    auto const corsikaCode = qgsjetII2corsika[codeInt - minQgsjetII];
     if (corsikaCode == Code::Unknown) {
       throw std::runtime_error(std::string("QGSJETII/CORSIKA conversion of pCodeInt=")
-                                   .append(std::to_string(pCodeInt))
+                                   .append(std::to_string(codeInt))
                                    .append(" impossible"));
     }
     return corsikaCode;
   }
 
-  QgsjetIICodeIntType constexpr convertToQgsjetIIRaw(Code pCode) {
-    return static_cast<QgsjetIICodeIntType>(convertToQgsjetII(pCode));
+  QgsjetIICodeIntType constexpr convertToQgsjetIIRaw(Code const code) {
+    return static_cast<QgsjetIICodeIntType>(convertToQgsjetII(code));
   }
 
-  QgsjetIIXSClass constexpr getQgsjetIIXSCode(Code pCode) {
-    // if (pCode == corsika::particles::Code::Nucleus)
-    // static_cast(QgsjetIIXSClassIntType>();
-    return corsika2qgsjetIIXStype[static_cast<CodeIntType>(pCode)];
+  QgsjetIIXSClass constexpr getQgsjetIIXSCode(Code const code) {
+    return is_nucleus(code) ? QgsjetIIXSClass::Baryons
+                            : corsika2qgsjetIIXStype[static_cast<CodeIntType>(code)];
   }
 
-  QgsjetIIXSClassIntType constexpr getQgsjetIIXSCodeRaw(Code pCode) {
-    return static_cast<QgsjetIIXSClassIntType>(getQgsjetIIXSCode(pCode));
+  QgsjetIIXSClassIntType constexpr getQgsjetIIXSCodeRaw(Code const code) {
+    return is_nucleus(code)
+               ? static_cast<QgsjetIIXSClassIntType>(QgsjetIIXSClass::Baryons)
+               : static_cast<QgsjetIIXSClassIntType>(getQgsjetIIXSCode(code));
   }
 
-  bool constexpr canInteract(Code pCode) {
-    return getQgsjetIIXSCode(pCode) != QgsjetIIXSClass::CannotInteract;
+  bool constexpr canInteract(Code const code) {
+    if (is_nucleus(code)) return true;
+    return getQgsjetIIXSCode(code) != QgsjetIIXSClass::CannotInteract;
   }
 
-  QgsjetIIHadronType constexpr getQgsjetIIHadronType(Code pCode) {
-    return corsika2qgsjetIIHadronType[static_cast<CodeIntType>(pCode)];
+  QgsjetIIHadronType constexpr getQgsjetIIHadronType(Code const code) {
+    if (is_nucleus(code)) return QgsjetIIHadronType::NucleusType;
+    return corsika2qgsjetIIHadronType[static_cast<CodeIntType>(code)];
   }
 
 } // namespace corsika::qgsjetII

corsika/modules/qgsjetII/QGSJetIIFragmentsStack.hpp

 /*
  * (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
@@ -50,16 +49,16 @@ namespace corsika::qgsjetII {
     }
   };
 
-  template <typename StackIteratorInterface>
-  class FragmentsInterface : public corsika::ParticleBase<StackIteratorInterface> {
+  template <typename TStackIterator>
+  class FragmentsInterface : public corsika::ParticleBase<TStackIterator> {
 
-    using corsika::ParticleBase<StackIteratorInterface>::getStackData;
-    using corsika::ParticleBase<StackIteratorInterface>::getIndex;
+    using corsika::ParticleBase<TStackIterator>::getStackData;
+    using corsika::ParticleBase<TStackIterator>::getIndex;
 
   public:
     void setParticleData(const int vSize) { setFragmentSize(vSize); }
 
-    void setParticleData(FragmentsInterface<StackIteratorInterface>& /*parent*/,
+    void setParticleData(FragmentsInterface<TStackIterator>& /*parent*/,
                          const int vSize) {
       setFragmentSize(vSize);
     }

corsika/modules/qgsjetII/QGSJetIIStack.hpp

 /*
  * (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
@@ -44,17 +43,17 @@ namespace corsika::qgsjetII {
     void decrementSize();
   };
 
-  template <typename StackIteratorInterface>
-  class ParticleInterface : public corsika::ParticleBase<StackIteratorInterface> {
+  template <typename TStackIterator>
+  class ParticleInterface : public corsika::ParticleBase<TStackIterator> {
 
-    using corsika::ParticleBase<StackIteratorInterface>::getStackData;
-    using corsika::ParticleBase<StackIteratorInterface>::getIndex;
+    using corsika::ParticleBase<TStackIterator>::getStackData;
+    using corsika::ParticleBase<TStackIterator>::getIndex;
 
   public:
     void setParticleData(const int vID, const HEPEnergyType vE, const MomentumVector& vP,
                          const HEPMassType);
-    void setParticleData(ParticleInterface<StackIteratorInterface>& /*parent*/,
-                         const int vID, const HEPEnergyType vE, const MomentumVector& vP,
+    void setParticleData(ParticleInterface<TStackIterator>& /*parent*/, const int vID,
+                         const HEPEnergyType vE, const MomentumVector& vP,
                          const HEPMassType);
 
     void setEnergy(const HEPEnergyType v);

corsika/modules/qgsjetII/Random.hpp

-/*
- * (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/framework/random/RNGManager.hpp>
-#include <random>
-
-namespace qgsjetII {
-
-  double rndm_interface() {
-    static corsika::default_prng_type& rng =
-        corsika::RNGManager::getInstance().getRandomStream("qgsjet");
-    std::uniform_real_distribution<double> dist;
-    return dist(rng);
-  }
-
-} // namespace qgsjetII

corsika/modules/qgsjetII/qgsjet-II-04.hpp

 /*
  * (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
@@ -42,7 +41,7 @@ extern struct {
 } qgarr55_;
 
 /**
-   Small helper class to provide a data-directory name in the format qgsjetII expects
+ * Small helper class to provide a data-directory name in the format qgsjetII expects.
  */
 class datadir {
 private:
@@ -60,33 +59,33 @@ void qgaini_(
     const char* datdir); // Note: there is a length limiation 132 from fortran-qgsjet here
 
 /**
-   additional initialization procedure per event
-
-   @param e0n  - interaction energy (per hadron/nucleon),
-   @param icp0 - hadron type (+-1 - pi+-, +-2 - p(p~), +-3 - n(n~), +-4 - K+-, +-5 -
-   K_l/s),
-   @param iap  - projectile mass number (1 - for a hadron),
-   @param iat  - target mass number
-*/
+ * Additional initialization procedure per event.
+ *
+ * @param e0n  - interaction energy (per hadron/nucleon),
+ * @param icp0 - hadron type (+-1 - pi+-, +-2 - p(p~), +-3 - n(n~), +-4 - K+-, +-5 -
+ *               K_l/s),
+ * @param iap  - projectile mass number (1 - for a hadron),
+ * @param iat  - target mass number
+ */
 void qgini_(const double& e0n, const int& icp0, const int& iap, const int& iat);
 
 /**
-   generate one event configuration
-*/
+ * Generate one event configuration.
+ */
 void qgconf_();
 
 /**
-   hadron-nucleus (hadron-nucleus) particle production cross section
-
-   @param e0n lab. energy per projectile nucleon (hadron)
-   @param icz hadron class (1 - pion, 2 - nucleon, 3 - kaon)
-   @param iap projectile mass number (1=<iap<=iapmax),
-   @param iat target mass number     (1=<iat<=iapmax)
+ * Hadron-nucleus (hadron-nucleus) particle production cross section.
+ *
+ * @param e0n lab. energy per projectile nucleon (hadron)
+ * @param icz hadron class (1 - pion, 2 - nucleon, 3 - kaon)
+ * @param iap0 projectile mass number (1=<iap<=iapmax),
+ * @param iat0 target mass number     (1=<iat<=iapmax)
  */
 double qgsect_(const double& e0n, const int& icz, const int& iap0, const int& iat0);
 
 /**
-   link to random number generation
+ * Link to random number generation.
  */
 double qgran_(int&);
 }