From 67e78cad72c03171ab5fb03f26211cf209220aaa Mon Sep 17 00:00:00 2001
From: ralfulrich <ralf.ulrich@kit.edu>
Date: Sun, 16 Sep 2018 23:50:44 +0200
Subject: [PATCH] - some majro updates to the stack-interface - revised version
 of process interface - added tests for stacks, processes, cascade

---
 Documentation/Examples/CMakeLists.txt         |   5 +-
 Documentation/Examples/stack_example.cc       |   6 +-
 .../Examples/staticsequence_example.cc        |  39 +++---
 Framework/CMakeLists.txt                      |   1 +
 Framework/Cascade/CMakeLists.txt              |  66 ++++++++++
 Framework/Cascade/Cascade.cc                  |  16 +--
 Framework/Cascade/Cascade.h                   |  57 +++++++--
 Framework/Cascade/testCascade.cc              | 103 +++++++++++++++
 Framework/Geometry/CMakeLists.txt             |   1 +
 Framework/Geometry/Helix.h                    |   2 -
 Framework/Geometry/LineTrajectory.h           |  15 +--
 Framework/ProcessSequence/CMakeLists.txt      |  41 +++---
 Framework/ProcessSequence/ProcessSequence.h   |  29 +++--
 .../ProcessSequence/testProcessSequence.cc    |  87 +++++++++++++
 Framework/StackInterface/CMakeLists.txt       |   8 ++
 Framework/StackInterface/ParticleBase.h       |  53 ++++++++
 Framework/StackInterface/Stack.h              |  63 +++++++---
 Framework/StackInterface/StackIterator.h      | 104 +++++++---------
 .../StackInterface/testStackInterface.cc      | 117 ++++++++++++++++++
 Framework/Units/PhysicalUnits.h               |   7 +-
 Framework/Units/testUnits.cc                  |   1 +
 Stack/CMakeLists.txt                          |   1 +
 Stack/DummyStack/CMakeLists.txt               |  29 +++++
 Stack/DummyStack/DummyStack.h                 |  61 +++++++++
 Stack/DummyStack/SuperStupidStack.h           |  97 +++++++++++++++
 Stack/SuperStupidStack/SuperStupidStack.h     |  61 +++++----
 26 files changed, 884 insertions(+), 186 deletions(-)
 create mode 100644 Framework/Cascade/CMakeLists.txt
 create mode 100644 Framework/Cascade/testCascade.cc
 create mode 100644 Framework/ProcessSequence/testProcessSequence.cc
 create mode 100644 Framework/StackInterface/ParticleBase.h
 create mode 100644 Framework/StackInterface/testStackInterface.cc
 create mode 100644 Stack/DummyStack/CMakeLists.txt
 create mode 100644 Stack/DummyStack/DummyStack.h
 create mode 100644 Stack/DummyStack/SuperStupidStack.h

diff --git a/Documentation/Examples/CMakeLists.txt b/Documentation/Examples/CMakeLists.txt
index d4efb2fc7..237e7a0b2 100644
--- a/Documentation/Examples/CMakeLists.txt
+++ b/Documentation/Examples/CMakeLists.txt
@@ -15,6 +15,9 @@ target_link_libraries (stack_example SuperStupidStack CORSIKAunits CORSIKAloggin
 install (TARGETS stack_example DESTINATION share/examples)
 
 add_executable (staticsequence_example staticsequence_example.cc)
-target_link_libraries (staticsequence_example SuperStupidStack CORSIKAprocesssequence CORSIKAunits CORSIKAlogging)
+target_link_libraries (staticsequence_example
+  CORSIKAprocesssequence
+  CORSIKAunits
+  CORSIKAlogging)
 install (TARGETS staticsequence_example DESTINATION share/examples)
 
diff --git a/Documentation/Examples/stack_example.cc b/Documentation/Examples/stack_example.cc
index d54f20434..e70b7e6c8 100644
--- a/Documentation/Examples/stack_example.cc
+++ b/Documentation/Examples/stack_example.cc
@@ -13,7 +13,7 @@ using namespace corsika::stack;
 void fill(corsika::stack::super_stupid::SuperStupidStack& s) {
   for (int i = 0; i < 11; ++i) {
     auto p = s.NewParticle();
-    p.SetId(corsika::particles::Code::Electron);
+    p.SetPID(corsika::particles::Code::Electron);
     p.SetEnergy(1.5_GeV * i);
   }
 }
@@ -21,9 +21,9 @@ void fill(corsika::stack::super_stupid::SuperStupidStack& s) {
 void read(corsika::stack::super_stupid::SuperStupidStack& s) {
   cout << "found Stack with " << s.GetSize() << " particles. " << endl;
   EnergyType Etot;
-  for (auto p : s) {
+  for (auto& p : s) {
     Etot += p.GetEnergy();
-    cout << "particle: " << p.GetId() << " with " << p.GetEnergy() / 1_GeV << " GeV"
+    cout << "particle: " << p.GetPID() << " with " << p.GetEnergy() / 1_GeV << " GeV"
          << endl;
   }
   cout << "Etot=" << Etot << " = " << Etot / 1_GeV << " GeV" << endl;
diff --git a/Documentation/Examples/staticsequence_example.cc b/Documentation/Examples/staticsequence_example.cc
index a28150c23..6b8406739 100644
--- a/Documentation/Examples/staticsequence_example.cc
+++ b/Documentation/Examples/staticsequence_example.cc
@@ -7,21 +7,21 @@
 using namespace std;
 using namespace corsika::process;
 
-class Process1 : public corsika::process::BaseProcess<Process1> {
+class Process1 : public BaseProcess<Process1> {
 public:
   Process1() {}
-  template <typename D>
-  void DoContinuous(D& d) const {
+  template <typename D, typename T, typename S>
+  void DoContinuous(D& d, T& t, S& s) const {
     for (int i = 0; i < 10; ++i) d.p[i] += 1;
   }
 };
 
-class Process2 : public corsika::process::BaseProcess<Process2> {
+class Process2 : public BaseProcess<Process2> {
 public:
   Process2() {}
 
-  template <typename D>
-  inline void DoContinuous(D& d) const {
+  template <typename D, typename T, typename S>
+  inline void DoContinuous(D& d, T& t, S& s) const {
     // for (int i=0; i<10; ++i) d.p[i] *= 2;
   }
 };
@@ -31,8 +31,8 @@ public:
   // Process3(const int v) :fV(v) {}
   Process3() {}
 
-  template <typename D>
-  inline void DoContinuous(D& d) const {
+  template <typename D, typename T, typename S>
+  inline void DoContinuous(D& d, T& t, S& s) const {
     // for (int i=0; i<10; ++i) d.p[i] += fV;
   }
 
@@ -44,8 +44,8 @@ class Process4 : public BaseProcess<Process4> {
 public:
   // Process4(const int v) : fV(v) {}
   Process4() {}
-  template <typename D>
-  inline void DoContinuous(D& d) const {
+  template <typename D, typename T, typename S>
+  inline void DoContinuous(D& d, T& t, S& s) const {
     // for (int i=0; i<10; ++i) d.p[i] /= fV;
   }
 
@@ -53,13 +53,13 @@ private:
   // int fV;
 };
 
-class Data {
-public:
-  std::array<double, 10> p{{0.}};
+struct DummyData {
+  double p[10];
 };
+struct DummyStack {};
+struct DummyTrajectory {};
 
 void modular() {
-  Data d0;
 
   Process1 m1;
   Process2 m2;
@@ -68,13 +68,14 @@ void modular() {
 
   const auto sequence = m1 + m2 + m3 + m4;
 
-  const int n = 100000000;
-  for (int i = 0; i < n; ++i) { sequence.DoContinuous(d0); }
+  DummyData p;
+  DummyTrajectory t;
+  DummyStack s;
 
-  double s = 0;
-  for (int i = 0; i < 10; ++i) { s += d0.p[i]; }
+  const int n = 100000000;
+  for (int i = 0; i < n; ++i) { sequence.DoContinuous(p, t, s); }
 
-  cout << scientific << " v=" << s << " n=" << n << endl;
+  cout << " done (nothing...) " << endl;
 }
 
 int main() {
diff --git a/Framework/CMakeLists.txt b/Framework/CMakeLists.txt
index 2fdb7c85a..836d67110 100644
--- a/Framework/CMakeLists.txt
+++ b/Framework/CMakeLists.txt
@@ -5,3 +5,4 @@ add_subdirectory (Particles)
 add_subdirectory (Logging)
 add_subdirectory (StackInterface)
 add_subdirectory (ProcessSequence)
+add_subdirectory (Cascade)
diff --git a/Framework/Cascade/CMakeLists.txt b/Framework/Cascade/CMakeLists.txt
new file mode 100644
index 000000000..63fe2cb85
--- /dev/null
+++ b/Framework/Cascade/CMakeLists.txt
@@ -0,0 +1,66 @@
+
+# namespace of library -> location of header files
+set (
+  CORSIKAcascade_NAMESPACE
+  corsika/cascade
+  )
+
+# header files of this library
+set (
+  CORSIKAcascade_HEADERS
+  Cascade.h
+  )
+
+#set (
+#  CORSIKAcascade_SOURCES
+#  Cascade.cc
+#  )
+
+#add_library (CORSIKAcascade STATIC ${CORSIKAcascade_SOURCES})
+add_library (CORSIKAcascade INTERFACE)
+
+CORSIKA_COPY_HEADERS_TO_NAMESPACE (CORSIKAcascade ${CORSIKAcascade_NAMESPACE} ${CORSIKAcascade_HEADERS})
+
+#target_link_libraries (
+#  CORSIKAcascade
+#  CORSIKAparticles
+#  CORSIKAunits
+#  CORSIKAthirdparty # for catch2
+#  )
+
+# include directive for upstream code
+target_include_directories (
+  CORSIKAcascade
+  INTERFACE
+  $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
+  $<INSTALL_INTERFACE:include/>
+  )
+
+# install library
+install (
+  FILES ${CORSIKAcascade_HEADERS}
+  DESTINATION include/${CORSIKAcascade_NAMESPACE}
+  )
+
+# ----------------
+# code unit testing
+add_executable (
+  testCascade
+  testCascade.cc
+  )
+
+target_link_libraries (
+  testCascade
+  CORSIKAcascade
+  CORSIKAparticles
+  CORSIKAgeometry
+  CORSIKAprocesssequence
+  SuperStupidStack
+  CORSIKAunits
+  CORSIKAthirdparty # for catch2
+  )
+
+add_test (
+  NAME testCascade
+  COMMAND testCascade
+  )
diff --git a/Framework/Cascade/Cascade.cc b/Framework/Cascade/Cascade.cc
index f4496f536..1bd6b5b2b 100644
--- a/Framework/Cascade/Cascade.cc
+++ b/Framework/Cascade/Cascade.cc
@@ -1,20 +1,20 @@
+#include <corsika/cascade/Cascade.h>
 
+using namespace corsika::cascade;
 
-namespace cascade;
-
-template <typename Sequence, typename Trajectory>
-void Cascade::Cascade() {
-  kkk;
-  kk;
+template <typename ProcessList, typename Particle, typename Trajectory, typename Stack>
+Cascade<ProcessList, Particle, Trajectory, Stack>::Cascade() {
+  //  kkk;
+  //  kk;
 }
 
-template <typename Sequence, typename Trajectory>
+template <typename ProcessList, typename Particle, typename Trajectory, typename Stack>
 void Cascade::Init() {
   fStack.Init();
   fProcesseList.Init();
 }
 
-template <typename Sequence, typename Trajectory>
+template <typename ProcessList, typename Particle, typename Trajectory, typename Stack>
 void Cascade::Run() {
   if (!fStack.IsEmpty()) {
     if (!fStack.IsEmpty()) {
diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h
index a8af260b5..8c29a7798 100644
--- a/Framework/Cascade/Cascade.h
+++ b/Framework/Cascade/Cascade.h
@@ -1,23 +1,62 @@
 #ifndef _include_Cascade_h_
 #define _include_Cascade_h_
 
-namespace cascade {
+#include <corsika/geometry/LineTrajectory.h> // to be removed
+#include <corsika/geometry/Point.h>          // to be removed
+#include <corsika/units/PhysicalUnits.h>
 
-  template <typename Processes, typename Trajectory, typename Stack>
+using namespace corsika::units;
+
+namespace corsika::cascade {
+
+  template <typename Trajectory, typename ProcessList, typename Stack>
   class Cascade {
 
+    typedef typename Stack::ParticleType Particle;
+
   public:
-    Cascade();
+    Cascade(ProcessList& pl, Stack& stack)
+        : fProcesseList(pl)
+        , fStack(stack) {}
+
+    void Init() {
+      fStack.Init();
+      fProcesseList.Init();
+    }
+
+    void Run() {
+      while (!fStack.IsEmpty()) {
+        while (!fStack.IsEmpty()) {
+          Particle& p = *fStack.GetNextParticle();
+          Step(p);
+        }
+        // do cascade equations, which can put new particles on Stack,
+        // thus, the double loop
+        // DoCascadeEquations(); //
+      }
+    }
 
-    void Init();
-    void Run();
-    void Step(Particle& particle);
+    void Step(Particle& particle) {
+      double nextStep = fProcesseList.MinStepLength(particle);
+      corsika::geometry::CoordinateSystem root;
+      Trajectory trajectory(
+          corsika::geometry::Point(root, {0_m, 0_m, 0_m}),
+          corsika::geometry::Vector<corsika::units::SpeedType::dimension_type>(
+              root, 0 * 1_m / second, 0 * 1_m / second, 1 * 1_m / second));
+      fProcesseList.DoContinuous(particle, trajectory, fStack);
+      // if (particle.IsMarkedToBeDeleted())
+      {
+        // std::cout << "DELETET THISSKSKJD!" << std::endl;
+        // fStack.Delete(particle);
+      }
+      fProcesseList.DoDiscrete(particle, fStack);
+    }
 
   private:
-    Stack fStack;
-    Processes fProcesseList;
+    Stack& fStack;
+    ProcessList& fProcesseList;
   };
 
-} // namespace cascade
+} // namespace corsika::cascade
 
 #endif
diff --git a/Framework/Cascade/testCascade.cc b/Framework/Cascade/testCascade.cc
new file mode 100644
index 000000000..3f9f1b0fc
--- /dev/null
+++ b/Framework/Cascade/testCascade.cc
@@ -0,0 +1,103 @@
+#include <corsika/cascade/Cascade.h>
+#include <corsika/geometry/LineTrajectory.h>
+#include <corsika/process/ProcessSequence.h>
+#include <corsika/stack/super_stupid/SuperStupidStack.h>
+
+#define CATCH_CONFIG_MAIN // This tells Catch to provide a main() - only do this in one
+                          // cpp file
+#include <catch2/catch.hpp>
+
+using namespace corsika::process;
+using namespace corsika::units;
+
+#include <iostream>
+using namespace std;
+
+class ProcessSplit : public corsika::process::BaseProcess<ProcessSplit> {
+public:
+  ProcessSplit() {}
+
+  template <typename Particle>
+  double MinStepLength(Particle&) const {
+    return 0;
+  }
+
+  template <typename Particle, typename Trajectory, typename Stack>
+  void DoContinuous(Particle& p, Trajectory& t, Stack& s) const {}
+
+  template <typename Particle, typename Stack>
+  void DoDiscrete(Particle& p, Stack& s) const {
+    EnergyType E = p.GetEnergy();
+    if (E < 1_GeV) {
+      p.Delete();
+      fCount++;
+    } else {
+      p.SetEnergy(E / 2);
+      s.NewParticle().SetEnergy(E / 2);
+    }
+  }
+
+  void Init() { fCount = 0; }
+
+  int GetCount() { return fCount; }
+
+private:
+  mutable int fCount = 0;
+};
+
+class ProcessReport : public corsika::process::BaseProcess<ProcessReport> {
+  bool fReport = false;
+
+public:
+  ProcessReport(bool v)
+      : fReport(v) {}
+
+  template <typename Particle>
+  double MinStepLength(Particle&) const {
+    return 0;
+  }
+
+  template <typename Particle, typename Trajectory, typename Stack>
+  void DoContinuous(Particle& p, Trajectory& t, Stack& s) const {
+    if (!fReport) return;
+    static int fCount = 0;
+    std::cout << "generation  " << fCount << std::endl;
+    int i = 0;
+    for (auto& iterP : s) {
+      EnergyType E = iterP.GetEnergy();
+      std::cout << " particle data: " << i++ << ", id=" << iterP.GetPID()
+                << ", E=" << double(E / 1_GeV) << " GeV "
+                << " | " << std::endl;
+    }
+    fCount++;
+  }
+
+  template <typename Particle, typename Stack>
+  void DoDiscrete(Particle& p, Stack& s) const {}
+  void Init() {}
+};
+
+TEST_CASE("Cascade", "[Cascade]") {
+
+  ProcessReport p0(false);
+  ProcessSplit p1;
+  const auto sequence = p0 + p1;
+  corsika::stack::super_stupid::SuperStupidStack stack;
+
+  corsika::cascade::Cascade<corsika::geometry::LineTrajectory, decltype(sequence),
+                            decltype(stack)>
+      EAS(sequence, stack);
+
+  SECTION("sectionTwo") {
+    for (int i = 0; i < 5; ++i) {
+      stack.Clear();
+      auto particle = stack.NewParticle();
+      EnergyType E0 = 100_GeV * pow(10, i);
+      particle.SetEnergy(E0);
+      EAS.Init();
+      EAS.Run();
+
+      cout << "E0=" << E0 / 1_GeV << "GeV, count=" << p1.GetCount() << endl;
+    }
+  }
+}
diff --git a/Framework/Geometry/CMakeLists.txt b/Framework/Geometry/CMakeLists.txt
index c55295c89..36fab64cc 100644
--- a/Framework/Geometry/CMakeLists.txt
+++ b/Framework/Geometry/CMakeLists.txt
@@ -13,6 +13,7 @@ set (
   Helix.h
   BaseVector.h
   QuantityVector.h
+  LineTrajectory.h
   )
 
 set (
diff --git a/Framework/Geometry/Helix.h b/Framework/Geometry/Helix.h
index c392356d4..ca97fa8bf 100644
--- a/Framework/Geometry/Helix.h
+++ b/Framework/Geometry/Helix.h
@@ -17,12 +17,10 @@ namespace corsika::geometry {
   class Helix // TODO: inherit from to-be-implemented "Trajectory"
   {
     using SpeedVec = Vector<SpeedType::dimension_type>;
-
     Point const r0;
     FrequencyType const omegaC;
     SpeedVec const vPar;
     SpeedVec vPerp, uPerp;
-
     LengthType const radius;
 
   public:
diff --git a/Framework/Geometry/LineTrajectory.h b/Framework/Geometry/LineTrajectory.h
index 4d56a598f..fc550d4bd 100644
--- a/Framework/Geometry/LineTrajectory.h
+++ b/Framework/Geometry/LineTrajectory.h
@@ -1,26 +1,27 @@
 #ifndef _include_LINETRAJECTORY_H
 #define _include_LINETRAJECTORY_H
 
-#include <Units/PhysicalUnits.h>
-#include <corsika/Point.h>
-#include <corsika/Vector.h>
+#include <corsika/geometry/Point.h>
+#include <corsika/geometry/Vector.h>
+#include <corsika/units/PhysicalUnits.h>
 
-namesapce corsika {
+namespace corsika::geometry {
 
   class LineTrajectory // TODO: inherit from Trajectory
   {
-    using SpeedVec = Vector<Speed::dimension_type>;
+    using SpeedVec = Vector<corsika::units::SpeedType::dimension_type>;
 
     Point const r0;
     SpeedVec const v0;
 
+  public:
     LineTrajectory(Point const& pR0, SpeedVec const& pV0)
         : r0(r0)
         , v0(pV0) {}
 
-    auto GetPosition(Time t) const { return r0 + v0 * t; }
+    auto GetPosition(corsika::units::TimeType t) const { return r0 + v0 * t; }
   };
 
-} // end namesapce
+} // namespace corsika::geometry
 
 #endif
diff --git a/Framework/ProcessSequence/CMakeLists.txt b/Framework/ProcessSequence/CMakeLists.txt
index dd8475845..ac853c837 100644
--- a/Framework/ProcessSequence/CMakeLists.txt
+++ b/Framework/ProcessSequence/CMakeLists.txt
@@ -1,13 +1,13 @@
 
 add_library (CORSIKAprocesssequence INTERFACE)
 
-# namespace of library -> location of header files
+#namespace of library->location of header files
 set (
   CORSIKAprocesssequence_NAMESPACE
   corsika/process
   )
 
-# header files of this library
+#header files of this library
 set (
   CORSIKAprocesssequence_HEADERS
   ProcessSequence.h
@@ -15,7 +15,7 @@ set (
 
 CORSIKA_COPY_HEADERS_TO_NAMESPACE (CORSIKAprocesssequence ${CORSIKAprocesssequence_NAMESPACE} ${CORSIKAprocesssequence_HEADERS})
 
-# include directive for upstream code
+#include directive for upstream code
 target_include_directories (
   CORSIKAprocesssequence
   INTERFACE
@@ -23,29 +23,26 @@ target_include_directories (
   $<INSTALL_INTERFACE:include/>
   )
 
-# install library
+#install library
 install (
   FILES ${CORSIKAprocesssequence_HEADERS}
   DESTINATION include/${CORSIKAprocesssequence_NAMESPACE}
   )
 
-# ----------------
-# code unit testing
-#add_executable (
-#  testLogging
-#  testLogging.cc
-#  )
-
-#target_link_libraries (
-#  testLogging
-#  CORSIKAprocesssequence
-#  CORSIKAthirdparty # for catch2
-#  )
-
-#add_test (
-#  NAME testLogging
-#  COMMAND testLogging
-#  )
+#-- -- -- -- -- -- -- --
+#code unit testing
+add_executable (
+  testProcessSequence
+  testProcessSequence.cc
+  )
 
+target_link_libraries (
+  testProcessSequence
+  CORSIKAprocesssequence
+  CORSIKAthirdparty # for catch2
+  )
 
-  
+add_test (
+  NAME testProcessSequence
+  COMMAND testProcessSequence
+  )
diff --git a/Framework/ProcessSequence/ProcessSequence.h b/Framework/ProcessSequence/ProcessSequence.h
index 621b7bb04..75ab8dad3 100644
--- a/Framework/ProcessSequence/ProcessSequence.h
+++ b/Framework/ProcessSequence/ProcessSequence.h
@@ -1,13 +1,14 @@
 #ifndef _include_ProcessSequence_h_
 #define _include_ProcessSequence_h_
 
+#include <cmath>
 #include <iostream>
 #include <typeinfo>
 
 namespace corsika::process {
 
   /**
-     /class BaseProcess
+     \class BaseProcess
 
      The structural base type of a process object in a
      ProcessSequence. Both, the ProcessSequence and all its elements
@@ -17,6 +18,7 @@ namespace corsika::process {
 
   template <typename derived>
   struct BaseProcess {
+    derived& GetRef() { return static_cast<derived&>(*this); }
     const derived& GetRef() const { return static_cast<const derived&>(*this); }
   };
 
@@ -41,28 +43,35 @@ namespace corsika::process {
         : A(in_A)
         , B(in_B) {}
 
-    template <typename D>
-    inline void DoContinuous(D& d) const {
-      A.DoContinuous(d);
-      B.DoContinuous(d);
+    template <typename Particle, typename Trajectory, typename Stack>
+    inline void DoContinuous(Particle& p, Trajectory& t, Stack& s) const {
+      A.DoContinuous(p, t, s);
+      B.DoContinuous(p, t, s);
     } // add trajectory
 
     template <typename D>
     inline double MinStepLength(D& d) const {
-      return min(A.MinStepLength(d), B.MinStepLength(d));
+      return std::min(A.MinStepLength(d), B.MinStepLength(d));
     }
 
     // template<typename D>
     // inline Trajectory Transport(D& d, double& length) const { A.Transport(d, length);
     // B.Transport(d, length); }
 
-    template <typename D>
-    inline void DoDiscrete(D& d) const {
-      A.DoDiscrete(d);
-      B.DoDiscrete(d);
+    template <typename Particle, typename Stack>
+    void DoDiscrete(Particle& p, Stack& s) const {
+      A.DoDiscrete(p, s);
+      B.DoDiscrete(p, s);
+    }
+
+    /// TODO the const_cast is not nice, think about the constness here
+    inline void Init() const {
+      const_cast<T1*>(&A)->Init();
+      const_cast<T2*>(&B)->Init();
     }
   };
 
+  /// the + operator that assembles more BaseProcess objects into a ProcessSequence
   template <typename T1, typename T2>
   inline const ProcessSequence<T1, T2> operator+(const BaseProcess<T1>& A,
                                                  const BaseProcess<T2>& B) {
diff --git a/Framework/ProcessSequence/testProcessSequence.cc b/Framework/ProcessSequence/testProcessSequence.cc
new file mode 100644
index 000000000..d057facff
--- /dev/null
+++ b/Framework/ProcessSequence/testProcessSequence.cc
@@ -0,0 +1,87 @@
+#define CATCH_CONFIG_MAIN // This tells Catch to provide a main() - only do this in one
+                          // cpp file
+#include <catch2/catch.hpp>
+
+#include <array>
+#include <iomanip>
+#include <iostream>
+
+#include <corsika/process/ProcessSequence.h>
+
+using namespace std;
+using namespace corsika::process;
+
+class Process1 : public BaseProcess<Process1> {
+public:
+  Process1() {}
+  void Init() {} // cout << "Process1::Init" << endl; }
+  template <typename D, typename T, typename S>
+  void DoContinuous(D& d, T& t, S& s) const {
+    for (int i = 0; i < 10; ++i) d.p[i] += 1 + i;
+  }
+};
+
+class Process2 : public BaseProcess<Process2> {
+public:
+  Process2() {}
+  void Init() {} // cout << "Process2::Init" << endl; }
+  template <typename D, typename T, typename S>
+  inline void DoContinuous(D& d, T& t, S& s) const {
+    for (int i = 0; i < 10; ++i) d.p[i] *= 0.7;
+  }
+};
+
+class Process3 : public BaseProcess<Process3> {
+public:
+  Process3() {}
+  void Init() {} // cout << "Process3::Init" << endl; }
+  template <typename D, typename T, typename S>
+  inline void DoContinuous(D& d, T& t, S& s) const {
+    for (int i = 0; i < 10; ++i) d.p[i] += 0.933;
+  }
+};
+
+class Process4 : public BaseProcess<Process4> {
+public:
+  Process4() {}
+  void Init() {} // cout << "Process4::Init" << endl; }
+  template <typename D, typename T, typename S>
+  inline void DoContinuous(D& d, T& t, S& s) const {
+    for (int i = 0; i < 10; ++i) d.p[i] /= 1.2;
+  }
+  // inline double MinStepLength(D& d) {
+  // void DoDiscrete(Particle& p, Stack& s) const {
+};
+
+struct DummyData {
+  double p[10] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
+};
+struct DummyStack {};
+struct DummyTrajectory {};
+
+TEST_CASE("Cascade", "[Cascade]") {
+
+  SECTION("sectionTwo") {
+
+    Process1 m1;
+    Process2 m2;
+    Process3 m3;
+    Process4 m4;
+
+    const auto sequence = m1 + m2 + m3 + m4;
+
+    DummyData p;
+    DummyTrajectory t;
+    DummyStack s;
+
+    sequence.Init();
+
+    const int n = 100;
+    INFO("Running loop with n=" << n);
+    for (int i = 0; i < n; ++i) { sequence.DoContinuous(p, t, s); }
+
+    for (int i = 0; i < 10; i++) { INFO("data[" << i << "]=" << p.p[i]); }
+  }
+
+  SECTION("sectionThree") {}
+}
diff --git a/Framework/StackInterface/CMakeLists.txt b/Framework/StackInterface/CMakeLists.txt
index 45a3ebb4d..11a30200e 100644
--- a/Framework/StackInterface/CMakeLists.txt
+++ b/Framework/StackInterface/CMakeLists.txt
@@ -2,6 +2,7 @@ set (
   CORSIKAstackinterface_HEADERS
   Stack.h
   StackIterator.h
+  ParticleBase.h
   )
 
 set (
@@ -29,3 +30,10 @@ install (
   DESTINATION
   include/${CORSIKAstackinterface_NAMESPACE}
   )
+
+
+# code testing
+add_executable (testStackInterface testStackInterface.cc)
+target_link_libraries (testStackInterface CORSIKAstackinterface CORSIKAthirdparty) # for catch2
+add_test(NAME testStackInterface COMMAND testStackInterface)
+
diff --git a/Framework/StackInterface/ParticleBase.h b/Framework/StackInterface/ParticleBase.h
new file mode 100644
index 000000000..aa0bb3d11
--- /dev/null
+++ b/Framework/StackInterface/ParticleBase.h
@@ -0,0 +1,53 @@
+#ifndef _include_particleBase_h_
+#define _include_particleBase_h_
+
+class StackData; // forward decl
+
+namespace corsika::stack {
+
+  //  template <typename> class PI;// : public ParticleBase<StackIteratorInterface> {
+  // template <typename, template <typename> typename> class Stack; // forward decl
+
+  /**
+   \class ParticleBase
+
+   The base class to define the readout of particle properties from a
+   particle stack. Every stack must implement this readout via the
+   ParticleBase class.
+  */
+
+  template <typename StackIterator>
+  class ParticleBase {
+
+    //    friend class Stack<StackData, PI>;  // for access to GetIterator
+  public:
+    ParticleBase() {}
+
+  private:
+    ParticleBase(ParticleBase&);
+    // ParticleBase& operation=(ParticleBase& p);
+
+  public:
+    /// delete this particle on the stack. The corresponding iterator
+    /// will be invalidated by this operation
+    void Delete() { GetIterator().GetStack().Delete(GetIterator()); }
+
+    //  protected: // todo should be proteced, but don't now how to 'friend Stack'
+    /// Function to provide CRTP access to inheriting class (type)
+    StackIterator& GetIterator() { return static_cast<StackIterator&>(*this); }
+    const StackIterator& GetIterator() const {
+      return static_cast<const StackIterator&>(*this);
+    }
+
+  protected:
+    /// access to underling stack data
+    auto& GetStackData() { return GetIterator().GetStackData(); }
+    const auto& GetStackData() const { return GetIterator().GetStackData(); }
+
+    /// return the index number of the underlying iterator object
+    int GetIndex() const { return GetIterator().GetIndex(); }
+  };
+
+}; // namespace corsika::stack
+
+#endif
diff --git a/Framework/StackInterface/Stack.h b/Framework/StackInterface/Stack.h
index 42187c1ae..d493ef5ef 100644
--- a/Framework/StackInterface/Stack.h
+++ b/Framework/StackInterface/Stack.h
@@ -9,46 +9,71 @@
 
 namespace corsika::stack {
 
+  template <typename>
+  class PI; // forward decl
+
   /**
      Interface definition of a Stack object. The Stack implements the
      std-type begin/end function to allow integration in normal for
      loops etc.
    */
 
-  template <typename DataImpl, typename Particle>
-  class Stack : public DataImpl {
+  template <typename StackData, template <typename> typename PI>
+  class Stack : public StackData {
+
+  public:
+    typedef Stack<StackData, PI> StackType;
+    typedef StackIteratorInterface<StackData, PI> StackIterator;
+    typedef const StackIterator ConstStackIterator;
+    typedef typename StackIterator::ParticleInterfaceType ParticleType;
+    friend class StackIteratorInterface<StackData, PI>;
 
   public:
-    using DataImpl::GetCapacity;
-    using DataImpl::GetSize;
+    using StackData::GetCapacity;
+    using StackData::GetSize;
 
-    using DataImpl::Clear;
-    using DataImpl::Copy;
+    using StackData::Clear;
+    using StackData::Copy;
 
-    using DataImpl::DecrementSize;
-    using DataImpl::IncrementSize;
+    using StackData::DecrementSize;
+    using StackData::IncrementSize;
 
-  public:
-    typedef Particle iterator;
-    typedef const Particle const_iterator;
+    using StackData::Init;
 
+  public:
     /// these are functions required by std containers and std loops
-    iterator begin() { return iterator(*this, 0); }
-    iterator end() { return iterator(*this, GetSize()); }
-    iterator last() { return iterator(*this, GetSize() - 1); }
+    StackIterator begin() { return StackIterator(*this, 0); }
+    StackIterator end() { return StackIterator(*this, GetSize()); }
+    StackIterator last() { return StackIterator(*this, GetSize() - 1); }
 
     /// these are functions required by std containers and std loops
-    const_iterator cbegin() const { return const_iterator(*this, 0); }
-    const_iterator cend() const { return const_iterator(*this, GetSize()); }
-    const_iterator clast() const { return const_iterator(*this, GetSize() - 1); }
+    ConstStackIterator cbegin() const { return ConstStackIterator(*this, 0); }
+    ConstStackIterator cend() const { return ConstStackIterator(*this, GetSize()); }
+    ConstStackIterator clast() const { return ConstStackIterator(*this, GetSize() - 1); }
 
     /// increase stack size, create new particle at end of stack
-    iterator NewParticle() {
+    StackIterator NewParticle() {
       IncrementSize();
-      return iterator(*this, GetSize() - 1);
+      return StackIterator(*this, GetSize() - 1);
+    }
+    /// delete this particle
+    void Delete(StackIterator& p) {
+      if (GetSize() == 0) { /*error*/
+      }
+      if (p.GetIndex() < GetSize() - 1) Copy(GetSize() - 1, p.GetIndex());
+      DeleteLast();
+      // p.SetInvalid();
     }
+    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
+    bool IsEmpty() { return GetSize() == 0; }
+    StackIterator GetNextParticle() { return last(); }
+
+  protected:
+    StackData& GetStackData() { return static_cast<StackData&>(*this); }
+    const StackData& GetStackData() const { return static_cast<const StackData&>(*this); }
   };
 
 } // namespace corsika::stack
diff --git a/Framework/StackInterface/StackIterator.h b/Framework/StackInterface/StackIterator.h
index 4a951f00a..8137e9b53 100644
--- a/Framework/StackInterface/StackIterator.h
+++ b/Framework/StackInterface/StackIterator.h
@@ -1,14 +1,17 @@
 #ifndef _include_StackIterator_h__
 #define _include_StackIterator_h__
 
+#include <corsika/stack/ParticleBase.h>
+
 #include <iomanip>
 #include <iostream>
 
+class StackData; // forward decl
+
 namespace corsika::stack {
 
-  // forward decl.
-  template <class Stack, class Particle>
-  class StackIteratorInfo;
+  template <typename StackData, template <typename> typename ParticleInterface>
+  class Stack; // forward decl
 
   /**
      @class StackIterator
@@ -34,82 +37,65 @@ namespace corsika::stack {
      for a particle entry at any index fIndex.
   */
 
-  template <typename Stack, typename Particle>
-  class StackIterator : public Particle {
-    friend Stack;
-    friend Particle;
-    friend StackIteratorInfo<Stack, Particle>;
+  template <typename StackData, template <typename> typename ParticleInterface>
+  class StackIteratorInterface
+      : public ParticleInterface<StackIteratorInterface<StackData, ParticleInterface> > {
 
-  private:
-    int fIndex;
+    typedef Stack<StackData, ParticleInterface> StackType;
+    typedef ParticleInterface<StackIteratorInterface<StackData, ParticleInterface> >
+        ParticleInterfaceType;
 
-    //#warning stacks should not be copied because of this:
-    Stack* fData;
+    // friend class ParticleInterface<StackIterator<StackData>>; // to access GetStackData
+    friend class Stack<StackData, ParticleInterface>;  // for access to GetIndex
+    friend class ParticleBase<StackIteratorInterface>; // for access to GetStackData
+
+  private:
+    int fIndex = 0;
+    StackType* fData = 0; // todo is this problematic, when stacks are copied?
 
   public:
-    StackIterator()
-        : fData(0)
-        , fIndex(0) {}
-    StackIterator(Stack& data, const int index)
+    // StackIterator() : fData(0), fIndex(0) { }
+    StackIteratorInterface(StackType& data, const int index)
         : fData(&data)
         , fIndex(index) {}
-    StackIterator(const StackIterator& mit)
+
+  private:
+    StackIteratorInterface(const StackIteratorInterface& mit)
         : fData(mit.fData)
         , fIndex(mit.fIndex) {}
+    StackIteratorInterface& operator=(const StackIteratorInterface& mit) {
+      fData = mit.fData;
+      fIndex = mit.fIndex;
+    }
 
-    StackIterator& operator++() {
+  public:
+    StackIteratorInterface& operator++() {
       ++fIndex;
       return *this;
     }
-    StackIterator operator++(int) {
-      StackIterator tmp(*this);
+    StackIteratorInterface operator++(int) {
+      StackIteratorInterface tmp(*this);
       ++fIndex;
       return tmp;
     }
-    bool operator==(const StackIterator& rhs) { return fIndex == rhs.fIndex; }
-    bool operator!=(const StackIterator& rhs) { return fIndex != rhs.fIndex; }
+    bool operator==(const StackIteratorInterface& rhs) { return fIndex == rhs.fIndex; }
+    bool operator!=(const StackIteratorInterface& rhs) { return fIndex != rhs.fIndex; }
 
-    StackIterator& operator*() { return *this; }
-    const StackIterator& operator*() const { return *this; }
+    ParticleInterfaceType& operator*() {
+      return static_cast<ParticleInterfaceType&>(*this);
+    }
+    const ParticleInterfaceType& operator*() const {
+      return static_cast<const ParticleInterfaceType&>(*this);
+    }
 
   protected:
     int GetIndex() const { return fIndex; }
-    Stack& GetStack() { return *fData; }
-    const Stack& GetStack() const { return *fData; }
-
-    // this is probably not needed rigth now:
-    // inline StackIterator<Stack,Particle>& BaseRef() { return
-    // static_cast<StackIterator<Stack, Particle>&>(*this); } inline const
-    // StackIterator<Stack,Particle>& BaseRef() const { return static_cast<const
-    // StackIterator<Stack, Particle>&>(*this); }
-  };
+    StackType& GetStack() { return *fData; }
+    const StackType& GetStack() const { return *fData; }
+    StackData& GetStackData() { return fData->GetStackData(); }
+    const StackData& GetStackData() const { return fData->GetStackData(); }
 
-  /**
-     @class StackIteratorInfo
-
-     This is the class where custom ...
-     Internal helper class for StackIterator. Document better...
-   */
-
-  template <typename _Stack, typename Particle>
-  class StackIteratorInfo {
-
-    friend Particle;
-
-  private:
-    StackIteratorInfo() {}
-
-  protected:
-    inline _Stack& GetStack() {
-      return static_cast<StackIterator<_Stack, Particle>*>(this)->GetStack();
-    }
-    inline int GetIndex() const {
-      return static_cast<const StackIterator<_Stack, Particle>*>(this)->GetIndex();
-    }
-    inline const _Stack& GetStack() const {
-      return static_cast<const StackIterator<_Stack, Particle>*>(this)->GetStack();
-    }
-  };
+  }; // end class StackIterator
 
 } // namespace corsika::stack
 
diff --git a/Framework/StackInterface/testStackInterface.cc b/Framework/StackInterface/testStackInterface.cc
new file mode 100644
index 000000000..e00ecc183
--- /dev/null
+++ b/Framework/StackInterface/testStackInterface.cc
@@ -0,0 +1,117 @@
+#include <corsika/stack/Stack.h>
+
+#include <iomanip>
+#include <iostream>
+#include <vector>
+
+#define CATCH_CONFIG_MAIN // This tells Catch to provide a main() - only do this in one
+                          // cpp file
+#include <catch2/catch.hpp>
+
+using namespace corsika::stack;
+using namespace std;
+
+// definition of stack-data object
+class StackOneData {
+
+public:
+  // these functions are needed for the Stack interface
+  void Init() {}
+  void Clear() { fData.clear(); }
+  int GetSize() const { return fData.size(); }
+  int GetCapacity() const { return fData.size(); }
+  void Copy(const int i1, const int i2) { fData[i2] = fData[i1]; }
+  void Swap(const int i1, const int i2) {
+    double tmp0 = fData[i1];
+    fData[i1] = fData[i2];
+    fData[i2] = tmp0;
+  }
+
+  // custom data access function
+  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() {
+    if (fData.size() > 0) { fData.pop_back(); }
+  }
+
+  // custom private data section
+private:
+  std::vector<double> fData;
+};
+
+// defintion of a stack-readout object, the iteractor dereference
+// operator will deliver access to these function
+template <typename StackIteratorInterface>
+class ParticleInterface : public ParticleBase<StackIteratorInterface> {
+  //  using ParticleBase<StackIteratorInterface>::Delete;
+  using ParticleBase<StackIteratorInterface>::GetStackData;
+  using ParticleBase<StackIteratorInterface>::GetIndex;
+
+public:
+  void SetData(const double v) { GetStackData().SetData(GetIndex(), v); }
+  double GetData() const { return GetStackData().GetData(GetIndex()); }
+};
+
+TEST_CASE("Stack", "[Stack]") {
+
+  SECTION("StackInterface") {
+
+    // construct a valid Stack object
+    typedef Stack<StackOneData, ParticleInterface> StackTest;
+    StackTest s;
+    s.Init();
+    s.Clear();
+    s.IncrementSize();
+    s.Copy(0, 0);
+    s.Swap(0, 0);
+    s.GetCapacity();
+    REQUIRE(s.GetSize() == 1);
+    s.DecrementSize();
+    REQUIRE(s.GetSize() == 0);
+  }
+
+  SECTION("write") {
+
+    // construct a valid Stack object
+    typedef Stack<StackOneData, ParticleInterface> StackTest;
+    StackTest s;
+  }
+
+  SECTION("read") {
+
+    typedef Stack<StackOneData, ParticleInterface> StackTest;
+    StackTest s;
+    s.NewParticle().SetData(9.9);
+    cout << "kk" << endl;
+    double v = 0;
+    for (auto& p : s) {
+      cout << typeid(p).name() << endl;
+      v += p.GetData();
+    }
+    cout << "k222k" << endl;
+
+    REQUIRE(v == 9.9);
+  }
+
+  SECTION("delete_stack") {
+
+    typedef Stack<StackOneData, ParticleInterface> StackTest;
+    StackTest s;
+    auto p = s.NewParticle();
+    p.SetData(9.9);
+    s.Delete(p);
+  }
+
+  SECTION("delete_particle") {
+
+    typedef Stack<StackOneData, ParticleInterface> StackTest;
+    StackTest s;
+    auto p = s.NewParticle();
+    p.SetData(9.9);
+    p.Delete();
+  }
+}
diff --git a/Framework/Units/PhysicalUnits.h b/Framework/Units/PhysicalUnits.h
index 9aeefb151..91f816163 100644
--- a/Framework/Units/PhysicalUnits.h
+++ b/Framework/Units/PhysicalUnits.h
@@ -17,8 +17,11 @@ namespace phys {
     namespace literals {
       QUANTITY_DEFINE_SCALING_LITERALS(eV, energy_d,
                                        magnitude(corsika::units::constants::eV))
-    }
-  } // namespace units
+
+      // phys::units::quantity<energy_d/mass_d> Joule2Kg = c2; // 1_Joule / 1_kg;
+
+    } // namespace literals
+  }   // namespace units
 } // namespace phys
 
 namespace corsika::units {
diff --git a/Framework/Units/testUnits.cc b/Framework/Units/testUnits.cc
index 080aabe8c..5ea86b37a 100644
--- a/Framework/Units/testUnits.cc
+++ b/Framework/Units/testUnits.cc
@@ -62,6 +62,7 @@ TEST_CASE("PhysicalUnits", "[Units]") {
   SECTION("Formulas") {
     const EnergyType E2 = 20_GeV * 2;
     REQUIRE(E2 == 40_GeV);
+    REQUIRE(E2 / 1_GeV == 40);
 
     const double lgE = log10(E2 / 1_GeV);
     REQUIRE(lgE == Approx(log10(40.)));
diff --git a/Stack/CMakeLists.txt b/Stack/CMakeLists.txt
index 4971f73f8..4cf17a616 100644
--- a/Stack/CMakeLists.txt
+++ b/Stack/CMakeLists.txt
@@ -1,2 +1,3 @@
 
+add_subdirectory (DummyStack)
 add_subdirectory (SuperStupidStack)
diff --git a/Stack/DummyStack/CMakeLists.txt b/Stack/DummyStack/CMakeLists.txt
new file mode 100644
index 000000000..7e29fe3fb
--- /dev/null
+++ b/Stack/DummyStack/CMakeLists.txt
@@ -0,0 +1,29 @@
+
+set (DummyStack_HEADERS DummyStack.h)
+set (DummyStack_NAMESPACE corsika/stack/dummy)
+
+add_library (DummyStack INTERFACE)
+
+CORSIKA_COPY_HEADERS_TO_NAMESPACE (DummyStack ${DummyStack_NAMESPACE} ${DummyStack_HEADERS})
+
+target_link_libraries (
+  DummyStack
+  INTERFACE
+  CORSIKAstackinterface
+  CORSIKAunits
+  CORSIKAparticles
+  )
+
+target_include_directories (
+  DummyStack
+  INTERFACE
+  $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
+  $<INSTALL_INTERFACE:include>
+  )
+
+install (
+  FILES
+  ${DummyStack_HEADERS}
+  DESTINATION
+  include/${DummyStack_NAMESPACE}
+  )
diff --git a/Stack/DummyStack/DummyStack.h b/Stack/DummyStack/DummyStack.h
new file mode 100644
index 000000000..8dd528a7e
--- /dev/null
+++ b/Stack/DummyStack/DummyStack.h
@@ -0,0 +1,61 @@
+#ifndef _include_dummystack_h_
+#define _include_dummystack_h_
+
+#include <corsika/particles/ParticleProperties.h>
+#include <corsika/stack/Stack.h>
+#include <corsika/units/PhysicalUnits.h>
+
+#include <string>
+#include <vector>
+
+namespace corsika::stack {
+
+  namespace dummy {
+
+    /**
+     * Example of a particle object on the stack.
+     */
+
+    template <typename _Stack>
+    class ParticleRead : public StackIteratorInfo<_Stack, ParticleRead<_Stack> > {
+
+      using StackIteratorInfo<_Stack, ParticleRead>::GetIndex;
+      using StackIteratorInfo<_Stack, ParticleRead>::GetStack;
+
+    public:
+    };
+
+    /**
+     *
+     * Memory implementation of the most simple (stupid) particle stack object.
+     */
+
+    class DummyStackImpl {
+
+    public:
+      void Init() {}
+
+      void Clear() {}
+
+      int GetSize() const { return 0; }
+      int GetCapacity() const { return 0; }
+
+      /**
+       *   Function to copy particle at location i2 in stack to i1
+       */
+      void Copy(const int i1, const int i2) {}
+
+    protected:
+      void IncrementSize() {}
+      void DecrementSize() {}
+
+    }; // end class DummyStackImpl
+
+    typedef StackIterator<DummyStackImpl, ParticleRead<DummyStackImpl> > Particle;
+    typedef Stack<DummyStackImpl, Particle> DummyStack;
+
+  } // namespace dummy
+
+} // namespace corsika::stack
+
+#endif
diff --git a/Stack/DummyStack/SuperStupidStack.h b/Stack/DummyStack/SuperStupidStack.h
new file mode 100644
index 000000000..a4e931b39
--- /dev/null
+++ b/Stack/DummyStack/SuperStupidStack.h
@@ -0,0 +1,97 @@
+#ifndef _include_superstupidstack_h_
+#define _include_superstupidstack_h_
+
+#include <corsika/particles/ParticleProperties.h>
+#include <corsika/stack/Stack.h>
+#include <corsika/units/PhysicalUnits.h>
+
+#include <string>
+#include <vector>
+
+namespace corsika::stack {
+
+  namespace super_stupid {
+
+    using corsika::particles::Code;
+    using corsika::units::EnergyType;
+    using corsika::units::operator""_GeV; // literals;
+
+    /**
+     * Example of a particle object on the stack.
+     */
+
+    template <typename _Stack>
+    class ParticleRead : public StackIteratorInfo<_Stack, ParticleRead<_Stack> > {
+
+      using StackIteratorInfo<_Stack, ParticleRead>::GetIndex;
+      using StackIteratorInfo<_Stack, ParticleRead>::GetStack;
+
+    public:
+      void SetId(const Code id) { GetStack().SetId(GetIndex(), id); }
+      void SetEnergy(const EnergyType& e) { GetStack().SetEnergy(GetIndex(), e); }
+
+      Code GetId() const { return GetStack().GetId(GetIndex()); }
+      const EnergyType& GetEnergy() const { return GetStack().GetEnergy(GetIndex()); }
+    };
+
+    /**
+     *
+     * Memory implementation of the most simple (stupid) particle stack object.
+     */
+
+    class SuperStupidStackImpl {
+
+    public:
+      void Init() {}
+
+      void Clear() {
+        fDataE.clear();
+        fDataId.clear();
+      }
+
+      int GetSize() const { return fDataId.size(); }
+      int GetCapacity() const { return fDataId.size(); }
+
+      void SetId(const int i, const Code id) { fDataId[i] = id; }
+      void SetEnergy(const int i, const EnergyType& e) { fDataE[i] = e; }
+
+      const Code GetId(const int i) const { return fDataId[i]; }
+      const EnergyType& GetEnergy(const int i) const { return fDataE[i]; }
+
+      /**
+       *   Function to copy particle at location i2 in stack to i1
+       */
+      void Copy(const int i1, const int i2) {
+        fDataE[i2] = fDataE[i1];
+        fDataId[i2] = fDataId[i1];
+      }
+
+    protected:
+      void IncrementSize() {
+        fDataE.push_back(0_GeV);
+        fDataId.push_back(Code::unknown);
+      }
+      void DecrementSize() {
+        if (fDataE.size() > 0) {
+          fDataE.pop_back();
+          fDataId.pop_back();
+        }
+      }
+
+    private:
+      /// the actual memory to store particle data
+
+      std::vector<Code> fDataId;
+      std::vector<EnergyType> fDataE;
+
+    }; // end class SuperStupidStackImpl
+
+    typedef StackIterator<SuperStupidStackImpl, ParticleRead<SuperStupidStackImpl> >
+        Particle;
+    typedef Stack<SuperStupidStackImpl, Particle> SuperStupidStack;
+
+  } // namespace super_stupid
+
+} // namespace corsika::stack
+
+#endif
diff --git a/Stack/SuperStupidStack/SuperStupidStack.h b/Stack/SuperStupidStack/SuperStupidStack.h
index 1f02e9087..2ab8ea790 100644
--- a/Stack/SuperStupidStack/SuperStupidStack.h
+++ b/Stack/SuperStupidStack/SuperStupidStack.h
@@ -1,13 +1,12 @@
 #ifndef _include_superstupidstack_h_
 #define _include_superstupidstack_h_
 
-#include <string>
-#include <vector>
-
 #include <corsika/particles/ParticleProperties.h>
 #include <corsika/stack/Stack.h>
 #include <corsika/units/PhysicalUnits.h>
 
+#include <vector>
+
 namespace corsika::stack {
 
   namespace super_stupid {
@@ -20,18 +19,18 @@ namespace corsika::stack {
      * Example of a particle object on the stack.
      */
 
-    template <typename _Stack>
-    class ParticleRead : public StackIteratorInfo<_Stack, ParticleRead<_Stack> > {
+    template <typename StackIteratorInterface>
+    class ParticleInterface : public ParticleBase<StackIteratorInterface> {
 
-      using StackIteratorInfo<_Stack, ParticleRead>::GetIndex;
-      using StackIteratorInfo<_Stack, ParticleRead>::GetStack;
+      using ParticleBase<StackIteratorInterface>::GetStackData;
+      using ParticleBase<StackIteratorInterface>::GetIndex;
 
     public:
-      void SetId(const Code id) { GetStack().SetId(GetIndex(), id); }
-      void SetEnergy(const EnergyType& e) { GetStack().SetEnergy(GetIndex(), e); }
+      void SetPID(const Code id) { GetStackData().SetPID(GetIndex(), id); }
+      void SetEnergy(const EnergyType& e) { GetStackData().SetEnergy(GetIndex(), e); }
 
-      Code GetId() const { return GetStack().GetId(GetIndex()); }
-      const EnergyType& GetEnergy() const { return GetStack().GetEnergy(GetIndex()); }
+      Code GetPID() const { return GetStackData().GetPID(GetIndex()); }
+      EnergyType GetEnergy() const { return GetStackData().GetEnergy(GetIndex()); }
     };
 
     /**
@@ -42,51 +41,63 @@ namespace corsika::stack {
     class SuperStupidStackImpl {
 
     public:
+      void Init() {}
+
       void Clear() {
         fDataE.clear();
-        fDataId.clear();
+        fDataPID.clear();
       }
 
-      int GetSize() const { return fDataId.size(); }
-      int GetCapacity() const { return fDataId.size(); }
+      int GetSize() const { return fDataPID.size(); }
+      int GetCapacity() const { return fDataPID.size(); }
 
-      void SetId(const int i, const Code id) { fDataId[i] = id; }
-      void SetEnergy(const int i, const EnergyType& e) { fDataE[i] = e; }
+      void SetPID(const int i, const Code id) { fDataPID[i] = id; }
+      void SetEnergy(const int i, const EnergyType e) { fDataE[i] = e; }
 
-      const Code GetId(const int i) const { return fDataId[i]; }
-      const EnergyType& GetEnergy(const int i) const { return fDataE[i]; }
+      Code GetPID(const int i) const { return fDataPID[i]; }
+      EnergyType GetEnergy(const int i) const { return fDataE[i]; }
 
       /**
        *   Function to copy particle at location i2 in stack to i1
        */
       void Copy(const int i1, const int i2) {
         fDataE[i2] = fDataE[i1];
-        fDataId[i2] = fDataId[i1];
+        fDataPID[i2] = fDataPID[i1];
+      }
+
+      /**
+       *   Function to copy particle at location i2 in stack to i1
+       */
+      void Swap(const int i1, const int i2) {
+        EnergyType tE = fDataE[i2];
+        Code tC = fDataPID[i2];
+        fDataE[i2] = fDataE[i1];
+        fDataPID[i2] = fDataPID[i1];
+        fDataE[i1] = tE;
+        fDataPID[i1] = tC;
       }
 
     protected:
       void IncrementSize() {
         fDataE.push_back(0_GeV);
-        fDataId.push_back(Code::unknown);
+        fDataPID.push_back(Code::unknown);
       }
       void DecrementSize() {
         if (fDataE.size() > 0) {
           fDataE.pop_back();
-          fDataId.pop_back();
+          fDataPID.pop_back();
         }
       }
 
     private:
       /// the actual memory to store particle data
 
-      std::vector<Code> fDataId;
+      std::vector<Code> fDataPID;
       std::vector<EnergyType> fDataE;
 
     }; // end class SuperStupidStackImpl
 
-    typedef StackIterator<SuperStupidStackImpl, ParticleRead<SuperStupidStackImpl> >
-        Particle;
-    typedef Stack<SuperStupidStackImpl, Particle> SuperStupidStack;
+    typedef Stack<SuperStupidStackImpl, ParticleInterface> SuperStupidStack;
 
   } // namespace super_stupid
 
-- 
GitLab