corsika/framework/process/ContinuousProcess.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
@@ -18,47 +17,46 @@
 namespace corsika {
 
   /**
-     @ingroup Processes
-     @{
-
-     Processes with continuous effects along a particle Trajectory
-
-     Create a new ContinuousProcess, e.g. for XYModel, via
-     @code{.cpp}
-     class XYModel : public ContinuousProcess<XYModel> {};
-     @endcode
-
-     and provide two necessary interface methods:
-     @code{.cpp}
-     template <typename TParticle, typename TTrack>
-     LengthType getMaxStepLength(TParticle const& p, TTrack const& track) const;
-     @endcode
-
-     which allows any ContinuousProcess to tell to CORSIKA a maximum
-     allowed step length. Such step-length limitation, if it turns out
-     to be smaller/sooner than any other limit (decay length,
-     interaction length, other continuous processes, geometry, etc.)
-     will lead to a limited step length.
-
-     @code{.cpp}
-     template <typename TParticle, typename TTrack>
-     ProcessReturn doContinuous(TParticle& p, TTrack const& t, bool const stepLimit)
-     const;
-     @endcode
-
-     which applied any continuous effects on Particle p along
-     Trajectory t. The particle in all typical scenarios will be
-     altered by a doContinuous. The flag stepLimit will be true if the
-     preious evaluation of getMaxStepLength resulted in this
-     particular ContinuousProcess to be responsible for the step
-     length limit on the current track t. This information can be
-     expoited and avoid e.g. any uncessary calculations.
-
-     Particle and Track are the valid classes to
-     access particles and track (Trajectory) data on the Stack. Those two methods
-     do not need to be templated, they could use the types
-     e.g. corsika::setup::Stack::particle_type -- but by the cost of
-     loosing all flexibility otherwise provided.
+    * @ingroup Processes
+    * @{
+    * Processes with continuous effects along a particle Trajectory.
+    *
+    * Create a new ContinuousProcess, e.g. for XYModel, via:
+    * @code{.cpp}
+    * class XYModel : public ContinuousProcess<XYModel> {};
+    * @endcode
+    *
+    * and provide two necessary interface methods:
+    * @code{.cpp}
+    * template <typename TParticle, typename TTrack>
+    * LengthType getMaxStepLength(TParticle const& p, TTrack const& track) const;
+    * @endcode
+
+    * which allows any ContinuousProcess to tell to CORSIKA a maximum
+    * allowed step length. Such step-length limitation, if it turns out
+    * to be smaller/sooner than any other limit (decay length,
+    * interaction length, other continuous processes, geometry, etc.)
+    * will lead to a limited step length.
+
+    * @code{.cpp}
+    * template <typename TParticle, typename TTrack>
+    * ProcessReturn doContinuous(TParticle& p, TTrack const& t, bool const stepLimit)
+    * const;
+    * @endcode
+
+    * which applied any continuous effects on Particle p along
+    * Trajectory t. The particle in all typical scenarios will be
+    * altered by a doContinuous. The flag stepLimit will be true if the
+    * preious evaluation of getMaxStepLength resulted in this
+    * particular ContinuousProcess to be responsible for the step
+    * length limit on the current track t. This information can be
+    * expoited and avoid e.g. any uncessary calculations.
+
+    * Particle and Track are the valid classes to
+    * access particles and track (Trajectory) data on the Stack. Those two methods
+    * do not need to be templated, they could use the types
+    * e.g. corsika::setup::Stack::particle_type -- but by the cost of
+    * loosing all flexibility otherwise provided.
 
    */
 
@@ -68,8 +66,8 @@ namespace corsika {
   };
 
   /**
-   * ProcessTraits specialization to flag ContinuousProcess objects
-   **/
+   * ProcessTraits specialization to flag ContinuousProcess objects.
+   */
   template <typename TProcess>
   struct is_continuous_process<
       TProcess, std::enable_if_t<

corsika/framework/process/ContinuousProcessIndex.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.
+ * 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
@@ -11,26 +10,24 @@
 namespace corsika {
 
   /**
-     @ingroup Processes
-
-     To index individual processes (continuous processes) inside a
-     ProcessSequence.
-
-   **/
+   * @ingroup Processes
+   * To index individual processes (continuous processes) inside a
+   * ProcessSequence.
+   */
 
   class ContinuousProcessIndex {
   public:
     ContinuousProcessIndex()
-        : id_(-1) {} // default
-    ContinuousProcessIndex(int const id)
+        : id_(nullptr) {} // default
+    ContinuousProcessIndex(void const* id)
         : id_(id) {}
-    void setIndex(int const id) { id_ = id; }
-    int getIndex() const { return id_; }
+    void setIndex(void const* id) { id_ = id; }
+    void const* getIndex() const { return id_; }
     bool operator==(ContinuousProcessIndex const v) const { return id_ == v.id_; }
-    bool operator!=(ContinuousProcessIndex const v) const { return id_ != v.id_; }
+    bool operator!=(ContinuousProcessIndex const v) const { return !(*this == v); }
 
   private:
-    int id_;
+    void const* id_;
   };
 
 } // namespace corsika

corsika/framework/process/ContinuousProcessStepLength.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.
+ * 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,12 +13,10 @@
 namespace corsika {
 
   /**
-     @ingroup Processes
-
-     To store step length in LengthType and unique index in ProcessSequence of shortest
-     step ContinuousProcess.
-
-   **/
+   * @ingroup Processes
+   * To store step length in LengthType and unique index in ProcessSequence of shortest
+   * step ContinuousProcess.
+   */
 
   class ContinuousProcessStepLength {
   public:

corsika/framework/process/DecayProcess.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
@@ -36,7 +35,7 @@ namespace corsika {
      @endcode
 
      Where, of course, SecondaryView and Particle are the valid
-     classes to access particles on the Stack. Those two methods do
+     classes to access particles on the Stack. In user code those two methods do
      not need to be templated, they could use the types
      e.g. corsika::setup::Stack::particle_type -- but by the cost of
      loosing all flexibility otherwise provided.
@@ -50,7 +49,7 @@ namespace corsika {
   template <typename TDerived>
   struct DecayProcess : BaseProcess<TDerived> {
   public:
-    using BaseProcess<TDerived>::ref;
+    using BaseProcess<TDerived>::getRef;
 
     template <typename TParticle>
     InverseTimeType getInverseLifetime(TParticle const& particle) {
@@ -61,7 +60,7 @@ namespace corsika {
                     "getInteractionLength(TParticle const&)\" required for "
                     "InteractionProcess<TDerived>. ");
 
-      return 1. / ref().getLifetime(particle);
+      return 1. / getRef().getLifetime(particle);
     }
   };
 

corsika/framework/process/DynamicInteractionProcess.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.
+ */
+
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <corsika/framework/geometry/FourVector.hpp>
+#include <corsika/framework/core/Logging.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/process/InteractionProcess.hpp>
+#include <corsika/framework/process/ProcessTraits.hpp>
+
+#include <boost/type_index.hpp>
+
+#include <memory>
+
+namespace corsika {
+
+  /**
+   * This class allows selecting/using different InteractionProcesses at runtime without
+   * recompiling the process sequence. The implementation is based on the "type-erasure"
+   * technique.
+   *
+   * @tparam TStack the stack type; has to match the one used by Cascade together with
+   * the ProcessSequence containing the DynamicInteractionProcess.
+   */
+  template <typename TStack>
+  class DynamicInteractionProcess
+      : InteractionProcess<DynamicInteractionProcess<TStack>> {
+  public:
+    using stack_view_type = typename TStack::stack_view_type;
+
+    DynamicInteractionProcess() = default;
+
+    /**
+     * Create new DynamicInteractionProcess. Calls to this instance will be forwared
+     * to the process referred to by obj. The newly created DynamicInteractionProcess
+     * shares ownership of the underlying process, so that you do not need to worry
+     * about it going out of scope.
+     */
+    template <typename TInteractionProcess>
+    DynamicInteractionProcess(std::shared_ptr<TInteractionProcess> obj)
+        : concept_{std::make_unique<ConcreteModel<TInteractionProcess>>(std::move(obj))} {
+      CORSIKA_LOG_DEBUG("creating DynamicInteractionProcess from {}",
+                        boost::typeindex::type_id<TInteractionProcess>().pretty_name());
+
+      // poor man's check while we don't have C++20 concepts
+      static_assert(is_interaction_process_v<TInteractionProcess>);
+
+      // since we only forward interaction process API, all other types of corsika
+      // processes are prohibited
+      static_assert(!is_continuous_process_v<TInteractionProcess>);
+      static_assert(!is_decay_process_v<TInteractionProcess>);
+      static_assert(!is_stack_process_v<TInteractionProcess>);
+      static_assert(!is_cascade_equations_process_v<TInteractionProcess>);
+      static_assert(!is_secondaries_process_v<TInteractionProcess>);
+      static_assert(!is_boundary_process_v<TInteractionProcess>);
+      static_assert(!is_cascade_equations_process_v<TInteractionProcess>);
+    }
+
+    /// forwards arguments to doInteraction() of wrapped instance
+    void doInteraction(stack_view_type& view, Code projectileId, Code targetId,
+                       FourMomentum const& projectileP4, FourMomentum const& targetP4) {
+      return concept_->doInteraction(view, projectileId, targetId, projectileP4,
+                                     targetP4);
+    }
+
+    /// forwards arguments to getCrossSection() of wrapped instance
+    CrossSectionType getCrossSection(Code projectileId, Code targetId,
+                                     FourMomentum const& projectileP4,
+                                     FourMomentum const& targetP4) const {
+      return concept_->getCrossSection(projectileId, targetId, projectileP4, targetP4);
+    }
+
+  private:
+    struct IInteractionModel {
+      virtual ~IInteractionModel() = default;
+      virtual void doInteraction(stack_view_type&, Code, Code, FourMomentum const&,
+                                 FourMomentum const&) = 0;
+      virtual CrossSectionType getCrossSection(Code, Code, FourMomentum const&,
+                                               FourMomentum const&) const = 0;
+    };
+
+    template <typename TModel>
+    struct ConcreteModel final : IInteractionModel {
+      ConcreteModel(std::shared_ptr<TModel> obj) noexcept
+          : model_{std::move(obj)} {}
+
+      void doInteraction(stack_view_type& view, Code projectileId, Code targetId,
+                         FourMomentum const& projectileP4,
+                         FourMomentum const& targetP4) override {
+        return model_->doInteraction(view, projectileId, targetId, projectileP4,
+                                     targetP4);
+      }
+
+      CrossSectionType getCrossSection(Code projectileId, Code targetId,
+                                       FourMomentum const& projectileP4,
+                                       FourMomentum const& targetP4) const override {
+        return model_->getCrossSection(projectileId, targetId, projectileP4, targetP4);
+      }
+
+    private:
+      std::shared_ptr<TModel> model_;
+    };
+
+    std::unique_ptr<IInteractionModel> concept_;
+  };
+} // namespace corsika

corsika/framework/process/InteractionCounter.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
 
 #include <corsika/framework/process/InteractionHistogram.hpp>
 #include <corsika/framework/process/InteractionProcess.hpp>
+#include <corsika/framework/geometry/FourVector.hpp>
 
 namespace corsika {
 
-  /*!
-    @ingroup Processes
-    @{
-
+  /**
+   * @ingroup Processes
+   * @{
+   *
    * Wrapper around an InteractionProcess that fills histograms of the number
    * of calls to `doInteraction()` binned in projectile energy (both in
-   * lab and center-of-mass frame) and species
+   * lab and center-of-mass frame) and species.
    *
-   * Use by wrapping a normal InteractionProcess
+   * Use by wrapping a normal InteractionProcess:
    * @code{.cpp}
    * InteractionProcess collision1;
    * InteractionClounter<collision1> counted_collision1;
    * @endcode
-   *
    */
+
   template <class TCountedProcess>
   class InteractionCounter
       : public InteractionProcess<InteractionCounter<TCountedProcess>> {
@@ -35,17 +35,24 @@ namespace corsika {
   public:
     InteractionCounter(TCountedProcess& process);
 
-    //! wrapper around internall process doInteraction
+    /**
+     * Wrapper around internal process doInteraction.
+     */
     template <typename TSecondaryView>
-    void doInteraction(TSecondaryView& view);
+    void doInteraction(TSecondaryView& view, Code const, Code const, FourMomentum const&,
+                       FourMomentum const&);
 
-    ///! returns internal process getInteractionLength
-    template <typename TParticle>
-    GrammageType getInteractionLength(TParticle const& particle) const;
+    /**
+     * Wrapper around internal process getCrossSection.
+     */
+    CrossSectionType getCrossSection(Code const, Code const, FourMomentum const&,
+                                     FourMomentum const&) const;
 
-    /** returns the filles histograms
-        @return InteractionHistogram, which contains the histogram data
-    */
+    /**
+     * returns the filles histograms.
+     *
+     * @return InteractionHistogram, which contains the histogram data
+     */
     InteractionHistogram const& getHistogram() const;
 
   private:

corsika/framework/process/InteractionHistogram.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
@@ -57,11 +56,11 @@ namespace corsika {
        @param projectile_id corsika::Code of particle
        @param lab_energy Energy in lab. frame
        @param mass_target Mass of target particle
-       @param A if projectile_id is corsika::Nucleus : Mass of nucleus
-       @param Z if projectile_id is corsika::Nucleus : Charge of nucleus
+       @param A if projectile_id is Nucleus : Mass of nucleus
+       @param Z if projectile_id is Nucleus : Charge of nucleus
     */
-    void fill(Code projectile_id, HEPEnergyType lab_energy, HEPEnergyType mass_target,
-              int A = 0, int Z = 0);
+    void fill(Code const projectile_id, HEPEnergyType const lab_energy,
+              HEPEnergyType const mass_target);
 
     //! return histogram in c.m.s. frame
     hist_type const& CMSHist() const { return inthist_cms_; }

corsika/framework/process/InteractionLengthModifier.hpp

+/*
+ * (c) Copyright 2021 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 <functional>
+
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/process/InteractionProcess.hpp>
+
+namespace corsika {
+  /*!
+    @ingroup Processes
+    @{
+
+   * Wrapper around an InteractionProcess that allows modifying the interaction
+   * length returned by the underlying process in a user-defined way.   *
+   */
+  template <class TUnderlyingProcess>
+  class InteractionLengthModifier
+      : public InteractionProcess<InteractionLengthModifier<TUnderlyingProcess>> {
+
+    //! identity function as default modifier
+    static auto constexpr non_modifying_functor = [](GrammageType original, corsika::Code,
+                                                     HEPEnergyType) { return original; };
+
+  public:
+    //! The signature of the modifying functor. Arguments are original int. length, PID,
+    //! energy
+    using functor_signature = GrammageType(GrammageType, corsika::Code, HEPEnergyType);
+
+    /**
+     * Create wrapper around InteractionProcess. Note that the passed process object
+     * itself may no longer be used, only through this class.
+     */
+    InteractionLengthModifier(TUnderlyingProcess&& process,
+                              std::function<functor_signature> modifier);
+
+    //! wrapper around internal process doInteraction
+    template <typename TSecondaryView>
+    void doInteraction(TSecondaryView& view);
+
+    ///! returns underlying process getInteractionLength modified
+    template <typename TParticle>
+    GrammageType getInteractionLength(TParticle const& particle);
+
+    ///! obtain reference to wrapped process
+    TUnderlyingProcess const& getProcess() const;
+    TUnderlyingProcess& getProcess();
+
+  private:
+    TUnderlyingProcess process_;
+    std::function<functor_signature> const modifier_{non_modifying_functor};
+  };
+
+  //! @}
+
+} // namespace corsika
+
+#include <corsika/detail/framework/process/InteractionLengthModifier.inl>

corsika/framework/process/InteractionProcess.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
 
 #include <corsika/framework/process/BaseProcess.hpp>
 #include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/media/NuclearComposition.hpp>
 
 #include <corsika/detail/framework/process/InteractionProcess.hpp> // for extra traits, method/interface checking
 
 namespace corsika {
 
   /**
-     @ingroup Processes
-     @{
-
-     Process describing the interaction of particles
-
-     Create a new InteractionProcess, e.g. for XYModel, via
-     @code
-     class XYModel : public InteractionProcess<XYModel> {};
-     @endcode
-
-     and provide the two necessary interface methods
-     @code
-     template <typename TSecondaryView>
-     void XYModel::doInteraction(TSecondaryView&);
-
-     template <typename TParticle>
-     GrammageType XYModel::getInteractionLength(TParticle const&)
-     @endcode
-
-     Where, of course, SecondaryView and Particle are the valid
-     classes to access particles on the Stack. Those two methods do
-     not need to be templated, they could use the types
-     e.g. corsika::setup::Stack::particle_type -- but by the cost of
-     loosing all flexibility otherwise provided.
-
-     SecondaryView allows to retrieve the properties of the projectile
-     particles, AND to store new particles (secondaries) which then
-     subsequently can be processes by SecondariesProcess. This is how
-     the output of interactions can be studied right away.
-
+   * @ingroup Processes
+   * @{
+   *
+   * Process describing the interaction of particles.
+   *
+   * Create a new InteractionProcess, e.g. for XYModel, via:
+   * @code
+   * class XYModel : public InteractionProcess<XYModel> {};
+   * @endcode
+   *
+   * and provide the two necessary interface methods:
+   * @code
+   * template <typename TSecondaryView>
+   * void XYModel::doInteraction(TSecondaryView&);
+   *
+   * template <typename TParticle>
+   * GrammageType XYModel::getInteractionLength(TParticle const&)
+   * @endcode
+   *
+   * Where, of course, SecondaryView and Particle are the valid
+   * classes to access particles on the Stack. In user code, those two methods do
+   * not need to be templated, they could use the types
+   * e.g. corsika::setup::Stack::particle_type -- but by the cost of
+   * loosing all flexibility otherwise provided.
+   *
+   * SecondaryView allows to retrieve the properties of the projectile
+   * particles, AND to store new particles (secondaries) which then
+   * subsequently can be processes by SecondariesProcess. This is how
+   * the output of interactions can be studied right away.
    */
 
-  template <typename TDerived>
-  class InteractionProcess : public BaseProcess<TDerived> {
+  template <typename TModel>
+  class InteractionProcess : public BaseProcess<TModel> {
 
   public:
-    using BaseProcess<TDerived>::ref;
-
-    template <typename TParticle>
-    InverseGrammageType getInverseInteractionLength(TParticle const& particle) {
-
-      // interface checking on TProcess1
-      static_assert(
-          has_method_getInteractionLength_v<TDerived, GrammageType, TParticle const&>,
-          "TDerived has no method with correct signature \"GrammageType "
-          "getInteractionLength(TParticle const&)\" required for "
-          "InteractionProcess<TDerived>. ");
-
-      return 1. / ref().getInteractionLength(particle);
-    }
+    using BaseProcess<TModel>::getRef;
   };
 
   /**
-   * ProcessTraits specialization to flag InteractionProcess objects
-   **/
+   * ProcessTraits specialization to flag InteractionProcess objects.
+   */
   template <typename TProcess>
   struct is_interaction_process<
       TProcess, std::enable_if_t<

corsika/framework/process/NullModel.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,11 +13,11 @@
 namespace corsika {
 
   /**
-     @ingroup Processes
-     @{
-
-     Process that does nothing. It is not even derived from
-     BaseProcess. But it can be added to a ProcessSequence.
+   * @ingroup Processes
+   * @{
+   *
+   * Process that does nothing. It is not even derived from
+   * BaseProcess. But it can be added to a ProcessSequence.
    */
 
   class NullModel {
@@ -31,23 +30,23 @@ namespace corsika {
     static bool const is_switch_process_sequence = false;
 
     //! Default number of processes is zero, obviously
-    static unsigned int constexpr getNumberOfProcesses() { return 0; }
+    static size_t constexpr getNumberOfProcesses() { return 0; }
   };
 
   /**
-     is_process traits specialization to indicate compatibility with BaseProcess
-  */
+   * is_process traits specialization to indicate compatibility with BaseProcess.
+   */
   template <typename TNull>
   struct is_process<
       TNull, std::enable_if_t<std::is_base_of_v<NullModel, typename std::decay_t<TNull>>>>
       : std::true_type {};
 
   /**
-     count_processes traits specialization to increase process count by one.
+   * count_processes traits specialization to increase process count by one.
    */
   template <int N>
   struct count_processes<NullModel, N, void> {
-    static unsigned int constexpr count = N;
+    static size_t constexpr count = N;
   };
 
   //! @}

corsika/framework/process/ProcessReturn.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/framework/process/ProcessSequence.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
@@ -20,132 +19,138 @@
 #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>
 #include <corsika/framework/process/NullModel.hpp>
+#include <corsika/framework/core/Step.hpp>
 #include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/core/ParticleProperties.hpp>
+
+#include <corsika/framework/geometry/FourVector.hpp>
 
 namespace corsika {
 
+  class COMBoost;           // fwd-decl
+  class NuclearComposition; // fwd-decl
+
   /**
-     count_processes traits specialization to increase process count by
-     getNumberOfProcesses(). This is used to statically count processes in the sequence
-  **/
+   * count_processes traits specialization to increase process count by
+   * getNumberOfProcesses(). This is used to statically count processes in the sequence.
+   */
   template <typename TProcess, int N>
   struct count_processes<
       TProcess, N,
       typename std::enable_if_t<is_process_v<std::decay_t<TProcess>> &&
                                 std::decay_t<TProcess>::is_process_sequence>> {
-    static unsigned int constexpr count =
-        N + std::decay_t<TProcess>::getNumberOfProcesses();
+    static size_t constexpr count = N + std::decay_t<TProcess>::getNumberOfProcesses();
   };
 
   /**
-     @defgroup Processes Physics Processes and Modules
-
-     Physics processes in CORSIKA 8 are clustered in ProcessSequence and
-     SwitchProcessSequence containers. The former is a mere (ordered) collection, while
-     the latter has the option to switch between two alternative ProcessSequences.
-
-     Depending on the type of data to act on and on the allowed actions of processes there
-     are several interface options:
-     - InteractionProcess
-     - DecayProcess
-     - ContinuousProcess
-     - StackProcess
-     - SecondariesProcess
-     - BoundaryCrossingProcess
-
-     And all processes (including ProcessSequence and SwitchProcessSequence) are derived
-     from BaseProcess.
-
-     Processes of any type (e.g. p1, p2, p3,...) can be assembled into a ProcessSequence
-     using the `make_sequence` factory function.
-
-     @code{.cpp}
-       auto sequence1 = make_sequence(p1, p2, p3);
-       auto sequence2 = make_sequence(p4, p5, p6, p7);
-       auto sequence3 = make_sequence(sequence1, sequemce2, p8, p9);
-     @endcode
-
-     Note, if the order of processes
-     matters, the order of occurence
-     in the ProcessSequence determines
-     the executiion order.
-
-     SecondariesProcess alyways act on
-     new secondaries produced (i.e. in
-     InteractionProcess and
-     DecayProcess) in the scope of
-     their ProcessSequence. For
-     example if i1 and i2 are
-     InteractionProcesses and s1 is a
-     SecondariesProcess, then
-
-     @code{.cpp}
-       auto sequence = make_sequence(i1, make_sequence(i2, s1))
-     @endcode
-
-     will result in s1 acting only on
-     the particles produced by i2 and
-     not by i1. This can be very
-     useful, e.g. to fine tune thinning.
-
-     A special type of ProcessSequence
-     is SwitchProcessSequence, which
-     has two branches and a functor
-     that can select between these two
-     branches.
-
-     @code{.cpp}
-       auto sequence = make_switch(sequence1, sequence2, selector);
-     @endcode
-
-     where the only requirement to
-     `selector` is that it
-     provides a `SwitchResult operator()(Particle const& particle) const` method. Thus,
-     based on the dynamic properties
-     of `particle` the functor
-     can make its decision. This is
-     clearly important for switching
-     between low-energy and
-     high-energy models, but not
-     limited to this. The selection
-     can even be done with a lambda
-     function.
-
-
-     @class ProcessSequence
-     @ingroup Processes
-
-       Definition of a static process list/sequence
-
-       A compile time static list of processes. The compiler will
-       generate a new type based on template logic containing all the
-       elements provided by the user.
-
-       TProcess1 and TProcess2 must both be derived from BaseProcess,
-       and are both references if possible (lvalue), otherwise (rvalue)
-       they are just classes. This allows us to handle both, rvalue as
-       well as lvalue Processes in the ProcessSequence.
-
-       (For your potential interest,
-       the static version of the
-       ProcessSequence and all Process
-       types are based on the CRTP C++
-       design pattern)
-
-      Template parameters:
-        @tparam TProcess1 is of type BaseProcess, either a dedicatd process, or a
-     ProcessSequence
-        @tparam TProcess2 is of type BaseProcess, either a dedicatd process, or a
-     ProcessSequence
-      @tparam IndexFirstProcess to count and index each Process in the entire
-  process-chain. The offset is the starting value for this ProcessSequence
-      @tparam IndexOfProcess1 index of TProcess1 (counting of Process)
-      @tparam IndexOfProcess2 index of TProcess2 (counting of Process)
-  **/
+   * @defgroup Processes Physics Processes and Modules
+   *
+   * Physics processes in CORSIKA 8 are clustered in ProcessSequence and
+   * SwitchProcessSequence containers. The former is a mere (ordered) collection, while
+   * the latter has the option to switch between two alternative ProcessSequences.
+   *
+   * Depending on the type of data to act on and on the allowed actions of processes
+   * there are several interface options:
+   * - InteractionProcess
+   * - DecayProcess
+   * - ContinuousProcess
+   * - StackProcess
+   * - SecondariesProcess
+   * - BoundaryCrossingProcess
+   *
+   * And all processes (including ProcessSequence and SwitchProcessSequence) are derived
+   * from BaseProcess.
+   *
+   * Processes of any type (e.g. p1, p2, p3,...) can be assembled into a ProcessSequence
+   * using the `make_sequence` factory function.
+   *
+   * @code
+   *   auto sequence1 = make_sequence(p1, p2, p3);
+   *   auto sequence2 = make_sequence(p4, p5, p6, p7);
+   *   auto sequence3 = make_sequence(sequence1, sequemce2, p8, p9);
+   * @endcode
+   *
+   * Note, if the order of processes
+   * matters, the order of occurence
+   * in the ProcessSequence determines
+   * the executiion order.
+   *
+   * SecondariesProcess alyways act on
+   * new secondaries produced (i.e. in
+   * InteractionProcess and
+   * DecayProcess) in the scope of
+   * their ProcessSequence. For
+   * example if i1 and i2 are
+   * InteractionProcesses and s1 is a
+   * SecondariesProcess, then:
+   *
+   * @code{.cpp}
+   *   auto sequence = make_sequence(i1, make_sequence(i2, s1))
+   * @endcode
+   *
+   * will result in s1 acting only on
+   * the particles produced by i2 and
+   * not by i1. This can be very
+   * useful, e.g. to fine tune thinning.
+   *
+   * A special type of ProcessSequence
+   * is SwitchProcessSequence, which
+   * has two branches and a functor
+   * that can select between these two
+   * branches.
+   *
+   * @code{.cpp}
+   *   auto sequence = make_switch(sequence1, sequence2, selector);
+   * @endcode
+   *
+   * where the only requirement to
+   * `selector` is that it
+   * provides a `SwitchResult operator()(Particle const& particle) const` method. Thus,
+   * based on the dynamic properties
+   * of `particle` the functor
+   * can make its decision. This is
+   * clearly important for switching
+   * between low-energy and
+   * high-energy models, but not
+   * limited to this. The selection
+   * can even be done with a lambda
+   * function.
+   *
+   * @class ProcessSequence
+   * @ingroup Processes
+   *
+   *   Definition of a static process list/sequence.
+   *
+   *   A compile time static list of processes. The compiler will
+   *   generate a new type based on template logic containing all the
+   *   elements provided by the user.
+   *
+   *   TProcess1 and TProcess2 must both be derived from BaseProcess,
+   *   and are both references if possible (lvalue), otherwise (rvalue)
+   *   they are just classes. This allows us to handle both, rvalue as
+   *   well as lvalue Processes in the ProcessSequence.
+   *
+   *   (For your potential interest,
+   *   the static version of the
+   *   ProcessSequence and all Process
+   *   types are based on the CRTP C++
+   *   design pattern).
+   *
+   *  Template parameters:
+   *  @tparam TProcess1 is of type BaseProcess, either a dedicatd process, or a
+   *           ProcessSequence
+   *  @tparam TProcess2 is of type BaseProcess, either a dedicatd process, or a
+   *           ProcessSequence
+   *  @tparam IndexFirstProcess to count and index each Process in the entire
+   *          process-chain. The offset is the starting value for this ProcessSequence
+   *  @tparam IndexOfProcess1 index of TProcess1 (counting of Process)
+   *  @tparam IndexOfProcess2 index of TProcess2 (counting of Process)
+   */
 
   template <typename TProcess1, typename TProcess2 = NullModel,
             int ProcessIndexOffset = 0,
@@ -170,91 +175,155 @@ namespace corsika {
     static bool const is_process_sequence = true;
 
     /**
-      Only valid user constructor will create fully initialized object
-
-      ProcessSequence supports and encourages move semantics. You can
-      use object, l-value references or r-value references to
-      construct sequences.
-
-      @param in_A BaseProcess or switch/process list
-      @param in_B BaseProcess or switch/process list
-     **/
+     * Only valid user constructor will create fully initialized object.
+     *
+     * ProcessSequence supports and encourages move semantics. You can
+     * use object, l-value references or r-value references to
+     * construct sequences.
+     *
+     * @param in_A BaseProcess or switch/process list
+     * @param in_B BaseProcess or switch/process list
+     */
     ProcessSequence(TProcess1 in_A, TProcess2 in_B);
 
+    /**
+     * List of all BoundaryProcess.
+     *
+     * @tparam TParticle
+     * @param particle The particle.
+     * @param from Volume the particle is exiting.
+     * @param to Volume the particle is entering.
+     * @return ProcessReturn
+     */
     template <typename TParticle>
     ProcessReturn doBoundaryCrossing(TParticle& particle,
                                      typename TParticle::node_type const& from,
                                      typename TParticle::node_type const& to);
 
-    template <typename TParticle, typename TTrack>
-    ProcessReturn doContinuous(TParticle& particle, TTrack& vT,
+    template <typename TParticle>
+    ProcessReturn doContinuous(Step<TParticle>& step,
                                ContinuousProcessIndex const limitID);
 
+    /**
+     * Process all secondaries in TSecondaries.
+     *
+     * The seondaries produced by other processes and accessible via TSecondaries
+     * are processed by all SecondariesProcesse via a call here.
+     *
+     * @tparam TSecondaries
+     * @param vS
+     */
     template <typename TSecondaries>
     void doSecondaries(TSecondaries& vS);
 
     /**
-       The processes of type StackProcess do have an internal counter,
-       so they can be exectuted only each N steps. Often these are
-       "maintenacne processes" that do not need to run after each
-       single step of the simulations. In the CheckStep function it is
-       tested if either A_ or B_ are StackProcess and if they are due
-       for execution.
+     * The processes of type StackProcess do have an internal counter,
+     * so they can be exectuted only each N steps. Often these are
+     * "maintenacne processes" that do not need to run after each
+     * single step of the simulations. In the CheckStep function it is
+     * tested if either A_ or B_ are StackProcess and if they are due
+     * for execution.
      */
     bool checkStep();
 
     /**
-       Execute the StackProcess-es in the ProcessSequence
+     * Execute the StackProcess-es in the ProcessSequence.
      */
     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
+     * ContinuousProcess-es.
      *
      * The maximum allowed step length is the minimum of the allowed track lenght over all
      * ContinuousProcess-es in the ProcessSequence.
      *
+     * @tparam TParticle particle type.
+     * @tparam TTrack the trajectory type.
+     * @param particle The particle data object.
+     * @param track The track data object.
      * @return ContinuousProcessStepLength which contains the step length itself in
      *          LengthType, and a unique identifier of the related ContinuousProcess.
-     **/
+     */
 
     template <typename TParticle, typename TTrack>
-    ContinuousProcessStepLength getMaxStepLength(TParticle& particle, TTrack& vTrack);
-
-    template <typename TParticle>
-    GrammageType getInteractionLength(TParticle&& particle) {
-      return 1. / getInverseInteractionLength(particle);
-    }
+    ContinuousProcessStepLength getMaxStepLength(TParticle&& particle, TTrack&& vTrack);
 
+    /**
+     * @brief Calculates the cross section of a projectile with a target.
+     *
+     * @tparam TParticle
+     * @param projectile
+     * @param targetId
+     * @param targetP4
+     * @return CrossSectionType
+     */
     template <typename TParticle>
-    InverseGrammageType getInverseInteractionLength(TParticle&& particle);
-
-    template <typename TSecondaryView>
-    ProcessReturn selectInteraction(
-        TSecondaryView& view, [[maybe_unused]] InverseGrammageType lambda_inv_select,
-        [[maybe_unused]] InverseGrammageType lambda_inv_sum =
-            InverseGrammageType::zero());
+    CrossSectionType getCrossSection(TParticle const& projectile, Code const targetId,
+                                     FourMomentum const& targetP4) const;
 
     template <typename TParticle>
     TimeType getLifetime(TParticle& particle) {
       return 1. / getInverseLifetime(particle);
     }
 
+    /**
+     * Selects one concrete InteractionProcess and samples a target nucleus from
+     * the material.
+     *
+     * The selectInteraction method statically loops over all active InteractionProcess
+     * and calculates the material-weighted cross section for all of them. In an iterative
+     * way those cross sections are summed up. The random number cx_select, uniformely
+     * drawn from the cross section before energy losses, is used to discriminate the
+     * selected sub-process here. If the cross section after the step smaller than it was
+     * before, there is a non-zero probability that the particle survives and no
+     * interaction takes place. This method becomes imprecise when cross section rise with
+     * falling energies.
+     *
+     * If a sub-process was selected, the target nucleus is selected from the material
+     * (weighted with cross section). The interaction is then executed.
+     *
+     * @tparam TSecondaryView Object type as storage for new secondary particles.
+     * @tparam TRNG Object type to produce random numbers.
+     * @param view Object to store new secondary particles.
+     * @param projectileP4 The four momentum of the projectile.
+     * @param composition The environment/material composition.
+     * @param rng Random number object.
+     * @param cx_select Drawn random numer, uniform between [0, cx_initial]
+     * @param cx_sum For interal use, to sum up cross section contributions.
+     * @return ProcessReturn
+     */
+    template <typename TSecondaryView, typename TRNG>
+    inline ProcessReturn selectInteraction(
+        TSecondaryView&& view, FourMomentum const& projectileP4,
+        NuclearComposition const& composition, TRNG&& rng,
+        CrossSectionType const cx_select,
+        CrossSectionType cx_sum = CrossSectionType::zero());
+
     template <typename TParticle>
     InverseTimeType getInverseLifetime(TParticle&& particle);
 
     // select decay process
     template <typename TSecondaryView>
-    ProcessReturn selectDecay(
-        TSecondaryView& view, [[maybe_unused]] InverseTimeType decay_inv_select,
-        [[maybe_unused]] InverseTimeType decay_inv_sum = InverseTimeType::zero());
+    ProcessReturn selectDecay(TSecondaryView&& view, InverseTimeType decay_inv_select,
+                              InverseTimeType decay_inv_sum = InverseTimeType::zero());
 
     /**
      * static counter to uniquely index (count) all ContinuousProcess in switch sequence.
-     **/
-    static unsigned int constexpr getNumberOfProcesses() { return numberOfProcesses_; }
+     */
+    static size_t constexpr getNumberOfProcesses() { return numberOfProcesses_; }
 
 #ifdef CORSIKA_UNIT_TESTING
     TProcess1 getProcess1() const { return A_; }
@@ -262,40 +331,40 @@ namespace corsika {
 #endif
 
   private:
-    TProcess1 A_; /// process/list A, this is a reference, if possible
-    TProcess2 B_; /// process/list B, this is a reference, if possible
+    TProcess1 A_; //! process/list A, this is a reference, if possible
+    TProcess2 B_; //! process/list B, this is a reference, if possible
 
-    static unsigned int constexpr numberOfProcesses_ = IndexOfProcess1; // static counter
+    static size_t constexpr numberOfProcesses_ = IndexOfProcess1; // static counter
   };
 
   /**
-    @fn make_sequence
-    @ingroup Processes
-
-    Factory function to create a ProcessSequence
-
-    to construct ProcessSequences in a flexible and dynamic way the
-    `sequence` factory functions are provided
-
-    Any objects of type
-     - BaseProcess
-     - ContinuousProcess and
-     - InteractionProcess/DecayProcess
-     - StackProcess
-     - SecondariesProcess
-     - BoundaryCrossingProcess
-
-    can be assembled into a ProcessSequence, all
-    combinatorics are allowed.
-
-    The sequence function checks that all its arguments are all of
-    types derived from BaseProcess. Also the ProcessSequence itself
-    is derived from type BaseProcess
-
-    @tparam TProcesses parameter pack with objects of type BaseProcess
-    @tparam TProcess1 another BaseProcess
-    @param vA needs to derive from BaseProcess
-    @param vB paramter-pack, needs to derive BaseProcess
+   * @fn make_sequence
+   * @ingroup Processes
+   *
+   * Factory function to create a ProcessSequence.
+   *
+   * to construct ProcessSequences in a flexible and dynamic way the
+   * `sequence` factory functions are provided.
+   *
+   * Any objects of type
+   *  - BaseProcess
+   *  - ContinuousProcess and
+   *  - InteractionProcess/DecayProcess
+   *  - StackProcess
+   *  - SecondariesProcess
+   *  - BoundaryCrossingProcess
+   *
+   * can be assembled into a ProcessSequence, all
+   * combinatorics are allowed.
+   *
+   * The sequence function checks that all its arguments are all of
+   * types derived from BaseProcess. Also the ProcessSequence itself
+   * is derived from type BaseProcess.
+   *
+   * @tparam TProcesses parameter pack with objects of type BaseProcess
+   * @tparam TProcess1 another BaseProcess
+   * @param vA needs to derive from BaseProcess
+   * @param vB paramter-pack, needs to derive BaseProcess
    */
   template <typename... TProcesses, typename TProcess1>
   ProcessSequence<TProcess1, decltype(make_sequence(std::declval<TProcesses>()...))>
@@ -306,34 +375,34 @@ namespace corsika {
   }
 
   /**
-    @fn make_sequence
-    @ingroup Processes
-
-    Factory function to create ProcessSequence
-
-    specialization for two input objects (no paramter pack in vB).
-
-    @tparam TProcess1 another BaseProcess
-    @tparam TProcess2 another BaseProcess
-    @param vA needs to derive from BaseProcess
-    @param vB needs to derive BaseProcess
-  */
+   * @fn make_sequence
+   * @ingroup Processes
+   *
+   * Factory function to create ProcessSequence.
+   *
+   * specialization for two input objects (no paramter pack in vB).
+   *
+   * @tparam TProcess1 another BaseProcess
+   * @tparam TProcess2 another BaseProcess
+   * @param vA needs to derive from BaseProcess
+   * @param vB needs to derive BaseProcess
+   */
   template <typename TProcess1, typename TProcess2>
   ProcessSequence<TProcess1, TProcess2> make_sequence(TProcess1&& vA, TProcess2&& vB) {
     return ProcessSequence<TProcess1, TProcess2>(vA, vB);
   }
 
   /**
-    @fn make_sequence
-    @ingroup Processes
-
-    Factory function to create ProcessSequence from a single BaseProcess
-
-    also allow a single Process in ProcessSequence, accompany by
-    `NullModel`
-
-    @tparam TProcess1 another BaseProcess
-    @param vA needs to derive from BaseProcess
+   * @fn make_sequence
+   * @ingroup Processes
+   *
+   * Factory function to create ProcessSequence from a single BaseProcess.
+   *
+   * also allow a single Process in ProcessSequence, accompany by
+   * `NullModel`.
+   *
+   * @tparam TProcess1 another BaseProcess
+   * @param vA needs to derive from BaseProcess
    */
   template <typename TProcess>
   ProcessSequence<TProcess, NullModel> make_sequence(TProcess&& vA) {

corsika/framework/process/ProcessTraits.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
@@ -13,11 +12,12 @@
  */
 
 #include <type_traits>
+#include <cstddef>
 
 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,40 +99,11 @@ 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;
+    static size_t constexpr count = N;
   };
 
-  namespace detail {
-
-    /**
-       Helper traits class (partial) for static compile time checking.
-
-       Note, this is a poor replacement for C++20 concepts... they are
-       eagerly awaited!
-
-       It defines the default body of a generic test function returning
-       std::false_type.
-
-       In addition it defines the pattern for class-method matching with a
-       return type TReturn and function arguments TArgs... . Right now
-       both method signatures, "const" and "not const", are matched.
-     */
-    template <typename TReturn, typename... TArgs>
-    struct has_method_signature {
-
-      // the non-const version
-      template <class T>
-      static std::true_type testSignature(TReturn (T::*)(TArgs...));
-
-      // the const version
-      template <class T>
-      static std::true_type testSignature(TReturn (T::*)(TArgs...) const);
-    };
-
-  } // namespace detail
-
 } // namespace corsika

corsika/framework/process/SecondariesProcess.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
@@ -16,24 +15,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 +41,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

 /*
  * (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
@@ -16,32 +15,38 @@
 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 (@f$n_{step}@f$) 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.
+   *
+   * The number of *steps* during the cascade processing after
+   * which a StackProcess is going to be executed is determined from `iStep_` and `nStep_`
+   * using the modulo: @f$ !(iStep\_ \; \% \; nStep\_) @f$. And `iStep_` is increased
+   * for each evaluation (step).
    */
 
   template <typename TDerived>
@@ -56,7 +61,7 @@ namespace corsika {
     static bool const is_stack_process = true;
 
     //! return the current Cascade step counter
-    int getStep() const { return iStep_; }
+    auto getStep() const { return iStep_; }
 
     //! check if current step is where StackProcess should be executed, this also
     //! increases the internal step counter implicitly
@@ -64,10 +69,8 @@ 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 Execution control and counters.
+     * @{
      */
     unsigned int nStep_ = 0;
     unsigned long int iStep_ = 0;
@@ -75,8 +78,8 @@ namespace corsika {
   };
 
   /**
-   * ProcessTraits specialization to flag StackProcess objects
-   **/
+   * ProcessTraits specialization to flag StackProcess objects.
+   */
   template <typename TProcess>
   struct is_stack_process<
       TProcess,

corsika/framework/process/SwitchProcessSequence.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
@@ -17,12 +16,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>
@@ -115,16 +117,16 @@ namespace corsika {
     static bool const is_switch_process_sequence = true;
 
     /**
-      Only valid user constructor will create fully initialized object
-
-      SwitchProcessSequence supports and encourages move semantics. You can
-      use object, l-value references or r-value references to
-      construct sequences.
-
-      @param sel functor to switch between branch A and B
-      @param in_A process branch A
-      @param in_A process branch B
-     **/
+     * Only valid user constructor will create fully initialized object.
+     *
+     * SwitchProcessSequence supports and encourages move semantics. You can
+     * use object, l-value references or r-value references to
+     * construct sequences.
+     *
+     * @param sel functor to switch between branch A and B
+     * @param in_A process branch A
+     * @param in_B process branch B
+     */
     SwitchProcessSequence(TCondition sel, TSequence in_A, USequence in_B)
         : select_(sel)
         , A_(in_A)
@@ -135,8 +137,8 @@ namespace corsika {
                                      typename TParticle::node_type const& from,
                                      typename TParticle::node_type const& to);
 
-    template <typename TParticle, typename TTrack>
-    ProcessReturn doContinuous(TParticle& particle, TTrack& vT,
+    template <typename TParticle>
+    ProcessReturn doContinuous(Step<TParticle>& particle,
                                ContinuousProcessIndex const limitId);
 
     template <typename TSecondaries>
@@ -146,18 +148,15 @@ namespace corsika {
     ContinuousProcessStepLength getMaxStepLength(TParticle& particle, TTrack& vTrack);
 
     template <typename TParticle>
-    GrammageType getInteractionLength(TParticle&& particle) {
-      return 1. / getInverseInteractionLength(particle);
-    }
-
-    template <typename TParticle>
-    InverseGrammageType getInverseInteractionLength(TParticle&& particle);
+    CrossSectionType getCrossSection(TParticle const& projectile, Code const targetId,
+                                     FourMomentum const& targetP4) const;
 
-    template <typename TSecondaryView>
-    ProcessReturn selectInteraction(
-        TSecondaryView& view, [[maybe_unused]] InverseGrammageType lambda_inv_select,
-        [[maybe_unused]] InverseGrammageType lambda_inv_sum =
-            InverseGrammageType::zero());
+    template <typename TSecondaryView, typename TRNG>
+    ProcessReturn selectInteraction(TSecondaryView& view,
+                                    FourMomentum const& projectileP4,
+                                    NuclearComposition const& composition, TRNG& rng,
+                                    CrossSectionType const cx_select,
+                                    CrossSectionType cx_sum = CrossSectionType::zero());
 
     template <typename TParticle>
     TimeType getLifetime(TParticle&& particle) {
@@ -175,8 +174,8 @@ namespace corsika {
 
     /**
      * static counter to uniquely index (count) all ContinuousProcess in switch sequence.
-     **/
-    static unsigned int constexpr getNumberOfProcesses() { return numberOfProcesses_; }
+     */
+    static size_t constexpr getNumberOfProcesses() { return numberOfProcesses_; }
 
 #ifdef CORSIKA_UNIT_TESTING
     TCondition getCondition() const { return select_; }

corsika/framework/random/ExponentialDistribution.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,16 @@
 
 namespace corsika {
 
+  /**
+   * Describes a random distribution with \f[ \beta e^{-\beta X} \f] for a physical
+   * quantity of type Quantity.
+   *
+   * @tparam Quantity is the type of the physical quantity.
+   */
+
   template <typename Quantity>
   class ExponentialDistribution {
+    static_assert(is_quantity_v<Quantity>, "usable only with Quantity types");
 
     typedef typename Quantity::value_type real_type;
     typedef std::exponential_distribution<real_type> distribution_type;
@@ -24,7 +31,7 @@ namespace corsika {
 
     ExponentialDistribution() = delete;
 
-    ExponentialDistribution(value_type const& beta)
+    ExponentialDistribution(value_type beta)
         : beta_(beta) {}
 
     ExponentialDistribution(ExponentialDistribution<value_type> const& other)
@@ -38,30 +45,25 @@ namespace corsika {
     }
 
     /**
-     * @fn value_type getBeta()const
-     * @brief Get parameter of exponential distribution \f[ \beta e^{-X}\f]
-     * @pre
-     * @post
+     * @fn value_type getBeta() const
+     * @brief Get parameter of exponential distribution \f[ \beta e^{-\beta X}\f].
+     *
      * @return  value_type
      */
     value_type getBeta() const { return beta_; }
 
     /**
      * @fn void setBeta(value_type)
-     * @brief Set parameter of exponential distribution \f[ \beta e^{-X}\f]
+     * @brief Set parameter of exponential distribution \f[ \beta e^{-\beta X}\f].
      *
-     * @pre
-     * @post
      * @param vBeta
      */
     void setBeta(value_type const& beta) { beta_ = beta; }
 
     /**
      * @fn value_type operator ()(Generator&)
-     * @brief Generate a random number distributed like \f[ \beta e^{-X}\f]
+     * @brief Generate a random number distributed like \f[ \beta e^{-\beta X}\f].
      *
-     * @pre
-     * @post
      * @tparam Generator
      * @param g
      * @return

corsika/framework/random/PowerLawDistribution.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 <corsika/framework/core/PhysicalUnits.hpp>
+#include <random>
+
+namespace corsika {
+
+  /**
+   * Describes a random distribution with \f[ p(x) \sim x^{\alpha} \f] with spectral index
+   * \f[ \alpha \f] for a physical quantity \f[x\f] of type Quantity within an interval
+   * \f[ [x_0, x_1) \f].
+   *
+   * @tparam Quantity is the type of the physical quantity.
+   */
+
+  template <typename Quantity>
+  class PowerLawDistribution {
+    static_assert(is_quantity_v<Quantity>, "usable only with Quantity types");
+
+    using real_type = typename Quantity::value_type;
+    using uniform_distribution_type = std::uniform_real_distribution<real_type>;
+
+  public:
+    using value_type = Quantity;
+
+    PowerLawDistribution(real_type index /** spectral index */,
+                         value_type low /** lower bound */,
+                         value_type high /** higher bound */)
+        : low_{low}
+        , high_{high}
+        , index_{index} {}
+
+    PowerLawDistribution(PowerLawDistribution<value_type> const& other) = default;
+    PowerLawDistribution<value_type>& operator=(
+        PowerLawDistribution<value_type> const& other) = default;
+
+    template <class Generator>
+    value_type operator()(Generator& g) {
+      real_type const u = dist_(g);
+      if (index_ == static_cast<real_type>(-1)) {
+        return low_ * pow(high_ / low_, u);
+      } else {
+        return value_type{phys::units::detail::magnitude_tag,
+                          pow(pow(low_.magnitude(), index_ + 1) * (1 + u) -
+                                  pow(high_.magnitude(), index_ + 1),
+                              1 / (index_ + 1))};
+      }
+    }
+
+  private:
+    uniform_distribution_type dist_{};
+    value_type low_, high_;
+    real_type index_;
+  };
+
+} // namespace corsika

corsika/framework/random/RNGManager.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,27 +14,27 @@
 
 #include <corsika/framework/utility/Singleton.hpp>
 #include <corsika/framework/core/Logging.hpp>
-
+#include <corsika/detail/framework/random/random_iterator/Stream.hpp>
 /*!
  * With this class modules can register streams of random numbers.
  */
 
 namespace corsika {
 
-  class RNGManager : public corsika::Singleton<RNGManager> {
+  template <typename CBPRNG = random_iterator::philox>
+  class RNGManager : public corsika::Singleton<RNGManager<CBPRNG>> {
 
-    friend class corsika::Singleton<RNGManager>;
+    friend class corsika::Singleton<RNGManager<CBPRNG>>;
 
   public:
-    typedef std::mt19937_64 prng_type;
+    typedef CBPRNG prng_type;
     typedef std::uint64_t seed_type;
     typedef std::string string_type;
     typedef std::map<std::string, prng_type> streams_type;
-    typedef std::map<std::string, std::seed_seq> seeds_type;
 
-    RNGManager(RNGManager const&) = delete; // since it is a singleton
+    RNGManager(RNGManager<CBPRNG> const&) = delete; // since it is a singleton
 
-    RNGManager& operator=(RNGManager const&) = delete;
+    RNGManager<CBPRNG>& operator=(RNGManager<CBPRNG> const&) = delete;
 
     /*!
      * This function is to be called by a module requiring a random-number
@@ -43,35 +42,24 @@ namespace corsika {
      *
      * \throws sth. when stream \a pModuleName is already registered
      */
-    void registerRandomStream(string_type const& streamName);
+    inline void registerRandomStream(string_type const& streamName);
 
     /*!
      * returns the pre-stored stream of given name \a pStreamName if
      * available
      */
-    prng_type& getRandomStream(string_type const& streamName);
+    inline prng_type& getRandomStream(string_type const& streamName);
 
     /*!
      * Check whether a stream has been registered.
      */
-    bool isRegistered(string_type const& streamName) const;
+    inline bool isRegistered(string_type const& streamName) const;
 
     /*!
      * dumps the names and states of all registered random-number streams
      * into a std::stringstream.
      */
-    std::stringstream dumpState() const;
-
-    /**
-     * Set explicit seeds for all currently registered streams. The actual seed values
-     * are incremented from \a vSeed.
-     */
-    void seedAll(seed_type seed);
-
-    /**
-     * Set seeds for all currently registered streams.
-     */
-    void seedAll(void); //!< seed all currently registered streams with "real" randomness
+    inline std::stringstream dumpState() const;
 
     /**
      * @fn const streams_type getRngs&()const
@@ -83,16 +71,6 @@ namespace corsika {
      */
     const streams_type& getRngs() const { return rngs_; }
 
-    /**
-     * @fn const seeds_type getSeeds&()const
-     * @brief Constant access to the seeds.
-     *
-     * @pre
-     * @post
-     * @return RNGManager::seeds_type
-     */
-    const seeds_type& getSeeds() const { return seeds_; }
-
     /**
      * @fn streams_type Rngs()
      * @brief Non-constant access to the streams.
@@ -101,27 +79,25 @@ namespace corsika {
      * @post
      * @return RNGManager::streams_type&
      */
-    streams_type Rngs() { return rngs_; }
+    streams_type& Rngs() { return rngs_; }
 
-    /**
-     * @fn seeds_type Seeds&()
-     * @brief Non-constant access to seeds.
-     *
-     * @pre
-     * @post
-     * @return RNGManager::seeds_type&
-     */
-    seeds_type& getSeeds() { return seeds_; }
+    seed_type getSeed() const { return seed_; }
+
+    void setSeed(seed_type seed) {
+      seed_ = seed;
+      // update the rgn states
+      for (auto& [streamName, rng] : rngs_) rng.setSeed(seed_);
+    }
 
   protected:
     RNGManager() = default;
 
   private:
     streams_type rngs_;
-    seeds_type seeds_;
+    seed_type seed_;
   };
 
-  typedef typename RNGManager::prng_type default_prng_type;
+  typedef typename RNGManager<>::prng_type default_prng_type;
 
 } // namespace corsika
 

corsika/framework/random/UniformRealDistribution.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,6 +14,7 @@ namespace corsika {
 
   template <typename Quantity>
   class UniformRealDistribution {
+    static_assert(is_quantity_v<Quantity>, "usable only with Quantity types");
 
     typedef typename Quantity::value_type real_type;
     typedef std::uniform_real_distribution<real_type> distribution_type;
@@ -24,86 +24,71 @@ namespace corsika {
 
     UniformRealDistribution() = delete;
 
-    UniformRealDistribution(Quantity const& b)
-        : min_{value_type(phys::units::detail::magnitude_tag, 0)}
-        , max_(b) {}
+    UniformRealDistribution(value_type b)
+        : min_{value_type::zero()}
+        , diff_{b} {}
 
-    UniformRealDistribution(value_type const& pmin, value_type const& pmax)
+    UniformRealDistribution(value_type pmin, value_type pmax)
         : min_(pmin)
-        , max_(pmax) {}
+        , diff_{pmax - pmin} {}
 
     UniformRealDistribution(UniformRealDistribution<value_type> const& other)
-        : min_(other.getMin())
-        , max_(other.getMax()) {}
+        : min_{other.min_}
+        , diff_{other.diff_} {}
 
     UniformRealDistribution<value_type>& operator=(
         UniformRealDistribution<value_type> const& other) {
       if (this == &other) return *this;
-      min_ = other.getMin();
-      max_ = other.getMax();
+      min_ = other.min_;
+      diff_ = other.diff_;
       return *this;
     }
 
     /**
-     * @fn quantity_type getMax()const
-     * @brief Get the upper limit.
+     * Get the upper limit.
      *
-     * @pre
-     * @post
-     * @return quantity_type
+     * @return value_type
      */
-    value_type getMax() const { return max_; }
+    value_type getMax() const { return min_ + diff_; }
 
     /**
-     * @fn void setMax(quantity_type)
-     * @brief Set the upper limit.
+     * Set the upper bound.
      *
-     * @pre
-     * @post
-     * @param vMax
+     * @param pmax: new upper bound
      */
-    void setMax(value_type const& pmax) { max_ = pmax; }
+    void setMax(value_type const pmax) { diff_ = pmax - min_; }
 
     /**
-     * @fn quantity_type getMin()const
-     * @brief Get the lower limit.
+     * Get the lower bound.
      *
-     * @pre
-     * @post
-     * @return
+     * @return The minimum possible value.
      */
     value_type getMin() const { return min_; }
 
     /**
-     * @fn void setMin(quantity_type)
-     * @brief  Set the lower limit.
+     * Set the lower bound.
      *
-     * @pre
-     * @post
-     * @param vMin
+     * @param pmin is the new lower bound.
      */
-    void setMin(value_type const& pmin) { min_ = pmin; }
+    void setMin(value_type const pmin) { min_ = pmin; }
 
     /**
-     * @fn quantity_type operator ()(Generator&)
-     * @brief Generate a random numberin the range [min, max]
+     * Generate a random number in the range [min, max).
      *
-     * @pre
-     * @post
-     * @tparam Generator
+     * @tparam TGenerator
      * @param g
-     * @return quantity_type
+     * @return value_type
      */
-    template <class Generator>
-    value_type operator()(Generator& g) {
-      return min_ + dist_(g) * (max_ - min_);
+    template <class TGenerator>
+    value_type operator()(TGenerator& g) {
+      return min_ + dist_(g) * diff_;
     }
 
   private:
     distribution_type dist_{real_type(0.), real_type(1.)};
 
     value_type min_;
-    value_type max_;
+    value_type diff_;
   };
 
 } // namespace corsika