Framework/ProcessSequence/DecayProcess.h

-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#pragma once
-
-#include <corsika/process/BaseProcess.h>
-#include <corsika/process/ProcessReturn.h> // for convenience
-#include <corsika/setup/SetupTrajectory.h>
-#include <corsika/units/PhysicalUnits.h>
-
-namespace corsika::process {
-
-  /**
-     \class DecayProcess
-
-     The structural base type of a process object in a
-     ProcessSequence. Both, the ProcessSequence and all its elements
-     are of type DecayProcess<T>
-
-   */
-
-  template <typename TDerived>
-  struct DecayProcess : BaseProcess<TDerived> {
-
-    using BaseProcess<TDerived>::GetRef;
-
-    /// here starts the interface-definition part
-    // -> enforce TDerived to implement DoDecay...
-    template <typename TParticle>
-    EProcessReturn DoDecay(TParticle&);
-
-    template <typename TParticle>
-    corsika::units::si::TimeType GetLifetime(TParticle& p);
-
-    template <typename TParticle>
-    corsika::units::si::InverseTimeType GetInverseLifetime(TParticle& vP) {
-      return 1. / GetRef().GetLifetime(vP);
-    }
-  };
-
-} // namespace corsika::process

Framework/ProcessSequence/InteractionProcess.h

-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#pragma once
-
-#include <corsika/process/BaseProcess.h>
-#include <corsika/process/ProcessReturn.h> // for convenience
-#include <corsika/setup/SetupTrajectory.h>
-#include <corsika/units/PhysicalUnits.h>
-
-namespace corsika::process {
-
-  /**
-     \class InteractionProcess
-
-     The structural base type of a process object in a
-     ProcessSequence. Both, the ProcessSequence and all its elements
-     are of type InteractionProcess<T>
-
-   */
-
-  template <typename TDerived>
-  struct InteractionProcess : public BaseProcess<TDerived> {
-
-    using BaseProcess<TDerived>::GetRef;
-
-    /// here starts the interface-definition part
-    // -> enforce TDerived to implement DoInteraction...
-    template <typename TParticle>
-    EProcessReturn DoInteraction(TParticle&);
-
-    template <typename TParticle>
-    corsika::units::si::GrammageType GetInteractionLength(TParticle& p);
-
-    template <typename TParticle>
-    corsika::units::si::InverseGrammageType GetInverseInteractionLength(TParticle& p) {
-      return 1. / GetRef().GetInteractionLength(p);
-    }
-  };
-
-} // namespace corsika::process

Framework/ProcessSequence/ProcessReturn.h

-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#pragma once
-
-namespace corsika::process {
-
-  /**
-     since in a process sequence many status updates can accumulate
-     for a single particle, this enum should define only bit-flags
-     that can be accumulated easily with "|="
-   */
-
-  enum class EProcessReturn : int {
-    eOk = (1 << 0),
-    eParticleAbsorbed = (1 << 2),
-    eInteracted = (1 << 3),
-    eDecayed = (1 << 4),
-  };
-
-  inline EProcessReturn operator|(EProcessReturn a, EProcessReturn b) {
-    return static_cast<EProcessReturn>(static_cast<int>(a) | static_cast<int>(b));
-  }
-
-  inline EProcessReturn& operator|=(EProcessReturn& a, const EProcessReturn b) {
-    return a = a | b;
-  }
-
-  inline EProcessReturn operator&(const EProcessReturn a, const EProcessReturn b) {
-    return static_cast<EProcessReturn>(static_cast<int>(a) & static_cast<int>(b));
-  }
-
-  inline bool operator==(const EProcessReturn a, const EProcessReturn b) {
-    return (static_cast<int>(a) & static_cast<int>(b)) != 0;
-  }
-
-  inline bool isAbsorbed(const EProcessReturn a) {
-    return static_cast<int>(a & EProcessReturn::eParticleAbsorbed);
-  }
-
-} // namespace corsika::process

Framework/ProcessSequence/ProcessSequence.h

-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#pragma once
-
-#include <corsika/process/BaseProcess.h>
-#include <corsika/process/BoundaryCrossingProcess.h>
-#include <corsika/process/ContinuousProcess.h>
-#include <corsika/process/DecayProcess.h>
-#include <corsika/process/InteractionProcess.h>
-#include <corsika/process/ProcessReturn.h>
-#include <corsika/process/SecondariesProcess.h>
-#include <corsika/process/StackProcess.h>
-#include <corsika/units/PhysicalUnits.h>
-
-#include <cmath>
-#include <limits>
-#include <type_traits>
-
-namespace corsika::process {
-
-  /**
-     \class ProcessSequence
-
-     A compile time static list of processes. The compiler will
-     generate a new type based on template logic containing all the
-     elements.
-
-     \comment Using CRTP pattern,
-     https://en.wikipedia.org/wiki/Curiously_recurring_template_pattern
-   */
-
-  // this is a marker to track which BaseProcess is also a ProcessSequence
-  template <typename TClass>
-  struct is_process_sequence : std::false_type {};
-
-  template <typename TClass>
-  bool constexpr is_process_sequence_v = is_process_sequence<TClass>::value;
-
-  // we also need a marker to identify SwitchProcess
-  namespace switch_process {
-    template <typename TProcess1, typename TProcess2>
-    class SwitchProcess; // fwd-decl.
-  }
-
-  // to detect SwitchProcesses inside the ProcessSequence
-  template <typename TClass>
-  struct is_switch_process : std::false_type {};
-
-  template <typename TClass>
-  bool constexpr is_switch_process_v = is_switch_process<TClass>::value;
-
-  template <typename Process1, typename Process2>
-  struct is_switch_process<switch_process::SwitchProcess<Process1, Process2>>
-      : std::true_type {};
-
-  /**
-     T1 and T2 are both references if possible (lvalue), otherwise
-     (rvalue) they are just classes. This allows us to handle both,
-     rvalue as well as lvalue Processes in the ProcessSequence.
-   */
-  template <typename T1, typename T2>
-  class ProcessSequence : public BaseProcess<ProcessSequence<T1, T2>> {
-
-    using T1type = typename std::decay<T1>::type;
-    using T2type = typename std::decay<T2>::type;
-
-    static bool constexpr t1ProcSeq = is_process_sequence_v<T1type>;
-    static bool constexpr t2ProcSeq = is_process_sequence_v<T2type>;
-
-    static bool constexpr t1SwitchProc = is_switch_process_v<T1type>;
-    static bool constexpr t2SwitchProc = is_switch_process_v<T2type>;
-
-  public:
-    T1 A; // this is a reference, if possible
-    T2 B; // this is a reference, if possible
-
-    ProcessSequence(T1 in_A, T2 in_B)
-        : A(in_A)
-        , B(in_B) {}
-
-    template <typename Particle, typename VTNType>
-    EProcessReturn DoBoundaryCrossing(Particle& p, VTNType const& from,
-                                      VTNType const& to) {
-      EProcessReturn ret = EProcessReturn::eOk;
-
-      if constexpr (std::is_base_of_v<BoundaryCrossingProcess<T1type>, T1type> ||
-                    t1ProcSeq) {
-        ret |= A.DoBoundaryCrossing(p, from, to);
-      }
-
-      if constexpr (std::is_base_of_v<BoundaryCrossingProcess<T2type>, T2type> ||
-                    t2ProcSeq) {
-        ret |= B.DoBoundaryCrossing(p, from, to);
-      }
-
-      return ret;
-    }
-
-    template <typename TParticle, typename TTrack>
-    EProcessReturn DoContinuous(TParticle& vP, TTrack& vT) {
-      EProcessReturn ret = EProcessReturn::eOk;
-      if constexpr (std::is_base_of_v<ContinuousProcess<T1type>, T1type> || t1ProcSeq) {
-        if (!isAbsorbed(ret)) { ret |= A.DoContinuous(vP, vT); }
-      }
-      if constexpr (std::is_base_of_v<ContinuousProcess<T2type>, T2type> || t2ProcSeq) {
-        if (!isAbsorbed(ret)) { ret |= B.DoContinuous(vP, vT); }
-      }
-      return ret;
-    }
-
-    template <typename TSecondaries>
-    EProcessReturn DoSecondaries(TSecondaries& vS) {
-      EProcessReturn ret = EProcessReturn::eOk;
-      if constexpr (std::is_base_of_v<SecondariesProcess<T1type>, T1type> || t1ProcSeq) {
-        ret |= A.DoSecondaries(vS);
-      }
-      if constexpr (std::is_base_of_v<SecondariesProcess<T2type>, T2type> || t2ProcSeq) {
-        ret |= B.DoSecondaries(vS);
-      }
-      return ret;
-    }
-
-    /**
-       The processes of type StackProcess do have an internal counter,
-       so they can be exectuted only each N steps. Often these are
-       "maintenacne processes" that do not need to run after each
-       single step of the simulations. In the CheckStep function it is
-       tested if either A or B are StackProcess and if they are due
-       for execution.
-     */
-    bool CheckStep() {
-      bool ret = false;
-      if constexpr (std::is_base_of_v<StackProcess<T1type>, T1type> || t1ProcSeq) {
-        ret |= A.CheckStep();
-      }
-      if constexpr (std::is_base_of_v<StackProcess<T2type>, T2type> || t2ProcSeq) {
-        ret |= B.CheckStep();
-      }
-      return ret;
-    }
-
-    /**
-       Execute the StackProcess-es in the ProcessSequence
-     */
-    template <typename TStack>
-    EProcessReturn DoStack(TStack& vS) {
-      EProcessReturn ret = EProcessReturn::eOk;
-      if constexpr (std::is_base_of_v<StackProcess<T1type>, T1type> || t1ProcSeq) {
-        if (A.CheckStep()) { ret |= A.DoStack(vS); }
-      }
-      if constexpr (std::is_base_of_v<StackProcess<T2type>, T2type> || t2ProcSeq) {
-        if (B.CheckStep()) { ret |= B.DoStack(vS); }
-      }
-      return ret;
-    }
-
-    template <typename TParticle, typename TTrack>
-    corsika::units::si::LengthType MaxStepLength(TParticle& vP, TTrack& vTrack) {
-      corsika::units::si::LengthType
-          max_length = // if no other process in the sequence implements it
-          std::numeric_limits<double>::infinity() * corsika::units::si::meter;
-
-      if constexpr (std::is_base_of_v<ContinuousProcess<T1type>, T1type> || t1ProcSeq) {
-        corsika::units::si::LengthType const len = A.MaxStepLength(vP, vTrack);
-        max_length = std::min(max_length, len);
-      }
-      if constexpr (std::is_base_of_v<ContinuousProcess<T2type>, T2type> || t2ProcSeq) {
-        corsika::units::si::LengthType const len = B.MaxStepLength(vP, vTrack);
-        max_length = std::min(max_length, len);
-      }
-      return max_length;
-    }
-
-    template <typename TParticle>
-    corsika::units::si::GrammageType GetTotalInteractionLength(TParticle& vP) {
-      return 1. / GetInverseInteractionLength(vP);
-    }
-
-    template <typename TParticle>
-    corsika::units::si::InverseGrammageType GetTotalInverseInteractionLength(
-        TParticle& vP) {
-      return GetInverseInteractionLength(vP);
-    }
-
-    template <typename TParticle>
-    corsika::units::si::InverseGrammageType GetInverseInteractionLength(TParticle& vP) {
-      using namespace corsika::units::si;
-
-      InverseGrammageType tot = 0 * meter * meter / gram;
-
-      if constexpr (std::is_base_of_v<InteractionProcess<T1type>, T1type> || t1ProcSeq ||
-                    t1SwitchProc) {
-        tot += A.GetInverseInteractionLength(vP);
-      }
-      if constexpr (std::is_base_of_v<InteractionProcess<T2type>, T2type> || t2ProcSeq ||
-                    t2SwitchProc) {
-        tot += B.GetInverseInteractionLength(vP);
-      }
-      return tot;
-    }
-
-    template <typename TParticle, typename TSecondaries>
-    EProcessReturn SelectInteraction(
-        TParticle& vP, TSecondaries& vS,
-        [[maybe_unused]] corsika::units::si::InverseGrammageType lambda_select,
-        corsika::units::si::InverseGrammageType& lambda_inv_count) {
-
-      if constexpr (t1ProcSeq || t1SwitchProc) {
-        // if A is a process sequence --> check inside
-        const EProcessReturn ret =
-            A.SelectInteraction(vP, vS, lambda_select, lambda_inv_count);
-        // if A did succeed, stop routine
-        if (ret != EProcessReturn::eOk) { return ret; }
-      } else if constexpr (std::is_base_of_v<InteractionProcess<T1type>, T1type>) {
-        // if this is not a ContinuousProcess --> evaluate probability
-        lambda_inv_count += A.GetInverseInteractionLength(vP);
-        // check if we should execute THIS process and then EXIT
-        if (lambda_select < lambda_inv_count) {
-          A.DoInteraction(vS);
-          return EProcessReturn::eInteracted;
-        }
-      } // end branch A
-
-      if constexpr (t2ProcSeq || t2SwitchProc) {
-        // if A is a process sequence --> check inside
-        const EProcessReturn ret =
-            B.SelectInteraction(vP, vS, lambda_select, lambda_inv_count);
-        // if A did succeed, stop routine
-        if (ret != EProcessReturn::eOk) { return ret; }
-      } else if constexpr (std::is_base_of_v<InteractionProcess<T2type>, T2type>) {
-        // if this is not a ContinuousProcess --> evaluate probability
-        lambda_inv_count += B.GetInverseInteractionLength(vP);
-        // check if we should execute THIS process and then EXIT
-        if (lambda_select < lambda_inv_count) {
-          B.DoInteraction(vS);
-          return EProcessReturn::eInteracted;
-        }
-      } // end branch A
-      return EProcessReturn::eOk;
-    }
-
-    template <typename TParticle>
-    corsika::units::si::TimeType GetTotalLifetime(TParticle& p) {
-      return 1. / GetInverseLifetime(p);
-    }
-
-    template <typename TParticle>
-    corsika::units::si::InverseTimeType GetTotalInverseLifetime(TParticle& p) {
-      return GetInverseLifetime(p);
-    }
-
-    template <typename TParticle>
-    corsika::units::si::InverseTimeType GetInverseLifetime(TParticle& p) {
-      using namespace corsika::units::si;
-
-      corsika::units::si::InverseTimeType tot = 0 / second;
-
-      if constexpr (std::is_base_of_v<DecayProcess<T1type>, T1type> || t1ProcSeq) {
-        tot += A.GetInverseLifetime(p);
-      }
-      if constexpr (std::is_base_of_v<DecayProcess<T2type>, T2type> || t2ProcSeq) {
-        tot += B.GetInverseLifetime(p);
-      }
-      return tot;
-    }
-
-    // select decay process
-    template <typename TParticle, typename TSecondaries>
-    EProcessReturn SelectDecay(
-        TParticle& vP, TSecondaries& vS,
-        [[maybe_unused]] corsika::units::si::InverseTimeType decay_select,
-        corsika::units::si::InverseTimeType& decay_inv_count) {
-      if constexpr (t1ProcSeq) {
-        // if A is a process sequence --> check inside
-        const EProcessReturn ret = A.SelectDecay(vP, vS, decay_select, decay_inv_count);
-        // if A did succeed, stop routine
-        if (ret != EProcessReturn::eOk) { return ret; }
-      } else if constexpr (std::is_base_of_v<DecayProcess<T1type>, T1type>) {
-        // if this is not a ContinuousProcess --> evaluate probability
-        decay_inv_count += A.GetInverseLifetime(vP);
-        // check if we should execute THIS process and then EXIT
-        if (decay_select < decay_inv_count) { // more pedagogical: rndm_select <
-                                              // decay_inv_count / decay_inv_tot
-          A.DoDecay(vS);
-          return EProcessReturn::eDecayed;
-        }
-      } // end branch A
-
-      if constexpr (t2ProcSeq) {
-        // if A is a process sequence --> check inside
-        const EProcessReturn ret = B.SelectDecay(vP, vS, decay_select, decay_inv_count);
-        // if A did succeed, stop routine
-        if (ret != EProcessReturn::eOk) { return ret; }
-      } else if constexpr (std::is_base_of_v<DecayProcess<T2type>, T2type>) {
-        // if this is not a ContinuousProcess --> evaluate probability
-        decay_inv_count += B.GetInverseLifetime(vP);
-        // check if we should execute THIS process and then EXIT
-        if (decay_select < decay_inv_count) {
-          B.DoDecay(vS);
-          return EProcessReturn::eDecayed;
-        }
-      } // end branch B
-      return EProcessReturn::eOk;
-    }
-  };
-
-  // the << operator assembles many BaseProcess, ContinuousProcess, and
-  // Interaction/DecayProcess objects into a ProcessSequence, all combinatorics
-  // must be allowed, this is why we define a macro to define all
-  // combinations here:
-
-  // enable the << operator to construct ProcessSequence from two
-  // Processes, only if poth Processes derive from BaseProcesses
-
-  template <typename TProcess1, typename TProcess2>
-  inline typename std::enable_if<
-      std::is_base_of<BaseProcess<typename std::decay<TProcess1>::type>,
-                      typename std::decay<TProcess1>::type>::value &&
-          std::is_base_of<BaseProcess<typename std::decay<TProcess2>::type>,
-                          typename std::decay<TProcess2>::type>::value,
-      ProcessSequence<TProcess1, TProcess2>>::type
-  operator<<(TProcess1&& vA, TProcess2&& vB) {
-    return ProcessSequence<TProcess1, TProcess2>(vA, vB);
-  }
-
-  /// marker to identify objectas ProcessSequence
-  template <typename TProcess1, typename TProcess2>
-  struct is_process_sequence<corsika::process::ProcessSequence<TProcess1, TProcess2>>
-      : std::true_type {};
-
-} // namespace corsika::process

Framework/ProcessSequence/ProcessSignature.h

-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#pragma once
-
-#define FORCE_SIGNATURE(nameTrait, nameMethod, signatureMethod)                \
-  template <typename U>                                                        \
-  class nameTrait {                                                            \
-  private:                                                                     \
-    template <typename T, T>                                                   \
-    struct helper;                                                             \
-    template <typename T>                                                      \
-    static std::uint8_t check(helper<signatureMethod, &nameMethod>*);          \
-    template <typename T>                                                      \
-    static std::uint16_t check(...);                                           \
-                                                                               \
-  public:                                                                      \
-    static constexpr bool value = sizeof(check<U>(0)) == sizeof(std::uint8_t); \
-  }
-
-// FORCE_SIGNATURE(thisMustBeDefined, T::thisMustBeDefined, int(*)(void));

Framework/ProcessSequence/SecondariesProcess.h

-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#pragma once
-
-#include <corsika/process/BaseProcess.h>
-#include <corsika/process/ProcessReturn.h> // for convenience
-#include <corsika/setup/SetupTrajectory.h>
-#include <corsika/units/PhysicalUnits.h>
-
-namespace corsika::process {
-
-  /**
-     \class SecondariesProcess
-
-     The structural base type of a process object in a
-     ProcessSequence. Both, the ProcessSequence and all its elements
-     are of type SecondariesProcess<T>
-
-   */
-
-  template <typename TDerived>
-  struct SecondariesProcess : public BaseProcess<TDerived> {
-
-    /// here starts the interface-definition part
-    // -> enforce TDerived to implement DoSecondaries...
-    template <typename TSecondaries>
-    inline EProcessReturn DoSecondaries(TSecondaries&);
-  };
-
-} // namespace corsika::process

Framework/ProcessSequence/StackProcess.h

-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#pragma once
-
-#include <corsika/process/BaseProcess.h>
-#include <corsika/process/ProcessReturn.h> // for convenience
-#include <corsika/setup/SetupTrajectory.h>
-#include <corsika/units/PhysicalUnits.h>
-
-namespace corsika::process {
-
-  /**
-     \class StackProcess
-
-     The structural base type of a process object in a
-     ProcessSequence. Both, the ProcessSequence and all its elements
-     are of type StackProcess<T>
-
-   */
-
-  template <typename TDerived>
-  class StackProcess : public BaseProcess<TDerived> {
-
-  public:
-    StackProcess() = delete;
-    StackProcess(const unsigned int nStep)
-        : fNStep(nStep) {}
-
-    /// here starts the interface-definition part
-    // -> enforce TDerived to implement DoStack...
-    template <typename TStack>
-    inline EProcessReturn DoStack(TStack&);
-
-    int GetStep() const { return fIStep; }
-    bool CheckStep() { return !((++fIStep) % fNStep); }
-
-  private:
-    /**
-       @name The number of "steps" during the cascade processing after
-       which this StackProcess is going to be executed. The logic is
-       "fIStep modulo fNStep"
-       @{
-     */
-    unsigned int fNStep = 0;
-    unsigned long int fIStep = 0;
-    //! @}
-  };
-
-  // overwrite the default trait class, to mark BaseProcess<T> as useful process
-  template <class T>
-  std::true_type is_process_impl(const StackProcess<T>* impl);
-
-} // namespace corsika::process

Framework/ProcessSequence/testProcessSequence.cc

-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#include <catch2/catch.hpp>
-
-#include <array>
-#include <iomanip>
-#include <iostream>
-
-#include <corsika/process/ProcessSequence.h>
-#include <corsika/process/switch_process/SwitchProcess.h>
-
-using namespace corsika;
-using namespace corsika::units::si;
-using namespace corsika::process;
-using namespace std;
-
-static const int nData = 10;
-
-int globalCount = 0;
-
-class ContinuousProcess1 : public ContinuousProcess<ContinuousProcess1> {
-  int fV = 0;
-
-public:
-  ContinuousProcess1(const int v)
-      : fV(v) {
-
-    cout << "globalCount: " << globalCount << ", fV: " << fV << std::endl;
-    globalCount++;
-  }
-
-  template <typename D, typename T>
-  inline EProcessReturn DoContinuous(D& d, T&) const {
-    cout << "ContinuousProcess1::DoContinuous" << endl;
-    for (int i = 0; i < nData; ++i) d.p[i] += 0.933;
-    return EProcessReturn::eOk;
-  }
-};
-
-class ContinuousProcess2 : public ContinuousProcess<ContinuousProcess2> {
-  int fV = 0;
-
-public:
-  ContinuousProcess2(const int v)
-      : fV(v) {
-    cout << "globalCount: " << globalCount << ", fV: " << fV << std::endl;
-    globalCount++;
-  }
-
-  template <typename D, typename T>
-  inline EProcessReturn DoContinuous(D& d, T&) const {
-    cout << "ContinuousProcess2::DoContinuous" << endl;
-    for (int i = 0; i < nData; ++i) d.p[i] += 0.933;
-    return EProcessReturn::eOk;
-  }
-};
-
-class Process1 : public InteractionProcess<Process1> {
-public:
-  Process1(const int v)
-      : fV(v) {
-    cout << "globalCount: " << globalCount << ", fV: " << fV << std::endl;
-    globalCount++;
-  }
-
-  template <typename D, typename S>
-  inline EProcessReturn DoInteraction(D& d, S&) const {
-    for (int i = 0; i < nData; ++i) d.p[i] += 1 + i;
-    return EProcessReturn::eOk;
-  }
-
-private:
-  int fV;
-};
-
-class Process2 : public InteractionProcess<Process2> {
-  int fV = 0;
-
-public:
-  Process2(const int v)
-      : fV(v) {
-    cout << "globalCount: " << globalCount << ", fV: " << fV << std::endl;
-    globalCount++;
-  }
-
-  template <typename Particle>
-  inline EProcessReturn DoInteraction(Particle&) const {
-    cout << "Process2::DoInteraction" << endl;
-    return EProcessReturn::eOk;
-  }
-  template <typename Particle>
-  GrammageType GetInteractionLength(Particle&) const {
-    cout << "Process2::GetInteractionLength" << endl;
-    return 3_g / (1_cm * 1_cm);
-  }
-};
-
-class Process3 : public InteractionProcess<Process3> {
-  int fV = 0;
-
-public:
-  Process3(const int v)
-      : fV(v) {
-    cout << "globalCount: " << globalCount << ", fV: " << fV << std::endl;
-    globalCount++;
-  }
-
-  template <typename Particle>
-  inline EProcessReturn DoInteraction(Particle&) const {
-    cout << "Process3::DoInteraction" << endl;
-    return EProcessReturn::eOk;
-  }
-  template <typename Particle>
-  GrammageType GetInteractionLength(Particle&) const {
-    cout << "Process3::GetInteractionLength" << endl;
-    return 1_g / (1_cm * 1_cm);
-  }
-};
-
-class Process4 : public BaseProcess<Process4> {
-  int fV = 0;
-
-public:
-  Process4(const int v)
-      : fV(v) {
-    cout << "globalCount: " << globalCount << ", fV: " << fV << std::endl;
-    globalCount++;
-  }
-
-  template <typename D, typename T>
-  inline EProcessReturn DoContinuous(D& d, T&) const {
-    for (int i = 0; i < nData; ++i) { d.p[i] /= 1.2; }
-    return EProcessReturn::eOk;
-  }
-  template <typename Particle>
-  EProcessReturn DoInteraction(Particle&) const {
-    return EProcessReturn::eOk;
-  }
-};
-
-class Decay1 : public DecayProcess<Decay1> {
-
-public:
-  Decay1(const int v) {
-    cout << "Decay1()" << endl;
-    globalCount++;
-  }
-
-  template <typename Particle>
-  TimeType GetLifetime(Particle&) const {
-    return 1_s;
-  }
-  template <typename Particle>
-  EProcessReturn DoDecay(Particle&) const {
-    return EProcessReturn::eOk;
-  }
-};
-
-class Stack1 : public StackProcess<Stack1> {
-  int fCount = 0;
-
-public:
-  Stack1(const int n)
-      : StackProcess(n) {}
-  template <typename TStack>
-  EProcessReturn DoStack(TStack&) {
-    fCount++;
-    return EProcessReturn::eOk;
-  }
-  int GetCount() const { return fCount; }
-};
-
-struct DummyStack {};
-struct DummyData {
-  double p[nData] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
-};
-struct DummyTrajectory {};
-
-TEST_CASE("Process Sequence", "[Process Sequence]") {
-
-  SECTION("Check construction") {
-    globalCount = 0;
-    Process1 m1(0);
-    CHECK(globalCount == 1);
-    Process2 m2(1);
-    CHECK(globalCount == 2);
-    Process3 m3(2);
-    CHECK(globalCount == 3);
-    Process4 m4(3);
-    CHECK(globalCount == 4);
-
-    [[maybe_unused]] auto sequence = m1 << m2 << m3 << m4;
-  }
-
-  SECTION("interaction length") {
-    globalCount = 0;
-    ContinuousProcess1 cp1(0);
-    Process2 m2(1);
-    Process3 m3(2);
-
-    DummyData particle;
-
-    auto sequence2 = cp1 << m2 << m3;
-    GrammageType const tot = sequence2.GetTotalInteractionLength(particle);
-    InverseGrammageType const tot_inv =
-        sequence2.GetTotalInverseInteractionLength(particle);
-    cout << "lambda_tot=" << tot << "; lambda_tot_inv=" << tot_inv << endl;
-    globalCount = 0;
-  }
-
-  SECTION("lifetime") {
-    globalCount = 0;
-    ContinuousProcess1 cp1(0);
-    Process2 m2(1);
-    Process3 m3(2);
-    Decay1 d3(3);
-
-    DummyData particle;
-
-    auto sequence2 = cp1 << m2 << m3 << d3;
-    TimeType const tot = sequence2.GetTotalLifetime(particle);
-    InverseTimeType const tot_inv = sequence2.GetTotalInverseLifetime(particle);
-    cout << "lambda_tot=" << tot << "; lambda_tot_inv=" << tot_inv << endl;
-
-    globalCount = 0;
-  }
-
-  SECTION("sectionTwo") {
-    globalCount = 0;
-    ContinuousProcess1 cp1(0);
-    ContinuousProcess2 cp2(1);
-    Process2 m2(2);
-    Process3 m3(3);
-
-    auto sequence2 = cp1 << m2 << m3 << cp2;
-
-    DummyData particle;
-    DummyTrajectory track;
-
-    cout << "-->init sequence2" << endl;
-    globalCount = 0;
-    cout << "-->docont" << endl;
-
-    sequence2.DoContinuous(particle, track);
-    cout << "-->dodisc" << endl;
-    cout << "-->done" << endl;
-
-    const int nLoop = 5;
-    cout << "Running loop with n=" << nLoop << endl;
-    for (int i = 0; i < nLoop; ++i) { sequence2.DoContinuous(particle, track); }
-    for (int i = 0; i < nData; i++) {
-      cout << "data[" << i << "]=" << particle.p[i] << endl;
-    }
-    cout << "done" << endl;
-  }
-
-  SECTION("StackProcess") {
-
-    globalCount = 0;
-    Stack1 s1(1);
-    Stack1 s2(2);
-
-    auto sequence = s1 << s2;
-
-    DummyStack stack;
-
-    const int nLoop = 20;
-    for (int i = 0; i < nLoop; ++i) { sequence.DoStack(stack); }
-
-    CHECK(s1.GetCount() == 20);
-    CHECK(s2.GetCount() == 10);
-  }
-}
-
-/*
-  Note: there is a fine-grained dedicated test-suite for SwitchProcess
-  in Processes/SwitchProcess/testSwtichProcess
- */
-TEST_CASE("SwitchProcess") {
-  globalCount = 0;
-  Process1 p1(0);
-  Process2 p2(1);
-  switch_process::SwitchProcess s(p1, p2, 10_GeV);
-  CHECK(is_switch_process_v<decltype(s)>);
-}

Framework/Random/CMakeLists.txt

-set (
-  CORSIKArandom_SOURCES
-  RNGManager.cc
-  )
-
-set (
-  CORSIKArandom_HEADERS
-  RNGManager.h
-  UniformRealDistribution.h
-  ExponentialDistribution.h
-  )
-
-set (
-  CORSIKArandom_NAMESPACE
-  corsika/random
-  )
-
-add_library (CORSIKArandom STATIC ${CORSIKArandom_SOURCES})
-CORSIKA_COPY_HEADERS_TO_NAMESPACE (CORSIKArandom ${CORSIKArandom_NAMESPACE} ${CORSIKArandom_HEADERS})
-
-target_link_libraries (
-  CORSIKArandom
-  INTERFACE
-  CORSIKAutilities
-  CORSIKAunits
-  )
-
-target_include_directories (
-  CORSIKArandom
-  PUBLIC
-  $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
-  $<INSTALL_INTERFACE:include/>
-  )
-
-# target dependencies on other libraries (also the header onlys)
-# none
-
-install (
-  TARGETS CORSIKArandom
-  LIBRARY DESTINATION lib
-  ARCHIVE DESTINATION lib
-  PUBLIC_HEADER DESTINATION include/${CORSIKArandom_NAMESPACE}
-  )
-
-
-# --------------------
-# code unit testing
-CORSIKA_ADD_TEST(testRandom)
-target_link_libraries (
-  testRandom
-  CORSIKArandom
-  CORSIKAtesting
-  )

Framework/Random/ExponentialDistribution.h

-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#pragma once
-
-#include <corsika/units/PhysicalUnits.h>
-#include <random>
-
-namespace corsika::random {
-
-  template <class TQuantity>
-  class ExponentialDistribution {
-    using RealType = typename TQuantity::value_type;
-    std::exponential_distribution<RealType> dist{1.};
-
-    TQuantity const fBeta;
-
-  public:
-    ExponentialDistribution(TQuantity beta)
-        : fBeta(beta) {}
-
-    template <class Generator>
-    TQuantity operator()(Generator& g) {
-      return fBeta * dist(g);
-    }
-  };
-
-} // namespace corsika::random

Framework/Random/RNGManager.cc

-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#include <corsika/random/RNGManager.h>
-#include <sstream>
-
-void corsika::random::RNGManager::RegisterRandomStream(std::string const& pStreamName) {
-  corsika::random::RNG rng;
-
-  if (auto const& it = seeds.find(pStreamName); it != seeds.end()) {
-    rng.seed(it->second);
-  }
-
-  rngs[pStreamName] = std::move(rng);
-}
-
-corsika::random::RNG& corsika::random::RNGManager::GetRandomStream(
-    std::string 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.");
-  }
-}
-
-bool corsika::random::RNGManager::IsRegistered(std::string const& pStreamName) const {
-  return rngs.count(pStreamName) > 0;
-}
-
-std::stringstream corsika::random::RNGManager::dumpState() const {
-  std::stringstream buffer;
-  for (auto const& [streamName, rng] : rngs) {
-    buffer << '"' << streamName << "\" = \"" << rng << '"' << std::endl;
-  }
-
-  return buffer;
-}
-
-void corsika::random::RNGManager::SeedAll(uint64_t vSeed) {
-  for (auto& entry : rngs) { entry.second.seed(vSeed++); }
-}
-
-void corsika::random::RNGManager::SeedAll() {
-  std::random_device rd;
-
-  for (auto& entry : rngs) {
-    std::seed_seq sseq{rd(), rd(), rd(), rd(), rd(), rd()};
-    entry.second.seed(sseq);
-  }
-}

Framework/Random/RNGManager.h

-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#pragma once
-
-#include <corsika/utl/Singleton.h>
-
-#include <map>
-#include <random>
-#include <string>
-
-/*!
- * With this class modules can register streams of random numbers.
- */
-
-namespace corsika::random {
-
-  using RNG = std::mt19937; //!< the actual RNG type that will be used
-
-  /*!
-   * Manage random number generators.
-   */
-  class RNGManager final : public corsika::utl::Singleton<RNGManager> {
-
-    friend class corsika::utl::Singleton<RNGManager>;
-
-    std::map<std::string, RNG> rngs;
-    std::map<std::string, std::seed_seq> seeds;
-
-  protected:
-    RNGManager() {}
-
-  public:
-    /*!
-     * This function is to be called by a module requiring a random-number
-     * stream during its initialization.
-     *
-     * \throws sth. when stream \a pModuleName is already registered
-     */
-    void RegisterRandomStream(std::string const& pStreamName);
-
-    /*!
-     * returns the pre-stored stream of given name \a pStreamName if
-     * available
-     */
-    RNG& GetRandomStream(std::string const& pStreamName);
-
-    /*!
-     * Check whether a stream has been registered.
-     */
-    bool IsRegistered(std::string const& pStreamName) const;
-
-    /*!
-     * dumps the names and states of all registered random-number streams
-     * into a std::stringstream.
-     */
-    std::stringstream dumpState() const;
-
-    /**
-     * Set explicit seeds for all currently registered streams. The actual seed values
-     * are incremented from \a vSeed.
-     */
-    void SeedAll(uint64_t vSeed);
-
-    void SeedAll(); //!< seed all currently registered streams with "real" randomness
-  };
-
-} // namespace corsika::random

Framework/Random/UniformRealDistribution.h

-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#pragma once
-
-#include <corsika/units/PhysicalUnits.h>
-#include <random>
-
-namespace corsika::random {
-
-  template <class TQuantity>
-  class UniformRealDistribution {
-    using RealType = typename TQuantity::value_type;
-    std::uniform_real_distribution<RealType> dist{RealType(0.), RealType(1.)};
-
-    TQuantity const a, b;
-
-  public:
-    UniformRealDistribution(TQuantity b)
-        : a{TQuantity(phys::units::detail::magnitude_tag, 0)}
-        , b(b) {}
-    UniformRealDistribution(TQuantity a, TQuantity b)
-        : a(a)
-        , b(b) {}
-
-    template <class Generator>
-    TQuantity operator()(Generator& g) {
-      return a + dist(g) * (b - a);
-    }
-  };
-
-} // namespace corsika::random

Framework/StackInterface/CMakeLists.txt

-set (
-  CORSIKAstackinterface_HEADERS
-  ParticleBase.h
-  StackIteratorInterface.h
-  Stack.h
-  SecondaryView.h
-  CombinedStack.h
-  )
-
-set (
-  CORSIKAstackinterface_NAMESPACE
-  corsika/stack
-  )
-
-add_library (
-  CORSIKAstackinterface
-  INTERFACE
-  )
-
-CORSIKA_COPY_HEADERS_TO_NAMESPACE (
-  CORSIKAstackinterface ${CORSIKAstackinterface_NAMESPACE} ${CORSIKAstackinterface_HEADERS}
-  )
-
-target_include_directories (
-  CORSIKAstackinterface
-  INTERFACE
-  $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
-  $<INSTALL_INTERFACE:include>
-  )
-
-install (
-  FILES ${CORSIKAstackinterface_HEADERS}
-  DESTINATION include/${CORSIKAstackinterface_NAMESPACE}
-  )
-
-#code testing
-CORSIKA_ADD_TEST(testStackInterface)
-target_link_libraries (testStackInterface CORSIKAstackinterface CORSIKAtesting)
-
-CORSIKA_ADD_TEST(testSecondaryView)
-target_link_libraries (testSecondaryView CORSIKAstackinterface CORSIKAtesting)
-
-CORSIKA_ADD_TEST(testCombinedStack)
-target_link_libraries (testCombinedStack CORSIKAstackinterface CORSIKAtesting)

Framework/StackInterface/SecondaryView.h

-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#pragma once
-
-#include <corsika/stack/Stack.h>
-
-#include <stdexcept>
-#include <vector>
-
-namespace corsika::stack {
-
-  /**
-   * @class SecondaryView
-   *
-   * SecondaryView can only be constructed by giving a valid
-   * Projectile particle, following calls to AddSecondary will
-   * populate the original Stack, but will be directly accessible via
-   * the SecondaryView, e.g.
-
-     This allows to write code like
-     \verbatim
-     auto projectileInput = mainStack.GetNextParticle();
-     const unsigned int nMain = mainStack.GetSize();
-     SecondaryView<StackData, ParticleInterface> mainStackView(projectileInput);
-     mainStackView.AddSecondary(...data...);
-     mainStackView.AddSecondary(...data...);
-     mainStackView.AddSecondary(...data...);
-     mainStackView.AddSecondary(...data...);
-     assert(mainStackView.GetSize() == 4);
-     assert(mainStack.GetSize() = nMain+4);
-     \endverbatim
-
-     All operations possible on a Stack object are also possible on a
-     SecondaryView object. This means you can add, delete, copy, swap,
-     iterate, etc.
-
-     *Further information about implementation (for developers):* All
-     data is stored in the original stack privided at construction
-     time. The secondary particle (view) indices are stored in an
-     extra std::vector of SecondaryView class 'fIndices' referring to
-     the original stack slot indices. The index of the primary
-     projectle particle is also explicitly stored in
-     'fProjectileIndex'. StackIterator indices
-     'i = StackIterator::GetIndex()' are referring to those numbers,
-     where 'i==0' refers to the 'fProjectileIndex', and
-     'StackIterator::GetIndex()>0' to 'fIndices[i-1]', see function
-     GetIndexFromIterator.
-   */
-
-  template <typename StackDataType, template <typename> typename ParticleInterface>
-  class SecondaryView : public Stack<StackDataType&, ParticleInterface> {
-
-    using ViewType = SecondaryView<StackDataType, ParticleInterface>;
-
-  private:
-    /**
-     * Helper type for inside this class
-     */
-    using InnerStackType = Stack<StackDataType&, ParticleInterface>;
-
-    /**
-     * @name We need this "special" types with non-reference StackData for
-     * the constructor of the SecondaryView class
-     * @{
-     */
-    using InnerStackTypeValue = Stack<StackDataType, ParticleInterface>;
-    using StackIteratorValue =
-        StackIteratorInterface<typename std::remove_reference<StackDataType>::type,
-                               ParticleInterface, InnerStackTypeValue>;
-    /// @}
-
-  public:
-    using StackIterator =
-        StackIteratorInterface<typename std::remove_reference<StackDataType>::type,
-                               ParticleInterface, ViewType>;
-    using ConstStackIterator =
-        ConstStackIteratorInterface<typename std::remove_reference<StackDataType>::type,
-                                    ParticleInterface, ViewType>;
-
-    /**
-     * this is the full type of the declared ParticleInterface: typedef typename
-     */
-    using ParticleType = StackIterator;
-    using ParticleInterfaceType = typename StackIterator::ParticleInterfaceType;
-
-    friend class StackIteratorInterface<
-        typename std::remove_reference<StackDataType>::type, ParticleInterface, ViewType>;
-
-    friend class ConstStackIteratorInterface<
-        typename std::remove_reference<StackDataType>::type, ParticleInterface, ViewType>;
-
-  private:
-    /**
-     * This is not accessible, since we don't want to allow creating a
-     * new stack.
-     */
-    template <typename... Args>
-    SecondaryView(Args... args) = delete;
-
-  public:
-    /**
-       SecondaryView can only be constructed passing it a valid
-       StackIterator to another Stack object
-     **/
-    SecondaryView(StackIteratorValue& vI)
-        : Stack<StackDataType&, ParticleInterface>(vI.GetStackData())
-        , fProjectileIndex(vI.GetIndex()) {}
-
-    StackIterator GetProjectile() {
-      // NOTE: 0 is special marker here for PROJECTILE, see GetIndexFromIterator
-      return StackIterator(*this, 0);
-    }
-
-    template <typename... Args>
-    auto AddSecondary(const Args... v) {
-      StackIterator proj = GetProjectile();
-      return AddSecondary(proj, v...);
-    }
-
-    template <typename... Args>
-    auto AddSecondary(StackIterator& proj, const Args... v) {
-      // make space on stack
-      InnerStackType::GetStackData().IncrementSize();
-      // get current number of secondaries on stack
-      const unsigned int idSec = GetSize();
-      // determine index on (inner) stack where new particle will be located
-      const unsigned int index = InnerStackType::GetStackData().GetSize() - 1;
-      fIndices.push_back(index);
-      // NOTE: "+1" is since "0" is special marker here for PROJECTILE, see
-      // GetIndexFromIterator
-      return StackIterator(*this, idSec + 1, proj, v...);
-    }
-
-    /**
-     * overwrite Stack::GetSize to return actual number of secondaries
-     */
-    unsigned int GetSize() const { return fIndices.size(); }
-
-    /**
-     * @name These are functions required by std containers and std loops
-     * The Stack-versions must be overwritten, since here we need the correct
-     * SecondaryView::GetSize
-     * @{
-     */
-    // NOTE: the "+1" is since "0" is special marker here for PROJECTILE, see
-    // GetIndexFromIterator
-    auto begin() { return StackIterator(*this, 0 + 1); }
-    auto end() { return StackIterator(*this, GetSize() + 1); }
-    auto last() { return StackIterator(*this, GetSize() - 1 + 1); }
-
-    auto begin() const { return ConstStackIterator(*this, 0 + 1); }
-    auto end() const { return ConstStackIterator(*this, GetSize() + 1); }
-    auto last() const { return ConstStackIterator(*this, GetSize() - 1 + 1); }
-
-    auto cbegin() const { return ConstStackIterator(*this, 0 + 1); }
-    auto cend() const { return ConstStackIterator(*this, GetSize() + 1); }
-    auto clast() const { return ConstStackIterator(*this, GetSize() - 1 + 1); }
-    /// @}
-
-    /**
-     * need overwrite Stack::Delete, since we want to call
-     * SecondaryView::DeleteLast
-     *
-     * The particle is deleted on the underlying (internal) stack. The
-     * local references in SecondaryView in fIndices must be fixed,
-     * too.  The approach is to a) check if the particle 'p' is at the
-     * very end of the internal stack, b) if not: move it there by
-     * copying the last particle to the current particle location, c)
-     * remove the last particle.
-     *
-     */
-    void Delete(StackIterator p) {
-      if (IsEmpty()) { /* error */
-        throw std::runtime_error("Stack, cannot delete entry since size is zero");
-      }
-      const int innerSize = InnerStackType::GetSize();
-      const int innerIndex = GetIndexFromIterator(p.GetIndex());
-      if (innerIndex < innerSize - 1)
-        InnerStackType::GetStackData().Copy(innerSize - 1,
-                                            GetIndexFromIterator(p.GetIndex()));
-      DeleteLast();
-    }
-
-    /**
-     * need overwrite Stack::Delete, since we want to call SecondaryView::DeleteLast
-     */
-    void Delete(ParticleInterfaceType p) { Delete(p.GetIterator()); }
-
-    /**
-     * delete last particle on stack by decrementing stack size
-     */
-    void DeleteLast() {
-      fIndices.pop_back();
-      InnerStackType::GetStackData().DecrementSize();
-    }
-
-    /**
-     * return next particle from stack, need to overwrtie Stack::GetNextParticle to get
-     * right reference
-     */
-    StackIterator GetNextParticle() { return last(); }
-
-    /**
-     * check if there are no further particles on stack
-     */
-    bool IsEmpty() { return GetSize() == 0; }
-
-  protected:
-    /**
-     * We only want to 'see' secondaries indexed in fIndices. In this
-     * function the conversion form iterator-index to stack-index is
-     * performed.
-     */
-    unsigned int GetIndexFromIterator(const unsigned int vI) const {
-      if (vI == 0) return fProjectileIndex;
-      return fIndices[vI - 1];
-    }
-
-  private:
-    unsigned int fProjectileIndex;
-    std::vector<unsigned int> fIndices;
-  };
-
-  /*
-    See Issue 161
-
-    unfortunately clang does not support this in the same way (yet) as
-    gcc, so we have to distinguish here. If clang cataches up, we
-    could remove the #if here and elsewhere. The gcc code is much more
-    generic and universal.
-  */
-#if not defined(__clang__) && defined(__GNUC__) || defined(__GNUG__)
-  template <typename S, template <typename> typename _PIType = S::template PIType>
-  struct MakeView {
-    using type = corsika::stack::SecondaryView<typename S::StackImpl, _PIType>;
-  };
-#endif
-
-} // namespace corsika::stack

Framework/StackInterface/Stack.dox

-/**
-  @page Stack Description of particle stacks
-
-  In the CORSIKA 8 framework particle data is always stored in
-  particle stacks. A particle is, thus, always a reference (iterator)
-  to an entry on a stack, e.g.
-
-  \verbatim
-  ModelStack stack;
-  stack.begin(); // returns an iterator: StackIterator, ConstStackIterator
-
-  *stack.begin(); // return a reference to ParticleInterfaceType, which is the class provided by the user to read/write particle properties
-  
-  \endverbatim
-
-  All functionality and algorithms for stack data access is located in the namespace corsika::stack
-
-  The minimal example of how to use this is shown in stack_example.cc  
-
-*/
\ No newline at end of file

Framework/StackInterface/Stack.h

-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#pragma once
-
-#include <corsika/stack/StackIteratorInterface.h>
-// must be after StackIteratorInterface
-#include <corsika/stack/SecondaryView.h>
-#include <corsika/utl/MetaProgramming.h>
-
-#include <stdexcept>
-#include <type_traits>
-
-/**
-   All classes around management of particles on a stack.
- */
-
-namespace corsika::stack {
-
-  /**
-     This is just a forward declatation for the user-defined
-     ParticleInterface, which is one of the essential template
-     parameters for the Stack.
-
-     <b>Important:</b> ParticleInterface must inherit from ParticleBase !
-   */
-
-  template <typename>
-  class ParticleInterface;
-
-  /**
-     The Stack class provides (and connects) the main particle data storage machinery.
-
-     The StackDataType type is the user-provided bare data storage
-     object. This can be of any complexity, from a simple struct
-     (fortran common block), to a combination of different and
-     distributed data sources.
-
-     The user-provided ParticleInterface template type is the base
-     class type of the StackIteratorInterface class (CRTP) and must
-     provide all functions to read single particle data from the
-     StackDataType, given an 'unsigned int' index.
-
-     The Stack implements the
-     std-type begin/end function to allow integration in normal for
-     loops, ranges, etc.
-   */
-
-  template <typename StackDataType, template <typename> typename ParticleInterface>
-  class Stack {
-    using StackDataValueType = std::remove_reference_t<StackDataType>;
-
-    StackDataType fData; ///< this in general holds all the data and can be quite big
-
-  private:
-    Stack(Stack&) = delete; ///< since Stack can be very big, we don't want to copy it
-    Stack& operator=(Stack&) =
-        delete; ///< since Stack can be very big, we don't want to copy it
-
-  public:
-    /**
-     * if StackDataType is a reference member we *HAVE* to initialize
-     * it in the constructor, this is typically needed for SecondaryView
-     */
-    template <typename _ = StackDataType, typename = utl::enable_if<std::is_reference<_>>>
-    Stack(StackDataType vD)
-        : fData(vD) {}
-
-    /**
-     * This constructor takes any argument and passes it on to the
-     * StackDataType user class. If the user did not provide a suited
-     * constructor this will fail with an error message.
-     *
-     * Furthermore, this is disabled with enable_if for SecondaryView
-     * stacks, where the inner data container is always a reference
-     * and cannot be initialized here.
-     */
-    template <typename... Args, typename _ = StackDataType,
-              typename = utl::disable_if<std::is_reference<_>>>
-    Stack(Args... args)
-        : fData(args...) {}
-
-  public:
-    typedef StackDataType
-        StackImpl; ///< this is the type of the user-provided data structure
-
-    template <typename SI>
-    using PIType = ParticleInterface<SI>;
-
-    /**
-     * Via the StackIteratorInterface and ConstStackIteratorInterface
-     * specialization, the type of the StackIterator
-     * template class is declared for a particular stack data
-     * object. Using CRTP, this also determines the type of
-     * ParticleInterface template class simultaneously.
-     */
-    using StackIterator =
-        StackIteratorInterface<StackDataValueType, ParticleInterface, Stack>;
-    using ConstStackIterator =
-        ConstStackIteratorInterface<StackDataValueType, ParticleInterface, Stack>;
-
-    /**
-     * this is the full type of the user-declared ParticleInterface
-     */
-    using ParticleInterfaceType = typename StackIterator::ParticleInterfaceType;
-    /**
-     * In all programming context, the object to access, copy, and
-     * transport particle data is via the StackIterator
-     */
-    using ParticleType = StackIterator;
-
-    // friends are needed since they need access to protected members
-    friend class StackIteratorInterface<StackDataValueType, ParticleInterface, Stack>;
-    friend class ConstStackIteratorInterface<StackDataValueType, ParticleInterface,
-                                             Stack>;
-
-  public:
-    /**
-     * @name Most generic proxy methods for StackDataType fData
-     * @{
-     */
-    unsigned int GetCapacity() const { return fData.GetCapacity(); }
-    unsigned int GetSize() const { return fData.GetSize(); }
-
-    template <typename... Args>
-    auto Clear(Args... args) {
-      return fData.Clear(args...);
-    }
-    ///@}
-
-  public:
-    /**
-     * @name These are functions required by std containers and std loops
-     * @{
-     */
-    StackIterator begin() { return StackIterator(*this, 0); }
-    StackIterator end() { return StackIterator(*this, GetSize()); }
-    StackIterator last() { return StackIterator(*this, GetSize() - 1); }
-
-    ConstStackIterator begin() const { return ConstStackIterator(*this, 0); }
-    ConstStackIterator end() const { return ConstStackIterator(*this, GetSize()); }
-    ConstStackIterator last() const { return ConstStackIterator(*this, GetSize() - 1); }
-
-    ConstStackIterator cbegin() const { return ConstStackIterator(*this, 0); }
-    ConstStackIterator cend() const { return ConstStackIterator(*this, GetSize()); }
-    ConstStackIterator clast() const { return ConstStackIterator(*this, GetSize() - 1); }
-    /// @}
-
-    /**
-     * increase stack size, create new particle at end of stack
-     */
-    template <typename... Args>
-    StackIterator AddParticle(const Args... v) {
-      fData.IncrementSize();
-      return StackIterator(*this, GetSize() - 1, v...);
-    }
-
-    /**
-     * increase stack size, create new particle at end of stack, related to parent
-     * particle/projectile
-     */
-    template <typename... Args>
-    StackIterator AddSecondary(StackIterator& parent, const Args... v) {
-      fData.IncrementSize();
-      return StackIterator(*this, GetSize() - 1, parent, v...);
-    }
-
-    void Swap(StackIterator a, StackIterator b) {
-      fData.Swap(a.GetIndex(), b.GetIndex());
-    }
-    void Swap(ConstStackIterator a, ConstStackIterator b) {
-      fData.Swap(a.GetIndex(), b.GetIndex());
-    }
-    void Copy(StackIterator a, StackIterator b) {
-      fData.Copy(a.GetIndex(), b.GetIndex());
-    }
-    void Copy(ConstStackIterator a, StackIterator b) {
-      fData.Copy(a.GetIndex(), b.GetIndex());
-    }
-
-    /**
-     * delete this particle
-     */
-    void Delete(StackIterator p) {
-      if (GetSize() == 0) { /*error*/
-        throw std::runtime_error("Stack, cannot delete entry since size is zero");
-      }
-      if (p.GetIndex() < GetSize() - 1) fData.Copy(GetSize() - 1, p.GetIndex());
-      DeleteLast();
-    }
-    /**
-     * delete this particle
-     */
-    void Delete(ParticleInterfaceType p) { Delete(p.GetIterator()); }
-
-    /**
-     * delete last particle on stack by decrementing stack size
-     */
-    void DeleteLast() { fData.DecrementSize(); }
-
-    /**
-     * check if there are no further particles on stack
-     */
-    bool IsEmpty() { return GetSize() == 0; }
-
-    /**
-     * return next particle from stack
-     */
-    StackIterator GetNextParticle() { return last(); }
-
-  protected:
-    /**
-     * Function to perform eventual transformation from
-     * StackIterator::GetIndex() to index in data stored in
-     * StackDataType fData. By default (and in almost all cases) this
-     * should just be identiy. See class SecondaryView for an alternative implementation.
-     */
-    unsigned int GetIndexFromIterator(const unsigned int vI) const { return vI; }
-
-    /**
-     * @name Return reference to StackDataType object fData for data access
-     * @{
-     */
-    StackDataValueType& GetStackData() { return fData; }
-    const StackDataValueType& GetStackData() const { return fData; }
-    ///@}
-  };
-
-} // namespace corsika::stack

Framework/StackInterface/StackIteratorInterface.h

-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#pragma once
-
-#include <corsika/stack/ParticleBase.h>
-
-namespace corsika::stack {
-
-  template <typename StackDataType, template <typename> typename ParticleInterface>
-  class Stack; // forward decl
-
-  template <typename StackDataType, template <typename> typename ParticleInterface>
-  class SecondaryView; // forward decl
-
-  /**
-     @class StackIteratorInterface
-
-     The StackIteratorInterface is the main interface to iterator over
-     particles on a stack. At the same time StackIteratorInterface is a
-     Particle object by itself, thus there is no difference between
-     type and ref_type for convenience of the physicist.
-
-     This allows to write code like
-     \verbatim
-     for (auto& p : theStack) { p.SetEnergy(newEnergy); }
-     \endverbatim
-
-     The template argument Stack determines the type of Stack object
-     the data is stored in. A pointer to the Stack object is part of
-     the StackIteratorInterface. In addition to Stack the iterator only knows
-     the index fIndex in the Stack data.
-
-     The template argument `ParticleInterface` acts as a policy to provide
-     readout function of Particle data from the stack. The ParticleInterface
-     class must know how to retrieve information from the Stack data
-     for a particle entry at any index fIndex.
-
-     The ParticleInterface class must be written and provided by the
-     user, it contains methods like <code> auto GetData() const {
-     return GetStackData().GetData(GetIndex()); }</code>, where
-     StackIteratorInterface::GetStackData() return a reference to the
-     object storing the particle data of type StackDataType. And
-     StackIteratorInterface::GetIndex() provides the iterator index to
-     be readout. The StackDataType is another user-provided class to
-     store data and must implement functions compatible with
-     ParticleInterface, in this example StackDataType::GetData(const unsigned int
-     vIndex).
-
-     For two examples see stack_example.cc, or the
-     corsika::processes::sibyll::SibStack class
-  */
-
-  template <typename StackDataType, template <typename> typename ParticleInterface,
-            typename StackType = Stack<StackDataType, ParticleInterface>>
-  class StackIteratorInterface
-      : public ParticleInterface<
-            StackIteratorInterface<StackDataType, ParticleInterface, StackType>> {
-
-  public:
-    using ParticleInterfaceType =
-        ParticleInterface<corsika::stack::StackIteratorInterface<
-            StackDataType, ParticleInterface, StackType>>;
-
-    // friends are needed for access to protected methods
-    friend class Stack<StackDataType,
-                       ParticleInterface>; // for access to GetIndex for Stack
-    friend class Stack<StackDataType&, ParticleInterface>; // for access to GetIndex
-                                                           // SecondaryView : public Stack
-    friend class ParticleBase<StackIteratorInterface>; // for access to GetStackDataType
-    friend class SecondaryView<StackDataType,
-                               ParticleInterface>; // access for SecondaryView
-
-  private:
-    unsigned int fIndex = 0;
-    StackType* fData = 0; // info: Particles and StackIterators become invalid when parent
-                          // Stack is copied or deleted!
-
-    // it is not allowed to create a "dangling" stack iterator
-    StackIteratorInterface() = delete;
-
-  public:
-    StackIteratorInterface(StackIteratorInterface const& vR)
-        : fIndex(vR.fIndex)
-        , fData(vR.fData) {}
-
-    StackIteratorInterface& operator=(StackIteratorInterface const& vR) {
-      fIndex = vR.fIndex;
-      fData = vR.fData;
-      return *this;
-    }
-
-    /** iterator must always point to data, with an index:
-          @param data reference to the stack [rw]
-          @param index index on stack
-       */
-    StackIteratorInterface(StackType& data, const unsigned int index)
-        : fIndex(index)
-        , fData(&data) {}
-
-    /** constructor that also sets new values on particle data object
-        @param data reference to the stack [rw]
-        @param index index on stack
-        @param args variadic list of data to initialize stack entry, this must be
-       consistent with the definition of the user-provided
-       ParticleInterfaceType::SetParticleData(...) function
-     */
-    template <typename... Args>
-    StackIteratorInterface(StackType& data, const unsigned int index, const Args... args)
-        : fIndex(index)
-        , fData(&data) {
-      (**this).SetParticleData(args...);
-    }
-
-    /** constructor that also sets new values on particle data object, including reference
-        to parent particle
-        @param data reference to the stack [rw]
-        @param index index on stack
-        @param reference to parent particle [rw]. This can be used for thinning, particle
-       counting, history, etc.
-        @param args variadic list of data to initialize stack entry, this must be
-       consistent with the definition of the user-provided
-       ParticleInterfaceType::SetParticleData(...) function
-    */
-    template <typename... Args>
-    StackIteratorInterface(StackType& data, const unsigned int index,
-                           StackIteratorInterface& parent, const Args... args)
-        : fIndex(index)
-        , fData(&data) {
-      (**this).SetParticleData(*parent, args...);
-    }
-
-  public:
-    /** @name Iterator interface
-        @{
-    */
-    StackIteratorInterface& operator++() {
-      ++fIndex;
-      return *this;
-    }
-    StackIteratorInterface operator++(int) {
-      StackIteratorInterface tmp(*this);
-      ++fIndex;
-      return tmp;
-    }
-    StackIteratorInterface operator+(int delta) {
-      return StackIteratorInterface(*fData, fIndex + delta);
-    }
-    bool operator==(const StackIteratorInterface& rhs) { return fIndex == rhs.fIndex; }
-    bool operator!=(const StackIteratorInterface& rhs) { return fIndex != rhs.fIndex; }
-
-    /**
-     * Convert iterator to value type, where value type is the user-provided particle
-     * readout class
-     */
-    ParticleInterfaceType& operator*() {
-      return static_cast<ParticleInterfaceType&>(*this);
-    }
-    /**
-     * Convert iterator to const value type, where value type is the user-provided
-     * particle readout class
-     */
-    const ParticleInterfaceType& operator*() const {
-      return static_cast<const ParticleInterfaceType&>(*this);
-    }
-    ///@}
-
-  protected:
-    /**
-     * @name Stack data access
-     * @{
-     */
-    /// Get current particle index
-    inline unsigned int GetIndex() const { return fIndex; }
-    /// Get current particle Stack object
-    inline StackType& GetStack() { return *fData; }
-    /// Get current particle const Stack object
-    inline const StackType& GetStack() const { return *fData; }
-    /// Get current user particle StackDataType object
-    inline StackDataType& GetStackData() { return fData->GetStackData(); }
-    /// Get current const user particle StackDataType object
-    inline const StackDataType& GetStackData() const { return fData->GetStackData(); }
-    /// Get data index as mapped in Stack class
-    inline unsigned int GetIndexFromIterator() const {
-      return fData->GetIndexFromIterator(fIndex);
-    }
-    ///@}
-  }; // end class StackIterator
-
-  /**
-     @class ConstStackIteratorInterface
-
-     This is the iterator class for const-access to stack data
-   */
-
-  template <typename StackDataType, template <typename> typename ParticleInterface,
-            typename StackType = Stack<StackDataType, ParticleInterface>>
-  class ConstStackIteratorInterface
-      : public ParticleInterface<
-            ConstStackIteratorInterface<StackDataType, ParticleInterface, StackType>> {
-
-  public:
-    typedef ParticleInterface<
-        ConstStackIteratorInterface<StackDataType, ParticleInterface, StackType>>
-        ParticleInterfaceType;
-
-    friend class Stack<StackDataType, ParticleInterface>;   // for access to GetIndex
-    friend class ParticleBase<ConstStackIteratorInterface>; // for access to
-                                                            // GetStackDataType
-
-  private:
-    unsigned int fIndex = 0;
-    const StackType* fData = 0; // info: Particles and StackIterators become invalid when
-                                // parent Stack is copied or deleted!
-
-    // we don't want to allow dangling iterators to exist
-    ConstStackIteratorInterface() = delete;
-
-  public:
-    ConstStackIteratorInterface(const StackType& data, const unsigned int index)
-        : fIndex(index)
-        , fData(&data) {}
-
-    /**
-       @class ConstStackIteratorInterface
-
-       The const counterpart of StackIteratorInterface, which is used
-       for read-only iterator access on particle stack:
-
-       \verbatim
-       for (const auto& p : theStack) { E += p.GetEnergy(); }
-       \endverbatim
-
-       See documentation of StackIteratorInterface for more details.
-    */
-
-  public:
-    /** @name Iterator interface
-     */
-    ///@{
-    ConstStackIteratorInterface& operator++() {
-      ++fIndex;
-      return *this;
-    }
-    ConstStackIteratorInterface operator++(int) {
-      ConstStackIteratorInterface tmp(*this);
-      ++fIndex;
-      return tmp;
-    }
-    ConstStackIteratorInterface operator+(int delta) {
-      return ConstStackIteratorInterface(*fData, fIndex + delta);
-    }
-    bool operator==(const ConstStackIteratorInterface& rhs) {
-      return fIndex == rhs.fIndex;
-    }
-    bool operator!=(const ConstStackIteratorInterface& rhs) {
-      return fIndex != rhs.fIndex;
-    }
-
-    const ParticleInterfaceType& operator*() const {
-      return static_cast<const ParticleInterfaceType&>(*this);
-    }
-    ///@}
-
-  protected:
-    /** @name Stack data access
-        Only the const versions for read-only access
-     */
-    ///@{
-    inline unsigned int GetIndex() const { return fIndex; }
-    inline const StackType& GetStack() const { return *fData; }
-    inline const StackDataType& GetStackData() const { return fData->GetStackData(); }
-    /// Get data index as mapped in Stack class
-    inline unsigned int GetIndexFromIterator() const {
-      return fData->GetIndexFromIterator(fIndex);
-    }
-    ///@}
-  }; // end class ConstStackIterator
-
-} // namespace corsika::stack

Framework/StackInterface/testSecondaryView.cc

-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#include <corsika/stack/SecondaryView.h>
-#include <corsika/stack/Stack.h>
-
-#include <testTestStack.h> // for testing: simple stack. This is a
-// test-build, and inluce file is obtained from CMAKE_CURRENT_SOURCE_DIR
-
-#include <iomanip>
-#include <iostream>
-#include <vector>
-
-#include <catch2/catch.hpp>
-
-using namespace corsika;
-using namespace corsika::stack;
-using namespace std;
-
-typedef Stack<TestStackData, TestParticleInterface> StackTest;
-
-/*
-  See Issue 161
-
-  unfortunately clang does not support this in the same way (yet) as
-  gcc, so we have to distinguish here. If clang cataches up, we could
-  remove the clang branch here and also in corsika::Cascade. The gcc
-  code is much more generic and universal.
- */
-#if defined(__clang__)
-using StackTestView = SecondaryView<TestStackData, TestParticleInterface>;
-#elif defined(__GNUC__) || defined(__GNUG__)
-using StackTestView = MakeView<StackTest>::type;
-#endif
-
-using Particle = typename StackTest::ParticleType;
-
-TEST_CASE("SecondaryStack", "[stack]") {
-
-  // helper function for sum over stack data
-  auto sum = [](const StackTest& stack) {
-    double v = 0;
-    for (const auto& p : stack) v += p.GetData();
-    return v;
-  };
-
-  SECTION("secondary view") {
-    StackTest s;
-    REQUIRE(s.GetSize() == 0);
-    s.AddParticle(std::tuple{9.9});
-    s.AddParticle(std::tuple{8.8});
-    const double sumS = 9.9 + 8.8;
-
-    auto particle = s.GetNextParticle();
-
-    StackTestView view(particle);
-    REQUIRE(view.GetSize() == 0);
-
-    {
-      auto proj = view.GetProjectile();
-      REQUIRE(proj.GetData() == particle.GetData());
-      proj.AddSecondary(std::tuple{4.4});
-    }
-
-    view.AddSecondary(std::tuple{4.5});
-    view.AddSecondary(std::tuple{4.6});
-
-    REQUIRE(view.GetSize() == 3);
-    REQUIRE(s.GetSize() == 5);
-    REQUIRE(!view.IsEmpty());
-
-    auto sumView = [](const StackTestView& stack) {
-      double value = 0;
-      for (const auto& p : stack) { value += p.GetData(); }
-      return value;
-    };
-
-    REQUIRE(sum(s) == sumS + 4.4 + 4.5 + 4.6);
-    REQUIRE(sumView(view) == 4.4 + 4.5 + 4.6);
-
-    view.DeleteLast();
-    REQUIRE(view.GetSize() == 2);
-    REQUIRE(s.GetSize() == 4);
-
-    REQUIRE(sum(s) == sumS + 4.4 + 4.5);
-    REQUIRE(sumView(view) == 4.4 + 4.5);
-
-    auto pDel = view.GetNextParticle();
-    view.Delete(pDel);
-    REQUIRE(view.GetSize() == 1);
-    REQUIRE(s.GetSize() == 3);
-
-    REQUIRE(sum(s) == sumS + 4.4 + 4.5 - pDel.GetData());
-    REQUIRE(sumView(view) == 4.4 + 4.5 - pDel.GetData());
-
-    view.Delete(view.GetNextParticle());
-    REQUIRE(sum(s) == sumS);
-    REQUIRE(sumView(view) == 0);
-    REQUIRE(view.IsEmpty());
-
-    {
-      auto proj = view.GetProjectile();
-      REQUIRE(proj.GetData() == particle.GetData());
-    }
-  }
-
-  SECTION("secondary view, construct from ParticleType") {
-    StackTest s;
-    REQUIRE(s.GetSize() == 0);
-    s.AddParticle(std::tuple{9.9});
-    s.AddParticle(std::tuple{8.8});
-
-    auto iterator = s.GetNextParticle();
-    typename StackTest::ParticleType& particle = iterator; // as in corsika::Cascade
-
-    StackTestView view(particle);
-    REQUIRE(view.GetSize() == 0);
-
-    view.AddSecondary(std::tuple{4.4});
-
-    REQUIRE(view.GetSize() == 1);
-  }
-
-  SECTION("deletion") {
-    StackTest stack;
-    stack.AddParticle(std::tuple{-99.});
-    stack.AddParticle(std::tuple{0.});
-
-    {
-      auto particle = stack.GetNextParticle();
-      StackTestView view(particle);
-
-      auto proj = view.GetProjectile();
-      proj.AddSecondary(std::tuple{-2.});
-      proj.AddSecondary(std::tuple{-1.});
-      proj.AddSecondary(std::tuple{1.});
-      proj.AddSecondary(std::tuple{2.});
-
-      CHECK(stack.GetSize() == 6); // -99, 0, -2, -1, 1, 2
-      CHECK(view.GetSize() == 4);  // -2, -1, 1, 2
-
-      // now delete all negative entries, i.e. -1 and -2
-      auto p = view.begin();
-      while (p != view.end()) {
-        auto data = p.GetData();
-        if (data < 0) {
-          p.Delete();
-        } else {
-          ++p;
-        }
-      }
-      CHECK(stack.GetSize() == 4); // -99, 0, 2, 1 (order changes during deletion)
-      CHECK(view.GetSize() == 2);  // 2, 1
-    }
-
-    // repeat
-
-    {
-      auto particle = stack.GetNextParticle();
-      StackTestView view(particle);
-
-      // put -2,...,+2 on stack
-      auto proj = view.GetProjectile();
-      proj.AddSecondary(std::tuple{-2.});
-      proj.AddSecondary(std::tuple{-1.});
-      proj.AddSecondary(std::tuple{1.});
-      proj.AddSecondary(std::tuple{2.});
-      // stack should contain -99, 0, 2, 1, [-2, -1, 1, 2]
-
-      auto p = view.begin();
-      while (p != view.end()) {
-        auto data = p.GetData();
-        if (data < 0) {
-          p.Delete();
-        } else {
-          ++p;
-        }
-      }
-
-      // stack should contain -99, 0, 2, 1, [2, 1]
-      // view should contain 1, 2
-
-      CHECK(stack.GetSize() == 6);
-    }
-  }
-}

Framework/StackInterface/testStackInterface.cc

-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#include <corsika/stack/Stack.h>
-
-#include <testTestStack.h> // simple test-stack for testing. This is
-                           // for testing only: include from
-                           // CMAKE_CURRENT_SOURCE_DIR
-
-#include <iomanip>
-#include <iostream>
-#include <tuple>
-#include <vector>
-
-#include <catch2/catch.hpp>
-
-using namespace corsika;
-using namespace corsika::stack;
-using namespace std;
-
-typedef Stack<TestStackData, TestParticleInterface> StackTest;
-
-TEST_CASE("Stack", "[Stack]") {
-
-  // helper function for sum over stack data
-  auto sum = [](const StackTest& stack) {
-    double v = 0;
-    for (const auto& p : stack) v += p.GetData();
-    return v;
-  };
-
-  SECTION("StackInterface") {
-
-    // construct a valid Stack object
-    StackTest s;
-    s.Clear();
-    s.AddParticle(std::tuple{0.});
-    s.Copy(s.cbegin(), s.begin());
-    s.Swap(s.begin(), s.begin());
-    REQUIRE(s.GetSize() == 1);
-  }
-
-  SECTION("construct") {
-
-    // construct a valid, empty Stack object
-    StackTest s;
-  }
-
-  SECTION("write and read") {
-
-    StackTest s;
-    s.AddParticle(std::tuple{9.9});
-    const double v = sum(s);
-    REQUIRE(v == 9.9);
-  }
-
-  SECTION("delete from stack") {
-
-    StackTest s;
-    REQUIRE(s.GetSize() == 0);
-    StackTest::StackIterator p =
-        s.AddParticle(std::tuple{0.}); // valid way to access particle data
-    p.SetData(9.9);
-    REQUIRE(s.GetSize() == 1);
-    s.Delete(p);
-    REQUIRE(s.GetSize() == 0);
-  }
-
-  SECTION("delete particle") {
-
-    StackTest s;
-    REQUIRE(s.GetSize() == 0);
-    auto p = s.AddParticle(
-        std::tuple{9.9}); // also valid way to access particle data, identical to above
-    REQUIRE(s.GetSize() == 1);
-    p.Delete();
-    REQUIRE(s.GetSize() == 0);
-  }
-
-  SECTION("create secondaries") {
-
-    StackTest s;
-    REQUIRE(s.GetSize() == 0);
-    auto iter = s.AddParticle(std::tuple{9.9});
-    StackTest::ParticleInterfaceType& p =
-        *iter; // also this is valid to access particle data
-    REQUIRE(s.GetSize() == 1);
-    p.AddSecondary(std::tuple{4.4});
-    REQUIRE(s.GetSize() == 2);
-    /*p.AddSecondary(3.3, 2.2);
-    REQUIRE(s.GetSize() == 3);
-    double v = 0;
-    for (auto& p : s) { v += p.GetData(); }
-    REQUIRE(v == 9.9 + 4.4 + 3.3 + 2.2);*/
-  }
-
-  SECTION("get next particle") {
-    StackTest s;
-    REQUIRE(s.GetSize() == 0);
-    s.AddParticle(std::tuple{9.9});
-    s.AddParticle(std::tuple{8.8});
-    auto particle = s.GetNextParticle(); // first particle
-    REQUIRE(particle.GetData() == 8.8);
-
-    particle.Delete();
-    auto particle2 = s.GetNextParticle(); // first particle
-    REQUIRE(particle2.GetData() == 9.9);
-    particle2.Delete();
-
-    REQUIRE(s.GetSize() == 0);
-  }
-}