From 6405950506d99c02b75cf08d2b3bbfa05398eeba Mon Sep 17 00:00:00 2001
From: ralfulrich <ralf.ulrich@kit.edu>
Date: Wed, 13 Oct 2021 13:23:29 +0200
Subject: [PATCH] added CascadeEquationsProcess

---
 corsika/detail/framework/core/Cascade.inl     |  2 +-
 .../framework/process/ProcessSequence.inl     | 54 +++++++++++++++-
 .../framework/process/SecondariesProcess.hpp  | 17 +++--
 corsika/detail/modules/conex/CONEXhybrid.inl  |  5 +-
 corsika/framework/process/ProcessSequence.hpp |  7 ++
 corsika/framework/process/ProcessTraits.hpp   | 32 ++++++----
 .../framework/process/SecondariesProcess.hpp  | 40 ++++++------
 corsika/framework/process/StackProcess.hpp    | 64 +++++++++----------
 corsika/modules/conex/CONEXhybrid.hpp         |  3 +-
 tests/modules/testCONEX.cpp                   |  9 ++-
 10 files changed, 154 insertions(+), 79 deletions(-)

diff --git a/corsika/detail/framework/core/Cascade.inl b/corsika/detail/framework/core/Cascade.inl
index 07ec49546..2bb692d06 100644
--- a/corsika/detail/framework/core/Cascade.inl
+++ b/corsika/detail/framework/core/Cascade.inl
@@ -48,7 +48,7 @@ namespace corsika {
       }
       // do cascade equations, which can put new particles on Stack,
       // thus, the double loop
-      // doCascadeEquations();
+      sequence_.doCascadeEquations(stack_);
     }
   }
 
diff --git a/corsika/detail/framework/process/ProcessSequence.inl b/corsika/detail/framework/process/ProcessSequence.inl
index 5eb3d20ae..78a28da60 100644
--- a/corsika/detail/framework/process/ProcessSequence.inl
+++ b/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,55 @@ 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>
   template <typename TSecondaryView>
@@ -485,8 +535,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 <
diff --git a/corsika/detail/framework/process/SecondariesProcess.hpp b/corsika/detail/framework/process/SecondariesProcess.hpp
index ee76cce85..f00555814 100644
--- a/corsika/detail/framework/process/SecondariesProcess.hpp
+++ b/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;
diff --git a/corsika/detail/modules/conex/CONEXhybrid.inl b/corsika/detail/modules/conex/CONEXhybrid.inl
index 95b47ace8..47a072f8a 100644
--- a/corsika/detail/modules/conex/CONEXhybrid.inl
+++ b/corsika/detail/modules/conex/CONEXhybrid.inl
@@ -126,7 +126,7 @@ namespace corsika {
 
     // 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};
+    std::array<int, 3> ioseed{1, 1, 1};
 
     double xminp = injectionHeight_ / 1_m;
 
@@ -233,7 +233,8 @@ namespace corsika {
     return true;
   }
 
-  inline void CONEXhybrid::solveCE() {
+  template <typename TStack>
+  inline void CONEXhybrid::doCascadeEquations(TStack& stack) {
 
     ::conex::conexcascade_();
 
diff --git a/corsika/framework/process/ProcessSequence.hpp b/corsika/framework/process/ProcessSequence.hpp
index 05191c248..432b78a28 100644
--- a/corsika/framework/process/ProcessSequence.hpp
+++ b/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,12 @@ namespace corsika {
     template <typename TStack>
     void doStack(TStack& stack);
 
+    /**
+       Execute the CascadeEquationsProcess-es in the ProcessSequence
+     */
+    template <typename TStack>
+    void doCascadeEquations(TStack& stack);
+
     /**
      * Calculate the maximum allowed length of the next tracking step, based on all
      * ContinuousProcess-es
diff --git a/corsika/framework/process/ProcessTraits.hpp b/corsika/framework/process/ProcessTraits.hpp
index cccd80aba..24910fece 100644
--- a/corsika/framework/process/ProcessTraits.hpp
+++ b/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;
diff --git a/corsika/framework/process/SecondariesProcess.hpp b/corsika/framework/process/SecondariesProcess.hpp
index 975b1e785..8fdc6b67d 100644
--- a/corsika/framework/process/SecondariesProcess.hpp
+++ b/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<
diff --git a/corsika/framework/process/StackProcess.hpp b/corsika/framework/process/StackProcess.hpp
index 773265f49..34a9c50bc 100644
--- a/corsika/framework/process/StackProcess.hpp
+++ b/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,
diff --git a/corsika/modules/conex/CONEXhybrid.hpp b/corsika/modules/conex/CONEXhybrid.hpp
index aae3e0f99..cbc73b7af 100644
--- a/corsika/modules/conex/CONEXhybrid.hpp
+++ b/corsika/modules/conex/CONEXhybrid.hpp
@@ -46,7 +46,8 @@ namespace corsika {
     /**
      * Cascade equations are solved basoned on the data in the tables
      */
-    void solveCE();
+    template <typename TStack>
+    void doCascadeEquations(TStack& stack);
 
     /**
      * Internal function to fill particle data inside CONEX
diff --git a/tests/modules/testCONEX.cpp b/tests/modules/testCONEX.cpp
index 5a943cd01..baf808aed 100644
--- a/tests/modules/testCONEX.cpp
+++ b/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);
@@ -102,7 +104,7 @@ TEST_CASE("CONEXSourceCut") {
 
   CONEXhybrid conex(center, showerAxis, t, injectionHeight, E0, get_PDG(Code::Proton));
   conex.init();
-  
+
   HEPEnergyType const Eem{1_PeV};
   auto const momentum = showerAxis.getDirection() * Eem;
 
@@ -123,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>
-- 
GitLab