Merge branch '106-need-stackview-object' into 'master'

Resolve "Need "StackView" object"

Closes #106

See merge request AirShowerPhysics/corsika!71

COAST/COASTStack.h

@@ -92,11 +92,12 @@ namespace corsika::coast {
     const crs::CParticle* fParticle2 = 0;
 
   public:
-    void Init() {}
-    void SetParticle(const crs::CParticle* v1, const crs::CParticle* v2) {
+    COASTStackImpl(const crs::CParticle* v1, const crs::CParticle* v2) {
       fParticle1 = v1;
       fParticle2 = v2;
     }
+
+    void Init() {}
     void Clear() {}
 
     // there is one particle only
@@ -177,11 +178,10 @@ namespace corsika::coast {
      */
     void Swap(const int, const int) {}
 
-  protected:
-    // size cannot be increased
+    // size cannot be increased, do nothing
     void IncrementSize() {}
 
-    // size cannot be decremented
+    // size cannot be decremented, do nothing
     void DecrementSize() {}
 
   }; // end class COASTStackImpl

COAST/COASTUserLib.cc

@@ -32,7 +32,6 @@ using namespace std;
 using namespace corsika;
 using namespace corsika::units::si;
 
-coast::COASTStack gCOASTStack;
 corsika::coast::COASTProcess gCorsikaProcess;
 
 /*
@@ -94,14 +93,14 @@ extern "C" void track_([[maybe_unused]] const crs::CParticle& pre,
     int    particleId;
     int    hadronicGeneration;
   */
-  gCOASTStack.SetParticle(&pre, &post);
-  const auto particle = gCOASTStack.GetNextParticle();
+  coast::COASTStack stack(&pre, &post);
+  const auto particle = stack.GetNextParticle();
   const geometry::CoordinateSystem& rootCS =
       geometry::RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
   geometry::Line const line(particle.GetPosition(rootCS), particle.GetVelocity(rootCS));
   const TimeType time = particle.GetTimeInterval();
   const geometry::Trajectory<geometry::Line> track(line, time);
-  gCorsikaProcess.DoContinuous(particle, track, gCOASTStack);
+  gCorsikaProcess.DoContinuous(particle, track, stack);
 }
 
 extern "C" void tabularizedatmosphere_([[maybe_unused]] const int& nPoints,

Documentation/Doxygen/Doxyfile.in

@@ -3,6 +3,7 @@ PROJECT_NUMBER         = 0.0.0
 
 GENERATE_HTML          = YES
 GENERATE_LATEX         = YES
+GENERATE_XML           = YES
 
 OUTPUT_DIRECTORY       = @CMAKE_CURRENT_BINARY_DIR@/
 INPUT                  = @PROJECT_SOURCE_DIR@ @PROJECT_BINARY_DIR@/Framework

Documentation/Examples/stack_example.cc

@@ -19,9 +19,10 @@
 #include <iomanip>
 #include <iostream>
 
+using namespace corsika;
 using namespace corsika::units::si;
 using namespace corsika::stack;
-using namespace corsika;
+using namespace corsika::geometry;
 using namespace std;
 
 void fill(corsika::stack::super_stupid::SuperStupidStack& s) {

Environment/FlatAtmosphere/FlatAtmosphere.h

@@ -8,3 +8,5 @@
  * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
  * the license.
  */
+
+

Framework/Cascade/Cascade.h

@@ -22,8 +22,33 @@
 #include <cmath>
 #include <iostream>
 
+/**
+ * The cascade namespace assembles all objects needed to simulate full particles cascades.
+ */
+
 namespace corsika::cascade {
 
+  /**
+   * The Cascade class is constructed from template arguments making
+   * it very versatile. Via the template arguments physics models are
+   * plugged into the cascade simulation.
+   *
+   * <b>Tracking</b> must be a class according to the
+   * TrackingInterface providing the functions: <code>void
+   * Init();</code> and <code>auto GetTrack(Particle const& p)</auto>,
+   * where the latter has a return type of <code>
+   * geometry::Trajectory<corsika::geometry::Line or Helix> </code>
+   *
+   * <b>ProcessList</b> must be a ProcessSequence.
+   *            TimeOfIntersection(corsika::geometry::Line const& line,
+   *
+   * <b>Stack</b> is the storage object for particle data, i.e. with
+   * Particle class type <code>Stack::ParticleType</code>
+   *
+   *
+
+   */
+
   template <typename Tracking, typename ProcessList, typename Stack>
   class Cascade {
     using Particle = typename Stack::ParticleType;
@@ -31,6 +56,10 @@ namespace corsika::cascade {
     Cascade() = delete;
 
   public:
+    /**
+     * Cascade class cannot be default constructed, but needs a valid
+     * list of physics processes for configuration at construct time.
+     */
     Cascade(corsika::environment::Environment const& env, Tracking& tr, ProcessList& pl,
             Stack& stack)
         : fEnvironment(env)
@@ -38,12 +67,20 @@ namespace corsika::cascade {
         , fProcessSequence(pl)
         , fStack(stack) {}
 
+    /**
+     * The Init function is called before the actual cascade simulations.
+     * All components of the Cascade simulation must be configured here.
+     */
     void Init() {
       fTracking.Init();
       fProcessSequence.Init();
       fStack.Init();
     }
 
+    /**
+     * The Run function is the main simulation loop, which processes
+     * particles from the Stack until the Stack is empty.
+     */
     void Run() {
       while (!fStack.IsEmpty()) {
         while (!fStack.IsEmpty()) {
@@ -57,6 +94,16 @@ namespace corsika::cascade {
     }
 
   private:
+    /**
+     * The Step function is executed for each particle from the
+     * stack. It will calcualte geometric transport of the particles,
+     * and apply continuous and stochastic processes to it, which may
+     * lead to energy losses, scattering, absorption, decays and the
+     * production of secondary particles.
+     *
+     * New particles produced in one step are subject to further
+     * processing, e.g. thinning, etc.
+     */
     void Step(Particle& particle) {
       using namespace corsika::units::si;
 

Framework/Logging/Logger.h

@@ -31,9 +31,6 @@
 #include <corsika/logging/NoSink.h>
 #include <corsika/logging/Sink.h>
 
-using namespace std;
-using namespace boost;
-
 namespace corsika::logging {
 
   /**

Framework/StackInterface/CMakeLists.txt

 set (
   CORSIKAstackinterface_HEADERS
   Stack.h
+  SecondaryView.h
   StackIteratorInterface.h
   ParticleBase.h
   )

Framework/StackInterface/ParticleBase.h

@@ -12,7 +12,7 @@
 #ifndef _include_particleBase_h_
 #define _include_particleBase_h_
 
-class StackData; // forward decl
+#include <type_traits>
 
 namespace corsika::stack {
 
@@ -89,19 +89,21 @@ namespace corsika::stack {
 
   protected:
     /**
-        @name Access to underlying stack data
+        @name Access to underlying stack data, these are service
+        function for user classes. User code can only rely on GetIndex
+        and GetStackData to retrieve data
         @{
     */
     auto& GetStackData() { return GetIterator().GetStackData(); }
     const auto& GetStackData() const { return GetIterator().GetStackData(); }
     auto& GetStack() { return GetIterator().GetStack(); }
     const auto& GetStack() const { return GetIterator().GetStack(); }
-    ///@}
 
     /**
      * return the index number of the underlying iterator object
      */
-    unsigned int GetIndex() const { return GetIterator().GetIndex(); }
+    unsigned int GetIndex() const { return GetIterator().GetIndexFromIterator(); }
+    ///@}
   };
 
 } // namespace corsika::stack

Framework/StackInterface/SecondaryView.h

+#ifndef _include_corsika_stack_secondaryview_h_
+#define _include_corsika_stack_secondaryview_h_
+
+#include <corsika/stack/Stack.h>
+
+#include <algorithm>
+#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.
+   */
+
+  /* INTERNAL NOTE FOR DEVELOPERS The secondary particle indices are
+     stored in a std::vector fIndices, the index of the primary projectle
+     particle is explicitly stored in fProjectileIndex. StackIterator
+     indices are refering to those numbers, where
+     StackIterator::GetIndex()==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 InnerStackTypeV = Stack<StackDataType, ParticleInterface>;
+    typedef StackIteratorInterface<typename std::remove_reference<StackDataType>::type,
+                                   ParticleInterface, InnerStackTypeV>
+        StackIteratorV;
+    /// @}
+
+  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 = 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);
+
+  public:
+    SecondaryView(StackIteratorV& vP)
+        : Stack<StackDataType&, ParticleInterface>(vP.GetStackData())
+        , fProjectileIndex(vP.GetIndex()) {}
+
+    auto GetProjectile() {
+      // NOTE: 0 is special marker here for PROJECTILE, see GetIndexFromIterator
+      return StackIterator(*this, 0);
+    }
+
+    template <typename... Args>
+    auto AddSecondary(const Args... v) {
+      InnerStackType::GetStackData().IncrementSize();
+      const unsigned int idSec = GetSize();
+      const unsigned int index = InnerStackType::GetStackData().GetSize() - 1;
+      fIndices.push_back(index);
+      StackIterator proj = GetProjectile();
+      // 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
+     */
+    void Delete(StackIterator p) {
+      if (IsEmpty()) { /* error */
+        throw std::runtime_error("Stack, cannot delete entry since size is zero");
+      }
+      if (p.GetIndex() < GetSize() - 1)
+        InnerStackType::GetStackData().Copy(GetSize() - 1, p.GetIndex());
+      DeleteLast();
+    }
+
+    /**
+     * need overwrite Stack::Delete, since we want to call SecondaryView::DeleteLast
+     */
+    void Delete(ParticleType 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:
+    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;
+  };
+
+} // namespace corsika::stack
+
+#endif

Framework/StackInterface/Stack.h

@@ -15,6 +15,25 @@
 #include <corsika/stack/StackIteratorInterface.h>
 
 #include <stdexcept>
+#include <type_traits>
+
+#include <corsika/stack/SecondaryView.h>
+
+// SFINAE test
+template <typename T>
+class HasGetIndexFromIterator {
+private:
+  typedef char YesType[1];
+  typedef char NoType[2];
+
+  template <typename C>
+  static YesType& test(decltype(&C::GetIndexFromIterator));
+  template <typename C>
+  static NoType& test(...);
+
+public:
+  enum { value = sizeof(test<T>(0)) == sizeof(YesType) };
+};
 
 /**
    All classes around management of particles on a stack.
@@ -36,7 +55,7 @@ namespace corsika::stack {
   /**
      The Stack class provides (and connects) the main particle data storage machinery.
 
-     The StackData type is the user-provided bare data storage
+     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.
@@ -44,23 +63,48 @@ namespace corsika::stack {
      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
-     StackData, given an 'unsigned int' index.
+     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 StackData, template <typename> typename ParticleInterface>
-  class Stack : public StackData {
+  template <typename StackDataType, template <typename> typename ParticleInterface>
+  class Stack {
+
+    using StackType = Stack<StackDataType, ParticleInterface>;
+
+    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:
+    //    template<typename = std::enable_if_t<std::is_reference<StackDataType>{}>>
+    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.
+     */
+    template <typename... Args,
+              typename = std::enable_if_t<!std::is_reference<StackDataType>{}>>
+    Stack(Args... args)
+        : fData(args...) {}
+    // , typename std::enable_if<!std::is_reference<StackDataType>::value,
+    // std::nullptr_t>::type = nullptr)
 
   public:
-    typedef StackData StackImpl; ///< this is the type of the user-provided data structure
+    typedef StackDataType
+        StackImpl; ///< this is the type of the user-provided data structure
+
     template <typename SI>
     using PIType = ParticleInterface<SI>;
-    // typedef ParticleInterface<StackIteratorInterface> StackParticleInterface;  ///<
-    // this is the type of the user-provided ParticleInterface typedef Stack<StackData,
-    // ParticleInterface> StackType;
 
     /**
      * Via the StackIteratorInterface and ConstStackIteratorInterface
@@ -69,33 +113,44 @@ namespace corsika::stack {
      * object. Using CRTP, this also determines the type of
      * ParticleInterface template class simultaneously.
      */
-    typedef StackIteratorInterface<StackData, ParticleInterface> StackIterator;
-    typedef ConstStackIteratorInterface<StackData, ParticleInterface> ConstStackIterator;
+    using StackIterator =
+        StackIteratorInterface<typename std::remove_reference<StackDataType>::type,
+                               ParticleInterface, StackType>;
+    using ConstStackIterator =
+        ConstStackIteratorInterface<typename std::remove_reference<StackDataType>::type,
+                                    ParticleInterface, StackType>;
+
     /**
      * this is the full type of the declared ParticleInterface: typedef typename
      */
     typedef typename StackIterator::ParticleInterfaceType ParticleType;
 
-    friend class StackIteratorInterface<StackData, ParticleInterface>;
-    friend class ConstStackIteratorInterface<StackData, ParticleInterface>;
+    friend class StackIteratorInterface<
+        typename std::remove_reference<StackDataType>::type, ParticleInterface,
+        StackType>;
 
-  protected:
-    using StackData::Copy;
-    using StackData::Swap;
+    friend class ConstStackIteratorInterface<
+        typename std::remove_reference<StackDataType>::type, ParticleInterface,
+        StackType>;
 
   public:
-    using StackData::GetCapacity;
-    using StackData::GetSize;
+    unsigned int GetCapacity() const { return fData.GetCapacity(); }
+    unsigned int GetSize() const { return fData.GetSize(); }
 
-    using StackData::Clear;
-
-    using StackData::DecrementSize;
-    using StackData::IncrementSize;
-
-    using StackData::Init;
+    template <typename... Args>
+    auto Init(Args... args) {
+      return fData.Init(args...);
+    }
+    template <typename... Args>
+    auto Clear(Args... args) {
+      return fData.Clear(args...);
+    }
 
   public:
-    /// these are functions required by std containers and std loops
+    /**
+     * @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); }
@@ -107,43 +162,81 @@ namespace corsika::stack {
     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
+    /**
+     * increase stack size, create new particle at end of stack
+     */
     template <typename... Args>
     StackIterator AddParticle(const Args... v) {
-      IncrementSize();
+      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) {
-      IncrementSize();
+      fData.IncrementSize();
       return StackIterator(*this, GetSize() - 1, parent, v...);
     }
-    void Swap(StackIterator a, StackIterator b) { Swap(a.GetIndex(), b.GetIndex()); }
+
+    void Swap(StackIterator a, StackIterator b) {
+      fData.Swap(a.GetIndex(), b.GetIndex());
+    }
     void Swap(ConstStackIterator a, ConstStackIterator b) {
-      Swap(a.GetIndex(), b.GetIndex());
+      fData.Swap(a.GetIndex(), b.GetIndex());
     }
-    void Copy(StackIterator a, StackIterator b) { Copy(a.GetIndex(), b.GetIndex()); }
-    void Copy(ConstStackIterator a, StackIterator b) { Copy(a.GetIndex(), b.GetIndex()); }
-    /// delete this particle
+    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) Copy(GetSize() - 1, p.GetIndex());
+      if (p.GetIndex() < GetSize() - 1) fData.Copy(GetSize() - 1, p.GetIndex());
       DeleteLast();
       // p.SetInvalid();
     }
+    /**
+     * delete this particle
+     */
     void Delete(ParticleType p) { Delete(p.GetIterator()); }
-    /// delete last particle on stack by decrementing stack size
-    void DeleteLast() { DecrementSize(); }
-    /// check if there are no further particles on stack
+
+    /**
+     * 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:
-    StackData& GetStackData() { return static_cast<StackData&>(*this); }
-    const StackData& GetStackData() const { return static_cast<const StackData&>(*this); }
+    // typename std::enable_if<HasGetIndexFromIterator<T>::value, unsigned int>::type
+    // typename std::enable_if<std::is_base_of<decltype(*this)>,
+    // SecondaryView<StackDataType, ParticleInterface>>::value, unsigned int>::type
+    unsigned int GetIndexFromIterator(const unsigned int vI) const { return vI; }
+
+    typename std::remove_reference<StackDataType>::type& GetStackData() { return fData; }
+    const typename std::remove_reference<StackDataType>::type& GetStackData() const {
+      return fData;
+    }
   };
 
 } // namespace corsika::stack

Framework/StackInterface/StackIteratorInterface.h

@@ -16,13 +16,14 @@
 
 #include <type_traits>
 
-class StackData; // forward decl
-
 namespace corsika::stack {
 
-  template <typename StackData, template <typename> typename ParticleInterface>
+  template <typename StackDataType, template <typename> typename ParticleInterface>
   class Stack; // forward decl
 
+  template <typename StackDataType, template <typename> typename ParticleInterface>
+  class SecondaryView; // forward decl
+
   /**
      @class StackIteratorInterface
 
@@ -47,22 +48,32 @@ namespace corsika::stack {
      for a particle entry at any index fIndex.
   */
 
-  template <typename StackData, template <typename> typename ParticleInterface>
+  template <typename StackDataType, template <typename> typename ParticleInterface,
+            typename StackType = Stack<StackDataType, ParticleInterface>>
   class StackIteratorInterface
-      : public ParticleInterface<StackIteratorInterface<StackData, ParticleInterface>> {
+      : public ParticleInterface<
+            StackIteratorInterface<StackDataType, ParticleInterface, StackType>> {
 
   public:
-    typedef Stack<StackData, ParticleInterface> StackType;
+    typedef StackDataType SD;
+    // typedef Stack<StackDataType, ParticleInterface> StackType;
     /*typedef
-        typename std::conditional<std::is_const<StackData>::value,
-                                  const Stack<const StackData, ParticleInterface>&,
-                                  Stack<StackData, ParticleInterface>&>::type StackType;*/
+        typename std::conditional<std::is_const<StackDataType>::value,
+                                  const Stack<const StackDataType, ParticleInterface>&,
+                                  Stack<StackDataType, ParticleInterface>&>::type
+       StackType;*/
 
-    typedef ParticleInterface<StackIteratorInterface<StackData, ParticleInterface>>
+    typedef ParticleInterface<
+        StackIteratorInterface<StackDataType, ParticleInterface, StackType>>
         ParticleInterfaceType;
 
-    friend class Stack<StackData, ParticleInterface>;  // for access to GetIndex
-    friend class ParticleBase<StackIteratorInterface>; // for access to GetStackData
+    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;
@@ -92,7 +103,8 @@ namespace corsika::stack {
     StackIteratorInterface(StackType& data, const unsigned int index, const Args... args)
         : fIndex(index)
         , fData(&data) {
-      (**this).SetParticleData(args...);
+      ParticleInterfaceType& p = **this;
+      p.SetParticleData(args...);
     }
 
     /** constructor that also sets new values on particle data object, including reference
@@ -110,7 +122,8 @@ namespace corsika::stack {
                            StackIteratorInterface& parent, const Args... args)
         : fIndex(index)
         , fData(&data) {
-      (**this).SetParticleData(*parent, args...);
+      ParticleInterfaceType& p = **this;
+      p.SetParticleData(*parent, args...);
     }
 
   public:
@@ -131,30 +144,42 @@ namespace corsika::stack {
     }
     bool operator==(const StackIteratorInterface& rhs) { return fIndex == rhs.fIndex; }
     bool operator!=(const StackIteratorInterface& rhs) { return fIndex != rhs.fIndex; }
-    /// Convert to value type
+
+    /**
+     * Convert iterator to value type, where value type is the user-provided particle
+     * readout class
+     */
     ParticleInterfaceType& operator*() {
       return static_cast<ParticleInterfaceType&>(*this);
     }
-    /// Convert to const value type
+    /**
+     * 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
+    /**
+     * @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 StackData object
-    inline StackData& GetStackData() { return fData->GetStackData(); }
-    /// Get current const user particle StackData object
-    inline const StackData& GetStackData() const { return fData->GetStackData(); }
+    /// 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
 
@@ -164,18 +189,20 @@ namespace corsika::stack {
      This is the iterator class for const-access to stack data
    */
 
-  template <typename StackData, template <typename> typename ParticleInterface>
+  template <typename StackDataType, template <typename> typename ParticleInterface,
+            typename StackType = Stack<StackDataType, ParticleInterface>>
   class ConstStackIteratorInterface
       : public ParticleInterface<
-            ConstStackIteratorInterface<StackData, ParticleInterface>> {
+            ConstStackIteratorInterface<StackDataType, ParticleInterface, StackType>> {
 
   public:
-    typedef Stack<StackData, ParticleInterface> StackType;
-    typedef ParticleInterface<ConstStackIteratorInterface<StackData, ParticleInterface>>
+    typedef ParticleInterface<
+        ConstStackIteratorInterface<StackDataType, ParticleInterface, StackType>>
         ParticleInterfaceType;
 
-    friend class Stack<StackData, ParticleInterface>;       // for access to GetIndex
-    friend class ParticleBase<ConstStackIteratorInterface>; // for access to GetStackData
+    friend class Stack<StackDataType, ParticleInterface>;   // for access to GetIndex
+    friend class ParticleBase<ConstStackIteratorInterface>; // for access to
+                                                            // GetStackDataType
 
   private:
     unsigned int fIndex = 0;
@@ -238,7 +265,11 @@ namespace corsika::stack {
     ///@{
     inline unsigned int GetIndex() const { return fIndex; }
     inline const StackType& GetStack() const { return *fData; }
-    inline const StackData& GetStackData() const { return fData->GetStackData(); }
+    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
 

Framework/StackInterface/testStackInterface.cc

@@ -9,8 +9,18 @@
  * the license.
  */
 
+#include <corsika/stack/SecondaryView.h>
 #include <corsika/stack/Stack.h>
 
+#include <boost/type_index.hpp>
+#include <type_traits>
+using boost::typeindex::type_id_with_cvr;
+
+template <typename T>
+struct ID {
+  typedef T type;
+};
+
 #include <iomanip>
 #include <iostream>
 #include <vector>
@@ -42,7 +52,6 @@ public:
   void SetData(const int i, const double v) { fData[i] = v; }
   double GetData(const int i) const { return fData[i]; }
 
-protected:
   // these functions are also needed by the Stack interface
   void IncrementSize() { fData.push_back(0.); }
   void DecrementSize() {
@@ -64,15 +73,26 @@ class TestParticleInterface : public ParticleBase<StackIteratorInterface> {
   using ParticleBase<StackIteratorInterface>::GetIterator;
 
 public:
-  // one version
-  StackIteratorInterface& AddSecondary(const double v) {
-    GetStack().AddParticle(v);
-    return GetIterator();
-  }
-  // another version
-  void AddSecondary(const double v, const double p) { GetStack().AddParticle(v + p); }
+  /* normally this function does not need to be specified here, it is
+   part of ParticleBase<StackIteratorInterface>.
+
+   However, when overloading this function again with a different
+   parameter set, as here, seems to require also the declaration of
+   the original function here...
+  */
 
+  // default version for particle-creation from input data
   void SetParticleData(const double v) { SetData(v); }
+  void SetParticleData(TestParticleInterface<StackIteratorInterface>& /*parent*/,
+                       const double v) {
+    SetData(v);
+  }
+  /// alternative set-particle data for non-standard construction from different inputs
+  void SetParticleData(const double v, const double p) { SetData(v + p); }
+  void SetParticleData(TestParticleInterface<StackIteratorInterface>& /*parent*/,
+                       const double v, const double p) {
+    SetData(v + p);
+  }
 
   void SetData(const double v) { GetStackData().SetData(GetIndex(), v); }
   double GetData() const { return GetStackData().GetData(GetIndex()); }
@@ -96,13 +116,10 @@ TEST_CASE("Stack", "[Stack]") {
     StackTest s;
     s.Init();
     s.Clear();
-    s.IncrementSize();
+    s.AddParticle(0.);
     s.Copy(s.cbegin(), s.begin());
     s.Swap(s.begin(), s.begin());
-    s.GetCapacity();
     REQUIRE(s.GetSize() == 1);
-    s.DecrementSize();
-    REQUIRE(s.GetSize() == 0);
   }
 
   SECTION("construct") {
@@ -156,4 +173,81 @@ TEST_CASE("Stack", "[Stack]") {
     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(9.9);
+    s.AddParticle(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);
+  }
+
+  SECTION("secondary view") {
+    StackTest s;
+    REQUIRE(s.GetSize() == 0);
+    s.AddParticle(9.9);
+    s.AddParticle(8.8);
+    const double sumS = 9.9 + 8.8;
+
+    auto particle = s.GetNextParticle();
+
+    typedef SecondaryView<TestStackData, TestParticleInterface> StackTestView;
+    StackTestView v(particle);
+    REQUIRE(v.GetSize() == 0);
+
+    {
+      auto proj = v.GetProjectile();
+      REQUIRE(proj.GetData() == particle.GetData());
+    }
+
+    v.AddSecondary(4.4);
+    v.AddSecondary(4.5);
+    v.AddSecondary(4.6);
+
+    REQUIRE(v.GetSize() == 3);
+    REQUIRE(s.GetSize() == 5);
+    REQUIRE(!v.IsEmpty());
+
+    auto sumView = [](const StackTestView& stack) {
+      double v = 0;
+      for (const auto& p : stack) { v += p.GetData(); }
+      return v;
+    };
+
+    REQUIRE(sum(s) == sumS + 4.4 + 4.5 + 4.6);
+    REQUIRE(sumView(v) == 4.4 + 4.5 + 4.6);
+
+    v.DeleteLast();
+    REQUIRE(v.GetSize() == 2);
+    REQUIRE(s.GetSize() == 4);
+
+    REQUIRE(sum(s) == sumS + 4.4 + 4.5);
+    REQUIRE(sumView(v) == 4.4 + 4.5);
+
+    auto pDel = v.GetNextParticle();
+    v.Delete(pDel);
+    REQUIRE(v.GetSize() == 1);
+    REQUIRE(s.GetSize() == 3);
+
+    REQUIRE(sum(s) == sumS + 4.4 + 4.5 - pDel.GetData());
+    REQUIRE(sumView(v) == 4.4 + 4.5 - pDel.GetData());
+
+    v.Delete(v.GetNextParticle());
+    REQUIRE(sum(s) == sumS);
+    REQUIRE(sumView(v) == 0);
+    REQUIRE(v.IsEmpty());
+
+    {
+      auto proj = v.GetProjectile();
+      REQUIRE(proj.GetData() == particle.GetData());
+    }
+  }
 }

Main/Environment.h

@@ -8,3 +8,5 @@
  * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
  * the license.
  */
+
+

Main/Stack.h

@@ -8,3 +8,5 @@
  * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
  * the license.
  */
+
+

Processes/Sibyll/Decay.h

@@ -90,7 +90,6 @@ namespace corsika::process {
 
         using std::cout;
         using std::endl;
-        // using namespace corsika::io;
         using namespace corsika::units::si;
 
         // name? also makes EM particles stable

Processes/Sibyll/Interaction.h

@@ -150,6 +150,8 @@ namespace corsika::process::sibyll {
 
           auto const [productionCrossSection, elaCrossSection, numberOfNucleons] =
               GetCrossSection(corsikaBeamId, targetId, ECoM);
+          [[maybe_unused]] auto elaCrossSectionCopy =
+              elaCrossSection; // ONLY TO AVOID COMPILER WARNING
 
           std::cout << "Interaction: "
                     << " IntLength: sibyll return (mb): "
@@ -278,8 +280,10 @@ namespace corsika::process::sibyll {
           const auto [sigProd, sigEla, nNuc] =
               GetCrossSection(corsikaBeamId, targetId, Ecm);
           cross_section_of_components[i] = sigProd;
-          int ideleteme = nNuc;  // to avoid not used warning in array binding
-          ideleteme = ideleteme; // to avoid not used warning in array binding
+          [[maybe_unused]] int ideleteme =
+              nNuc; // to avoid not used warning in array binding
+          [[maybe_unused]] auto sigElaCopy =
+              sigEla; // to avoid not used warning in array binding
         }
 
         const auto targetCode = currentNode->GetModelProperties().SampleTarget(

Processes/Sibyll/NuclearInteraction.h

@@ -176,6 +176,7 @@ namespace corsika::process::sibyll {
           // TODO: remove number of nucleons, avg target mass is available in environment
           auto const [productionCrossSection, elaCrossSection, numberOfNucleons] =
               GetCrossSection(corsikaBeamId, targetId, ECoM);
+          [[maybe_unused]] auto elaCrossSectionCopy = elaCrossSection;
 
           std::cout << "NuclearInteraction: "
                     << "IntLength: nuclib return (mb): "
@@ -322,6 +323,8 @@ namespace corsika::process::sibyll {
         cout << "beam id: " << targetId << endl;
         const auto [sigProd, sigEla, nNuc] = GetCrossSection(beamId, targetId, EcmNN);
         cross_section_of_components[i] = sigProd;
+        [[maybe_unused]] auto sigElaCopy = sigEla; // ONLY TO AVOID COMPILER WARNINGS
+        [[maybe_unused]] auto sigNucCopy = nNuc;   // ONLY TO AVOID COMPILER WARNINGS
       }
 
       const auto targetCode = currentNode->GetModelProperties().SampleTarget(
@@ -350,6 +353,7 @@ namespace corsika::process::sibyll {
       const auto protonId = corsika::particles::Proton::GetCode();
       const auto [prodCrossSection, elaCrossSection, dum] =
           fHadronicInteraction.GetCrossSection(protonId, protonId, EcmNN);
+      [[maybe_unused]] auto dumCopy = dum; // ONLY TO AVOID COMPILER WARNING
       const double sigProd = prodCrossSection / 1_mbarn;
       const double sigEla = elaCrossSection / 1_mbarn;
       // sample number of interactions (only input variables, output in common cnucms)

Processes/Sibyll/SibStack.h

@@ -80,7 +80,6 @@ namespace corsika::process::sibyll {
         std::swap(s_plist_.p[i][i1], s_plist_.p[i][i2]);
     }
 
-  protected:
     void IncrementSize() { s_plist_.np++; }
     void DecrementSize() {
       if (s_plist_.np > 0) { s_plist_.np--; }

Processes/Sibyll/sibyll2.3c.h

@@ -70,7 +70,6 @@ extern struct {
   int lun;
 } s_debug_;
 
-  
 // lund random generator setup
 // extern struct {int mrlu[6]; float rrlu[100]; }ludatr_;
 
@@ -79,7 +78,7 @@ void sibyll_(const int&, const int&, const double&);
 
 // subroutine to initiate sibyll
 void sibyll_ini_();
- 
+
 // subroutine to SET DECAYS
 void dec_ini_();
 
@@ -106,7 +105,5 @@ double get_sibyll_mass2(int&);
 
 // phojet random generator setup
 void pho_rndin_(int&, int&, int&, int&);
-
-  
 }
 #endif