Merge branch 'conex_multi_event' into 'master'

make conex module fit for multiple events. Added init function

See merge request AirShowerPhysics/corsika!385

corsika/detail/framework/core/Cascade.inl

@@ -31,10 +31,12 @@ namespace corsika {
     setNodes(); // put each particle on stack in correct environment volume
 
     while (!stack_.isEmpty()) {
+
+      sequence_.initCascadeEquations();
+
       while (!stack_.isEmpty()) {
         CORSIKA_LOG_TRACE("Stack: {}", stack_.asString());
         count_++;
-
         auto pNext = stack_.getNextParticle();
 
         CORSIKA_LOG_TRACE(
@@ -46,9 +48,10 @@ namespace corsika {
         step(pNext);
         sequence_.doStack(stack_);
       }
+
       // do cascade equations, which can put new particles on Stack,
       // thus, the double loop
-      // doCascadeEquations();
+      sequence_.doCascadeEquations(stack_);
     }
   }
 

corsika/detail/framework/process/CascadeEquationsProcess.hpp

+/*
+ * (c) Copyright 2021 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/process/ProcessTraits.hpp>
+#include <corsika/framework/utility/HasMethodSignature.hpp>
+
+/**
+ * @file CascadeEquationsProcess.hpp
+ */
+
+namespace corsika {
+
+  /**
+   * traits test for CascadeEquationsProcess::doCascadeEquations method.
+   */
+  template <class TProcess, typename TReturn, typename... TArg>
+  struct has_method_doCascadeEquations
+      : public detail::has_method_signature<TReturn, TArg...> {
+
+    //! method signature
+    using detail::has_method_signature<TReturn, TArg...>::testSignature;
+
+    //! the default value
+    template <class T>
+    static std::false_type test(...);
+
+    //! templated parameter option
+    template <class T>
+    static decltype(testSignature(&T::template doCascadeEquations<TArg...>)) test(
+        std::nullptr_t);
+
+    //! non templated parameter option
+    template <class T>
+    static decltype(testSignature(&T::doCascadeEquations)) test(std::nullptr_t);
+
+  public:
+    /**
+     *  @name traits results
+     *  @{
+     */
+    using type = decltype(test<std::decay_t<TProcess>>(nullptr));
+    static const bool value = type::value;
+    //! @}
+  };
+
+  /**
+   * value traits type.
+   */
+  template <class TProcess, typename TReturn, typename... TArg>
+  bool constexpr has_method_doCascadeEquations_v =
+      has_method_doCascadeEquations<TProcess, TReturn, TArg...>::value;
+
+} // namespace corsika

corsika/detail/framework/process/ProcessSequence.inl

@@ -19,6 +19,7 @@
 #include <corsika/framework/process/ProcessReturn.hpp>
 #include <corsika/framework/process/SecondariesProcess.hpp>
 #include <corsika/framework/process/StackProcess.hpp>
+#include <corsika/framework/process/CascadeEquationsProcess.hpp>
 
 #include <cmath>
 #include <limits>
@@ -331,6 +332,79 @@ namespace corsika {
     return tot;
   }
 
+  template <typename TProcess1, typename TProcess2, int IndexStart, int IndexProcess1,
+            int IndexProcess2>
+  template <typename TStack>
+  void ProcessSequence<TProcess1, TProcess2, IndexStart, IndexProcess1,
+                       IndexProcess2>::doCascadeEquations(TStack& stack) {
+
+    if constexpr (is_process_v<process1_type>) { // to protect from further compiler
+                                                 // errors if process1_type is invalid
+      if constexpr (process1_type::is_process_sequence &&
+                    !process1_type::is_switch_process_sequence) {
+
+        A_.doCascadeEquations(stack);
+
+      } else if constexpr (is_cascade_equations_process_v<process1_type>) {
+
+        // interface checking on TProcess1
+        static_assert(has_method_doCascadeEquations_v<TProcess1,
+                                                      void,     // return type
+                                                      TStack&>, // parameter
+                      "TDerived has no method with correct signature \"void "
+                      "doCascadeEquations(TStack&)\" required for "
+                      "CascadeEquationsProcess<TDerived>. ");
+
+        A_.doCascadeEquations(stack);
+      }
+    }
+
+    if constexpr (is_process_v<process2_type>) { // to protect from further compiler
+                                                 // errors if process2_type is invalid
+      if constexpr (process2_type::is_process_sequence &&
+                    !process2_type::is_switch_process_sequence) {
+
+        B_.doCascadeEquations(stack);
+
+      } else if constexpr (is_cascade_equations_process_v<process2_type>) {
+
+        // interface checking on TProcess2
+        static_assert(has_method_doCascadeEquations_v<TProcess2,
+                                                      void,     // return type
+                                                      TStack&>, // parameter
+                      "TDerived has no method with correct signature \"void "
+                      "doCascadeEquations(TStack&)\" required for "
+                      "CascadeEquationsProcess<TDerived>. ");
+
+        B_.doCascadeEquations(stack);
+      }
+    }
+  }
+
+  template <typename TProcess1, typename TProcess2, int IndexStart, int IndexProcess1,
+            int IndexProcess2>
+  void ProcessSequence<TProcess1, TProcess2, IndexStart, IndexProcess1,
+                       IndexProcess2>::initCascadeEquations() {
+
+    if constexpr (is_process_v<process1_type>) { // to protect from further compiler
+                                                 // errors if process1_type is invalid
+      if constexpr ((process1_type::is_process_sequence &&
+                     !process1_type::is_switch_process_sequence) ||
+                    is_cascade_equations_process_v<process1_type>) {
+        A_.initCascadeEquations();
+      }
+    }
+
+    if constexpr (is_process_v<process2_type>) { // to protect from further compiler
+                                                 // errors if process2_type is invalid
+      if constexpr ((process2_type::is_process_sequence &&
+                     !process2_type::is_switch_process_sequence) ||
+                    is_cascade_equations_process_v<process2_type>) {
+        B_.initCascadeEquations();
+      }
+    }
+  } // namespace corsika
+
   template <typename TProcess1, typename TProcess2, int IndexStart, int IndexProcess1,
             int IndexProcess2>
   template <typename TSecondaryView>
@@ -485,8 +559,8 @@ namespace corsika {
   }
 
   /**
-   * traits marker to identify objects containing any StackProcesses
-   **/
+   * traits marker to identify objects containing any StackProcesses.
+   */
   namespace detail {
     // need helper alias to achieve this:
     template <

corsika/detail/framework/process/SecondariesProcess.hpp

@@ -14,8 +14,8 @@
 namespace corsika {
 
   /**
-     traits test for SecondariesProcess::doSecondaries method
-  */
+   * traits test for SecondariesProcess::doSecondaries method.
+   */
   template <class TProcess, typename TReturn, typename... TArg>
   struct has_method_doSecondaries
       : public detail::has_method_signature<TReturn, TArg...> {
@@ -38,16 +38,19 @@ namespace corsika {
 
   public:
     /**
-        @name traits results
-        @{
-    */
+     *  @name traits results
+     *  @{
+     */
     using type = decltype(test<std::decay_t<TProcess>>(nullptr));
     static const bool value = type::value;
     //! @}
   };
 
-  //! @file DecayProcess.hpp
-  //! value traits type
+  /**
+   * @file SecondariesProcess.hpp
+   *
+   * @brief value traits type.
+   */
   template <class TProcess, typename TReturn, typename... TArg>
   bool constexpr has_method_doSecondaries_v =
       has_method_doSecondaries<TProcess, TReturn, TArg...>::value;

corsika/detail/modules/StackInspector.inl

@@ -61,28 +61,44 @@ namespace corsika {
     if (dE < dE_threshold_) return;
     double const progress = dE / E0_;
 
-    double const eta_seconds = elapsed_seconds.count() / progress;
     std::time_t const start_time = std::chrono::system_clock::to_time_t(StartTime_);
-    std::time_t const eta_time = std::chrono::system_clock::to_time_t(
-        StartTime_ + std::chrono::seconds((int)eta_seconds));
-
-    int const yday0 = std::localtime(&start_time)->tm_yday;
-    int const yday1 = std::localtime(&eta_time)->tm_yday;
-    int const dyday = yday1 - yday0;
-
-    CORSIKA_LOG_INFO(
-        "StackInspector: "
-        " time={}"
-        ", running={} seconds"
-        " ( {:.1f}%)"
-        ", nStep={}"
-        ", stackSize={}"
-        ", Estack={} GeV"
-        ", ETA={}{}",
-        std::put_time(std::localtime(&now_time), "%T"), elapsed_seconds.count(),
-        (progress * 100), getStep(), vS.getSize(), Etot / 1_GeV,
-        (dyday == 0 ? "" : fmt::format("+{}d ", dyday)),
-        std::put_time(std::localtime(&eta_time), "%T"));
+
+    if (progress > 0) {
+
+      double const eta_seconds = elapsed_seconds.count() / progress;
+      std::time_t const eta_time = std::chrono::system_clock::to_time_t(
+          StartTime_ + std::chrono::seconds((int)eta_seconds));
+
+      int const yday0 = std::localtime(&start_time)->tm_yday;
+      int const yday1 = std::localtime(&eta_time)->tm_yday;
+      int const dyday = yday1 - yday0;
+
+      CORSIKA_LOG_INFO(
+          "StackInspector: "
+          " time={}"
+          ", running={} seconds"
+          " ( {:.1f}%)"
+          ", nStep={}"
+          ", stackSize={}"
+          ", Estack={} GeV"
+          ", ETA={}{}",
+          std::put_time(std::localtime(&now_time), "%T"), elapsed_seconds.count(),
+          (progress * 100), getStep(), vS.getSize(), Etot / 1_GeV,
+          (dyday == 0 ? "" : fmt::format("+{}d ", dyday)),
+          std::put_time(std::localtime(&eta_time), "%T"));
+    } else {
+      CORSIKA_LOG_INFO(
+          "StackInspector: "
+          " time={}"
+          ", running={} seconds"
+          " ( {:.1f}%)"
+          ", nStep={}"
+          ", stackSize={}"
+          ", Estack={} GeV"
+          ", ETA={}{}",
+          std::put_time(std::localtime(&now_time), "%T"), elapsed_seconds.count(),
+          (progress * 100), getStep(), vS.getSize(), Etot / 1_GeV, "---", "---");
+    }
   }
 
 } // namespace corsika

corsika/detail/modules/conex/CONEXhybrid.inl

@@ -28,6 +28,9 @@ namespace corsika {
       : center_{center}
       , showerAxis_{showerAxis}
       , groundDist_{groundDist}
+      , injectionHeight_{injectionHeight}
+      , primaryEnergy_{primaryEnergy}
+      , primaryPDG_{primaryPDG}
       , showerCore_{showerAxis_.getStart() + showerAxis_.getDirection() * groundDist_}
       , conexObservationCS_{std::invoke([&]() {
         auto const& c8cs = center.getCoordinateSystem();
@@ -77,7 +80,9 @@ namespace corsika {
     CORSIKA_LOG_DEBUG("showerCore (C8): {}",
                       showerCore_.getCoordinates(center.getCoordinateSystem()));
 
-    int randomSeeds[3] = {1234, 0, 0}; // will be overwritten later??
+    int randomSeeds[3] = {1234, 0,
+                          0}; // SEEDS ARE NOT USED. All random numbers are obtained from
+                              // the CORSIKA 8 stream "conex" and "epos"!
     int heModel = eSibyll23;
 
     int nShower = 1; // large to avoid final stats.
@@ -92,8 +97,9 @@ namespace corsika {
                         particleListMode,
 #endif
                         configPath.c_str(), configPath.size());
+  }
 
-    double eprima = primaryEnergy / 1_GeV;
+  inline void CONEXhybrid::initCascadeEquations() {
 
     // set phi, theta
     Vector<length_d> ez{conexObservationCS_, {0._m, 0._m, -1_m}};
@@ -111,16 +117,17 @@ namespace corsika {
         "; phi (deg) = {}",
         theta, phi);
 
-    int ipart = static_cast<int>(primaryPDG);
-    auto rng = RNGManager<>::getInstance().getRandomStream("conex");
+    int ipart = static_cast<int>(primaryPDG_);
 
     double dimpact = 0.; // valid only if shower core is fixed on the observation plane;
                          // for skimming showers an offset is needed like in CONEX
 
-    std::array<int, 3> ioseed{static_cast<int>(rng()), static_cast<int>(rng()),
-                              static_cast<int>(rng())};
+    // SEEDS ARE NOT USED. All random numbers are obtained from
+    // the CORSIKA 8 stream "conex" and "epos"!
+    std::array<int, 3> ioseed{1, 1, 1};
 
-    double xminp = injectionHeight / 1_m;
+    double eprima = primaryEnergy_ / 1_GeV;
+    double xminp = injectionHeight_ / 1_m;
 
     ::conex::conexrun_(ipart, eprima, theta, phi, xminp, dimpact, ioseed.data());
   }
@@ -225,7 +232,8 @@ namespace corsika {
     return true;
   }
 
-  inline void CONEXhybrid::solveCE() {
+  template <typename TStack>
+  inline void CONEXhybrid::doCascadeEquations(TStack&) {
 
     ::conex::conexcascade_();
 

corsika/detail/modules/sibyll/Interaction.inl

@@ -206,10 +206,8 @@ namespace corsika::sibyll {
     // just for show:
     // boost projecticle
     [[maybe_unused]] auto const PprojCoM = boost.toCoM(PprojLab);
-
     // boost target
     [[maybe_unused]] auto const PtargCoM = boost.toCoM(PtargLab);
-
     CORSIKA_LOG_DEBUG(
         "Interaction: ebeam CoM: {} GeV "
         "Interaction: pbeam CoM: {} GeV ",
@@ -328,15 +326,16 @@ namespace corsika::sibyll {
         Ecm_final += psib.getEnergy();
       }
       CORSIKA_LOG_DEBUG(
-          "conservation (all GeV):"
-          "Ecm_initial(per nucleon)={}, Ecm_final(per nucleon)={}, "
-          "Elab_initial={}, Elab_final={}, "
-          "diff (%)={}, "
-          "E in nucleons={}, "
-          "Plab_initial={}, "
-          "Plab_final={} ",
+          "conservation (all GeV): "
+          "Ecm_initial(per nucleon)={:.2f}, Ecm_final(per nucleon)={:.2f}, "
+          "Elab_initial={:.2f}, Elab_final={:.2f}, "
+          "Elab-diff (%)={:.2f}, "
+          "m in target nucleons={:.2f}, "
+          "Plab_initial={:.2f}, "
+          "Plab_final={:.2f} ",
           Ecm / 1_GeV, Ecm_final * 2. / (get_nwounded() + 1) / 1_GeV, Etot / 1_GeV,
-          Elab_final / 1_GeV, (Elab_final / Etot / get_nwounded() - 1) * 100,
+          Elab_final / 1_GeV,
+          (Elab_final / (Etot + get_nwounded() * constants::nucleonMass) - 1) * 100,
           constants::nucleonMass * get_nwounded() / 1_GeV,
           (pProjectileLab / 1_GeV).getComponents(), (Plab_final / 1_GeV).getComponents());
     }

corsika/detail/modules/tracking/Intersect.inl

@@ -60,7 +60,8 @@ namespace corsika {
     for (auto const& node : volumeNode.getChildNodes()) {
 
       Intersections const time_intersections = TDerived::intersect(particle, *node);
-      CORSIKA_LOG_DEBUG("intersection times with child volume {}", fmt::ptr(node));
+      CORSIKA_LOG_DEBUG("search intersection times with child volume {}:",
+                        fmt::ptr(node));
       if (!time_intersections.hasIntersections()) { continue; }
       auto const t_entry = time_intersections.getEntry();
       auto const t_exit = time_intersections.getExit();

corsika/framework/process/CascadeEquationsProcess.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/process/BaseProcess.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+
+#include <corsika/detail/framework/process/CascadeEquationsProcess.hpp> // for extra traits, method/interface checking
+
+namespace corsika {
+
+  /**
+   * @ingroup Processes
+   * @{
+   *
+   * Processes executing cascade-equations calculations.
+   *
+   * Create a new CascadeEquationsProcess, e.g. for XYModel, via:
+   * @code{.cpp}
+   * class XYModel : public CascadeEquationsProcess<XYModel> {};
+   * @endcode
+   *
+   * and provide the necessary interface method:
+   * @code{.cpp}
+   * template <typename TStack>
+   * void doCascadeEquations(TStack& stack);
+   * @endcode
+   *
+   * Cascade equation processes may generate new particles on the stack. They also
+   * typically generate output.
+   */
+
+  template <typename TDerived>
+  class CascadeEquationsProcess : public BaseProcess<TDerived> {
+  public:
+  };
+
+  /**
+   * ProcessTraits specialization to flag CascadeEquationsProcess objects.
+   */
+  template <typename TProcess>
+  struct is_cascade_equations_process<
+      TProcess, std::enable_if_t<std::is_base_of_v<
+                    CascadeEquationsProcess<typename std::decay_t<TProcess>>,
+                    typename std::decay_t<TProcess>>>> : std::true_type {};
+
+  //! @}
+
+} // namespace corsika

corsika/framework/process/ProcessSequence.hpp

@@ -20,6 +20,7 @@
 #include <corsika/framework/process/ContinuousProcessStepLength.hpp>
 #include <corsika/framework/process/DecayProcess.hpp>
 #include <corsika/framework/process/InteractionProcess.hpp>
+#include <corsika/framework/process/CascadeEquationsProcess.hpp>
 #include <corsika/framework/process/ProcessReturn.hpp>
 #include <corsika/framework/process/SecondariesProcess.hpp>
 #include <corsika/framework/process/StackProcess.hpp>
@@ -209,6 +210,17 @@ namespace corsika {
     template <typename TStack>
     void doStack(TStack& stack);
 
+    /**
+     * Execute the CascadeEquationsProcess-es in the ProcessSequence.
+     */
+    template <typename TStack>
+    void doCascadeEquations(TStack& stack);
+
+    /**
+     * Init the CascadeEquationsProcess-es in the ProcessSequence.
+     */
+    void initCascadeEquations();
+
     /**
      * Calculate the maximum allowed length of the next tracking step, based on all
      * ContinuousProcess-es

corsika/framework/process/ProcessTraits.hpp

@@ -17,7 +17,7 @@
 namespace corsika {
 
   /**
-   * A traits marker to identify BaseProcess, thus any type of process
+   * A traits marker to identify BaseProcess, thus any type of process.
    */
   template <typename TProcess, typename TEnable = void>
   struct is_process : std::false_type {};
@@ -26,7 +26,7 @@ namespace corsika {
   bool constexpr is_process_v = is_process<TProcess>::value;
 
   /**
-   * A traits marker to identify ContinuousProcess
+   * A traits marker to identify ContinuousProcess.
    */
   template <typename TProcess, typename Enable = void>
   struct is_continuous_process : std::false_type {};
@@ -35,7 +35,7 @@ namespace corsika {
   bool constexpr is_continuous_process_v = is_continuous_process<TProcess>::value;
 
   /**
-   * A traits marker to identify DecayProcess
+   * A traits marker to identify DecayProcess.
    */
   template <typename TProcess, typename Enable = void>
   struct is_decay_process : std::false_type {};
@@ -44,7 +44,7 @@ namespace corsika {
   bool constexpr is_decay_process_v = is_decay_process<TProcess>::value;
 
   /**
-   * A traits marker to identify StackProcess
+   * A traits marker to identify StackProcess.
    */
   template <typename TProcess, typename Enable = void>
   struct is_stack_process : std::false_type {};
@@ -53,7 +53,17 @@ namespace corsika {
   bool constexpr is_stack_process_v = is_stack_process<TProcess>::value;
 
   /**
-   * A traits marker to identify SecondariesProcess
+   * A traits marker to identify CascadeEquationsProcess.
+   */
+  template <typename TProcess, typename Enable = void>
+  struct is_cascade_equations_process : std::false_type {};
+
+  template <typename TProcess>
+  bool constexpr is_cascade_equations_process_v =
+      is_cascade_equations_process<TProcess>::value;
+
+  /**
+   * A traits marker to identify SecondariesProcess.
    */
   template <typename TProcess, typename Enable = void>
   struct is_secondaries_process : std::false_type {};
@@ -62,7 +72,7 @@ namespace corsika {
   bool constexpr is_secondaries_process_v = is_secondaries_process<TProcess>::value;
 
   /**
-   * A traits marker to identify BoundaryProcess
+   * A traits marker to identify BoundaryProcess.
    */
   template <typename TProcess, typename Enable = void>
   struct is_boundary_process : std::false_type {};
@@ -71,7 +81,7 @@ namespace corsika {
   bool constexpr is_boundary_process_v = is_boundary_process<TProcess>::value;
 
   /**
-   * A traits marker to identify InteractionProcess
+   * A traits marker to identify InteractionProcess.
    */
   template <typename TProcess, typename Enable = void>
   struct is_interaction_process : std::false_type {};
@@ -80,8 +90,8 @@ namespace corsika {
   bool constexpr is_interaction_process_v = is_interaction_process<TProcess>::value;
 
   /**
-   * A traits marker to identify ProcessSequence that contain a StackProcess
-   **/
+   * A traits marker to identify ProcessSequence that contain a StackProcess.
+   */
   template <typename TClass>
   struct contains_stack_process : std::false_type {};
 
@@ -89,8 +99,8 @@ namespace corsika {
   bool constexpr contains_stack_process_v = contains_stack_process<TClass>::value;
 
   /**
-   * traits class to count any type of Process, general version
-   **/
+   * traits class to count any type of Process, general version.
+   */
   template <typename TProcess, int N = 0, typename Enable = void>
   struct count_processes {
     static unsigned int constexpr count = N;

corsika/framework/process/SecondariesProcess.hpp

@@ -16,24 +16,24 @@
 namespace corsika {
 
   /**
-     @ingroup Processes
-     @{
-
-     Processes acting on the secondaries produced by other processes.
-
-     Create a new SecondariesProcess, e.g. for XYModel, via
-     @code{.cpp}
-     class XYModel : public SecondariesProcess<XYModel> {};
-     @endcode
-
-     and provide the necessary interface method:
-     @code{.cpp}
-     template <typename TStackView>
-     void doSecondaries(TStackView& StackView);
-     @endcode
-
-     where StackView is an object that can store secondaries on a
-     stack and also iterate over these secondaries.
+   * @ingroup Processes
+   * @{
+   *
+   * Processes acting on the secondaries produced by other processes.
+   *
+   * Create a new SecondariesProcess, e.g. for XYModel, via
+   * @code{.cpp}
+   * class XYModel : public SecondariesProcess<XYModel> {};
+   * @endcode
+   *
+   * and provide the necessary interface method:
+   * @code{.cpp}
+   * template <typename TStackView>
+   * void doSecondaries(TStackView& StackView);
+   * @endcode
+   *
+   * where StackView is an object that can store secondaries on a
+   * stack and also iterate over these secondaries.
    */
 
   template <typename TDerived>
@@ -42,8 +42,8 @@ namespace corsika {
   };
 
   /**
-   * ProcessTraits specialization to flag SecondariesProcess objects
-   **/
+   * ProcessTraits specialization to flag SecondariesProcess objects.
+   */
   template <typename TProcess>
   struct is_secondaries_process<
       TProcess, std::enable_if_t<

corsika/framework/process/StackProcess.hpp

@@ -16,32 +16,32 @@
 namespace corsika {
 
   /**
-     @ingroup Processes
-     @{
-
-     Process to act on the entire particle stack
-
-     Create a new StackProcess, e.g. for XYModel, via
-     @code{.cpp}
-     class XYModel : public StackProcess<XYModel> {};
-     @endcode
-
-     and provide the necessary interface method
-     @code{.cpp}
-     template <typename TStack>
-     void XYModel::doStack(TStack&);
-     @endcode
-
-     Where, of course, Stack is the valid
-     class to access particles on the Stack. This methods does
-     not need to be templated, they could use the types
-     e.g. corsika::setup::Stack directly -- but by the cost of
-     loosing all flexibility otherwise provided.
-
-     A StackProcess has only one constructor `StackProcess::StackProcess(unsigned int
-     const nStep)` where nStep is the number of steps of the cascade stepping after which
-     the stack process should be run. Good values are on the order of 1000, which will not
-     compromise run time in the end, but provide all the benefits of the StackProcess.
+   * @ingroup Processes
+   * @{
+   *
+   * Process to act on the entire particle stack.
+   *
+   * Create a new StackProcess, e.g. for XYModel, via:
+   * @code{.cpp}
+   * class XYModel : public StackProcess<XYModel> {};
+   * @endcode
+   *
+   * and provide the necessary interface method:
+   * @code{.cpp}
+   * template <typename TStack>
+   * void XYModel::doStack(TStack&);
+   * @endcode
+   *
+   * Where, of course, Stack is the valid
+   * class to access particles on the Stack. This methods does
+   * not need to be templated, they could use the types
+   * e.g. corsika::setup::Stack directly -- but by the cost of
+   * loosing all flexibility otherwise provided.
+   *
+   * A StackProcess has only one constructor `StackProcess::StackProcess(unsigned int
+   * const nStep)` where nStep is the number of steps of the cascade stepping after which
+   * the stack process should be run. Good values are on the order of 1000, which will not
+   * compromise run time in the end, but provide all the benefits of the StackProcess.
    */
 
   template <typename TDerived>
@@ -64,10 +64,10 @@ namespace corsika {
 
   private:
     /**
-       @name The number of "steps" during the cascade processing after
-       which this StackProcess is going to be executed. The logic is
-       "iStep_ modulo nStep_"
-       @{
+     * @name The number of "steps" during the cascade processing after
+     * which this StackProcess is going to be executed. The logic is
+     * "iStep_ modulo nStep_"
+     * @{
      */
     unsigned int nStep_ = 0;
     unsigned long int iStep_ = 0;
@@ -75,8 +75,8 @@ namespace corsika {
   };
 
   /**
-   * ProcessTraits specialization to flag StackProcess objects
-   **/
+   * ProcessTraits specialization to flag StackProcess objects.
+   */
   template <typename TProcess>
   struct is_stack_process<
       TProcess,

corsika/modules/conex/CONEXhybrid.hpp

@@ -8,11 +8,13 @@
 
 #pragma once
 
+#include <corsika/framework/process/SecondariesProcess.hpp>
+#include <corsika/framework/process/CascadeEquationsProcess.hpp>
+
 #include <corsika/framework/core/PhysicalUnits.hpp>
 #include <corsika/framework/core/ParticleProperties.hpp>
 #include <corsika/framework/geometry/Point.hpp>
 #include <corsika/framework/geometry/Vector.hpp>
-#include <corsika/framework/process/SecondariesProcess.hpp>
 #include <corsika/media/ShowerAxis.hpp>
 
 #include <corsika/modules/conex/CONEX_f.hpp>
@@ -23,17 +25,37 @@ namespace corsika {
     LengthType constexpr earthRadius{6371315 * meter};
   } // namespace conex
 
-  class CONEXhybrid : public SecondariesProcess<CONEXhybrid> {
+  class CONEXhybrid : public CascadeEquationsProcess<CONEXhybrid>,
+                      public SecondariesProcess<CONEXhybrid> {
 
   public:
     CONEXhybrid(Point center, ShowerAxis const& showerAxis, LengthType groundDist,
                 LengthType injectionHeight, HEPEnergyType primaryEnergy, PDGCode pdg);
 
+    /**
+     * Main entry point to pass new particle data towards CONEX. If a
+     * particles is selected for CONEX, it is removed from the CORSIKA
+     * 8 stack.
+     */
     template <typename TStackView>
     void doSecondaries(TStackView&);
 
-    void solveCE();
+    /**
+     * init currently needs to be called to initializa a new
+     * event. All tables are cleared, etc.
+     */
+    void initCascadeEquations();
+
+    /**
+     * Cascade equations are solved basoned on the data in the tables
+     */
+    template <typename TStack>
+    void doCascadeEquations(TStack& stack);
 
+    /**
+     * Internal function to fill particle data inside CONEX
+     * tables. Only e.m. particles are selected right now.
+     */
     bool addParticle(Code pid, HEPEnergyType energy, HEPEnergyType mass,
                      Point const& position, Vector<dimensionless_d> const& direction,
                      TimeType t);
@@ -51,8 +73,11 @@ namespace corsika {
 
     Point const center_; //!< center of CONEX Earth
     ShowerAxis const& showerAxis_;
-    LengthType groundDist_;  //!< length from injection point to shower core
-    Point const showerCore_; //!< shower core
+    LengthType groundDist_;       //!< length from injection point to shower core
+    LengthType injectionHeight_;  //!< starting height of primary particle
+    HEPEnergyType primaryEnergy_; //!< primary particle energy
+    PDGCode primaryPDG_;          //!< primary particle PDG
+    Point const showerCore_;      //!< shower core
     CoordinateSystemPtr const conexObservationCS_; //!< CONEX observation frame
     DirectionVector const x_sf_,
         y_sf_; //!< unit vectors of CONEX shower frame, z_sf is shower axis direction

corsika/modules/conex/Random.hpp

@@ -17,6 +17,19 @@
  * This file is an integral part of the epos interface. It must be
  * linked to the executable linked to epos exactly once
  *
+ * Note, that the fortran random numbe interface functions are all
+ * defined in the epos corsika 8 interface:
+ *
+ * ranfst, ranfgt, rmmaqd, ranfini, ranfcv, rmmard, rangen, drangen
+ *
+ * All of them use the epos_random_interface registered as "epos" stream.
+ *
+ * Thus, the fortran part of CONEX will use the "epos" CORSIKA 8 random stream,
+ * only the CONEX c++ part will use the "conex" random stream.
+ *
+ * Since EPOS and CONEX use the same fortran symbols for access to
+ * random numbers this can only be changed by renaming inside the
+ * fortran part.
  */
 
 namespace conex {

examples/hybrid_MC.cpp

@@ -291,8 +291,6 @@ int main(int argc, char** argv) {
   cut.reset();
   eLoss.reset();
 
-  conex_model.solveCE();
-
   auto const hists = sibyllCounted.getHistogram() + sibyllNucCounted.getHistogram() +
                      urqmdCounted.getHistogram();
 

tests/modules/testCONEX.cpp

@@ -46,6 +46,8 @@ const std::string refDataDir = std::string(REFDATADIR); // from cmake
 template <typename T>
 using MExtraEnvirnoment = MediumPropertyModel<UniformMagneticField<T>>;
 
+struct DummyStack {};
+
 TEST_CASE("CONEXSourceCut") {
 
   logging::set_level(logging::level::info);
@@ -70,7 +72,7 @@ TEST_CASE("CONEXSourceCut") {
 
   builder.setNuclearComposition(
       {{Code::Nitrogen, Code::Oxygen},
-       {0.7847f, 1.f - 0.7847f}}); // values taken from AIRES manual, Ar removed for now
+       {0.7847, 1. - 0.7847}}); // values taken from AIRES manual, Ar removed for now
 
   builder.addExponentialLayer(1222.6562_g / (1_cm * 1_cm), 994186.38_cm, 4_km);
   builder.addExponentialLayer(1144.9069_g / (1_cm * 1_cm), 878153.55_cm, 10_km);
@@ -101,6 +103,7 @@ TEST_CASE("CONEXSourceCut") {
   [[maybe_unused]] corsika::sibyll::NuclearInteraction sibyllNuc(sibyll, env);
 
   CONEXhybrid conex(center, showerAxis, t, injectionHeight, E0, get_PDG(Code::Proton));
+  conex.initCascadeEquations();
 
   HEPEnergyType const Eem{1_PeV};
   auto const momentum = showerAxis.getDirection() * Eem;
@@ -122,7 +125,8 @@ TEST_CASE("CONEXSourceCut") {
   auto const momentumPhoton = showerAxis.getDirection() * 1_TeV;
   conex.addParticle(Code::Photon, 1_TeV, 0_eV, emPosition, momentumPhoton.normalized(),
                     0_s);
-  conex.solveCE();
+  DummyStack stack;
+  conex.doCascadeEquations(stack);
 }
 
 #include <algorithm>