CMakeLists.txt

@@ -253,6 +253,7 @@ set (CORSIKA_DATA_WITH_TEST ON) # we want to run the corsika-data unit test
 add_subdirectory (modules/data) # this is corsika-data (submodule)
 add_subdirectory (modules/pythia8)
 add_subdirectory (modules/sibyll)
+add_subdirectory (modules/sophia)
 add_subdirectory (modules/qgsjetII)
 add_subdirectory (modules/urqmd)
 add_subdirectory (modules/conex)

corsika/detail/framework/core/ParticleProperties.inl

@@ -147,7 +147,7 @@ namespace corsika {
     return get_mass(Code::Proton) * Z + (A - Z) * get_mass(Code::Neutron);
   }
 
-  inline std::string_view get_nucleus_name(Code const code) {
+  inline std::string get_nucleus_name(Code const code) {
     size_t const A = get_nucleus_A(code);
     size_t const Z = get_nucleus_Z(code);
     return fmt::format("Nucleus_A{}_Z{}", A, Z);

corsika/detail/framework/geometry/Path.inl

@@ -9,7 +9,9 @@
 #pragma once
 
 #include <deque>
+
 #include <corsika/framework/geometry/Point.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
 
 namespace corsika {
 
@@ -51,15 +53,21 @@ namespace corsika {
 
   inline LengthType Path::getLength() const { return length_; }
 
-  inline Point Path::getStart() const { return points_.front(); }
+  inline Point const& Path::getStart() const { return points_.front(); }
+
+  inline Point const& Path::getEnd() const { return points_.back(); }
+
+  inline Point const& Path::getPoint(std::size_t const index) const {
+    return points_.at(index);
+  }
 
-  inline Point Path::getEnd() const { return points_.back(); }
+  inline Path::const_iterator Path::begin() const { return points_.cbegin(); }
 
-  inline Point Path::getPoint(std::size_t const index) const { return points_.at(index); }
+  inline Path::const_iterator Path::end() const { return points_.cend(); }
 
-  inline auto Path::begin() { return points_.begin(); }
+  inline Path::iterator Path::begin() { return points_.begin(); }
 
-  inline auto Path::end() { return points_.end(); }
+  inline Path::iterator Path::end() { return points_.end(); }
 
   inline int Path::getNSegments() const { return points_.size() - 1; }
 

corsika/detail/framework/process/SwitchProcessSequence.inl

@@ -233,7 +233,7 @@ namespace corsika {
           return A_.getCrossSection(projectile, projectile.getPID(),
                                     {projectile.getEnergy(), projectile.getMomentum()});
         }
-      } else if (process1_type::is_process_sequence) {
+      } else if constexpr (process1_type::is_process_sequence) {
         return A_.getCrossSection(projectile, targetId, targetP4);
       }
 
@@ -252,7 +252,7 @@ namespace corsika {
         } else {
           return B_.getCrossSection(projectile, targetId, targetP4);
         }
-      } else if (process2_type::is_process_sequence) {
+      } else if constexpr (process2_type::is_process_sequence) {
         return B_.getCrossSection(projectile, targetId, targetP4);
       }
     }

corsika/detail/media/FlatExponential.inl

@@ -29,7 +29,7 @@ namespace corsika {
   template <typename T>
   inline MassDensityType FlatExponential<T>::getMassDensity(Point const& point) const {
     return BaseExponential<FlatExponential<T>>::getMassDensity(
-        (point - BaseExponential<FlatExponential<T>>::getAnchorPoint()).getNorm());
+        (point - BaseExponential<FlatExponential<T>>::getAnchorPoint()).dot(axis_));
   }
 
   template <typename T>

corsika/detail/modules/proposal/HadronicPhotonModel.inl

@@ -12,13 +12,16 @@
 #include <corsika/framework/core/EnergyMomentumOperations.hpp>
 
 #include <tuple>
+#include <random>
 
 namespace corsika::proposal {
 
-  template <typename THadronicModel>
-  inline HadronicPhotonModel<THadronicModel>::HadronicPhotonModel(
-      THadronicModel& _hadint, HEPEnergyType const& _heenthresholdNN)
-      : heHadronicInteraction_(_hadint)
+  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
@@ -27,8 +30,8 @@ namespace corsika::proposal {
     if (!heHadronicInteraction_.isValid(Code::Rho0, Code::Proton, sqrtS)) {
       CORSIKA_LOGGER_CRITICAL(
           logger_,
-          "Invalid energy threshold for hadron interaction model. theshold_lab= {} GeV, "
-          "theshold_com={} GeV",
+          "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!");
     }
@@ -37,65 +40,90 @@ namespace corsika::proposal {
         _heenthresholdNN / 1_GeV);
   }
 
-  template <typename THadronicModel>
+  template <typename THadronicLEModel, typename THadronicHEModel>
   template <typename TStackView>
-  inline ProcessReturn HadronicPhotonModel<THadronicModel>::doHadronicPhotonInteraction(
+  inline ProcessReturn
+  HadronicPhotonModel<THadronicLEModel, THadronicHEModel>::doHadronicPhotonInteraction(
       TStackView& view, CoordinateSystemPtr const& labCS, FourMomentum const& photonP4,
       Code const& targetId) {
-    if (photonP4.getTimeLikeComponent() > heHadronicModelThresholdLabNN_) {
-      CORSIKA_LOGGER_TRACE(
-          logger_, "HE photo-hadronic interaction! calling hadronic interaction model..");
 
-      //  copy from sibyll::NuclearInteractionModel
-      //  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;
-      Code const hadPhotonCode = Code::Rho0; // stand in for hadronic-photon
-      // target at rest
-      FourMomentum const targetP4(get_mass(targetId),
-                                  MomentumVector(labCS, {0_GeV, 0_GeV, 0_GeV}));
-      auto hadronicPhoton = photonStack.addParticle(std::make_tuple(
-          hadPhotonCode, 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);
+    //  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",
-                           hadPhotonCode, targetId,
-                           photonP4.getTimeLikeComponent() / 1_GeV);
-      // check if had. model can handle configuration
-      auto const sqrtSNN = (photonP4 + targetP4 / get_nucleus_A(targetId)).getNorm();
+    // 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
+    auto const sqrtSNN = (photonP4 + targetP4 / get_nucleus_A(targetId)).getNorm();
+    CORSIKA_LOGGER_DEBUG(logger_, "sqrtS={} GeV", sqrtSNN / 1_GeV);
+    if (photonP4.getTimeLikeComponent() > heHadronicModelThresholdLabNN_) {
+      CORSIKA_LOGGER_TRACE(logger_, "HE photo-hadronic interaction!");
       // 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(hadPhotonCode, targetId, sqrtSNN)) {
+      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!",
-            hadPhotonCode, targetId, get_nucleus_A(targetId), get_nucleus_Z(targetId),
+            Code::Rho0, targetId, get_nucleus_A(targetId), get_nucleus_Z(targetId),
             sqrtSNN / 1_GeV);
         return ProcessReturn::Ok;
       }
-      heHadronicInteraction_.doInteraction(photon_secondaries, hadPhotonCode, targetId,
+      heHadronicInteraction_.doInteraction(photon_secondaries, Code::Rho0, targetId,
                                            photonP4, targetP4);
-      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()));
-      }
     } else {
-      CORSIKA_LOGGER_TRACE(
-          logger_,
-          "LE photo-hadronic interaction! Production of secondaries not implemented..");
+      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);
+
+      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;
   }

corsika/detail/modules/proposal/InteractionModel.inl

@@ -18,16 +18,17 @@
 
 namespace corsika::proposal {
 
-  template <typename THadronicModel>
+  template <typename THadronicLEModel, typename THadronicHEModel>
   template <typename TEnvironment>
-  inline InteractionModel<THadronicModel>::InteractionModel(
-      TEnvironment const& _env, THadronicModel& _hadint,
+  inline InteractionModel<THadronicLEModel, THadronicHEModel>::InteractionModel(
+      TEnvironment const& _env, THadronicLEModel& _hadintLE, THadronicHEModel& _hadintHE,
       HEPEnergyType const& _enthreshold)
       : ProposalProcessBase(_env)
-      , HadronicPhotonModel<THadronicModel>(_hadint, _enthreshold) {}
+      , HadronicPhotonModel<THadronicLEModel, THadronicHEModel>(_hadintLE, _hadintHE,
+                                                                _enthreshold) {}
 
-  template <typename THadronicModel>
-  inline void InteractionModel<THadronicModel>::buildCalculator(
+  template <typename THadronicLEModel, typename THadronicHEModel>
+  inline void InteractionModel<THadronicLEModel, THadronicHEModel>::buildCalculator(
       Code code, NuclearComposition const& comp) {
     // search crosssection builder for given particle
     auto p_cross = cross.find(code);
@@ -52,9 +53,10 @@ namespace corsika::proposal {
         PROPOSAL::make_interaction(c, true, true));
   }
 
-  template <typename THadronicModel>
+  template <typename THadronicLEModel, typename THadronicHEModel>
   template <typename TStackView>
-  inline ProcessReturn InteractionModel<THadronicModel>::doInteraction(
+  inline ProcessReturn
+  InteractionModel<THadronicLEModel, THadronicHEModel>::doInteraction(
       TStackView& view, Code const projectileId, FourMomentum const& projectileP4) {
 
     auto const projectile = view.getProjectile();
@@ -138,9 +140,10 @@ namespace corsika::proposal {
     return ProcessReturn::Ok;
   }
 
-  template <typename THadronicModel>
+  template <typename THadronicLEModel, typename THadronicHEModel>
   template <typename TParticle>
-  inline CrossSectionType InteractionModel<THadronicModel>::getCrossSection(
+  inline CrossSectionType
+  InteractionModel<THadronicLEModel, THadronicHEModel>::getCrossSection(
       TParticle const& projectile, Code const projectileId,
       FourMomentum const& projectileP4) {
 

corsika/detail/modules/sophia/InteractionModel.inl

+/*
+ * (c) Copyright 2022 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/geometry/Point.hpp>
+
+#include <corsika/modules/sophia/ParticleConversion.hpp>
+#include <corsika/framework/utility/COMBoost.hpp>
+#include <corsika/modules/sophia/SophiaStack.hpp>
+#include <corsika/framework/core/EnergyMomentumOperations.hpp>
+
+#include <sophia.hpp>
+
+namespace corsika::sophia {
+
+  inline void InteractionModel::setVerbose(bool const flag) { sophia_listing_ = flag; }
+
+  inline InteractionModel::InteractionModel()
+      : sophia_listing_(false) {
+    // set all particles stable in SOPHIA
+    for (int i = 0; i < 49; ++i) so_csydec_.idb[i] = -abs(so_csydec_.idb[i]);
+  }
+
+  inline InteractionModel::~InteractionModel() {
+    CORSIKA_LOG_DEBUG("Sophia::Model n={}", count_);
+  }
+
+  inline bool constexpr InteractionModel::isValid(Code const projectileId,
+                                                  Code const targetId,
+                                                  HEPEnergyType const sqrtSnn) const {
+    if ((minEnergyCoM_ > sqrtSnn) || (sqrtSnn > maxEnergyCoM_)) { return false; }
+
+    if (!(targetId == Code::Proton || targetId == Code::Neutron ||
+          targetId == Code::Hydrogen))
+      return false;
+
+    if (projectileId != Code::Photon) return false;
+
+    return true;
+  }
+
+  template <typename TSecondaryView>
+  inline void InteractionModel::doInteraction(TSecondaryView& secondaries,
+                                              Code const projectileId,
+                                              Code const targetId,
+                                              FourMomentum const& projectileP4,
+                                              FourMomentum const& targetP4) {
+
+    CORSIKA_LOGGER_DEBUG(logger_, "projectile: Id={}, E={} GeV, p3={} GeV", projectileId,
+                         projectileP4.getTimeLikeComponent() / 1_GeV,
+                         projectileP4.getSpaceLikeComponents().getComponents() / 1_GeV);
+    CORSIKA_LOGGER_DEBUG(logger_, "target: Id={}, E={} GeV, p3={} GeV", targetId,
+                         targetP4.getTimeLikeComponent() / 1_GeV,
+                         targetP4.getSpaceLikeComponents().getComponents() / 1_GeV);
+
+    // sqrtS per target nucleon
+    HEPEnergyType const sqrtS = (projectileP4 + targetP4).getNorm();
+
+    CORSIKA_LOGGER_DEBUG(logger_, "sqrtS={}GeV", sqrtS / 1_GeV);
+
+    // accepts only photon-nucleon interactions
+    if (!isValid(projectileId, targetId, sqrtS)) {
+      CORSIKA_LOGGER_ERROR(logger_, "Invalid target/projectile/energy combination");
+      throw std::runtime_error("SOPHIA: Invalid target/projectile/energy combination");
+    }
+
+    COMBoost const boost(projectileP4, targetP4);
+
+    int nucleonSophiaCode = convertToSophiaRaw(targetId); // either proton or neutron
+    // initialize resonance spectrum
+    initial_(nucleonSophiaCode);
+    double Enucleon = targetP4.getTimeLikeComponent() / 1_GeV;
+    double Ephoton = projectileP4.getTimeLikeComponent() / 1_GeV;
+    double theta = 0.0; // set nucleon at rest in collision
+    int Imode = -1;     // overwritten inside SOPHIA
+    CORSIKA_LOGGER_DEBUG(logger_,
+                         "calling SOPHIA eventgen with L0={}, E0={}, eps={},theta={}",
+                         nucleonSophiaCode, Enucleon, Ephoton, theta);
+    count_++;
+    // call sophia
+    eventgen_(nucleonSophiaCode, Enucleon, Ephoton, theta, Imode);
+
+    if (sophia_listing_) {
+      int arg = 3;
+      print_event_(arg);
+    }
+
+    auto const& originalCS = boost.getOriginalCS();
+    // SOPHIA has photon along -z  and nucleon along +z (GZK calc..)
+    COMBoost const boostInternal(targetP4, projectileP4);
+    auto const& csPrime = boost.getRotatedCS();
+    CoordinateSystemPtr csPrimePrime =
+        make_rotation(csPrime, QuantityVector<length_d>{1_m, 0_m, 0_m}, M_PI);
+
+    SophiaStack ss;
+
+    MomentumVector P_final(originalCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
+    HEPEnergyType E_final = 0_GeV;
+    for (auto& psop : ss) {
+      // abort on particles that have decayed in SOPHIA. Should not happen!
+      if (psop.hasDecayed()) { // LCOV_EXCL_START
+        throw std::runtime_error("found particle that decayed in SOPHIA!");
+      } // LCOV_EXCL_STOP
+
+      auto momentumSophia = psop.getMomentum(csPrimePrime);
+      momentumSophia.rebase(csPrime);
+      auto const energySophia = psop.getEnergy();
+      auto const P4com = boostInternal.toCoM(FourVector{energySophia, momentumSophia});
+      auto const P4lab = boost.fromCoM(P4com);
+      SophiaCode const pidSophia = psop.getPID();
+      Code const pid = convertFromSophia(pidSophia);
+      auto momentum = P4lab.getSpaceLikeComponents();
+      momentum.rebase(originalCS);
+      HEPEnergyType const Ekin =
+          calculate_kinetic_energy(momentum.getNorm(), get_mass(pid));
+
+      CORSIKA_LOGGER_TRACE(logger_, "SOPHIA: pid={}, p={} GeV", pidSophia,
+                           momentumSophia.getComponents() / 1_GeV);
+
+      CORSIKA_LOGGER_TRACE(logger_, "CORSIKA: pid={}, p={} GeV", pid,
+                           momentum.getComponents() / 1_GeV);
+
+      auto pnew =
+          secondaries.addSecondary(std::make_tuple(pid, Ekin, momentum.normalized()));
+
+      P_final += pnew.getMomentum();
+      E_final += pnew.getEnergy();
+    }
+    CORSIKA_LOGGER_TRACE(logger_, "Efinal={} GeV,Pfinal={} GeV", E_final / 1_GeV,
+                         P_final.getComponents() / 1_GeV);
+  }
+
+} // namespace corsika::sophia
\ No newline at end of file

corsika/detail/modules/sophia/ParticleConversion.inl

+/*
+ * (c) Copyright 2022 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 <sophia.hpp>
+
+namespace corsika::sophia {
+
+  inline HEPMassType getSophiaMass(Code const pCode) {
+    if (is_nucleus(pCode)) throw std::runtime_error("Not defined for nuclei.");
+    auto sCode = convertToSophiaRaw(pCode);
+    if (sCode == 0)
+      throw std::runtime_error("getSophiaMass: unknown particle!");
+    else
+      return sqrt(get_sophia_mass2(sCode)) * 1_GeV;
+  }
+} // namespace corsika::sophia
\ No newline at end of file

corsika/framework/core/ParticleProperties.hpp

@@ -207,7 +207,7 @@ namespace corsika {
    * @param code
    * @return std::string_view
    */
-  inline std::string_view get_nucleus_name(Code const code);
+  inline std::string get_nucleus_name(Code const code);
 
   /**
    * @brief convert PDG code to CORSIKA 8 internal code.

corsika/framework/geometry/Path.hpp

@@ -9,6 +9,8 @@
 #pragma once
 
 #include <deque>
+
+#include <corsika/framework/core/PhysicalUnits.hpp>
 #include <corsika/framework/geometry/Point.hpp>
 
 namespace corsika {
@@ -18,8 +20,14 @@ namespace corsika {
    * points using N >= 1 straight-line segments.
    */
   class Path {
+
+  protected:
     std::deque<Point> points_;               ///< The points that make up this path.
     LengthType length_ = LengthType::zero(); ///< The length of the path.
+
+    using iterator = std::deque<Point>::iterator;
+    using const_iterator = std::deque<Point>::const_iterator;
+
   public:
     /**
      * Create a Path with a given starting Point.
@@ -49,27 +57,37 @@ namespace corsika {
     /**
      * Get the starting point of the path.
      */
-    inline Point getStart() const;
+    inline Point const& getStart() const;
 
     /**
      * Get the end point of the path.
      */
-    inline Point getEnd() const;
+    inline Point const& getEnd() const;
 
     /**
      * Get a specific point of the path.
      */
-    inline Point getPoint(std::size_t const index) const;
+    inline Point const& getPoint(std::size_t const index) const;
+
+    /**
+     * Return an iterator to the start of the Path.
+     */
+    inline const_iterator begin() const;
+
+    /**
+     * Return an iterator to the end of the Path.
+     */
+    inline const_iterator end() const;
 
     /**
      * Return an iterator to the start of the Path.
      */
-    inline auto begin();
+    inline iterator begin();
 
     /**
      * Return an iterator to the end of the Path.
      */
-    inline auto end();
+    inline iterator end();
 
     /**
      * Get the number of steps in the path.

corsika/framework/process/SwitchProcessSequence.hpp

@@ -17,12 +17,15 @@
 #include <corsika/framework/process/BoundaryCrossingProcess.hpp>
 #include <corsika/framework/process/ContinuousProcess.hpp>
 #include <corsika/framework/process/ContinuousProcessIndex.hpp>
+#include <corsika/framework/process/ContinuousProcessStepLength.hpp>
 #include <corsika/framework/process/DecayProcess.hpp>
 #include <corsika/framework/process/InteractionProcess.hpp>
 #include <corsika/framework/process/ProcessReturn.hpp>
 #include <corsika/framework/process/SecondariesProcess.hpp>
 #include <corsika/framework/process/StackProcess.hpp>
 #include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/core/Step.hpp>
+#include <corsika/framework/geometry/FourVector.hpp>
 
 #include <cmath>
 #include <limits>

corsika/framework/units/quantity.hpp

@@ -763,6 +763,8 @@ namespace phys {
 
     // The standard SI prefixes.
 
+    constexpr long double quetta = 1e+30L;
+    constexpr long double ronna = 1e+27L;
     constexpr long double yotta = 1e+24L;
     constexpr long double zetta = 1e+21L;
     constexpr long double exa = 1e+18L;
@@ -783,6 +785,8 @@ namespace phys {
     constexpr long double atto = 1e-18L;
     constexpr long double zepto = 1e-21L;
     constexpr long double yocto = 1e-24L;
+    constexpr long double ronto = 1e-27L;
+    constexpr long double quecto = 1e-30L;
 
     // Binary prefixes, pending adoption.
 
@@ -938,6 +942,8 @@ namespace phys {
   }
 
 #define QUANTITY_DEFINE_SCALING_LITERALS(pfx, dim, fact)    \
+  QUANTITY_DEFINE_SCALING_LITERAL(Q##pfx, dim, fact* quetta) \
+  QUANTITY_DEFINE_SCALING_LITERAL(R##pfx, dim, fact* ronna) \
   QUANTITY_DEFINE_SCALING_LITERAL(Y##pfx, dim, fact* yotta) \
   QUANTITY_DEFINE_SCALING_LITERAL(Z##pfx, dim, fact* zetta) \
   QUANTITY_DEFINE_SCALING_LITERAL(E##pfx, dim, fact* exa)   \
@@ -958,7 +964,9 @@ namespace phys {
   QUANTITY_DEFINE_SCALING_LITERAL(f##pfx, dim, fact* femto) \
   QUANTITY_DEFINE_SCALING_LITERAL(a##pfx, dim, fact* atto)  \
   QUANTITY_DEFINE_SCALING_LITERAL(z##pfx, dim, fact* zepto) \
-  QUANTITY_DEFINE_SCALING_LITERAL(y##pfx, dim, fact* yocto)
+  QUANTITY_DEFINE_SCALING_LITERAL(y##pfx, dim, fact* yocto) \
+  QUANTITY_DEFINE_SCALING_LITERAL(r##pfx, dim, fact* ronto) \
+  QUANTITY_DEFINE_SCALING_LITERAL(q##pfx, dim, fact* quecto)
 
 #define QUANTITY_DEFINE_LITERALS(pfx, dim) QUANTITY_DEFINE_SCALING_LITERALS(pfx, dim, 1)
 

corsika/framework/units/quantity_io.hpp

@@ -81,6 +81,10 @@ inline Rep prefix( std::string const prefix_ )
         { "da", deka },
         { "d", deci  },
         { "c", centi },
+        { "r", ronto },
+        { "q", quecto },
+        { "R", ronna },
+        { "Q", quetta },
     };
 
     auto pos = table.find( prefix_ );

corsika/modules/PROPOSAL.hpp

@@ -13,12 +13,15 @@
 
 namespace corsika::proposal {
 
-  template <typename THadronicModel>
-  class Interaction : public InteractionModel<THadronicModel>,
-                      public InteractionProcess<Interaction<THadronicModel>> {
+  template <typename THadronicLEModel, typename THadronicHEModel>
+  class Interaction
+      : public InteractionModel<THadronicLEModel, THadronicHEModel>,
+        public InteractionProcess<Interaction<THadronicLEModel, THadronicHEModel>> {
   public:
     template <typename TEnvironment>
-    Interaction(TEnvironment const& env, THadronicModel& model, HEPEnergyType const& thr)
-        : InteractionModel<THadronicModel>(env, model, thr) {}
+    Interaction(TEnvironment const& env, THadronicLEModel& modelLE,
+                THadronicHEModel& modelHE, HEPEnergyType const& thr)
+        : InteractionModel<THadronicLEModel, THadronicHEModel>(env, modelLE, modelHE,
+                                                               thr) {}
   };
 } // namespace corsika::proposal

corsika/modules/Random.hpp

@@ -17,6 +17,7 @@
   link this togehter, it will fail.
  */
 #include <corsika/modules/sibyll/Random.hpp>
+#include <corsika/modules/sophia/Random.hpp>
 #include <corsika/modules/epos/Random.hpp>
 #include <corsika/modules/urqmd/Random.hpp>
 #include <corsika/modules/qgsjetII/Random.hpp>

corsika/modules/Sophia.hpp

+/*
+ * (c) Copyright 2022 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/sophia/ParticleConversion.hpp>
+
+#include <corsika/modules/sophia/InteractionModel.hpp>
+
+#include <corsika/framework/process/InteractionProcess.hpp>

corsika/modules/proposal/HadronicPhotonModel.hpp

@@ -15,17 +15,18 @@
 namespace corsika::proposal {
 
   //! Implements the production of secondary hadrons for the hadronic interaction of real
-  //! and virtual photons. At high energies an external model
-  //! is needed that implements the doInteraction(TSecondaries& view, Code const
-  //! projectile, Code const target,FourMomentum const& projectileP4, FourMomentum const&
-  //! targetP4) routine. Low energy interactions are currently not implemented. The
+  //! 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 THadronicModel>
+  template <class THadronicLEModel, class THadronicHEModel>
   class HadronicPhotonModel {
   public:
-    HadronicPhotonModel(THadronicModel&, HEPEnergyType const&);
+    HadronicPhotonModel(THadronicLEModel&, THadronicHEModel&, HEPEnergyType const&);
     //!
     //! Calculate produce the hadronic secondaries in a hadronic photon interaction and
     //! store them on the particle stack.
@@ -36,7 +37,8 @@ namespace corsika::proposal {
 
   private:
     inline static auto logger_{get_logger("corsika_proposal_HadronicPhotonModel")};
-    THadronicModel& heHadronicInteraction_;
+    THadronicLEModel& leHadronicInteraction_;
+    THadronicHEModel& heHadronicInteraction_;
     //! threshold for high energy hadronic interaction model. Lab. energy per nucleon
     HEPEnergyType const heHadronicModelThresholdLabNN_;
   };

corsika/modules/proposal/InteractionModel.hpp

@@ -33,9 +33,10 @@ namespace corsika::proposal {
   //! @tparam THadronicModel
   //!
 
-  template <class THadronicModel>
-  class InteractionModel : public ProposalProcessBase,
-                           public HadronicPhotonModel<THadronicModel> {
+  template <class THadronicLEModel, class THadronicHEModel>
+  class InteractionModel
+      : public ProposalProcessBase,
+        public HadronicPhotonModel<THadronicLEModel, THadronicHEModel> {
 
     enum { eSECONDARIES, eINTERACTION };
     using calculator_t = std::tuple<std::unique_ptr<PROPOSAL::SecondariesCalculator>,
@@ -57,7 +58,8 @@ namespace corsika::proposal {
     //! compositions and stochastic description limited by the particle cut.
     //!
     template <typename TEnvironment>
-    InteractionModel(TEnvironment const& env, THadronicModel&, HEPEnergyType const&);
+    InteractionModel(TEnvironment const& env, THadronicLEModel&, THadronicHEModel&,
+                     HEPEnergyType const&);
 
     //!
     //! Calculate the rates for the different targets and interactions. Sample a

corsika/modules/radio/propagators/SignalPath.hpp

@@ -18,7 +18,7 @@ namespace corsika {
    *
    * This is basically a container class
    */
-  struct SignalPath final : private Path {
+  struct SignalPath final : public Path {
 
     // TODO: discuss if we need average refractivity or average refractive index
     TimeType const propagation_time_;       ///< The total propagation time.
@@ -28,8 +28,6 @@ namespace corsika {
         refractive_index_destination_; ///< The refractive index at the destination point.
     Vector<dimensionless_d> const emit_;    ///< The (unit-length) emission vector.
     Vector<dimensionless_d> const receive_; ///< The (unit-length) receive vector.
-    std::deque<Point> const
-        points_; ///< A collection of points that make up the geometrical path.
     LengthType const
         R_distance_; ///< The distance from the point of emission to an observer. TODO:
                      ///< optical path, not geometrical! (probably)