corsika/detail/framework/process/ContinuousProcess.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 <corsika/framework/process/ProcessTraits.hpp>
+#include <corsika/framework/utility/HasMethodSignature.hpp>
+
+namespace corsika {
+
+  /**
+     traits test for ContinuousProcess::doContinuous method
+  */
+  template <class TProcess, typename TReturn, typename TArg1, typename TArg2>
+  struct has_method_doContinuous
+      : public detail::has_method_signature<TReturn, TArg1, TArg2, bool> {
+
+    //! method signature
+    using detail::has_method_signature<TReturn, TArg1, TArg2, bool>::testSignature;
+
+    //! the default value
+    template <class T>
+    static std::false_type test(...);
+
+    //! templated parameter option
+    template <class T>
+    static decltype(testSignature(&T::template doContinuous<TArg1, TArg2>)) test(
+        std::nullptr_t);
+
+    //! non templated parameter option
+    template <class T>
+    static decltype(testSignature(&T::doContinuous)) test(std::nullptr_t);
+
+  public:
+    /**
+        @name traits results
+        @{
+    */
+    using type = decltype(test<std::decay_t<TProcess>>(nullptr));
+    static const bool value = type::value;
+    //! @}
+  };
+
+  //! @file ContinuousProcess.hpp
+  //! value traits type
+  template <class TProcess, typename TReturn, typename TArg1, typename TArg2>
+  bool constexpr has_method_doContinuous_v =
+      has_method_doContinuous<TProcess, TReturn, TArg1, TArg2>::value;
+
+  /**
+     traits test for ContinuousProcess::getMaxStepLength method
+  */
+
+  template <class TProcess, typename TReturn, typename... TArgs>
+  struct has_method_getMaxStepLength
+      : public detail::has_method_signature<TReturn, TArgs...> {
+
+    //! method signature
+    using detail::has_method_signature<TReturn, TArgs...>::testSignature;
+
+    //! the default value
+    template <class T>
+    static std::false_type test(...);
+
+    //! templated option
+    template <class T>
+    static decltype(testSignature(&T::template getMaxStepLength<TArgs...>)) test(
+        std::nullptr_t);
+
+    //! non templated option
+    template <class T>
+    static decltype(testSignature(&T::getMaxStepLength)) test(std::nullptr_t);
+
+  public:
+    /**
+        @name traits results
+        @{
+    */
+    using type = decltype(test<std::decay_t<TProcess>>(nullptr));
+    static const bool value = type::value;
+    //! @}
+  };
+
+  //! value traits type
+  template <class TProcess, typename TReturn, typename... TArgs>
+  bool constexpr has_method_getMaxStepLength_v =
+      has_method_getMaxStepLength<TProcess, TReturn, TArgs...>::value;
+
+} // namespace corsika

corsika/detail/framework/process/DecayProcess.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 <corsika/framework/process/ProcessTraits.hpp>
+#include <corsika/framework/utility/HasMethodSignature.hpp>
+
+namespace corsika {
+
+  /**
+     traits test for DecayProcess::doDecay method
+  */
+
+  template <class TProcess, typename TReturn, typename... TArgs>
+  struct has_method_doDecay : public detail::has_method_signature<TReturn, TArgs...> {
+
+    //! type of process to be studied
+    typedef std::decay_t<TProcess> process_type;
+
+    ///! method signature
+    using detail::has_method_signature<TReturn, TArgs...>::testSignature;
+
+    //! the default value
+    template <class T>
+    static std::false_type test(...);
+
+    //! signature of templated method
+    template <class T>
+    static decltype(testSignature(&T::template doDecay<TArgs...>)) test(std::nullptr_t);
+
+    //! signature of non-templated method
+    template <class T>
+    static decltype(testSignature(&T::doDecay)) test(std::nullptr_t);
+
+  public:
+    /**
+        @name traits results
+        @{
+    */
+    using type = decltype(test<process_type>(nullptr));
+    static const bool value = type::value;
+    //! @}
+  };
+
+  //! @file DecayProcess.hpp
+  //! value traits type
+  template <class TProcess, typename TReturn, typename... TArgs>
+  bool constexpr has_method_doDecay_v =
+      has_method_doDecay<TProcess, TReturn, TArgs...>::value;
+
+  /**
+     traits test for DecayProcess::getLifetime method
+  */
+
+  template <class TProcess, typename TReturn, typename... TArgs>
+  struct has_method_getLifetime : public detail::has_method_signature<TReturn, TArgs...> {
+
+    //! the moethod signature
+    using detail::has_method_signature<TReturn, TArgs...>::testSignature;
+
+    //! the default value
+    template <class T>
+    static std::false_type test(...);
+
+    //! signature of templated method
+    template <class T>
+    static decltype(testSignature(&T::template getLifetime<TArgs...>)) test(
+        std::nullptr_t);
+
+    //! signature of non-templated method
+    template <class T>
+    static decltype(testSignature(&T::getLifetime)) test(std::nullptr_t);
+
+  public:
+    /**
+        @name traits results
+        @{
+    */
+    using type = decltype(test<std::decay_t<TProcess>>(nullptr));
+    static const bool value = type::value;
+    //! @}
+  };
+
+  //! value traits type
+  template <class TProcess, typename TReturn, typename... TArgs>
+  bool constexpr has_method_getLifetime_v =
+      has_method_getLifetime<TProcess, TReturn, TArgs...>::value;
+
+} // namespace corsika

corsika/detail/framework/process/InteractionCounter.inl

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/process/InteractionHistogram.hpp>
+
+namespace corsika {
+
+  template <class TCountedProcess>
+  inline InteractionCounter<TCountedProcess>::InteractionCounter(TCountedProcess& process)
+      : process_(process) {}
+
+  template <class TCountedProcess>
+  template <typename TSecondaryView>
+  inline void InteractionCounter<TCountedProcess>::doInteraction(
+      TSecondaryView& view, Code const projectileId, Code const targetId,
+      FourMomentum const& projectileP4, FourMomentum const& targetP4) {
+    size_t const massNumber = is_nucleus(targetId) ? get_nucleus_A(targetId) : 1;
+    auto const massTarget = massNumber * constants::nucleonMass;
+    histogram_.fill(projectileId, projectileP4.getTimeLikeComponent(), massTarget);
+    process_.doInteraction(view, projectileId, targetId, projectileP4, targetP4);
+  }
+
+  template <class TCountedProcess>
+  inline CrossSectionType InteractionCounter<TCountedProcess>::getCrossSection(
+      Code const projectileId, Code const targetId, FourMomentum const& projectileP4,
+      FourMomentum const& targetP4) const {
+    return process_.getCrossSection(projectileId, targetId, projectileP4, targetP4);
+  }
+
+  template <class TCountedProcess>
+  inline InteractionHistogram const& InteractionCounter<TCountedProcess>::getHistogram()
+      const {
+    return histogram_;
+  }
+
+} // namespace corsika

corsika/detail/framework/process/InteractionHistogram.hpp

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <boost/histogram.hpp>
+
+namespace corsika {
+
+  namespace detail {
+    inline auto hist_factory(unsigned int const bin_number, double const e_low,
+                             double const e_high) {
+      namespace bh = boost::histogram;
+      namespace bha = bh::axis;
+
+      auto h = bh::make_histogram(
+          bha::category<int, bh::use_default, bha::option::growth_t>{{2212, 2112},
+                                                                     "projectile PDG"},
+          bha::regular<float, bha::transform::log>{bin_number, (float)e_low,
+                                                   (float)e_high, "energy/eV"});
+      return h;
+    }
+  } // namespace detail
+} // namespace corsika

corsika/detail/framework/process/InteractionHistogram.inl

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <fstream>
+#include <string>
+
+#include <corsika/framework/process/InteractionHistogram.hpp>
+#include <corsika/detail/framework/process/InteractionHistogram.hpp> // for detail namespace
+
+namespace corsika {
+
+  inline InteractionHistogram::InteractionHistogram()
+      : inthist_cms_{detail::hist_factory(num_bins_cms, lower_edge_cms, upper_edge_cms)}
+      , inthist_lab_{detail::hist_factory(num_bins_lab, lower_edge_lab, upper_edge_lab)} {
+  }
+
+  inline void InteractionHistogram::fill(Code const projectile_id,
+                                         HEPEnergyType const lab_energy,
+                                         HEPEnergyType const mass_target) {
+    auto constexpr inv_eV = 1 / 1_eV;
+    auto const projectile_mass = get_mass(projectile_id);
+    auto const sqrtS = sqrt(projectile_mass * projectile_mass +
+                            mass_target * mass_target + 2 * lab_energy * mass_target);
+
+    CORSIKA_LOG_DEBUG("pM={}, tM={}, pid={}, Elab={}, sqrtS={}, pdg={} a={} z={}",
+                      projectile_mass / 1_GeV, mass_target / 1_GeV, projectile_id,
+                      lab_energy / 1_GeV, sqrtS / 1_GeV, get_PDG(projectile_id),
+                      get_nucleus_A(projectile_id), get_nucleus_Z(projectile_id));
+
+    inthist_cms_(static_cast<int>(get_PDG(projectile_id)), sqrtS * inv_eV);
+    inthist_lab_(static_cast<int>(get_PDG(projectile_id)), lab_energy * inv_eV);
+  }
+
+  inline InteractionHistogram& InteractionHistogram::operator+=(
+      InteractionHistogram const& other) {
+    inthist_lab_ += other.inthist_lab_;
+    inthist_cms_ += other.inthist_cms_;
+    return *this;
+  }
+
+  inline InteractionHistogram InteractionHistogram::operator+(
+      InteractionHistogram other) const {
+    other.inthist_lab_ += inthist_lab_;
+    other.inthist_cms_ += inthist_cms_;
+    return other;
+  }
+} // namespace corsika

corsika/detail/framework/process/InteractionLengthModifier.inl

+/*
+ * (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 <utility>
+
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+
+namespace corsika {
+
+  template <class TUnderlyingProcess>
+  inline InteractionLengthModifier<TUnderlyingProcess>::InteractionLengthModifier(
+      TUnderlyingProcess&& process,
+      std::function<InteractionLengthModifier::functor_signature> modifier)
+      : process_{std::move(process)}
+      , modifier_{std::move(modifier)} {}
+
+  template <class TUnderlyingProcess>
+  template <typename TSecondaryView>
+  inline void InteractionLengthModifier<TUnderlyingProcess>::doInteraction(
+      TSecondaryView& view) {
+    process_.doInteraction(view);
+  }
+
+  template <class TUnderlyingProcess>
+  template <typename TParticle>
+  inline GrammageType InteractionLengthModifier<TUnderlyingProcess>::getInteractionLength(
+      TParticle const& particle) {
+    GrammageType const original = process_.getInteractionLength(particle);
+    Code const pid = particle.getPID();
+    HEPEnergyType const energy = particle.getEnergy();
+
+    return modifier_(original, pid, energy);
+  }
+
+  template <class TUnderlyingProcess>
+  inline TUnderlyingProcess const&
+  InteractionLengthModifier<TUnderlyingProcess>::getProcess() const {
+    return process_;
+  }
+
+  template <class TUnderlyingProcess>
+  inline TUnderlyingProcess& InteractionLengthModifier<TUnderlyingProcess>::getProcess() {
+    return process_;
+  }
+
+} // namespace corsika

corsika/detail/framework/process/InteractionProcess.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 <corsika/framework/process/ProcessTraits.hpp>
+#include <corsika/framework/utility/HasMethodSignature.hpp>
+
+namespace corsika {
+
+  /**
+   *  traits test for InteractionProcess::doInteraction method.
+   */
+
+  template <class TProcess, typename TReturn, typename TTemplate, typename... TArgs>
+  struct has_method_doInteract : public detail::has_method_signature<TReturn, TArgs...> {
+
+    ///! method signature
+    using detail::has_method_signature<TReturn, TArgs...>::testSignature;
+
+    //! the default value
+    template <class T>
+    static std::false_type test(...);
+
+    //! signature of templated method
+    template <class T>
+    static decltype(testSignature(&T::template doInteraction<TTemplate>)) test(
+        std::nullptr_t);
+
+    //! signature of non-templated method
+    template <class T>
+    static decltype(testSignature(&T::doInteraction)) test(std::nullptr_t);
+
+  public:
+    /**
+     *  @name traits results
+     *  @{
+     */
+    using type = decltype(test<std::decay_t<TProcess>>(nullptr));
+    static const bool value = type::value;
+    //! @}
+  };
+
+  //! value traits type
+  template <class TProcess, typename TReturn, typename TTemplate, typename... TArgs>
+  bool constexpr has_method_doInteract_v =
+      has_method_doInteract<TProcess, TReturn, TTemplate, TArgs...>::value;
+
+  /**
+   *  traits test for TEMPLATED InteractionProcess::getCrossSection method (PROPOSAL).
+   */
+
+  template <class TProcess, typename TReturn, typename TTemplate, typename... TArgs>
+  struct has_method_getCrossSectionTemplate
+      : public detail::has_method_signature<TReturn, TArgs...> {
+
+    ///! method signature
+    using detail::has_method_signature<TReturn, TArgs...>::testSignature;
+
+    //! the default value
+    template <class T>
+    static std::false_type test(...);
+
+    //! templated parameter option
+    template <class T>
+    static decltype(testSignature(&T::template getCrossSection<TTemplate>)) test(
+        std::nullptr_t);
+
+    //! non templated parameter option
+    template <class T>
+    static decltype(testSignature(&T::getCrossSection)) test(std::nullptr_t);
+
+  public:
+    /**
+     *  @name traits results
+     *  @{
+     */
+    using type = decltype(test<std::decay_t<TProcess>>(nullptr));
+    static const bool value = type::value;
+    //! @}
+  };
+
+  //! value traits type shortcut
+  template <class TProcess, typename TReturn, typename TTemplate, typename... TArgs>
+  bool constexpr has_method_getCrossSectionTemplate_v =
+      has_method_getCrossSectionTemplate<TProcess, TReturn, TTemplate, TArgs...>::value;
+
+  /**
+   *  traits test for InteractionProcess::getCrossSection method.
+   */
+
+  template <class TProcess, typename TReturn, typename... TArgs>
+  struct has_method_getCrossSection
+      : public detail::has_method_signature<TReturn, TArgs...> {
+
+    ///! method signature
+    using detail::has_method_signature<TReturn, TArgs...>::testSignature;
+
+    //! the default value
+    template <class T>
+    static std::false_type test(...);
+
+    //! templated parameter option
+    template <class T>
+    static decltype(testSignature(&T::template getCrossSection<TArgs...>)) test(
+        std::nullptr_t);
+
+    //! non templated parameter option
+    template <class T>
+    static decltype(testSignature(&T::getCrossSection)) test(std::nullptr_t);
+
+  public:
+    /**
+     *  @name traits results
+     * @{
+     */
+    using type = decltype(test<std::decay_t<TProcess>>(nullptr));
+    static const bool value = type::value;
+    //! @}
+  };
+
+  //! value traits type shortcut
+  template <class TProcess, typename TReturn, typename... TArgs>
+  bool constexpr has_method_getCrossSection_v =
+      has_method_getCrossSection<TProcess, TReturn, TArgs...>::value;
+
+} // namespace corsika

corsika/detail/framework/process/ProcessSequence.inl

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/core/PhysicalUnits.hpp>
+
+#include <corsika/framework/process/BaseProcess.hpp>
+#include <corsika/framework/process/BoundaryCrossingProcess.hpp>
+#include <corsika/framework/process/ContinuousProcess.hpp>
+#include <corsika/framework/process/ContinuousProcessStepLength.hpp>
+#include <corsika/framework/process/ContinuousProcessIndex.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/process/CascadeEquationsProcess.hpp>
+#include <corsika/framework/core/Step.hpp>
+
+#include <cmath>
+#include <limits>
+#include <type_traits>
+
+namespace corsika {
+
+  template <typename TProcess1, typename TProcess2, int IndexStart, int IndexProcess1,
+            int IndexProcess2>
+  inline ProcessSequence<TProcess1, TProcess2, IndexStart, IndexProcess1,
+                         IndexProcess2>::ProcessSequence(TProcess1 in_A, TProcess2 in_B)
+      : A_(in_A)
+      , B_(in_B) {
+
+    // make sure only BaseProcess types TProcess1/2 are passed
+
+    static_assert(is_process_v<TProcess1>,
+                  "can only use process derived from BaseProcess in "
+                  "ProcessSequence, for Process 1");
+    static_assert(is_process_v<TProcess2>,
+                  "can only use process derived from BaseProcess in "
+                  "ProcessSequence, for Process 2");
+  }
+
+  template <typename TProcess1, typename TProcess2, int IndexStart, int IndexProcess1,
+            int IndexProcess2>
+  template <typename TParticle>
+  inline ProcessReturn ProcessSequence<
+      TProcess1, TProcess2, IndexStart, IndexProcess1,
+      IndexProcess2>::doBoundaryCrossing(TParticle& particle,
+                                         typename TParticle::node_type const& from,
+                                         typename TParticle::node_type const& to) {
+
+    ProcessReturn ret = ProcessReturn::Ok;
+
+    if constexpr (is_process_v<process1_type>) { // to protect from further compiler
+                                                 // errors if process1_type is invalid
+      if constexpr (is_boundary_process_v<process1_type> ||
+                    process1_type::is_process_sequence) {
+
+        // interface checking on TProcess1
+        static_assert(
+            has_method_doBoundaryCrossing_v<TProcess1, ProcessReturn, TParticle>,
+            "TDerived has no method with correct signature \"ProcessReturn "
+            "doBoundaryCrossing(TParticle&, VolumeNode const&, VolumeNode const&)\" "
+            "required for "
+            "BoundaryCrossingProcess<TDerived>. ");
+
+        ret |= A_.doBoundaryCrossing(particle, from, to);
+      }
+    }
+
+    if constexpr (is_process_v<process2_type>) { // to protect from further compiler
+                                                 // errors if process2_type is invalid
+      if constexpr (is_boundary_process_v<process2_type> ||
+                    process2_type::is_process_sequence) {
+
+        // interface checking on TProcess2
+        static_assert(
+            has_method_doBoundaryCrossing_v<TProcess2, ProcessReturn, TParticle>,
+            "TDerived has no method with correct signature \"ProcessReturn "
+            "doBoundaryCrossing(TParticle&, VolumeNode const&, VolumeNode const&)\" "
+            "required for "
+            "BoundaryCrossingProcess<TDerived>. ");
+
+        ret |= B_.doBoundaryCrossing(particle, from, to);
+      }
+    }
+
+    return ret;
+  }
+
+  template <typename TProcess1, typename TProcess2, int IndexStart, int IndexProcess1,
+            int IndexProcess2>
+  template <typename TParticle>
+  inline ProcessReturn
+  ProcessSequence<TProcess1, TProcess2, IndexStart, IndexProcess1, IndexProcess2>::
+      doContinuous(Step<TParticle>& step,
+                   [[maybe_unused]] ContinuousProcessIndex const limitId) {
+    ProcessReturn ret = ProcessReturn::Ok;
+
+    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) {
+
+        ret |= A_.doContinuous(step, limitId);
+      } else if constexpr (is_continuous_process_v<process1_type>) {
+
+        // interface checking on TProcess1
+        //~ static_assert(
+        //~ has_method_doContinuous_v<TProcess1, ProcessReturn, TParticle&, TTrack&> ||
+        //~ has_method_doContinuous_v<TProcess1, ProcessReturn, TParticle&,
+        //~ TTrack const&> ||
+        //~ has_method_doContinuous_v<TProcess1, ProcessReturn, TParticle const&,
+        //~ TTrack const&>,
+        //~ "TDerived has no method with correct signature \"ProcessReturn "
+        //~ "doContinuous(TParticle [const]&,TTrack [const]&,bool)\" required for "
+        //~ "ContinuousProcess<TDerived>. ");
+
+        ret |= A_.doContinuous(
+            step, limitId == ContinuousProcessIndex(
+                                 static_cast<void const*>(std::addressof(A_))));
+      }
+    }
+
+    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) {
+        ret |= B_.doContinuous(step, limitId);
+      } else if constexpr (is_continuous_process_v<process2_type>) {
+
+        // interface checking on TProcess2
+        //~ static_assert(
+        //~ has_method_doContinuous_v<TProcess2, ProcessReturn, TParticle&, TTrack&> ||
+        //~ has_method_doContinuous_v<TProcess2, ProcessReturn, TParticle&,
+        //~ TTrack const&> ||
+        //~ has_method_doContinuous_v<TProcess2, ProcessReturn, TParticle const&,
+        //~ TTrack const&>,
+        //~ "TDerived has no method with correct signature \"ProcessReturn "
+        //~ "doContinuous(TParticle [const]&,TTrack [const]&,bool)\" required for "
+        //~ "ContinuousProcess<TDerived>. ");
+
+        ret |= B_.doContinuous(
+            step, limitId == ContinuousProcessIndex(
+                                 static_cast<void const*>(std::addressof(B_))));
+      }
+    }
+
+    return ret;
+  }
+
+  template <typename TProcess1, typename TProcess2, int IndexStart, int IndexProcess1,
+            int IndexProcess2>
+  template <typename TSecondaries>
+  inline void
+  ProcessSequence<TProcess1, TProcess2, IndexStart, IndexProcess1,
+                  IndexProcess2>::doSecondaries([[maybe_unused]] TSecondaries& vS) {
+
+    if constexpr (is_process_v<process1_type>) { // to protect from further compiler
+                                                 // errors if process1_type is invalid
+      if constexpr (is_secondaries_process_v<process1_type> ||
+                    process1_type::is_process_sequence) {
+
+        // interface checking on TProcess1
+        static_assert(
+            has_method_doSecondaries_v<TProcess1, void, TSecondaries&> ||
+                has_method_doSecondaries_v<TProcess1, void, TSecondaries const&>,
+            "TDerived has no method with correct signature \"void "
+            "doSecondaries(TStackView [const]&)\" required for "
+            "SecondariesProcessProcess<TDerived>. ");
+
+        A_.doSecondaries(vS);
+      }
+    }
+
+    if constexpr (is_process_v<process2_type>) { // to protect from further compiler
+                                                 // errors if process2_type is invalid
+      if constexpr (is_secondaries_process_v<process2_type> ||
+                    process2_type::is_process_sequence) {
+
+        // interface checking on TProcess2
+        static_assert(
+            has_method_doSecondaries_v<TProcess2, void, TSecondaries&> ||
+                has_method_doSecondaries_v<TProcess2, void, TSecondaries const&>,
+            "TDerived has no method with correct signature \"void "
+            "doSecondaries(TStackView [const]&)\" required for "
+            "SecondariesProcessProcess<TDerived>. ");
+
+        B_.doSecondaries(vS);
+      }
+    }
+  }
+
+  template <typename TProcess1, typename TProcess2, int IndexStart, int IndexProcess1,
+            int IndexProcess2>
+  inline bool ProcessSequence<TProcess1, TProcess2, IndexStart, IndexProcess1,
+                              IndexProcess2>::checkStep() {
+    bool ret = false;
+    if constexpr (is_process_v<process1_type>) { // to protect from further compiler
+                                                 // errors if process1_type is invalid
+      if constexpr (is_stack_process_v<process1_type> ||
+                    (process1_type::is_process_sequence &&
+                     !process1_type::is_switch_process_sequence)) {
+        ret |= A_.checkStep();
+      }
+    }
+    if constexpr (is_process_v<process2_type>) { // to protect from further compiler
+                                                 // errors if process2_type is invalid
+      if constexpr (is_stack_process_v<process2_type> ||
+                    (process2_type::is_process_sequence &&
+                     !process2_type::is_switch_process_sequence)) {
+        ret |= B_.checkStep();
+      }
+    }
+    return ret;
+  }
+
+  template <typename TProcess1, typename TProcess2, int IndexStart, int IndexProcess1,
+            int IndexProcess2>
+  template <typename TStack>
+  inline void ProcessSequence<TProcess1, TProcess2, IndexStart, IndexProcess1,
+                              IndexProcess2>::doStack(TStack& stack) {
+    if constexpr (is_process_v<process1_type>) { // to protect from further compiler
+                                                 // errors if process1_type is invalid
+      if constexpr (is_stack_process_v<process1_type> ||
+                    (process1_type::is_process_sequence &&
+                     !process1_type::is_switch_process_sequence)) {
+
+        // interface checking on TProcess1
+        static_assert(has_method_doStack_v<TProcess1, void, TStack&> ||
+                          has_method_doStack_v<TProcess1, void, TStack const&>,
+                      "TDerived has no method with correct signature \"void "
+                      "doStack(TStack [const]&)\" required for "
+                      "StackProcess<TDerived>. ");
+
+        if (A_.checkStep()) { A_.doStack(stack); }
+      }
+    }
+    if constexpr (is_process_v<process2_type>) { // to protect from further compiler
+                                                 // errors if process2_type is invalid
+      if constexpr (is_stack_process_v<process2_type> ||
+                    (process2_type::is_process_sequence &&
+                     !process2_type::is_switch_process_sequence)) {
+
+        // interface checking on TProcess1
+        static_assert(has_method_doStack_v<TProcess2, void, TStack&> ||
+                          has_method_doStack_v<TProcess2, void, TStack const&>,
+                      "TDerived has no method with correct signature \"void "
+                      "doStack(TStack [const]&)\" required for "
+                      "StackProcess<TDerived>. ");
+
+        if (B_.checkStep()) { B_.doStack(stack); }
+      }
+    }
+  }
+
+  template <typename TProcess1, typename TProcess2, int IndexStart, int IndexProcess1,
+            int IndexProcess2>
+  template <typename TParticle, typename TTrack>
+  inline ContinuousProcessStepLength
+  ProcessSequence<TProcess1, TProcess2, IndexStart, IndexProcess1,
+                  IndexProcess2>::getMaxStepLength(TParticle&& particle,
+                                                   TTrack&& vTrack) {
+    // if no other process in the sequence implements it
+    ContinuousProcessStepLength max_length(std::numeric_limits<double>::infinity() *
+                                           meter);
+
+    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) {
+        ContinuousProcessStepLength const step = A_.getMaxStepLength(particle, vTrack);
+        max_length = std::min(max_length, step);
+      } else if constexpr (is_continuous_process_v<process1_type>) {
+
+        // interface checking on TProcess1
+        static_assert(has_method_getMaxStepLength_v<TProcess1, LengthType,
+                                                    TParticle const&, TTrack const&>,
+                      "TDerived has no method with correct signature \"LengthType "
+                      "getMaxStepLength(TParticle const&, TTrack const&)\" required for "
+                      "ContinuousProcess<TDerived>. ");
+
+        ContinuousProcessStepLength const step(
+            A_.getMaxStepLength(particle, vTrack),
+            ContinuousProcessIndex(static_cast<void const*>(std::addressof(A_))));
+        max_length = std::min(max_length, step);
+      }
+    }
+
+    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) {
+        ContinuousProcessStepLength const step = B_.getMaxStepLength(particle, vTrack);
+        max_length = std::min(max_length, step);
+      } else if constexpr (is_continuous_process_v<process2_type>) {
+
+        // interface checking on TProcess2
+        static_assert(has_method_getMaxStepLength_v<TProcess2, LengthType,
+                                                    TParticle const&, TTrack const&>,
+                      "TDerived has no method with correct signature \"LengthType "
+                      "getMaxStepLength(TParticle const&, TTrack const&)\" required for "
+                      "ContinuousProcess<TDerived>. ");
+
+        ContinuousProcessStepLength const step(
+            B_.getMaxStepLength(particle, vTrack),
+            ContinuousProcessIndex(static_cast<void const*>(std::addressof(B_))));
+        max_length = std::min(max_length, step);
+      }
+    }
+    return max_length;
+  }
+
+  template <typename TProcess1, typename TProcess2, int IndexStart, int IndexProcess1,
+            int IndexProcess2>
+  template <typename TParticle>
+  inline CrossSectionType
+  ProcessSequence<TProcess1, TProcess2, IndexStart, IndexProcess1, IndexProcess2>::
+      getCrossSection([[maybe_unused]] TParticle const& projectile,
+                      [[maybe_unused]] Code const targetId,
+                      [[maybe_unused]] FourMomentum const& targetP4) const {
+
+    CrossSectionType tot = CrossSectionType::zero();
+
+    if constexpr (is_process_v<process1_type>) { // to protect from further compiler
+                                                 // errors if process1_type is invalid
+      if constexpr (is_interaction_process_v<process1_type>) {
+
+        bool constexpr has_signature_cx1 =
+            has_method_getCrossSection_v<TProcess1,        // process object
+                                         CrossSectionType, // return type
+                                         Code, Code,       // parameters
+                                         FourMomentum const&, FourMomentum const&>;
+
+        if constexpr (has_signature_cx1) {
+          tot += A_.getCrossSection(projectile.getPID(), targetId,
+                                    {projectile.getEnergy(), projectile.getMomentum()},
+                                    targetP4);
+        } else { // for PROPOSAL
+          tot += A_.getCrossSection(projectile, projectile.getPID(),
+                                    {projectile.getEnergy(), projectile.getMomentum()});
+        }
+
+      } else if constexpr (process1_type::is_process_sequence) {
+        tot += A_.getCrossSection(projectile, targetId, targetP4);
+      }
+    }
+    if constexpr (is_process_v<process2_type>) { // to protect from further compiler
+                                                 // errors if process2_type is invalid
+      if constexpr (is_interaction_process_v<process2_type>) {
+
+        bool constexpr has_signature_cx1 =
+            has_method_getCrossSection_v<TProcess2,        // process object
+                                         CrossSectionType, // return type
+                                         Code, Code,       // parameters
+                                         FourMomentum const&, FourMomentum const&>;
+
+        if constexpr (has_signature_cx1) {
+          tot += B_.getCrossSection(projectile.getPID(), targetId,
+                                    {projectile.getEnergy(), projectile.getMomentum()},
+                                    targetP4);
+        } else { // for PROPOSAL
+          tot += B_.getCrossSection(projectile, projectile.getPID(),
+                                    {projectile.getEnergy(), projectile.getMomentum()});
+        }
+
+      } else if constexpr (process2_type::is_process_sequence) {
+        tot += B_.getCrossSection(projectile, targetId, targetP4);
+      }
+    }
+    return tot;
+  }
+
+  template <typename TProcess1, typename TProcess2, int IndexStart, int IndexProcess1,
+            int IndexProcess2>
+  template <typename TStack>
+  void ProcessSequence<TProcess1, TProcess2, IndexStart, IndexProcess1,
+                       IndexProcess2>::doCascadeEquations(TStack& stack) {
+
+    if constexpr (is_process_v<process1_type>) { // to protect from further compiler
+                                                 // errors if process1_type is invalid
+      if constexpr (process1_type::is_process_sequence &&
+                    !process1_type::is_switch_process_sequence) {
+
+        A_.doCascadeEquations(stack);
+
+      } else if constexpr (is_cascade_equations_process_v<process1_type>) {
+
+        // interface checking on TProcess1
+        static_assert(has_method_doCascadeEquations_v<TProcess1,
+                                                      void,     // return type
+                                                      TStack&>, // parameter
+                      "TDerived has no method with correct signature \"void "
+                      "doCascadeEquations(TStack&)\" required for "
+                      "CascadeEquationsProcess<TDerived>. ");
+
+        A_.doCascadeEquations(stack);
+      }
+    }
+
+    if constexpr (is_process_v<process2_type>) { // to protect from further compiler
+                                                 // errors if process2_type is invalid
+      if constexpr (process2_type::is_process_sequence &&
+                    !process2_type::is_switch_process_sequence) {
+
+        B_.doCascadeEquations(stack);
+
+      } else if constexpr (is_cascade_equations_process_v<process2_type>) {
+
+        // interface checking on TProcess2
+        static_assert(has_method_doCascadeEquations_v<TProcess2,
+                                                      void,     // return type
+                                                      TStack&>, // parameter
+                      "TDerived has no method with correct signature \"void "
+                      "doCascadeEquations(TStack&)\" required for "
+                      "CascadeEquationsProcess<TDerived>. ");
+
+        B_.doCascadeEquations(stack);
+      }
+    }
+  }
+
+  template <typename TProcess1, typename TProcess2, int IndexStart, int IndexProcess1,
+            int IndexProcess2>
+  void ProcessSequence<TProcess1, TProcess2, IndexStart, IndexProcess1,
+                       IndexProcess2>::initCascadeEquations() {
+
+    if constexpr (is_process_v<process1_type>) { // to protect from further compiler
+                                                 // errors if process1_type is invalid
+      if constexpr ((process1_type::is_process_sequence &&
+                     !process1_type::is_switch_process_sequence) ||
+                    is_cascade_equations_process_v<process1_type>) {
+        A_.initCascadeEquations();
+      }
+    }
+
+    if constexpr (is_process_v<process2_type>) { // to protect from further compiler
+                                                 // errors if process2_type is invalid
+      if constexpr ((process2_type::is_process_sequence &&
+                     !process2_type::is_switch_process_sequence) ||
+                    is_cascade_equations_process_v<process2_type>) {
+        B_.initCascadeEquations();
+      }
+    }
+  } // namespace corsika
+
+  template <typename TProcess1, typename TProcess2, int IndexStart, int IndexProcess1,
+            int IndexProcess2>
+  template <typename TSecondaryView, typename TRNG>
+  inline ProcessReturn
+  ProcessSequence<TProcess1, TProcess2, IndexStart, IndexProcess1, IndexProcess2>::
+      selectInteraction(TSecondaryView&& view, FourMomentum const& projectileP4,
+                        [[maybe_unused]] NuclearComposition const& composition,
+                        [[maybe_unused]] TRNG&& rng,
+                        [[maybe_unused]] CrossSectionType const cx_select,
+                        [[maybe_unused]] CrossSectionType cx_sum) {
+
+    // TODO: add check for cx_select > cx_tot
+
+    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) {
+        // if A is a process sequence --> check inside
+        ProcessReturn const ret =
+            A_.selectInteraction(view, projectileP4, composition, rng, cx_select, cx_sum);
+        // if A_ did succeed, stop routine. Not checking other static branch B_.
+        if (ret != ProcessReturn::Ok) { return ret; }
+      } else if constexpr (is_interaction_process_v<process1_type>) {
+
+        auto const& projectile = view.parent();
+        Code const projectileId = projectile.getPID();
+
+        // get cross section vector for all material components
+        // for selected process A
+
+        bool constexpr has_signature_cx1 =
+            has_method_getCrossSection_v<TProcess1,        // process object
+                                         CrossSectionType, // return type
+                                         Code, Code,       // parameters
+                                         FourMomentum const&, FourMomentum const&>;
+        bool constexpr has_signature_cx2 = // needed for PROPOSAL interface
+            has_method_getCrossSectionTemplate_v<
+                TProcess1,                   // process object
+                CrossSectionType,            // return type
+                decltype(projectile) const&, // template argument
+                decltype(projectile) const&, // parameters
+                Code, FourMomentum const&>;
+
+        static_assert((has_signature_cx1 || has_signature_cx2),
+                      "TProcess1 has no method with correct signature \"CrossSectionType "
+                      "getCrossSection(Code, Code, FourMomentum const&, FourMomentum "
+                      "const&)\" required by "
+                      "InteractionProcess<TProcess1>. ");
+
+        std::vector<CrossSectionType> weightedCrossSections;
+        if constexpr (has_signature_cx1) {
+          /*std::vector<CrossSectionType> const*/ weightedCrossSections =
+              composition.getWeighted([=](Code const targetId) -> CrossSectionType {
+                FourMomentum const targetP4(
+                    get_mass(targetId),
+                    MomentumVector(projectile.getMomentum().getCoordinateSystem(),
+                                   {0_GeV, 0_GeV, 0_GeV}));
+                return A_.getCrossSection(projectileId, targetId, projectileP4, targetP4);
+              });
+
+          cx_sum +=
+              std::accumulate(weightedCrossSections.cbegin(),
+                              weightedCrossSections.cend(), CrossSectionType::zero());
+
+        } else { // this is for PROPOSAL
+          cx_sum += A_.template getCrossSection(projectile, projectileId, projectileP4);
+        }
+
+        // check if we should execute THIS process and then EXIT
+        if (cx_select < cx_sum) {
+
+          if constexpr (has_signature_cx1) {
+            // now also sample targetId from weighted cross sections
+            Code const targetId = composition.sampleTarget(weightedCrossSections, rng);
+            FourMomentum const targetP4(
+                get_mass(targetId),
+                MomentumVector(projectile.getMomentum().getCoordinateSystem(),
+                               {0_GeV, 0_GeV, 0_GeV}));
+
+            // interface checking on TProcess1
+            static_assert(
+                has_method_doInteract_v<TProcess1,       // process object
+                                        void,            // return type
+                                        TSecondaryView,  // template argument
+                                        TSecondaryView&, // method parameters
+                                        Code, Code, FourMomentum const&,
+                                        FourMomentum const&>,
+                "TProcess1 has no method with correct signature \"void "
+                "doInteraction<TSecondaryView>(TSecondaryView&, "
+                "Code, Code, FourMomentum const&, FourMomentum const&)\" required for "
+                "InteractionProcess<TProcess1>. ");
+
+            A_.template doInteraction(view, projectileId, targetId, projectileP4,
+                                      targetP4);
+
+          } else { // this is for PROPOSAL
+            A_.template doInteraction(view, projectileId, projectileP4);
+          }
+
+          return ProcessReturn::Interacted;
+        }
+      }
+    } // end branch A
+
+    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) {
+        // if B_ is a process sequence --> check inside
+        return B_.selectInteraction(view, projectileP4, composition, rng, cx_select,
+                                    cx_sum);
+      } else if constexpr (is_interaction_process_v<process2_type>) {
+
+        auto const& projectile = view.parent();
+        Code const projectileId = projectile.getPID();
+
+        // get cross section vector for all material components, for selected process B
+        bool constexpr has_signature_cx1 =
+            has_method_getCrossSection_v<TProcess2,        // process object
+                                         CrossSectionType, // return type
+                                         Code, Code,       // parameters
+                                         FourMomentum const&, FourMomentum const&>;
+        bool constexpr has_signature_cx2 = // needed for PROPOSAL interface
+            has_method_getCrossSectionTemplate_v<
+                TProcess2,                    // process object
+                CrossSectionType,             // return type
+                decltype(*projectile) const&, // template argument
+                decltype(*projectile) const&, // parameters
+                Code,                         // parameters
+                FourMomentum const&>;
+        static_assert((has_signature_cx1 || has_signature_cx2),
+                      "TProcess2 has no method with correct signature \"CrossSectionType "
+                      "getCrossSection(Code, Code, FourMomentum const&, FourMomentum "
+                      "const&)\" required by "
+                      "InteractionProcess<TProcess1>. ");
+
+        std::vector<CrossSectionType> weightedCrossSections;
+        if constexpr (has_signature_cx1) {
+          /* std::vector<CrossSectionType> const*/ weightedCrossSections =
+              composition.getWeighted([=](Code const targetId) -> CrossSectionType {
+                FourMomentum const targetP4(
+                    get_mass(targetId),
+                    MomentumVector(projectile.getMomentum().getCoordinateSystem(),
+                                   {0_GeV, 0_GeV, 0_GeV}));
+                return B_.getCrossSection(projectileId, targetId, projectileP4, targetP4);
+              });
+
+          cx_sum +=
+              std::accumulate(weightedCrossSections.begin(), weightedCrossSections.end(),
+                              CrossSectionType::zero());
+        } else { // this is for PROPOSAL
+          cx_sum += B_.template getCrossSection(projectile, projectileId, projectileP4);
+        }
+
+        // check if we should execute THIS process and then EXIT
+        if (cx_select < cx_sum) {
+
+          if constexpr (has_signature_cx1) {
+
+            // now also sample targetId from weighted cross sections
+            Code const targetId = composition.sampleTarget(weightedCrossSections, rng);
+            FourMomentum const targetP4(
+                get_mass(targetId),
+                MomentumVector(projectile.getMomentum().getCoordinateSystem(),
+                               {0_GeV, 0_GeV, 0_GeV}));
+
+            // interface checking on TProcess2
+            static_assert(
+                has_method_doInteract_v<TProcess2,       // process object
+                                        void,            // return type
+                                        TSecondaryView,  // template argument
+                                        TSecondaryView&, // method parameters
+                                        Code, Code, FourMomentum const&,
+                                        FourMomentum const&>,
+                "TProcess1 has no method with correct signature \"void "
+                "doInteraction<TSecondaryView>(TSecondaryView&, "
+                "Code, Code, FourMomentum const&, FourMomentum const&)\" required for "
+                "InteractionProcess<TProcess2>. ");
+
+            B_.doInteraction(view, projectileId, targetId, projectileP4, targetP4);
+          } else { // this is for PROPOSAL
+            B_.doInteraction(view, projectileId, projectileP4);
+          }
+          return ProcessReturn::Interacted;
+        }
+      }
+    } // end branch B_
+    return ProcessReturn::Ok;
+  }
+
+  template <typename TProcess1, typename TProcess2, int IndexStart, int IndexProcess1,
+            int IndexProcess2>
+  template <typename TParticle>
+  inline InverseTimeType
+  ProcessSequence<TProcess1, TProcess2, IndexStart, IndexProcess1,
+                  IndexProcess2>::getInverseLifetime(TParticle&& particle) {
+
+    InverseTimeType tot = 0 / second; // default value
+
+    if constexpr (is_process_v<process1_type>) { // to protect from further compiler
+                                                 // errors if process1_type is invalid
+      if constexpr (is_decay_process_v<process1_type> ||
+                    process1_type::is_process_sequence) {
+        tot += A_.getInverseLifetime(particle);
+      }
+    }
+    if constexpr (is_process_v<process2_type>) { // to protect from further compiler
+                                                 // errors if process2_type is invalid
+      if constexpr (is_decay_process_v<process2_type> ||
+                    process2_type::is_process_sequence) {
+        tot += B_.getInverseLifetime(particle);
+      }
+    }
+    return tot;
+  }
+
+  template <typename TProcess1, typename TProcess2, int IndexStart, int IndexProcess1,
+            int IndexProcess2>
+  // select decay process
+  template <typename TSecondaryView>
+  inline ProcessReturn ProcessSequence<
+      TProcess1, TProcess2, IndexStart, IndexProcess1,
+      IndexProcess2>::selectDecay(TSecondaryView&& view,
+                                  [[maybe_unused]] InverseTimeType decay_inv_select,
+                                  [[maybe_unused]] InverseTimeType decay_inv_sum) {
+
+    // TODO: add check for decay_inv_select>decay_inv_tot
+
+    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) {
+        // if A_ is a process sequence --> check inside
+        ProcessReturn const ret = A_.selectDecay(view, decay_inv_select, decay_inv_sum);
+        // if A_ did succeed, stop routine here (not checking other static branch B_)
+        if (ret != ProcessReturn::Ok) { return ret; }
+      } else if constexpr (is_decay_process_v<process1_type>) {
+        // if this is not a ContinuousProcess --> evaluate probability
+        decay_inv_sum += A_.getInverseLifetime(view.parent());
+        // check if we should execute THIS process and then EXIT
+        if (decay_inv_select < decay_inv_sum) { // more pedagogical: rndm_select <
+                                                // decay_inv_sum / decay_inv_tot
+          // interface checking on TProcess1
+          static_assert(has_method_doDecay_v<TProcess1, void, TSecondaryView&>,
+                        "TDerived has no method with correct signature \"void "
+                        "doDecay(TSecondaryView&)\" required for "
+                        "DecayProcess<TDerived>. ");
+
+          A_.doDecay(view);
+          return ProcessReturn::Decayed;
+        }
+      } // end branch A_
+    }
+
+    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) {
+        // if B_ is a process sequence --> check inside
+        return B_.selectDecay(view, decay_inv_select, decay_inv_sum);
+      } else if constexpr (is_decay_process_v<process2_type>) {
+        // if this is not a ContinuousProcess --> evaluate probability
+        decay_inv_sum += B_.getInverseLifetime(view.parent());
+        // check if we should execute THIS process and then EXIT
+        if (decay_inv_select < decay_inv_sum) {
+
+          // interface checking on TProcess1
+          static_assert(has_method_doDecay_v<TProcess2, void, TSecondaryView&>,
+                        "TDerived has no method with correct signature \"void "
+                        "doDecay(TSecondaryView&)\" required for "
+                        "DecayProcess<TDerived>. ");
+
+          B_.doDecay(view);
+          return ProcessReturn::Decayed;
+        }
+      } // end branch B_
+    }
+
+    return ProcessReturn::Ok;
+  }
+
+  /**
+   * traits marker to identify objects containing any StackProcesses.
+   */
+  namespace detail {
+    // need helper alias to achieve this:
+    template <
+        typename TProcess1, typename TProcess2,
+        typename = typename std::enable_if_t<
+            contains_stack_process_v<TProcess1> || is_stack_process_v<TProcess1> ||
+                contains_stack_process_v<TProcess2> || is_stack_process_v<TProcess2>,
+            int>>
+    using enable_if_stack = ProcessSequence<TProcess1, TProcess2>;
+  } // namespace detail
+
+  template <typename TProcess1, typename TProcess2>
+  struct contains_stack_process<detail::enable_if_stack<TProcess1, TProcess2>>
+      : std::true_type {};
+
+} // namespace corsika

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

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

corsika/detail/framework/process/SwitchProcessSequence.inl

+/*
+ * (c) Copyright 2018 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/process/BaseProcess.hpp>
+#include <corsika/framework/process/ProcessTraits.hpp>
+#include <corsika/framework/process/BoundaryCrossingProcess.hpp>
+#include <corsika/framework/process/ContinuousProcess.hpp>
+#include <corsika/framework/process/ContinuousProcessStepLength.hpp>
+#include <corsika/framework/process/ContinuousProcessIndex.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 <cmath>
+#include <limits>
+#include <type_traits>
+
+namespace corsika {
+
+  template <typename TCondition, typename TSequence, typename USequence, int IndexStart,
+            int IndexProcess1, int IndexProcess2>
+  template <typename TParticle>
+  inline ProcessReturn SwitchProcessSequence<
+      TCondition, TSequence, USequence, IndexStart, IndexProcess1,
+      IndexProcess2>::doBoundaryCrossing(TParticle& particle,
+                                         typename TParticle::node_type const& from,
+                                         typename TParticle::node_type const& to) {
+    if (select_(particle)) {
+      if constexpr (is_boundary_process_v<process1_type> ||
+                    process1_type::is_process_sequence) {
+
+        // interface checking on TSequence
+        if constexpr (is_boundary_process_v<process1_type>) {
+
+          static_assert(
+              has_method_doBoundaryCrossing_v<TSequence, ProcessReturn, TParticle&>,
+              "TDerived has no method with correct signature \"ProcessReturn "
+              "doBoundaryCrossing(TParticle&, VolumeNode const&, VolumeNode const&)\" "
+              "required for "
+              "BoundaryCrossingProcess<TDerived>. ");
+        }
+
+        return A_.doBoundaryCrossing(particle, from, to);
+      }
+    } else {
+
+      if constexpr (is_boundary_process_v<process2_type> ||
+                    process2_type::is_process_sequence) {
+
+        // interface checking on USequence
+        if constexpr (is_boundary_process_v<process2_type>) {
+
+          static_assert(
+              has_method_doBoundaryCrossing_v<USequence, ProcessReturn, TParticle>,
+              "TDerived has no method with correct signature \"ProcessReturn "
+              "doBoundaryCrossing(TParticle&, VolumeNode const&, VolumeNode const&)\" "
+              "required for "
+              "BoundaryCrossingProcess<TDerived>. ");
+        }
+
+        return B_.doBoundaryCrossing(particle, from, to);
+      }
+    }
+    return ProcessReturn::Ok;
+  }
+
+  template <typename TCondition, typename TSequence, typename USequence, int IndexStart,
+            int IndexProcess1, int IndexProcess2>
+  template <typename TParticle>
+  inline ProcessReturn SwitchProcessSequence<TCondition, TSequence, USequence, IndexStart,
+                                             IndexProcess1, IndexProcess2>::
+      doContinuous(Step<TParticle>& step,
+                   [[maybe_unused]] ContinuousProcessIndex const idLimit) {
+    if (select_(step.getParticlePre())) {
+      if constexpr (process1_type::is_process_sequence) {
+        return A_.doContinuous(step, idLimit);
+      }
+      if constexpr (is_continuous_process_v<process1_type>) {
+
+        //        static_assert(
+        //            has_method_doContinuous_v<TSequence, ProcessReturn, TParticle&,
+        //            TTrack&> ||
+        //                has_method_doContinuous_v<TSequence, ProcessReturn, TParticle&,
+        //                                          TTrack const&> ||
+        //                has_method_doContinuous_v<TSequence, ProcessReturn, TParticle
+        //                const&,
+        //                                          TTrack const&>,
+        //            "TDerived has no method with correct signature \"ProcessReturn "
+        //            "doContinuous(TParticle[const]&,TTrack[const]&,bool)\" required for
+        //            " "ContinuousProcess<TDerived>. ");
+
+        return A_.doContinuous(
+            step, idLimit == ContinuousProcessIndex(
+                                 static_cast<void const*>(std::addressof(A_))));
+      }
+    } else {
+      if constexpr (process2_type::is_process_sequence) {
+        return B_.doContinuous(step, idLimit);
+      }
+      if constexpr (is_continuous_process_v<process2_type>) {
+
+        // interface checking on USequence
+        //        static_assert(
+        //            has_method_doContinuous_v<USequence, ProcessReturn, TParticle&,
+        //            TTrack&> ||
+        //                has_method_doContinuous_v<USequence, ProcessReturn, TParticle&,
+        //                                          TTrack const&> ||
+        //                has_method_doContinuous_v<USequence, ProcessReturn, TParticle
+        //                const&,
+        //                                          TTrack const&>,
+        //            "TDerived has no method with correct signature \"ProcessReturn "
+        //            "doContinuous(TParticle [const]&,TTrack[const]&,bool)\" required for
+        //            " "ContinuousProcess<TDerived>. ");
+
+        return B_.doContinuous(
+            step, idLimit == ContinuousProcessIndex(
+                                 static_cast<void const*>(std::addressof(B_))));
+      }
+    }
+    return ProcessReturn::Ok;
+  }
+
+  template <typename TCondition, typename TSequence, typename USequence, int IndexStart,
+            int IndexProcess1, int IndexProcess2>
+  template <typename TSecondaries>
+  inline void
+  SwitchProcessSequence<TCondition, TSequence, USequence, IndexStart, IndexProcess1,
+                        IndexProcess2>::doSecondaries(TSecondaries& vS) {
+    const auto& particle = vS.parent();
+    if (select_(particle)) {
+      if constexpr (is_secondaries_process_v<process1_type> ||
+                    process1_type::is_process_sequence) {
+
+        // interface checking on TSequence
+        static_assert(
+            has_method_doSecondaries_v<TSequence, void, TSecondaries&> ||
+                has_method_doSecondaries_v<TSequence, void, TSecondaries const&>,
+            "TDerived has no method with correct signature \"void "
+            "doSecondaries(TStackView [const]&)\" required for "
+            "SecondariesProcessProcess<TDerived>. ");
+
+        A_.doSecondaries(vS);
+      }
+    } else {
+      if constexpr (is_secondaries_process_v<process2_type> ||
+                    process2_type::is_process_sequence) {
+
+        // interface checking on USequence
+        static_assert(
+            has_method_doSecondaries_v<USequence, void, TSecondaries&> ||
+                has_method_doSecondaries_v<USequence, void, TSecondaries const&>,
+            "TDerived has no method with correct signature \"void "
+            "doSecondaries(TStackView [const]&)\" required for "
+            "SecondariesProcessProcess<TDerived>. ");
+
+        B_.doSecondaries(vS);
+      }
+    }
+  }
+
+  template <typename TCondition, typename TSequence, typename USequence, int IndexStart,
+            int IndexProcess1, int IndexProcess2>
+  template <typename TParticle, typename TTrack>
+  inline ContinuousProcessStepLength
+  SwitchProcessSequence<TCondition, TSequence, USequence, IndexStart, IndexProcess1,
+                        IndexProcess2>::getMaxStepLength(TParticle& particle,
+                                                         TTrack& vTrack) {
+    if (select_(particle)) {
+      if constexpr (process1_type::is_process_sequence) {
+        return A_.getMaxStepLength(particle, vTrack);
+      }
+      if constexpr (is_continuous_process_v<process1_type>) {
+
+        // interface checking on TSequence
+        static_assert(has_method_getMaxStepLength_v<TSequence, LengthType,
+                                                    TParticle const&, TTrack const&>,
+                      "TDerived has no method with correct signature \"LengthType "
+                      "getMaxStepLength(TParticle const&, TTrack const&)\" required for "
+                      "ContinuousProcess<TDerived>. ");
+
+        return ContinuousProcessStepLength(
+            A_.getMaxStepLength(particle, vTrack),
+            ContinuousProcessIndex(static_cast<void const*>(std::addressof(A_))));
+      }
+    } else {
+      if constexpr (process2_type::is_process_sequence) {
+        return B_.getMaxStepLength(particle, vTrack);
+      }
+      if constexpr (is_continuous_process_v<process2_type>) {
+
+        // interface checking on USequence
+        static_assert(has_method_getMaxStepLength_v<USequence, LengthType,
+                                                    TParticle const&, TTrack const&>,
+                      "TDerived has no method with correct signature \"LengthType "
+                      "getMaxStepLength(TParticle const&, TTrack const&)\" required for "
+                      "ContinuousProcess<TDerived>. ");
+
+        return ContinuousProcessStepLength(
+            B_.getMaxStepLength(particle, vTrack),
+            ContinuousProcessIndex(static_cast<void const*>(std::addressof(B_))));
+      }
+    }
+
+    // if no other process in the sequence implements it
+    return ContinuousProcessStepLength(std::numeric_limits<double>::infinity() * meter);
+  }
+
+  template <typename TCondition, typename TSequence, typename USequence, int IndexStart,
+            int IndexProcess1, int IndexProcess2>
+  template <typename TParticle>
+  CrossSectionType SwitchProcessSequence<TCondition, TSequence, USequence, IndexStart,
+                                         IndexProcess1, IndexProcess2>::
+      getCrossSection(TParticle const& projectile, [[maybe_unused]] Code const targetId,
+                      [[maybe_unused]] FourMomentum const& targetP4) const {
+
+    if (select_(projectile)) {
+      if constexpr (is_interaction_process_v<process1_type>) {
+        bool constexpr has_signature_cx1 =
+            has_method_getCrossSection_v<TSequence,        // process object
+                                         CrossSectionType, // return type
+                                         Code, Code,       // parameters
+                                         FourMomentum const&, FourMomentum const&>;
+        if constexpr (has_signature_cx1) {
+
+          return A_.getCrossSection(projectile.getPID(), targetId,
+                                    {projectile.getEnergy(), projectile.getMomentum()},
+                                    targetP4);
+        } else {
+          return A_.getCrossSection(projectile, projectile.getPID(),
+                                    {projectile.getEnergy(), projectile.getMomentum()});
+        }
+      } else if constexpr (process1_type::is_process_sequence) {
+        return A_.getCrossSection(projectile, targetId, targetP4);
+      }
+
+    } else {
+      if constexpr (is_interaction_process_v<process2_type>) {
+        bool constexpr has_signature_cx1 =
+            has_method_getCrossSection_v<USequence,        // process object
+                                         CrossSectionType, // return type
+                                         Code, Code,       // parameters
+                                         FourMomentum const&, FourMomentum const&>;
+        if constexpr (has_signature_cx1) {
+
+          return B_.getCrossSection(projectile.getPID(), targetId,
+                                    {projectile.getEnergy(), projectile.getMomentum()},
+                                    targetP4);
+        } else {
+          return B_.getCrossSection(projectile, targetId, targetP4);
+        }
+      } else if constexpr (process2_type::is_process_sequence) {
+        return B_.getCrossSection(projectile, targetId, targetP4);
+      }
+    }
+    return CrossSectionType::zero(); // default value
+  }
+
+  template <typename TCondition, typename TSequence, typename USequence, int IndexStart,
+            int IndexProcess1, int IndexProcess2>
+  template <typename TSecondaryView, typename TRNG>
+  inline ProcessReturn SwitchProcessSequence<
+      TCondition, TSequence, USequence, IndexStart, IndexProcess1,
+      IndexProcess2>::selectInteraction(TSecondaryView& view,
+                                        FourMomentum const& projectileP4,
+                                        NuclearComposition const& composition, TRNG& rng,
+                                        [[maybe_unused]] CrossSectionType const cx_select,
+                                        [[maybe_unused]] CrossSectionType cx_sum) {
+
+    if (select_(view.parent())) {
+      if constexpr (process1_type::is_process_sequence) {
+        // if A_ is a process sequence --> check inside
+        return A_.selectInteraction(view, projectileP4, composition, rng, cx_select,
+                                    cx_sum);
+      } else if constexpr (is_interaction_process_v<process1_type>) {
+
+        auto const& projectile = view.parent();
+        Code const projectileId = projectile.getPID();
+
+        // get cross section vector for all material components
+        // for selected process A
+        bool constexpr has_signature_cx1 =
+            has_method_getCrossSection_v<TSequence,        // process object
+                                         CrossSectionType, // return type
+                                         Code, Code,       // parameters
+                                         FourMomentum const&, FourMomentum const&>;
+        bool constexpr has_signature_cx2 = // needed for PROPOSAL interface
+            has_method_getCrossSectionTemplate_v<
+                TSequence,                   // process object
+                CrossSectionType,            // return type
+                decltype(projectile) const&, // template argument
+                decltype(projectile) const&, // parameters
+                Code, FourMomentum const&>;
+
+        static_assert((has_signature_cx1 || has_signature_cx2),
+                      "TSequence has no method with correct signature \"CrossSectionType "
+                      "getCrossSection(Code, Code, FourMomentum const&, FourMomentum "
+                      "const&)\" required by "
+                      "InteractionProcess<TSequence>. ");
+
+        std::vector<CrossSectionType> weightedCrossSections;
+        if constexpr (has_signature_cx1) {
+          /*std::vector<CrossSectionType> const*/ weightedCrossSections =
+              composition.getWeighted([=](Code const targetId) -> CrossSectionType {
+                FourMomentum const targetP4(
+                    get_mass(targetId),
+                    MomentumVector(projectile.getMomentum().getCoordinateSystem(),
+                                   {0_GeV, 0_GeV, 0_GeV}));
+                return A_.getCrossSection(projectileId, targetId, projectileP4, targetP4);
+              });
+
+          cx_sum +=
+              std::accumulate(weightedCrossSections.cbegin(),
+                              weightedCrossSections.cend(), CrossSectionType::zero());
+        } else { // this is for PROPOSAL
+          cx_sum += A_.template getCrossSection(projectile, projectileId, projectileP4);
+        }
+
+        if (cx_select < cx_sum) {
+
+          if constexpr (has_signature_cx1) {
+
+            // now also sample targetId from weighted cross sections
+            Code const targetId = composition.sampleTarget(weightedCrossSections, rng);
+            FourMomentum const targetP4(
+                get_mass(targetId),
+                MomentumVector(projectile.getMomentum().getCoordinateSystem(),
+                               {0_GeV, 0_GeV, 0_GeV}));
+
+            // interface checking on TProcess1
+            static_assert(
+                has_method_doInteract_v<TSequence,       // process object
+                                        void,            // return type
+                                        TSecondaryView,  // template argument
+                                        TSecondaryView&, // method parameters
+                                        Code, Code, FourMomentum const&,
+                                        FourMomentum const&>,
+                "TSequence has no method with correct signature \"void "
+                "doInteraction<TSecondaryView>(TSecondaryView&, "
+                "Code, Code, FourMomentum const&, FourMomentum const&)\" required for "
+                "InteractionProcess<TSequence>. ");
+
+            A_.template doInteraction(view, projectileId, targetId, projectileP4,
+                                      targetP4);
+          } else { // this is for PROPOSAL
+            A_.template doInteraction(view, projectileId, projectileP4);
+          }
+
+          return ProcessReturn::Interacted;
+        } // end collision branch A
+      }
+
+    } else { // selection: end branch A, start branch B
+
+      if constexpr (process2_type::is_process_sequence) {
+        // if B_ is a process sequence --> check inside
+        return B_.selectInteraction(view, projectileP4, composition, rng, cx_select,
+                                    cx_sum);
+      } else if constexpr (is_interaction_process_v<process2_type>) {
+
+        auto const& projectile = view.parent();
+        Code const projectileId = projectile.getPID();
+
+        // get cross section vector for all material components, for selected process B
+        bool constexpr has_signature_cx1 =
+            has_method_getCrossSection_v<USequence,        // process object
+                                         CrossSectionType, // return type
+                                         Code, Code,       // parameters
+                                         FourMomentum const&, FourMomentum const&>;
+        bool constexpr has_signature_cx2 = // needed for PROPOSAL interface
+            has_method_getCrossSectionTemplate_v<
+                USequence,                   // process object
+                CrossSectionType,            // return type
+                decltype(projectile) const&, // template argument
+                decltype(projectile) const&, // parameters
+                Code, FourMomentum const&>;
+
+        static_assert((has_signature_cx1 || has_signature_cx2),
+                      "USequence has no method with correct signature \"CrossSectionType "
+                      "getCrossSection(Code, Code, FourMomentum const&, FourMomentum "
+                      "const&)\" required by "
+                      "InteractionProcess<USequence>. ");
+
+        std::vector<CrossSectionType> weightedCrossSections;
+        if constexpr (has_signature_cx1) {
+          /* std::vector<CrossSectionType> const*/ weightedCrossSections =
+              composition.getWeighted([=](Code const targetId) -> CrossSectionType {
+                FourMomentum const targetP4(
+                    get_mass(targetId),
+                    MomentumVector(projectile.getMomentum().getCoordinateSystem(),
+                                   {0_GeV, 0_GeV, 0_GeV}));
+                return B_.getCrossSection(projectileId, targetId, projectileP4, targetP4);
+              });
+
+          cx_sum +=
+              std::accumulate(weightedCrossSections.begin(), weightedCrossSections.end(),
+                              CrossSectionType::zero());
+        } else { // this is for PROPOSAL
+          cx_sum += B_.template getCrossSection(projectile, projectileId, projectileP4);
+        }
+
+        // check if we should execute THIS process and then EXIT
+        if (cx_select < cx_sum) {
+
+          if constexpr (has_signature_cx1) {
+
+            // now also sample targetId from weighted cross sections
+            Code const targetId = composition.sampleTarget(weightedCrossSections, rng);
+            FourMomentum const targetP4(
+                get_mass(targetId),
+                MomentumVector(projectile.getMomentum().getCoordinateSystem(),
+                               {0_GeV, 0_GeV, 0_GeV}));
+
+            // interface checking on TProcess2
+            static_assert(
+                has_method_doInteract_v<USequence,       // process object
+                                        void,            // return type
+                                        TSecondaryView,  // template argument
+                                        TSecondaryView&, // method parameters
+                                        Code, Code, FourMomentum const&,
+                                        FourMomentum const&>,
+                "USequence has no method with correct signature \"void "
+                "doInteraction<TSecondaryView>(TSecondaryView&, "
+                "Code, Code, FourMomentum const&, FourMomentum const&)\" required for "
+                "InteractionProcess<USequence>. ");
+
+            B_.doInteraction(view, projectileId, targetId, projectileP4, targetP4);
+          } else { // this is for PROPOSAL
+            B_.doInteraction(view, projectileId, projectileP4);
+          }
+
+          return ProcessReturn::Interacted;
+        } // end collision in branch B
+      }
+    } // end branch B_
+
+    return ProcessReturn::Ok;
+  } // namespace corsika
+
+  template <typename TCondition, typename TSequence, typename USequence, int IndexStart,
+            int IndexProcess1, int IndexProcess2>
+  template <typename TParticle>
+  inline InverseTimeType
+  SwitchProcessSequence<TCondition, TSequence, USequence, IndexStart, IndexProcess1,
+                        IndexProcess2>::getInverseLifetime(TParticle&& particle) {
+
+    if (select_(particle)) {
+      if constexpr (is_decay_process_v<process1_type> ||
+                    process1_type::is_process_sequence) {
+        return A_.getInverseLifetime(particle);
+      }
+
+    } else {
+
+      if constexpr (is_decay_process_v<process2_type> ||
+                    process2_type::is_process_sequence) {
+        return B_.getInverseLifetime(particle);
+      }
+    }
+    return 0 / second; // default value
+  }
+
+  template <typename TCondition, typename TSequence, typename USequence, int IndexStart,
+            int IndexProcess1, int IndexProcess2>
+  // select decay process
+  template <typename TSecondaryView>
+  inline ProcessReturn SwitchProcessSequence<
+      TCondition, TSequence, USequence, IndexStart, IndexProcess1,
+      IndexProcess2>::selectDecay(TSecondaryView& view,
+                                  [[maybe_unused]] InverseTimeType decay_inv_select,
+                                  [[maybe_unused]] InverseTimeType decay_inv_sum) {
+    if (select_(view.parent())) {
+      if constexpr (process1_type::is_process_sequence) {
+        // if A_ is a process sequence --> check inside
+        ProcessReturn const ret = A_.selectDecay(view, decay_inv_select, decay_inv_sum);
+        // if A_ did succeed, stop routine here (not checking other static branch B_)
+        if (ret != ProcessReturn::Ok) { return ret; }
+      } else if constexpr (is_decay_process_v<process1_type>) {
+        // if this is not a ContinuousProcess --> evaluate probability
+        decay_inv_sum += A_.getInverseLifetime(view.parent());
+        // check if we should execute THIS process and then EXIT
+        if (decay_inv_select < decay_inv_sum) {
+          // more pedagogical: rndm_select < decay_inv_sum / decay_inv_tot
+
+          // interface checking on TSequence
+          static_assert(has_method_doDecay_v<TSequence, void, TSecondaryView&>,
+                        "TDerived has no method with correct signature \"void "
+                        "doDecay(TSecondaryView&)\" required for "
+                        "DecayProcess<TDerived>. ");
+
+          A_.doDecay(view);
+          return ProcessReturn::Decayed;
+        }
+      } // end branch A_
+
+    } else {
+
+      if constexpr (process2_type::is_process_sequence) {
+        // if B_ is a process sequence --> check inside
+        return B_.selectDecay(view, decay_inv_select, decay_inv_sum);
+      } else if constexpr (is_decay_process_v<process2_type>) {
+        // if this is not a ContinuousProcess --> evaluate probability
+        decay_inv_sum += B_.getInverseLifetime(view.parent());
+        // check if we should execute THIS process and then EXIT
+        if (decay_inv_select < decay_inv_sum) {
+
+          // interface checking on TSequence
+          static_assert(has_method_doDecay_v<USequence, void, TSecondaryView&>,
+                        "TDerived has no method with correct signature \"void "
+                        "doDecay(TSecondaryView&)\" required for "
+                        "DecayProcess<TDerived>. ");
+
+          B_.doDecay(view);
+          return ProcessReturn::Decayed;
+        }
+      } // end branch B_
+    }
+    return ProcessReturn::Ok;
+  }
+
+} // namespace corsika

corsika/detail/framework/random/RNGManager.inl

+/*
+ * (c) Copyright 2018 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 <iterator>
+#include <random>
+#include <sstream>
+#include <tuple>
+#include <utility>
+
+namespace corsika {
+
+  template <typename CBPRNG>
+  inline void RNGManager<CBPRNG>::registerRandomStream(string_type const& pStreamName) {
+
+    auto const& it = rngs_.find(pStreamName);
+
+    if (it == rngs_.end()) // key not in container, so create one and initialize the value
+      rngs_.emplace(std::piecewise_construct, std::forward_as_tuple(pStreamName.c_str()),
+                    std::forward_as_tuple(seed_, uint32_t(rngs_.size())));
+  }
+
+  template <typename CBPRNG>
+  inline typename RNGManager<CBPRNG>::prng_type& RNGManager<CBPRNG>::getRandomStream(
+      string_type const& pStreamName) {
+
+    if (isRegistered(pStreamName)) {
+      return rngs_.at(pStreamName);
+    } else { // this stream name is not in the map
+      throw std::runtime_error("'" + pStreamName + "' is not a registered stream.");
+    }
+  }
+
+  template <typename CBPRNG>
+  inline bool RNGManager<CBPRNG>::isRegistered(string_type const& pStreamName) const {
+    return rngs_.count(pStreamName) > 0;
+  }
+
+  template <typename CBPRNG>
+  inline std::stringstream RNGManager<CBPRNG>::dumpState() const {
+    std::stringstream buffer;
+    for (auto const& [streamName, rng] : rngs_) {
+      buffer << '"' << streamName << "\" = \"" << rng << '"' << std::endl;
+    }
+
+    return buffer;
+  }
+
+} // namespace corsika

corsika/detail/framework/random/random_iterator/RandomIterator.h

+/*----------------------------------------------------------------------------
+*
+* Copyright (C) 2021 - 2024 Antonio Augusto Alves Junior
+*
+* This file is part of RandomIterator.
+
+* Redistribution and use in source and binary forms, with or without modification,
+* are permitted provided that the following conditions are met:
+*
+* 1. Redistributions of source code must retain the above copyright notice,
+* this list of conditions and the following disclaimer.
+*
+* 2. Redistributions in binary form must reproduce the above copyright notice,
+* this list of conditions and the following disclaimer in the documentation and/or
+* other materials provided with the distribution.
+*
+* 3. Neither the name of the copyright holder nor the names of its contributors
+* may be used to endorse or promote products derived from this software without
+* specific prior written permission.
+*
+* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS “AS IS” AND
+* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+* WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
+* IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
+* INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
+* BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+* LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+* HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
+* STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+* ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
+* OF THE POSSIBILITY OF SUCH DAMAGE.
+*
+* ----------------------------------------------------------------------------*/
+
+/*
+ * RandomIterator.h
+ *
+ *  Created on: 21/02/2021
+ *      Author: Antonio Augusto Alves Junior
+ */
+
+#pragma once
+
+//  This is the only RandomIterator header that is guaranteed to
+//  change with every RandomIterator release.
+//
+//  RandomIterator_VERSION % 100 is the sub-minor version
+//  RandomIterator_VERSION / 100 % 1000 is the minor version
+//  RandomIterator_VERSION / 100000 is the major version
+//
+//  Because this header does not #include <RandomIterator/detail/Config.h>,
+//  it is the only RandomIterator header that does not cause
+//  RandomIterator_HOST_SYSTEM and RandomIterator_DEVICE_SYSTEM to be defined.
+
+/*! \def RandomIterator_VERSION
+ *  \brief The preprocessor macro \p RandomIterator_VERSION encodes the version
+ *         number of the RandomIterator.
+ *
+ *         <tt>RandomIterator_VERSION % 100</tt> is the patch version.
+ *         <tt>RandomIterator_VERSION / 100 % 1000</tt> is the feature version.
+ *         <tt>RandomIterator_VERSION / 100000</tt> is the major version.
+ */
+#define RandomIterator_VERSION 100001
+
+
+/*! \def RandomIterator_MAJOR_VERSION
+ *  \brief The preprocessor macro \p RandomIterator_MAJOR_VERSION encodes the
+ *         major version number of RandomIterator.
+ */
+#define RandomIterator_MAJOR_VERSION     (RandomIterator_VERSION / 100000)
+
+/*! \def RandomIterator_MINOR_VERSION
+ *  \brief The preprocessor macro \p RandomIterator_MINOR_VERSION encodes the
+ *         minor version number of RandomIterator.
+ */
+#define RandomIterator_MINOR_VERSION     (RandomIterator_VERSION / 100 % 1000)
+
+/*! \def RandomIterator_PATCH_NUMBER
+ *  \brief The preprocessor macro \p RandomIterator_PATCH_NUMBER encodes the
+ *         patch number of the RandomIterator library.
+ */
+#define RandomIterator_PATCH_NUMBER  0
+
+
+
+// Declare these namespaces here for the purpose of Doxygenating them
+
+/*!
+ *  \brief \p random_iterator is the top-level namespace which contains all RandomIterator
+ *         functions and types.
+ */
+namespace random_iterator{ }
+
+

corsika/detail/framework/random/random_iterator/Stream.hpp

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+/*
+ * Stream.hpp
+ *
+ *  Created on: 23/02/2021
+ *      Author: Antonio Augusto Alves Junior
+ */
+
+#pragma once
+
+#include <stdint.h>
+
+#include "detail/tbb/iterators.h"
+#include "detail/Engine.hpp"
+#include "detail/functors/EngineCaller.hpp"
+
+namespace random_iterator {
+
+  /// philox engine
+  typedef detail::philox philox;
+  /// threefry  engine
+  typedef detail::threefry threefry;
+  /// squares3_128  engine
+  typedef detail::squares3_128 squares3_128;
+  /// squares4_128  engine
+  typedef detail::squares4_128 squares4_128;
+  /// squares3_64  engine
+  typedef detail::squares3_64 squares3_64;
+  /// squares4_64  engine
+  typedef detail::squares4_64 squares4_64;
+
+#if RANDOM_ITERATOR_R123_USE_AES_NI
+  /// ars  engine
+  typedef detail::ars ars;
+#endif
+
+  /**
+   * \class Stream
+   *
+   * \brief representation for streams of pseudorandom numbers
+   *
+   * Objects of this class represent collections of up to \f$ {2}^{32} \f$  streams of
+   * random numbers for each seed. Each of such streams has a length of \f$ {2}^{64} \f$.
+   * If the chosen distribution generates data with 64bits, like ``double``, then each
+   * stream can handle up to 128 ExaByte of data. Stream objects are thread-safe and
+   * lazy-evaluated. Output is produced only when requested. It is user responsibility to
+   * retrieve and store, if needed, the generated output.
+   *
+   * \tparam DistibutionType STL compliant distribution type. For example:
+   * ```std::uniform_real_distribution```, ```std::exponential_distribution```.
+   *
+   * \tparam EngineType RNG type.
+   *
+   * There are seven options:
+   *
+   *  * random_iterator::philox
+   *  * random_iterator::threefry
+   *  * random_iterator::ars
+   *  * random_iterator::squares3_128
+   *  * random_iterator::squares4_128
+   *  * random_iterator::squares3_64
+   *  * random_iterator::squares4_64
+   *
+   *  Usage:
+   *
+   * 	1) Direct class instantiation:
+   *
+   *  \code
+   *  \\this code will print the 10 doubles randomly distributed in the range [0.0, 1.0].
+   *
+   *  \\instantiate real uniform distribution from C++ standard libray
+   *  std::uniform_real_distribution<double> udist(0.0, 1.0);
+   *
+   *  \\
+   *  Stream<std::uniform_real_distribution<double>, random_iterator::philox>
+   * uniform_stream(udist, 0x1a2b3c4d, 1 );
+   *
+   *  for( size_t i=0; i<10; ++i )
+   *   std::cout << uniform_stream[i] << std::endl;
+   *  \endcode
+   *
+   * 	2) Using ```random_iterator::make_stream(...)```:
+   *
+   * 	\code
+   *  \\this code will print the 10 doubles randomly distributed in the range [0.0, 1.0].
+   *
+   *  \\instantiate real uniform distribution from C++ standard libray
+   *  std::uniform_real_distribution<double> udist(0.0, 1.0);
+   *
+   *  \\instantiate the generator
+   *  random_iterator::philox   rng(0x1a2b3c4d);
+   *
+   *  \\create stream
+   *  auto uniform_stream = random_iterator::make_stream(udist, rng, 1 );
+   *
+   *  for( size_t i=0; i<10; ++i )
+   *   std::cout << uniform_stream[i] << std::endl;
+   * 	\endcode
+   */
+  template <typename DistibutionType, typename EngineType>
+  class Stream {
+    typedef detail::EngineCaller<DistibutionType, EngineType> caller_type; //!
+    typedef random_iterator_tbb::counting_iterator<typename EngineType::advance_type>
+        counting_iterator; //!
+
+  public:
+    typedef EngineType engine_type;                          //! Engine type
+    typedef typename engine_type::state_type state_type;     //! Type of RNG state
+    typedef typename engine_type::seed_type seed_type;       //! Type of RNG seed
+    typedef typename engine_type::advance_type advance_type; //! Type of RNG displacement
+    typedef typename engine_type::init_type init_type;       //! Type of RNG initializer
+
+    typedef DistibutionType distribution_type; //! type of the STL compliant distribution
+    typedef typename distribution_type::result_type result_type; //! type of the result
+
+    typedef random_iterator_tbb::transform_iterator<caller_type, counting_iterator>
+        iterator; //! type of stream iterator
+
+    Stream() = delete;
+
+    /**
+     * Constructor.
+     *
+     * @param distribution Distribution with stl compliant interface.
+     * @param seed Seed for pseudorandom number generation.
+     * @param stream Pseudorandom number stream.
+     */
+    Stream(distribution_type const& distribution, seed_type seed, uint32_t stream = 0)
+        : distribution_(distribution)
+        , engine_(seed, stream)
+        , seed_(seed)
+        , stream_(stream) {}
+
+    /**
+     * Copy constructor
+     *
+     * @param other
+     */
+    Stream(Stream<DistibutionType, EngineType> const& other)
+        : distribution_(other.getDistribution())
+        , engine_(other.getEngine())
+        , seed_(other.getSeed())
+        , stream_(other.getStream()) {}
+
+    /**
+     * Returns an iterator to the first element of the stream.
+     * @return  iterator to the first element of the stream.
+     */
+    inline iterator begin() const {
+
+      counting_iterator first(0);
+      caller_type caller(distribution_, seed_, stream_);
+
+      return iterator(first, caller);
+    }
+
+    /**
+     * Returns an iterator to the past last element of the stream.
+     * @return iterator to the past last element of the stream.
+     */
+    inline iterator end() const {
+
+      counting_iterator last(std::numeric_limits<advance_type>::max());
+      caller_type caller(distribution_, seed_, stream_);
+
+      return iterator(last, caller);
+    }
+
+    /**
+     * Returns the element at specified location n. No bounds checking is performed.
+     *
+     * @param n position in the sequence.
+     * @return element at location n
+     */
+    inline result_type operator[](advance_type n) const {
+
+      caller_type caller(distribution_, seed_, stream_);
+
+      return caller(n);
+    }
+
+    /**
+     *  At each call, this function will produce a pseudorandom number and advance the
+     * engine state.
+     *
+     * @return pseudorandom number distributed according with `DistributionType`
+     */
+    inline result_type operator()(void) { return distribution_(engine_); }
+
+    /**
+     * Get the distribution
+     *
+     * @return Copy of the object's distribution.
+     */
+    inline distribution_type getDistribution() const { return distribution_; }
+
+    /**
+     * Set the distribution
+     * @param dist
+     *
+     */
+    inline void setDistribution(distribution_type dist) { distribution_ = dist; }
+
+    /**
+     * Get the associated engine
+     *
+     * @return Copy of the object's engine
+     */
+    inline engine_type getEngine() const { return engine_; }
+
+    /**
+     * Set the object's engine
+     *
+     * @param engine
+     */
+    inline void setEngine(engine_type engine) { engine_ = engine; }
+
+    /**
+     * Get the seed
+     *
+     * @return Copy of the object's seed
+     */
+    inline const seed_type& getSeed() const { return seed_; }
+
+    /**
+     * Set the seed
+     *
+     * @param seed
+     */
+    inline void setSeed(const seed_type& seed) { seed_ = seed; }
+
+    /**
+     * Get the object's stream number
+     *
+     * @return stream number
+     */
+    inline uint32_t getStream() const { return stream_; }
+
+    /**
+     * Set the object's stream number
+     *
+     * @param stream
+     */
+    inline void setStream(uint32_t stream) { stream_ = stream; }
+
+    friend inline std::ostream& operator<<(
+        std::ostream& os, const Stream<DistibutionType, EngineType>& be) {
+      return os << "engine: " << be.getEngine() << " stream: " << be.getStream()
+                << " distribution: " << be.getDistribution();
+    }
+
+  private:
+    distribution_type distribution_;
+    engine_type engine_;
+    seed_type seed_;
+    uint32_t stream_;
+  };
+
+  /**
+   * Stream
+   *
+   * @param dist STL compliant distribution instance.
+   * @param eng Engine instance
+   * @param stream Stream number, in the range [0, 2^{32} - 1]
+   * @return Stream object.
+   */
+  template <typename DistibutionType, typename EngineType>
+  Stream<DistibutionType, EngineType> make_stream(DistibutionType const& dist,
+                                                  EngineType eng, uint32_t stream) {
+
+    return Stream<DistibutionType, EngineType>(dist, eng.getSeed(), stream);
+  }
+
+} // namespace random_iterator

corsika/detail/framework/random/random_iterator/detail/Engine.hpp

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+/*
+ * Engine.hpp
+ *
+ *  Created on: 22 de fev. de 2021
+ *      Author: Antonio Augusto Alves Junior
+ */
+
+#pragma once
+
+#include <limits>
+
+#include "EngineTraits.hpp"
+#include "SplitMix.hpp"
+#include "Squares3_128.hpp"
+#include "Squares4_128.hpp"
+#include "Squares3_64.hpp"
+#include "Squares4_64.hpp"
+
+namespace random_iterator {
+
+  namespace detail {
+
+    template <typename Engine123>
+    class Engine {
+
+      typedef unsigned trigger_type;
+
+    public:
+      typedef Engine123 engine_type;
+      typedef typename detail::random_traits<engine_type>::state_type state_type;
+      typedef typename detail::random_traits<engine_type>::seed_type seed_type;
+      typedef typename detail::random_traits<engine_type>::advance_type advance_type;
+      typedef typename detail::random_traits<engine_type>::init_type init_type;
+      typedef typename detail::random_traits<engine_type>::result_type result_type;
+
+      static const unsigned arity = detail::random_traits<engine_type>::arity;
+
+      Engine() = delete;
+
+      Engine(result_type seed)
+          : engine_(engine_type{})
+          , cache_(state_type{})
+          , state_(state_type{})
+          , seed_(seed_type{})
+          , trigger_(arity) {
+        init_type temp{{}};
+
+        for (unsigned i = 0; i < temp.size(); ++i) {
+          temp[i] = detail::splitmix<result_type>(seed);
+        }
+
+        seed_ = temp;
+      }
+
+      Engine(result_type seed, uint32_t stream)
+          : engine_(engine_type{})
+          , cache_(state_type{})
+          , state_(state_type{})
+          , seed_(seed_type{})
+          , trigger_(arity) {
+        init_type temp{{}};
+
+        for (unsigned i = 0; i < temp.size(); ++i) {
+          temp[i] = detail::splitmix<result_type>(seed);
+        }
+
+        state_[arity - 1] = stream;
+        state_[arity - 1] = state_[arity - 1] << 32;
+
+        seed_ = temp;
+      }
+
+      Engine(init_type seed)
+          : engine_(engine_type{})
+          , cache_(state_type{})
+          , state_(state_type{})
+          , seed_(seed)
+          , trigger_(arity) {}
+
+      Engine(init_type seed, uint32_t stream)
+          : engine_(engine_type{})
+          , cache_(state_type{})
+          , state_(state_type{})
+          , seed_(seed)
+          , trigger_(arity) {
+        state_[arity - 1] = stream;
+        state_[arity - 1] << 32;
+      }
+
+      Engine(Engine<Engine123> const& other)
+          : engine_(engine_type{})
+          , cache_(other.getCache())
+          , state_(other.getState())
+          , seed_(other.getSeed())
+          , trigger_(other.getTrigger()) {}
+
+      inline Engine<Engine123>& operator=(Engine<Engine123> const& other) {
+
+        if (this == &other) return *this;
+
+        engine_ = engine_type{};
+        cache_ = other.getCache();
+        state_ = other.getState();
+        seed_ = other.getSeed();
+        trigger_ = other.getTrigger();
+        return *this;
+      }
+
+      inline result_type operator()(void) {
+        result_type result = 0;
+
+        if (trigger_ == arity) {
+
+          trigger_ = 0;
+          cache_ = engine_(state_.incr(), seed_);
+          result = cache_[trigger_];
+          ++trigger_;
+        } else {
+          result = cache_[trigger_];
+          ++trigger_;
+        }
+
+        return result;
+      }
+
+      inline void discard(advance_type n) {
+
+        state_.incr(n);
+        trigger_ = arity;
+      }
+
+      inline void reset(void) {
+        trigger_ = arity;
+        state_ = state_type{};
+        cache_ = state_type{};
+      }
+
+      inline const seed_type& getSeed() const { return seed_; }
+
+      inline void setSeed(seed_type seed) { seed_ = seed; }
+
+      inline void setSeed(result_type seed) {
+
+        init_type temp{{}};
+        for (unsigned i = 0; i < temp.size(); ++i) {
+          temp[i] = detail::splitmix<result_type>(seed);
+        }
+        seed_ = temp;
+      }
+
+      inline const state_type& getState() const { return state_; }
+
+      inline void setState(const state_type& state) { state_ = state; }
+
+      inline trigger_type getTrigger() const { return trigger_; }
+
+      inline void setTrigger(trigger_type trigger) { trigger_ = trigger; }
+      static constexpr result_type min() { return 0; }
+
+      static constexpr result_type max() {
+        return std::numeric_limits<result_type>::max();
+      }
+
+      friend inline std::ostream& operator<<(std::ostream& os,
+                                             const Engine<Engine123>& be) {
+        return os << "state: " << be.getState() << " seed: " << be.getSeed()
+                  << " trigger: " << be.getTrigger();
+      }
+
+    private:
+      inline state_type getCache() const { return cache_; }
+
+      engine_type engine_;
+      state_type cache_;
+      state_type state_;
+      seed_type seed_;
+      trigger_type trigger_;
+    };
+
+    typedef Engine<random_iterator_r123::Philox2x64> philox;
+    typedef Engine<random_iterator_r123::Threefry2x64> threefry;
+    typedef detail::Squares3_64 squares3_64;
+    typedef detail::Squares4_64 squares4_64;
+    typedef detail::Squares3_128 squares3_128;
+    typedef detail::Squares4_128 squares4_128;
+
+#if RANDOM_ITERATOR_R123_USE_AES_NI
+    typedef Engine<random_iterator_r123::ARS2x64> ars;
+#endif
+
+  } // namespace detail
+
+} // namespace random_iterator

corsika/detail/framework/random/random_iterator/detail/EngineTraits.hpp

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+/*
+ * EngineTraits.hpp
+ *
+ *  Created on: 23 de fev. de 2021
+ *      Author: Antonio Augusto Alves Junior
+ */
+
+#pragma once
+
+#include <stdint.h>
+#include "Random123/array.h"
+#include "Random123/philox.h"
+#include "Random123/threefry.h"
+#include "Random123/ars.h"
+#include "Random123/ReinterpretCtr.hpp"
+
+namespace random_iterator {
+
+  namespace detail {
+
+    template <typename Engine>
+    struct random_traits;
+    /*
+     * random_traits<T>::state_type { counter, state}
+     * random_traits<T>::advance_type;
+     * random_traits<T>::init_type;
+     * random_traits<T>::result_type;
+     */
+
+    // philox
+    template <>
+    struct random_traits<random_iterator_r123::Philox2x64> {
+      typedef typename random_iterator_r123::Philox2x64::ctr_type state_type;
+      typedef typename random_iterator_r123::Philox2x64::key_type seed_type;
+      typedef typename random_iterator_r123::Philox2x64::ukey_type init_type;
+      typedef uint64_t advance_type;
+      typedef state_type::value_type result_type;
+
+      enum { arity = 2 };
+    };
+
+    template <>
+    struct random_traits<random_iterator_r123::Philox4x64> {
+      typedef typename random_iterator_r123::Philox4x64::ctr_type state_type;
+      typedef typename random_iterator_r123::Philox4x64::key_type seed_type;
+      typedef typename random_iterator_r123::Philox4x64::ukey_type init_type;
+      typedef uint64_t advance_type;
+      typedef state_type::value_type result_type;
+
+      enum { arity = 4 };
+    };
+    //
+    template <>
+    struct random_traits<random_iterator_r123::Threefry2x64> {
+      typedef typename random_iterator_r123::Threefry2x64::ctr_type state_type;
+      typedef typename random_iterator_r123::Threefry2x64::key_type seed_type;
+      typedef typename random_iterator_r123::Threefry2x64::ukey_type init_type;
+      typedef uint64_t advance_type;
+      typedef state_type::value_type result_type;
+
+      enum { arity = 2 };
+    };
+
+    //
+    template <>
+    struct random_traits<random_iterator_r123::Threefry4x64> {
+      typedef typename random_iterator_r123::Threefry4x64::ctr_type state_type;
+      typedef typename random_iterator_r123::Threefry4x64::key_type seed_type;
+      typedef typename random_iterator_r123::Threefry4x64::ukey_type init_type;
+      typedef uint64_t advance_type;
+      typedef state_type::value_type result_type;
+
+      enum { arity = 4 };
+    };
+
+#if RANDOM_ITERATOR_R123_USE_AES_NI
+    template <>
+    struct random_traits<random_iterator_r123::ARS4x32> {
+      typedef typename random_iterator_r123::ARS4x32::ctr_type state_type;
+      typedef typename random_iterator_r123::ARS4x32::key_type seed_type;
+      typedef typename random_iterator_r123::ARS4x32::ukey_type init_type;
+      typedef uint64_t advance_type;
+      typedef state_type::value_type result_type;
+
+      enum { arity = 4 };
+    };
+
+    template <>
+    struct random_traits<random_iterator_r123::ARS2x64> {
+      typedef typename random_iterator_r123::ARS2x64::ctr_type state_type;
+      typedef typename random_iterator_r123::ARS2x64::key_type seed_type;
+      typedef typename random_iterator_r123::ARS2x64::ukey_type init_type;
+      typedef uint64_t advance_type;
+      typedef state_type::value_type result_type;
+
+      enum { arity = 2 };
+    };
+
+#endif
+
+  } // namespace detail
+
+} // namespace random_iterator

corsika/detail/framework/random/random_iterator/detail/Macros.h

+/*----------------------------------------------------------------------------
+*
+* Copyright (C) 2021 - 2024 Antonio Augusto Alves Junior
+*
+* This file is part of RandomIterator.
+
+* Redistribution and use in source and binary forms, with or without modification,
+* are permitted provided that the following conditions are met:
+*
+* 1. Redistributions of source code must retain the above copyright notice,
+* this list of conditions and the following disclaimer.
+*
+* 2. Redistributions in binary form must reproduce the above copyright notice,
+* this list of conditions and the following disclaimer in the documentation and/or
+* other materials provided with the distribution.
+*
+* 3. Neither the name of the copyright holder nor the names of its contributors
+* may be used to endorse or promote products derived from this software without
+* specific prior written permission.
+*
+* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS “AS IS” AND
+* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+* WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
+* IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
+* INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
+* BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+* LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+* HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
+* STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+* ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
+* OF THE POSSIBILITY OF SUCH DAMAGE.
+*
+* ----------------------------------------------------------------------------*/
+/*
+ * Macros.h
+ *
+ *  Created on: 04/03/2021
+ *      Author: Antonio Augusto Alves Junior
+ */
+
+#pragma once
+
+#if defined(__has_builtin)
+    //__has_builtin(__builtin_expect)
+	#if __has_builtin(__builtin_expect)
+		 #define   _LIKELY_(x) __builtin_expect(x, 1)
+		 #define _UNLIKELY_(x) __builtin_expect(x, 0)
+	#else
+		 #define   _LIKELY_(x) x
+		 #define _UNLIKELY_(x) x
+	#endif //__has_builtin(__builtin_expect)
+#else
+	#define   _LIKELY_(x) x
+	#define _UNLIKELY_(x) x
+#endif //defined(__has_builtin)
+
+
+//detect endianess
+#if defined(__BYTE_ORDER) && __BYTE_ORDER == __BIG_ENDIAN || \
+    defined(__BIG_ENDIAN__) ||                               \
+    defined(__ARMEB__) ||                                    \
+    defined(__THUMBEB__) ||                                  \
+    defined(__AARCH64EB__) ||                                \
+    defined(_MIBSEB) || defined(__MIBSEB) || defined(__MIBSEB__)
+#ifndef __BIG_ENDIAN__
+#define __BIG_ENDIAN__
+#endif
+#elif defined(__BYTE_ORDER) && __BYTE_ORDER == __LITTLE_ENDIAN || \
+    defined(__LITTLE_ENDIAN__) ||                                 \
+    defined(__ARMEL__) ||                                         \
+    defined(__THUMBEL__) ||                                       \
+    defined(__AARCH64EL__) ||                                     \
+    defined(_MIPSEL) || defined(__MIPSEL) || defined(__MIPSEL__) || \
+    defined(_WIN32) || defined(__i386__) || defined(__x86_64__) || \
+    defined(_X86_) || defined(_IA64_)
+#ifndef __LITTLE_ENDIAN__
+#define __LITTLE_ENDIAN__
+#endif
+#else
+#error "I don't know what architecture this is!"
+#endif

corsika/detail/framework/random/random_iterator/detail/Random123/MicroURNG.hpp

+/*
+Copyright 2010-2011, D. E. Shaw Research.
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are
+met:
+
+* Redistributions of source code must retain the above copyright
+  notice, this list of conditions, and the following disclaimer.
+
+* Redistributions in binary form must reproduce the above copyright
+  notice, this list of conditions, and the following disclaimer in the
+  documentation and/or other materials provided with the distribution.
+
+* Neither the name of D. E. Shaw Research nor the names of its
+  contributors may be used to endorse or promote products derived from
+  this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+*/
+#ifndef __MicroURNG_dot_hpp__
+#define __MicroURNG_dot_hpp__
+
+#include <stdexcept>
+#include <limits>
+
+namespace random_iterator_r123 {
+  /**
+      Given a CBRNG whose ctr_type has an unsigned integral value_type,
+      MicroURNG<CBRNG>(c, k) is a type that satisfies the
+      requirements of a C++11 Uniform Random Number Generator.
+
+      The intended purpose is for a MicroURNG to be passed
+      as an argument to a C++11 Distribution, e.g.,
+      std::normal_distribution.  See examples/MicroURNG.cpp.
+
+      The MicroURNG functor has a period of "only"
+
+         ctr_type.size()*2^32,
+
+      after which it will silently repeat.
+
+      The high 32 bits of the highest word in the counter c, passed to
+      the constructor must be zero.  MicroURNG uses these bits to
+      "count".
+
+      Older versions of the library permitted a second template
+      parameter by which the caller could control the number of
+      bits devoted to the URNG's internal counter.  This flexibility
+      has been disabled because URNGs created with different
+      numbers of counter bits could, conceivably "collide".
+
+  \code
+         typedef ?someCBRNG? RNG;
+         RNG::ctr_type c = ...; // under application control
+         RNG::key_type k = ...; //
+         std::normal_distribution<float> nd;
+         MicroURNG<RNG> urng(c, k);
+         for(???){
+           ...
+           nd(urng);  // may be called several hundred times with BITS=10
+           ...
+         }
+  \endcode
+  */
+
+  template <typename CBRNG>
+  class MicroURNG {
+    // According to C++11, a URNG requires only a result_type,
+    // operator()(), min() and max() methods.  Everything else
+    // (ctr_type, key_type, reset() method, etc.) is "value added"
+    // for the benefit of users that "know" that they're dealing with
+    // a MicroURNG.
+  public:
+    typedef CBRNG cbrng_type;
+    static const int BITS = 32;
+    typedef typename cbrng_type::ctr_type ctr_type;
+    typedef typename cbrng_type::key_type key_type;
+    typedef typename cbrng_type::ukey_type ukey_type;
+    typedef typename ctr_type::value_type result_type;
+
+    RANDOM_ITERATOR_R123_STATIC_ASSERT(std::numeric_limits<result_type>::digits >= BITS,
+                                       "The result_type must have at least 32 bits");
+
+    result_type operator()() {
+      if (last_elem == 0) {
+        // jam n into the high bits of c
+        const size_t W = std::numeric_limits<result_type>::digits;
+        ctr_type c = c0;
+        c[c0.size() - 1] |= n << (W - BITS);
+        rdata = b(c, k);
+        n++;
+        last_elem = rdata.size();
+      }
+      return rdata[--last_elem];
+    }
+    MicroURNG(cbrng_type _b, ctr_type _c0, ukey_type _uk)
+        : b(_b)
+        , c0(_c0)
+        , k(_uk)
+        , n(0)
+        , last_elem(0) {
+      chkhighbits();
+    }
+    MicroURNG(ctr_type _c0, ukey_type _uk)
+        : b()
+        , c0(_c0)
+        , k(_uk)
+        , n(0)
+        , last_elem(0) {
+      chkhighbits();
+    }
+
+    // _Min and _Max work around a bug in the library shipped with MacOS Xcode 4.5.2.
+    // See the commment in conventional/Engine.hpp.
+    const static result_type _Min = 0;
+    const static result_type _Max = ~((result_type)0);
+
+    static RANDOM_ITERATOR_R123_CONSTEXPR result_type min
+    RANDOM_ITERATOR_R123_NO_MACRO_SUBST() {
+      return _Min;
+    }
+    static RANDOM_ITERATOR_R123_CONSTEXPR result_type max
+    RANDOM_ITERATOR_R123_NO_MACRO_SUBST() {
+      return _Max;
+    }
+    // extra methods:
+    const ctr_type& counter() const { return c0; }
+    void reset(ctr_type _c0, ukey_type _uk) {
+      c0 = _c0;
+      chkhighbits();
+      k = _uk;
+      n = 0;
+      last_elem = 0;
+    }
+
+  private:
+    cbrng_type b;
+    ctr_type c0;
+    key_type k;
+    RANDOM_ITERATOR_R123_ULONG_LONG n;
+    size_t last_elem;
+    ctr_type rdata;
+    void chkhighbits() {
+      result_type r = c0[c0.size() - 1];
+      result_type mask = ((uint64_t)std::numeric_limits<result_type>::max
+                              RANDOM_ITERATOR_R123_NO_MACRO_SUBST()) >>
+                         BITS;
+      if ((r & mask) != r)
+        throw std::runtime_error("MicroURNG: c0, does not have high bits clear");
+    }
+  };
+} // namespace random_iterator_r123
+#endif

corsika/detail/framework/random/random_iterator/detail/Random123/ReinterpretCtr.hpp

+/*
+Copyright 2010-2011, D. E. Shaw Research.
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are
+met:
+
+* Redistributions of source code must retain the above copyright
+  notice, this list of conditions, and the following disclaimer.
+
+* Redistributions in binary form must reproduce the above copyright
+  notice, this list of conditions, and the following disclaimer in the
+  documentation and/or other materials provided with the distribution.
+
+* Neither the name of D. E. Shaw Research nor the names of its
+  contributors may be used to endorse or promote products derived from
+  this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+*/
+#ifndef __ReinterpretCtr_dot_hpp__
+#define __ReinterpretCtr_dot_hpp__
+
+#include "features/compilerfeatures.h"
+#include <cstring>
+
+namespace random_iterator_r123 {
+  /*!
+    ReinterpretCtr uses memcpy to map back and forth
+    between a CBRNG's ctr_type and the specified ToType.  For example,
+    after:
+
+      typedef ReinterpretCtr<r123array4x32, Philox2x64> G;
+
+    G is a bona fide CBRNG with ctr_type r123array4x32.
+
+    WARNING:  ReinterpretCtr is endian dependent.  The
+    values returned by G, declared as above,
+    will depend on the endianness of the machine on which it runs.
+   */
+
+  template <typename ToType, typename CBRNG>
+  struct ReinterpretCtr {
+    typedef ToType ctr_type;
+    typedef typename CBRNG::key_type key_type;
+    typedef typename CBRNG::ctr_type bctype;
+    typedef typename CBRNG::ukey_type ukey_type;
+    RANDOM_ITERATOR_R123_STATIC_ASSERT(
+        sizeof(ToType) == sizeof(bctype) && sizeof(typename bctype::value_type) != 16,
+        "ReinterpretCtr:  sizeof(ToType) is not the same as sizeof(CBRNG::ctr_type) or "
+        "CBRNG::ctr_type::value_type looks like it might be __m128i");
+    // It's amazingly difficult to safely do conversions with __m128i.
+    // If we use the operator() implementation below with a CBRNG
+    // whose ctr_type is r123array1xm128i, gcc4.6 optimizes away the
+    // memcpys, inlines the operator()(c,k), and produces assembly
+    // language that ends with an aesenclast instruction with a
+    // destination operand pointing to an unaligned memory address ...
+    // Segfault!  See:  http://gcc.gnu.org/bugzilla/show_bug.cgi?id=50444
+    // MSVC also produces code that crashes.  We suspect a
+    // similar mechanism but haven't done the debugging necessary to
+    // be sure.  We were able to 'fix' gcc4.6 by making bc a mutable
+    // data member rather than declaring it in the scope of
+    // operator().  That didn't fix the MSVC problems, though.
+    //
+    // Conclusion - don't touch __m128i, at least for now.  The
+    // easiest (but highly imprecise) way to do that is the static
+    // assertion above that rejects bctype::value_types of size 16. -
+    // Sep 2011.
+    ctr_type operator()(ctr_type c, key_type k) {
+      bctype bc;
+      std::memcpy(&bc, &c, sizeof(c));
+      CBRNG b;
+      bc = b(bc, k);
+      std::memcpy(&c, &bc, sizeof(bc));
+      return c;
+    }
+  };
+} // namespace random_iterator_r123
+#endif

corsika/detail/framework/random/random_iterator/detail/Random123/aes.h

+/*
+Copyright 2010-2011, D. E. Shaw Research.
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are
+met:
+
+* Redistributions of source code must retain the above copyright
+  notice, this list of conditions, and the following disclaimer.
+
+* Redistributions in binary form must reproduce the above copyright
+  notice, this list of conditions, and the following disclaimer in the
+  documentation and/or other materials provided with the distribution.
+
+* Neither the name of D. E. Shaw Research nor the names of its
+  contributors may be used to endorse or promote products derived from
+  this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+*/
+#ifndef __Random123_aes_dot_hpp__
+#define __Random123_aes_dot_hpp__
+
+#include "features/compilerfeatures.h"
+#include "array.h"
+
+/* Implement a bona fide AES block cipher.  It's minimally
+// checked against the test vector in FIPS-197 in ut_aes.cpp. */
+#if RANDOM_ITERATOR_R123_USE_AES_NI
+
+/** @ingroup AESNI */
+typedef struct r123array1xm128i aesni1xm128i_ctr_t;
+/** @ingroup AESNI */
+typedef struct r123array1xm128i aesni1xm128i_ukey_t;
+/** @ingroup AESNI */
+typedef struct r123array4x32 aesni4x32_ukey_t;
+/** @ingroup AESNI */
+enum r123_enum_aesni1xm128i { aesni1xm128i_rounds = 10 };
+
+/** \cond HIDDEN_FROM_DOXYGEN */
+RANDOM_ITERATOR_R123_STATIC_INLINE __m128i AES_128_ASSIST (__m128i temp1, __m128i temp2) { 
+    __m128i temp3; 
+    temp2 = _mm_shuffle_epi32 (temp2 ,0xff); 
+    temp3 = _mm_slli_si128 (temp1, 0x4);
+    temp1 = _mm_xor_si128 (temp1, temp3);
+    temp3 = _mm_slli_si128 (temp3, 0x4);
+    temp1 = _mm_xor_si128 (temp1, temp3);
+    temp3 = _mm_slli_si128 (temp3, 0x4);
+    temp1 = _mm_xor_si128 (temp1, temp3);
+    temp1 = _mm_xor_si128 (temp1, temp2); 
+    return temp1; 
+}
+
+RANDOM_ITERATOR_R123_STATIC_INLINE void aesni1xm128iexpand(aesni1xm128i_ukey_t uk, __m128i ret[11])
+{
+    __m128i rkey = uk.v[0].m;
+    __m128i tmp2;
+
+    ret[0] = rkey;
+    tmp2 = _mm_aeskeygenassist_si128(rkey, 0x1);
+    rkey = AES_128_ASSIST(rkey, tmp2);
+    ret[1] = rkey;
+
+    tmp2 = _mm_aeskeygenassist_si128(rkey, 0x2);
+    rkey = AES_128_ASSIST(rkey, tmp2);
+    ret[2] = rkey;
+
+    tmp2 = _mm_aeskeygenassist_si128(rkey, 0x4);
+    rkey = AES_128_ASSIST(rkey, tmp2);
+    ret[3] = rkey;
+
+    tmp2 = _mm_aeskeygenassist_si128(rkey, 0x8);
+    rkey = AES_128_ASSIST(rkey, tmp2);
+    ret[4] = rkey;
+
+    tmp2 = _mm_aeskeygenassist_si128(rkey, 0x10);
+    rkey = AES_128_ASSIST(rkey, tmp2);
+    ret[5] = rkey;
+
+    tmp2 = _mm_aeskeygenassist_si128(rkey, 0x20);
+    rkey = AES_128_ASSIST(rkey, tmp2);
+    ret[6] = rkey;
+
+    tmp2 = _mm_aeskeygenassist_si128(rkey, 0x40);
+    rkey = AES_128_ASSIST(rkey, tmp2);
+    ret[7] = rkey;
+
+    tmp2 = _mm_aeskeygenassist_si128(rkey, 0x80);
+    rkey = AES_128_ASSIST(rkey, tmp2);
+    ret[8] = rkey;
+
+    tmp2 = _mm_aeskeygenassist_si128(rkey, 0x1b);
+    rkey = AES_128_ASSIST(rkey, tmp2);
+    ret[9] = rkey;
+
+    tmp2 = _mm_aeskeygenassist_si128(rkey, 0x36);
+    rkey = AES_128_ASSIST(rkey, tmp2);
+    ret[10] = rkey;
+}
+/** \endcond */
+    
+#ifdef __cplusplus
+/** @ingroup AESNI */
+struct aesni1xm128i_key_t{ 
+    __m128i k[11]; 
+    aesni1xm128i_key_t(){
+        aesni1xm128i_ukey_t uk;
+        uk.v[0].m = _mm_setzero_si128();
+        aesni1xm128iexpand(uk, k);
+    }
+    aesni1xm128i_key_t(const aesni1xm128i_ukey_t& uk){
+        aesni1xm128iexpand(uk, k);
+    }
+    aesni1xm128i_key_t(const aesni4x32_ukey_t& uk){
+        aesni1xm128i_ukey_t uk128;
+        uk128.v[0].m = _mm_set_epi32(uk.v[3], uk.v[2], uk.v[1], uk.v[0]);
+        aesni1xm128iexpand(uk128, k);
+    }
+    aesni1xm128i_key_t& operator=(const aesni1xm128i_ukey_t& uk){
+        aesni1xm128iexpand(uk, k);
+        return *this;
+    }
+    aesni1xm128i_key_t& operator=(const aesni4x32_ukey_t& uk){
+        aesni1xm128i_ukey_t uk128;
+        uk128.v[0].m = _mm_set_epi32(uk.v[3], uk.v[2], uk.v[1], uk.v[0]);
+        aesni1xm128iexpand(uk128, k);
+        return *this;
+    }
+    bool operator==(const aesni1xm128i_key_t& rhs) const{
+        for(int i=0; i<11; ++i){
+            // Sigh... No r123m128i(__m128i) constructor!
+            r123m128i li; li.m = k[i];
+            r123m128i ri; ri.m = rhs.k[i];
+            if( li != ri ) return false;
+        }
+        return true;
+    }
+    bool operator!=(const aesni1xm128i_key_t& rhs) const{
+        return !(*this == rhs);
+    }
+    friend std::ostream& operator<<(std::ostream& os, const aesni1xm128i_key_t& v){
+        r123m128i ki;
+        for(int i=0; i<10; ++i){
+            ki.m = v.k[i];
+            os << ki << " ";
+        }
+        ki.m = v.k[10];
+        return os << ki;
+    }
+    friend std::istream& operator>>(std::istream& is, aesni1xm128i_key_t& v){
+        r123m128i ki;
+        for(int i=0; i<11; ++i){
+            is >> ki;
+            v.k[i] = ki;
+        }
+        return is;
+    }
+};
+#else
+typedef struct { 
+    __m128i k[11]; 
+}aesni1xm128i_key_t;
+
+/** @ingroup AESNI */
+RANDOM_ITERATOR_R123_STATIC_INLINE aesni1xm128i_key_t aesni1xm128ikeyinit(aesni1xm128i_ukey_t uk){
+    aesni1xm128i_key_t ret;
+    aesni1xm128iexpand(uk, ret.k);
+    return ret;
+}
+#endif
+
+/** @ingroup AESNI */
+RANDOM_ITERATOR_R123_STATIC_INLINE aesni1xm128i_ctr_t aesni1xm128i(aesni1xm128i_ctr_t in, aesni1xm128i_key_t k) {
+    __m128i x = _mm_xor_si128(k.k[0], in.v[0].m);
+    x = _mm_aesenc_si128(x, k.k[1]);
+    x = _mm_aesenc_si128(x, k.k[2]);
+    x = _mm_aesenc_si128(x, k.k[3]);
+    x = _mm_aesenc_si128(x, k.k[4]);
+    x = _mm_aesenc_si128(x, k.k[5]);
+    x = _mm_aesenc_si128(x, k.k[6]);
+    x = _mm_aesenc_si128(x, k.k[7]);
+    x = _mm_aesenc_si128(x, k.k[8]);
+    x = _mm_aesenc_si128(x, k.k[9]);
+    x = _mm_aesenclast_si128(x, k.k[10]);
+    {
+      aesni1xm128i_ctr_t ret;
+      ret.v[0].m = x;
+      return ret;
+    }
+}
+
+/** @ingroup AESNI */
+RANDOM_ITERATOR_R123_STATIC_INLINE aesni1xm128i_ctr_t aesni1xm128i_R(unsigned R, aesni1xm128i_ctr_t in, aesni1xm128i_key_t k){
+    RANDOM_ITERATOR_R123_ASSERT(R==10);
+    return aesni1xm128i(in, k);
+}
+
+
+/** @ingroup AESNI */
+typedef struct r123array4x32 aesni4x32_ctr_t;
+/** @ingroup AESNI */
+typedef aesni1xm128i_key_t aesni4x32_key_t;
+/** @ingroup AESNI */
+enum r123_enum_aesni4x32 { aesni4x32_rounds = 10 };
+/** @ingroup AESNI */
+RANDOM_ITERATOR_R123_STATIC_INLINE aesni4x32_key_t aesni4x32keyinit(aesni4x32_ukey_t uk){
+    aesni1xm128i_ukey_t uk128;
+    aesni4x32_key_t ret;
+    uk128.v[0].m = _mm_set_epi32(uk.v[3], uk.v[2], uk.v[1], uk.v[0]);
+    aesni1xm128iexpand(uk128, ret.k);
+    return ret;
+}
+
+/** @ingroup AESNI */
+/** The aesni4x32_R function provides a C API to the @ref AESNI "AESNI" CBRNG, allowing the number of rounds to be specified explicitly **/
+RANDOM_ITERATOR_R123_STATIC_INLINE aesni4x32_ctr_t aesni4x32_R(unsigned int Nrounds, aesni4x32_ctr_t c, aesni4x32_key_t k){
+    aesni1xm128i_ctr_t c128;
+    c128.v[0].m = _mm_set_epi32(c.v[3], c.v[2], c.v[1], c.v[0]);
+    c128 = aesni1xm128i_R(Nrounds, c128, k);
+    _mm_storeu_si128((__m128i*)&c.v[0], c128.v[0].m);
+    return c;
+}
+
+#define aesni4x32_rounds aesni1xm128i_rounds
+
+/** The aesni4x32 macro provides a C API to the @ref AESNI "AESNI" CBRNG, uses the default number of rounds i.e. \c aesni4x32_rounds **/
+/** @ingroup AESNI */
+#define aesni4x32(c,k) aesni4x32_R(aesni4x32_rounds, c, k)
+
+#ifdef __cplusplus
+namespace random_iterator_r123{
+/** 
+@defgroup AESNI ARS and AESNI Classes and Typedefs
+
+The ARS4x32, ARS1xm128i, AESNI4x32 and AESNI1xm128i classes export the member functions, typedefs and
+operator overloads required by a @ref CBRNG "CBRNG" class.
+
+ARS1xm128i and AESNI1xm128i are based on the AES block cipher and rely on the AES-NI hardware instructions
+available on some some new (2011) CPUs.
+
+The ARS1xm128i CBRNG and the use of AES for random number generation are described in 
+<a href="http://dl.acm.org/citation.cfm?doid=2063405"><i>Parallel Random Numbers:  As Easy as 1, 2, 3</i> </a>.
+Although it uses some cryptographic primitives, ARS1xm128i uses a cryptographically weak key schedule and is \b not suitable for cryptographic use.
+
+@class AESNI1xm128i
+@ingroup AESNI
+AESNI exports the member functions, typedefs and operator overloads required by a @ref CBRNG class.
+
+AESNI1xm128i uses the crypotgraphic AES round function, including the cryptographic key schedule.
+
+In contrast to the other CBRNGs in the Random123 library, the AESNI1xm128i_R::key_type is opaque
+and is \b not identical to the AESNI1xm128i_R::ukey_type.  Creating a key_type, using either the constructor
+or assignment operator, is significantly more time-consuming than running the bijection (hundreds
+of clock cycles vs. tens of clock cycles).
+
+AESNI1xm128i is only available when the feature-test macro RANDOM_ITERATOR_R123_USE_AES_NI is true, which
+should occur only when the compiler is configured to generate AES-NI instructions (or
+when defaults are overridden by compile-time, compiler-command-line options).
+
+As of September 2011, the authors know of no statistical flaws with AESNI1xm128i.  It
+would be an event of major cryptographic note if any such flaws were ever found.
+*/
+struct AESNI1xm128i{
+    typedef aesni1xm128i_ctr_t ctr_type;
+    typedef aesni1xm128i_ukey_t ukey_type;
+    typedef aesni1xm128i_key_t key_type;
+    static const unsigned int rounds=10;
+    ctr_type operator()(ctr_type ctr, key_type key) const{
+        return aesni1xm128i(ctr, key);
+    }
+};
+
+/* @class AESNI4x32 */
+struct AESNI4x32{
+    typedef aesni4x32_ctr_t ctr_type;
+    typedef aesni4x32_ukey_t ukey_type;
+    typedef aesni4x32_key_t key_type;
+    static const unsigned int rounds=10;
+    ctr_type operator()(ctr_type ctr, key_type key) const{
+        return aesni4x32(ctr, key);
+    }
+};
+
+/** @ingroup AESNI
+    @class AESNI1xm128i_R
+
+AESNI1xm128i_R is provided for completeness, but is only instantiable with ROUNDS=10, in
+which case it is identical to AESNI1xm128i */
+template <unsigned ROUNDS=10> 
+struct AESNI1xm128i_R : public AESNI1xm128i{
+    RANDOM_ITERATOR_R123_STATIC_ASSERT(ROUNDS==10, "AESNI1xm128i_R<R> is only valid with R=10");
+};
+
+/** @class AESNI4x32_R **/
+template <unsigned ROUNDS=10> 
+struct AESNI4x32_R : public AESNI4x32{
+    RANDOM_ITERATOR_R123_STATIC_ASSERT(ROUNDS==10, "AESNI4x32_R<R> is only valid with R=10");
+};
+} // namespace random_iterator_r123
+#endif /* __cplusplus */
+
+#endif /* RANDOM_ITERATOR_R123_USE_AES_NI */
+
+#if RANDOM_ITERATOR_R123_USE_AES_OPENSSL
+#include "string.h"
+#include <openssl/aes.h>
+typedef struct r123array16x8 aesopenssl16x8_ctr_t;
+typedef struct r123array16x8 aesopenssl16x8_ukey_t;
+#ifdef __cplusplus
+struct aesopenssl16x8_key_t{
+    AES_KEY k;
+    aesopenssl16x8_key_t(){
+        aesopenssl16x8_ukey_t ukey={{}};
+        AES_set_encrypt_key((const unsigned char *)&ukey.v[0], 128, &k);
+    }
+    aesopenssl16x8_key_t(const aesopenssl16x8_ukey_t& ukey){
+        AES_set_encrypt_key((const unsigned char *)&ukey.v[0], 128, &k);
+    }
+    aesopenssl16x8_key_t& operator=(const aesopenssl16x8_ukey_t& ukey){
+        AES_set_encrypt_key((const unsigned char *)&ukey.v[0], 128, &k);
+        return *this;
+    }
+    bool operator==(const aesopenssl16x8_key_t& rhs) const{
+        return (k.rounds == rhs.k.rounds) && 0==::memcmp(&k.rd_key[0], &rhs.k.rd_key[0], (k.rounds+1) * 4 * sizeof(uint32_t));
+    }
+    bool operator!=(const aesopenssl16x8_key_t& rhs) const{
+        return !(*this == rhs);
+    }
+    friend std::ostream& operator<<(std::ostream& os, const aesopenssl16x8_key_t& v){
+        os << v.k.rounds;
+        const unsigned int *p = &v.k.rd_key[0];
+        for(int i=0; i<(v.k.rounds+1); ++i){
+            os << " " << p[0] << " " << p[1] << " " << p[2] << " " << p[3];
+            p += 4;
+        }
+        return os;
+    }
+    friend std::istream& operator>>(std::istream& is, aesopenssl16x8_key_t& v){
+        is >> v.k.rounds;
+        unsigned int *p = &v.k.rd_key[0];
+        for(int i=0; i<(v.k.rounds+1); ++i){
+            is >> p[0] >> p[1] >> p[2] >> p[3];
+            p += 4;
+        }
+        return is;
+    }
+};
+#else
+typedef struct aesopenssl16x8_key_t{
+    AES_KEY k;
+}aesopenssl16x8_key_t;
+RANDOM_ITERATOR_R123_STATIC_INLINE struct aesopenssl16x8_key_t aesopenssl16x8keyinit(aesopenssl16x8_ukey_t uk){
+    aesopenssl16x8_key_t ret;
+    AES_set_encrypt_key((const unsigned char *)&uk.v[0], 128, &ret.k);
+    return ret;
+}
+#endif
+
+RANDOM_ITERATOR_R123_STATIC_INLINE RANDOM_ITERATOR_R123_FORCE_INLINE(aesopenssl16x8_ctr_t aesopenssl16x8_R(aesopenssl16x8_ctr_t ctr, aesopenssl16x8_key_t key));
+RANDOM_ITERATOR_R123_STATIC_INLINE
+aesopenssl16x8_ctr_t aesopenssl16x8_R(aesopenssl16x8_ctr_t ctr, aesopenssl16x8_key_t key){
+    aesopenssl16x8_ctr_t ret;
+    AES_encrypt((const unsigned char*)&ctr.v[0], (unsigned char *)&ret.v[0], &key.k);
+    return ret;
+}
+
+#define aesopenssl16x8_rounds aesni4x32_rounds
+#define aesopenssl16x8(c,k) aesopenssl16x8_R(aesopenssl16x8_rounds)
+
+#ifdef __cplusplus
+namespace random_iterator_r123{
+struct AESOpenSSL16x8{
+    typedef aesopenssl16x8_ctr_t ctr_type;
+    typedef aesopenssl16x8_key_t key_type;
+    typedef aesopenssl16x8_ukey_t ukey_type;
+    static const unsigned int rounds=10;
+    ctr_type operator()(const ctr_type& in, const key_type& k){
+        ctr_type out;
+        AES_encrypt((const unsigned char *)&in[0], (unsigned char *)&out[0], &k.k);
+        return out;
+    }
+};
+} // namespace random_iterator_r123
+#endif /* __cplusplus */
+#endif /* RANDOM_ITERATOR_R123_USE_AES_OPENSSL */
+
+#endif