- some majro updates to the stack-interface

- revised version of process interface
- added tests for stacks, processes, cascade

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)
 

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;

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() {

Framework/CMakeLists.txt

@@ -5,3 +5,4 @@ add_subdirectory (Particles)
 add_subdirectory (Logging)
 add_subdirectory (StackInterface)
 add_subdirectory (ProcessSequence)
+add_subdirectory (Cascade)

Framework/Cascade/CMakeLists.txt

+
+# 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
+  )

Framework/Cascade/Cascade.cc

+#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()) {

Framework/Cascade/Cascade.h

 #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

Framework/Cascade/testCascade.cc

+#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;
+    }
+  }
+}

Framework/Geometry/CMakeLists.txt

@@ -13,6 +13,7 @@ set (
   Helix.h
   BaseVector.h
   QuantityVector.h
+  LineTrajectory.h
   )
 
 set (

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:

Framework/Geometry/LineTrajectory.h

 #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

Framework/ProcessSequence/CMakeLists.txt

 
 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
+  )

Framework/ProcessSequence/ProcessSequence.h

 #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) {

Framework/ProcessSequence/testProcessSequence.cc

+#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") {}
+}

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)
+

Framework/StackInterface/ParticleBase.h

+#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

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

Framework/StackInterface/StackIterator.h

 #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
 

Framework/StackInterface/testStackInterface.cc

+#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();
+  }
+}

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 {