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/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/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::history {
-  template <typename T>
-  class HistorySecondaryView; // forward decl. for befriending
-}
-
-namespace corsika::stack {
-
-  template <typename TStackData, template <typename> typename TParticleInterface>
-  class Stack; // forward decl
-
-  template <typename TStackData, template <typename> typename TParticleInterface>
-  class SecondaryView; // forward decl
-
-  template <typename TStackData, template <typename> typename TParticleInterface,
-            typename StackType>
-  class ConstStackIteratorInterface; // 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 index_ in the Stack data.
-
-     The template argument `TParticleInterface` acts as a policy to provide
-     readout function of Particle data from the stack. The TParticleInterface
-     class must know how to retrieve information from the Stack data
-     for a particle entry at any index index_.
-
-     The TParticleInterface 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 TStackData. And
-     StackIteratorInterface::GetIndex() provides the iterator index to
-     be readout. The TStackData is another user-provided class to
-     store data and must implement functions compatible with
-     TParticleInterface, in this example TStackData::GetData(const unsigned int
-     vIndex).
-
-     For two examples see stack_example.cc, or the
-     corsika::processes::sibyll::SibStack class
-  */
-
-  template <typename TStackData, template <typename> typename TParticleInterface,
-            typename StackType = Stack<TStackData, TParticleInterface>>
-  class StackIteratorInterface
-      : public TParticleInterface<
-            StackIteratorInterface<TStackData, TParticleInterface, StackType>> {
-
-  public:
-    using ParticleInterfaceType =
-        TParticleInterface<corsika::stack::StackIteratorInterface<
-            TStackData, TParticleInterface, StackType>>;
-
-    // friends are needed for access to protected methods
-    friend class Stack<TStackData,
-                       TParticleInterface>; // for access to GetIndex for Stack
-    friend class Stack<TStackData&, TParticleInterface>; // for access to GetIndex
-                                                         // SecondaryView : public Stack
-    friend class ParticleBase<StackIteratorInterface>;   // for access to GetStackData
-    friend class SecondaryView<TStackData,
-                               TParticleInterface>; // access for SecondaryView
-
-    template <typename T>
-    friend class corsika::history::HistorySecondaryView;
-
-    friend class ConstStackIteratorInterface<TStackData, TParticleInterface, StackType>;
-
-  protected:
-    unsigned int index_ = 0;
-
-  private:
-    StackType* data_ = 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)
-        : index_(vR.index_)
-        , data_(vR.data_) {}
-
-    StackIteratorInterface& operator=(StackIteratorInterface const& vR) {
-      index_ = vR.index_;
-      data_ = vR.data_;
-      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)
-        : index_(index)
-        , data_(&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)
-        : index_(index)
-        , data_(&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)
-        : index_(index)
-        , data_(&data) {
-      (**this).SetParticleData(*parent, args...);
-    }
-
-    bool isDeleted() const { return GetStack().isDeleted(*this); }
-
-  public:
-    /** @name Iterator interface
-        @{
-    */
-    StackIteratorInterface& operator++() {
-      do {
-        ++index_;
-      } while (
-          GetStack().isDeleted(*this)); // this also check the allowed bounds of index_
-      return *this;
-    }
-    StackIteratorInterface operator++(int) {
-      StackIteratorInterface tmp(*this);
-      do {
-        ++index_;
-      } while (
-          GetStack().isDeleted(*this)); // this also check the allowed bounds of index_
-      return tmp;
-    }
-    StackIteratorInterface operator+(int delta) const {
-      return StackIteratorInterface(*data_, index_ + delta);
-    }
-    bool operator==(const StackIteratorInterface& rhs) const {
-      return index_ == rhs.index_;
-    }
-    bool operator!=(const StackIteratorInterface& rhs) const {
-      return index_ != rhs.index_;
-    }
-    bool operator==(
-        const ConstStackIteratorInterface<TStackData, TParticleInterface, StackType>& rhs)
-        const; // implement below
-    bool operator!=(
-        const ConstStackIteratorInterface<TStackData, TParticleInterface, StackType>& rhs)
-        const; // implement below
-
-    /**
-     * 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 index_; }
-    /// Get current particle Stack object
-    inline StackType& GetStack() { return *data_; }
-    /// Get current particle const Stack object
-    inline const StackType& GetStack() const { return *data_; }
-    /// Get current user particle TStackData object
-    inline TStackData& GetStackData() { return data_->GetStackData(); }
-    /// Get current const user particle TStackData object
-    inline const TStackData& GetStackData() const { return data_->GetStackData(); }
-    /// Get data index as mapped in Stack class
-    inline unsigned int GetIndexFromIterator() const {
-      return data_->GetIndexFromIterator(index_);
-    }
-    ///@}
-  }; // end class StackIterator
-
-  /**
-     @class ConstStackIteratorInterface
-
-     This is the iterator class for const-access to stack data
-   */
-
-  template <typename TStackData, template <typename> typename TParticleInterface,
-            typename StackType = Stack<TStackData, TParticleInterface>>
-  class ConstStackIteratorInterface
-      : public TParticleInterface<
-            ConstStackIteratorInterface<TStackData, TParticleInterface, StackType>> {
-
-  public:
-    typedef TParticleInterface<
-        ConstStackIteratorInterface<TStackData, TParticleInterface, StackType>>
-        ParticleInterfaceType;
-
-    // friends are needed for access to protected methods
-    friend class Stack<TStackData,
-                       TParticleInterface>; // for access to GetIndex for Stack
-    friend class Stack<TStackData&, TParticleInterface>; // for access to GetIndex
-                                                         // SecondaryView : public Stack
-    friend class ParticleBase<ConstStackIteratorInterface>; // for access to GetStackData
-    friend class SecondaryView<TStackData,
-                               TParticleInterface>; // access for SecondaryView
-
-    friend class StackIteratorInterface<TStackData, TParticleInterface, StackType>;
-
-  protected:
-    unsigned int index_ = 0;
-
-  private:
-    const StackType* data_ = 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)
-        : index_(index)
-        , data_(&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.
-    */
-
-    bool isDeleted() const { return GetStack().isDeleted(*this); }
-
-  public:
-    /** @name Iterator interface
-     */
-    ///@{
-    ConstStackIteratorInterface& operator++() {
-      do {
-        ++index_;
-      } while (
-          GetStack().isDeleted(*this)); // this also check the allowed bounds of index_
-      return *this;
-    }
-    ConstStackIteratorInterface operator++(int) {
-      ConstStackIteratorInterface tmp(*this);
-      do {
-        ++index_;
-      } while (
-          GetStack().isDeleted(*this)); // this also check the allowed bounds of index_
-      return tmp;
-    }
-    ConstStackIteratorInterface operator+(const int delta) const {
-      return ConstStackIteratorInterface(*data_, index_ + delta);
-    }
-    bool operator==(const ConstStackIteratorInterface& rhs) const {
-      return index_ == rhs.index_;
-    }
-    bool operator!=(const ConstStackIteratorInterface& rhs) const {
-      return index_ != rhs.index_;
-    }
-    bool operator==(const StackIteratorInterface<TStackData, TParticleInterface,
-                                                 StackType>& rhs) const {
-      return index_ == rhs.index_;
-    }
-    bool operator!=(const StackIteratorInterface<TStackData, TParticleInterface,
-                                                 StackType>& rhs) const {
-      return index_ != rhs.index_;
-    }
-
-    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 index_; }
-    inline const StackType& GetStack() const { return *data_; }
-    inline const TStackData& GetStackData() const { return data_->GetStackData(); }
-    /// Get data index as mapped in Stack class
-    inline unsigned int GetIndexFromIterator() const {
-      return data_->GetIndexFromIterator(index_);
-    }
-    ///@}
-  }; // end class ConstStackIterator
-
-  template <typename TStackData, template <typename> typename TParticleInterface,
-            typename StackType>
-  bool StackIteratorInterface<TStackData, TParticleInterface, StackType>::operator==(
-      const ConstStackIteratorInterface<TStackData, TParticleInterface, StackType>& rhs)
-      const {
-    return index_ == rhs.index_;
-  }
-
-  template <typename TStackData, template <typename> typename TParticleInterface,
-            typename StackType>
-  bool StackIteratorInterface<TStackData, TParticleInterface, StackType>::operator!=(
-      const ConstStackIteratorInterface<TStackData, TParticleInterface, StackType>& rhs)
-      const {
-    return index_ != rhs.index_;
-  }
-
-} // namespace corsika::stack

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.
- */
-
-#define protected public // to also test the internal state of objects
-
-#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());
-    CHECK(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);
-    CHECK(v == 9.9);
-  }
-
-  SECTION("delete from stack") {
-
-    StackTest s;
-    CHECK(s.getSize() == 0);
-    StackTest::StackIterator p =
-        s.AddParticle(std::tuple{0.}); // valid way to access particle data
-    p.SetData(9.9);
-    CHECK(s.getSize() == 1);
-    CHECK(s.getEntries() == 1);
-    s.Delete(p);
-    CHECK(s.getSize() == 1);
-    CHECK(s.getEntries() == 0);
-  }
-
-  SECTION("delete particle") {
-
-    StackTest s;
-    CHECK(s.getSize() == 0);
-    s.AddParticle(std::tuple{8.9});
-    s.AddParticle(std::tuple{7.9});
-    auto p = s.AddParticle(
-        std::tuple{9.9}); // also valid way to access particle data, identical to above
-
-    CHECK(s.getSize() == 3);
-    CHECK(s.getEntries() == 3);
-    CHECK(!s.IsEmpty());
-
-    p.Delete(); // mark for deletion: size=3, entries=2
-    CHECK(s.getSize() == 3);
-    CHECK(s.getEntries() == 2);
-    CHECK(!s.IsEmpty());
-
-    s.last().Delete(); // mark for deletion: size=3, entries=1
-    CHECK(s.getSize() == 3);
-    CHECK(s.getEntries() == 1);
-    CHECK(!s.IsEmpty());
-
-    /*
-       GetNextParticle will find two entries marked as "deleted" and
-       will purge this from the end of the stack: size = 1
-    */
-    s.GetNextParticle().Delete(); // mark for deletion: size=3, entries=0
-    CHECK(s.getSize() == 1);
-    CHECK(s.getEntries() == 0);
-    CHECK(s.IsEmpty());
-  }
-
-  SECTION("create secondaries") {
-
-    StackTest s;
-    CHECK(s.getSize() == 0);
-    auto iter = s.AddParticle(std::tuple{9.9});
-    StackTest::ParticleInterfaceType& p =
-        *iter; // also this is valid to access particle data
-    CHECK(s.getSize() == 1);
-    p.AddSecondary(std::tuple{4.4});
-    CHECK(s.getSize() == 2);
-    /*p.AddSecondary(3.3, 2.2);
-    CHECK(s.getSize() == 3);
-    double v = 0;
-    for (auto& p : s) { v += p.GetData(); }
-    CHECK(v == 9.9 + 4.4 + 3.3 + 2.2);*/
-  }
-
-  SECTION("get next particle") {
-    StackTest s;
-    CHECK(s.getSize() == 0);
-    CHECK(s.getEntries() == 0);
-    CHECK(s.IsEmpty());
-
-    s.AddParticle(std::tuple{9.9});
-    s.AddParticle(std::tuple{8.8});
-    CHECK(s.getSize() == 2);
-    CHECK(s.getEntries() == 2);
-    CHECK(!s.IsEmpty());
-
-    auto particle = s.GetNextParticle(); // first particle
-    CHECK(particle.GetData() == 8.8);
-
-    particle.Delete(); // only marks (last) particle as deleted
-    CHECK(s.getSize() == 2);
-    CHECK(s.getEntries() == 1);
-    CHECK(!s.IsEmpty());
-
-    /*
-      This following call to GetNextParticle will realize that the
-      current last particle on the stack was marked "deleted" and will
-      purge it: stack size is reduced by one.
-     */
-    auto particle2 = s.GetNextParticle(); // first particle
-    CHECK(particle2.GetData() == 9.9);
-    CHECK(s.getSize() == 1);
-    CHECK(s.getEntries() == 1);
-    CHECK(!s.IsEmpty());
-
-    particle2.Delete(); // also mark this particle as deleted
-
-    CHECK(s.getSize() == 1);
-    CHECK(s.getEntries() == 0);
-    CHECK(s.IsEmpty());
-  }
-}

Framework/Testing/CMakeLists.txt

-set (
-  TESTING_SOURCES
-  TestMain.cc
-  )
-
-set (
-  TESTING_HEADERS
-  )
-
-set (
-  TESTING_NAMESPACE
-  corsika/testing
-  )
-
-add_library (CORSIKAtesting STATIC ${TESTING_SOURCES})
-CORSIKA_COPY_HEADERS_TO_NAMESPACE (CORSIKAtesting ${TESTING_NAMESPACE} ${TESTING_HEADERS})
-
-set_target_properties (
-  CORSIKAtesting
-  PROPERTIES
-  VERSION ${PROJECT_VERSION}
-  SOVERSION 1
-  PUBLIC_HEADER "${TESTING_HEADERS}"
-  )
-
-# target dependencies on other libraries (also the header onlys)
-target_link_libraries (
-  CORSIKAtesting
-  CORSIKAthirdparty
-  )
-
-# target_include_directories (
-#  CORSIKAtesting
-#  SYSTEM
-#  PUBLIC    ${CATCH2_INCLUDE_DIR}
-#  INTERFACE ${CATCH2_INCLUDE_DIR}
-#  )
-
-target_include_directories (
-  CORSIKAtesting
-  INTERFACE
-  $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
-  $<INSTALL_INTERFACE:include/include>
-  )
-
-install (
-  TARGETS CORSIKAtesting
-  LIBRARY DESTINATION lib
-  ARCHIVE DESTINATION lib
-  PUBLIC_HEADER DESTINATION include/${TESTING_NAMESPACE}
-  )
-

Framework/Units/CMakeLists.txt

-add_library (CORSIKAunits INTERFACE)
-
-set (CORSIKAunits_NAMESPACE corsika/units)
-set (
-  CORSIKAunits_HEADERS
-  PhysicalUnits.h
-  PhysicalConstants.h
-  )
-
-CORSIKA_COPY_HEADERS_TO_NAMESPACE (CORSIKAunits ${CORSIKAunits_NAMESPACE} ${CORSIKAunits_HEADERS})
-
-target_include_directories (CORSIKAunits
-  INTERFACE
-  $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
-  $<BUILD_INTERFACE:${PROJECT_SOURCE_DIR}/ThirdParty>
-  $<INSTALL_INTERFACE:include>
-  )
-
-install (FILES ${CORSIKAunits_HEADERS} DESTINATION include/${CORSIKAunits_NAMESPACE})
-
-# code testing
-CORSIKA_ADD_TEST (testUnits)
-target_link_libraries (testUnits CORSIKAunits CORSIKAtesting)

Framework/Utilities/Bit.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
-
-/**
-  \author Hans Dembinski
-  \author Lukas Nellen
-  \author Darko Veberic
-  \date 27 Jan 2014
-
-  \version $Id: Bit.h 25126 2014-02-03 22:13:10Z darko $
-*/
-
-#include <exception>
-
-// #include <utl/AugerException.h>
-
-namespace corsika::utl {
-
-  namespace Bit {
-
-    template <typename T>
-    class Array {
-    public:
-      Array(T& target)
-          : fTarget(target) {}
-
-      class Bit {
-      public:
-        Bit(T& target, T mask)
-            : fTarget(target)
-            , fMask(mask) {}
-
-        operator bool() const { return fTarget & fMask; }
-
-        bool operator~() const { return !bool(*this); }
-
-        Bit& operator=(const bool value) {
-          if (value)
-            fTarget |= fMask;
-          else
-            fTarget &= ~fMask;
-          return *this;
-        }
-
-        Bit& Flip() { return *this = ~(*this); }
-
-      private:
-        T& fTarget;
-        T fMask;
-      };
-
-      Bit operator[](unsigned int position) { return Bit(fTarget, T(1) << position); }
-
-      Bit At(unsigned int position) {
-        if (position >= 8 * sizeof(T))
-          // throw std::exceptionOutOfBoundException("Running out of bits.");
-          throw std::exception("Running out of bits.");
-        return (*this)[position];
-      }
-
-      template <typename M>
-      Array& Mask(const M mask, const bool value) {
-        Bit(fTarget, mask) = value;
-        return *this;
-      }
-
-      template <typename M>
-      T Get(const M mask) {
-        return fTarget & T(mask);
-      }
-
-    private:
-      T& fTarget;
-    };
-
-  } // namespace Bit
-
-  // helper
-  template <typename T>
-  inline Bit::Array<T> AsBitArray(T& target) {
-    return Bit::Array<T>(target);
-  }
-
-} // namespace corsika::utl

Framework/Utilities/CMakeLists.txt

-#
-# cfenv feature test - select implementation to use
-#
-try_compile (HAS_FEENABLEEXCEPT "${CMAKE_CURRENT_BINARY_DIR}" "${CMAKE_CURRENT_SOURCE_DIR}/try_feenableexcept.cc")
-if (HAS_FEENABLEEXCEPT)
-  set (CORSIKA_FENV "CorsikaFenvDefault.cc")
-  set_property(DIRECTORY ${CMAKE_HOME_DIRECTORY} APPEND PROPERTY COMPILE_DEFINITIONS "HAS_FEENABLEEXCEPT")
-else ()
-  if (APPLE)
-    set (CORSIKA_FENV "CorsikaFenvOSX.cc")
-  else()
-    set (CORSIKA_FENV "CorsikaFenvFallback.cc")
-  endif()
-endif ()
-
-#
-# library setup
-#
-set (
-  UTILITIES_SOURCES  
-  COMBoost.cc
-  CorsikaData.cc
-  ${CORSIKA_FENV})
-
-set (
-  UTILITIES_HEADERS
-  CorsikaData.h
-  COMBoost.h
-  Bit.h
-  Singleton.h
-  sgn.h
-  CorsikaFenv.h
-  MetaProgramming.h
-  )
-
-set (
-  UTILITIES_DEPENDS
-  CORSIKAgeometry
-  CORSIKAunits
-  C8::ext::boost # so far only for MetaProgramming
-  C8::ext::eigen3 # for COMboost
-  )
-
-if (TARGET cnpy)
-  LIST (APPEND
-    UTILITIES_HEADERS
-    SaveBoostHistogram.hpp
-    )
-  LIST (APPEND
-    UTILITIES_DEPENDS
-    cnpy # for SaveBoostHistogram
-    )
-endif (TARGET cnpy)
-  
-set (
-  UTILITIES_NAMESPACE
-  corsika/utl
-  )
-
-add_library (CORSIKAutilities STATIC ${UTILITIES_SOURCES})
-CORSIKA_COPY_HEADERS_TO_NAMESPACE (CORSIKAutilities ${UTILITIES_NAMESPACE} ${UTILITIES_HEADERS})
-
-set_target_properties (
-  CORSIKAutilities
-  PROPERTIES
-  VERSION ${PROJECT_VERSION}
-  SOVERSION 1
-  PUBLIC_HEADER "${UTILITIES_HEADERS}"
-  )
-
-# target dependencies on other libraries (also the header onlys)
-target_link_libraries (
-  CORSIKAutilities
-  ${UTILITIES_DEPENDS}
-  )
-
-target_include_directories (
-  CORSIKAutilities
-  PUBLIC
-  $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
-  $<INSTALL_INTERFACE:include/include>
-  )
-
-install (
-  TARGETS CORSIKAutilities
-  LIBRARY DESTINATION lib
-  ARCHIVE DESTINATION lib
-  PUBLIC_HEADER DESTINATION include/${UTILITIES_NAMESPACE}
-  )
-
-
-# --------------------
-# code unit testing
-CORSIKA_ADD_TEST(testCOMBoost)
-target_link_libraries (
-  testCOMBoost
-  CORSIKAutilities
-  CORSIKAtesting
-  )
-
-CORSIKA_ADD_TEST(testCorsikaFenv)
-target_link_libraries (
-  testCorsikaFenv
-  CORSIKAutilities
-  CORSIKAtesting
-)
-
-if (TARGET cnpy)
-  CORSIKA_ADD_TEST(testSaveBoostHistogram)
-  target_link_libraries (
-    testSaveBoostHistogram
-    CORSIKAutilities
-    CORSIKAtesting
-    )
-endif (TARGET cnpy)