transformed Stack to keep deleted particle, until they are pureged

Documentation/Examples/boundary_example.cc

@@ -81,6 +81,9 @@ private:
 // The example main program for a particle cascade
 //
 int main() {
+
+  std::cout << "boundary_example" << std::endl;
+
   feenableexcept(FE_INVALID);
   // initialize random number sequence(s)
   random::RNGManager::GetInstance().RegisterRandomStream("cascade");

Documentation/Examples/cascade_example.cc

@@ -56,6 +56,8 @@ using namespace corsika::units::si;
 //
 int main() {
 
+  std::cout << "cascade_example" << std::endl;
+
   const LengthType height_atmosphere = 112.8_km;
 
   feenableexcept(FE_INVALID);

Documentation/Examples/cascade_proton_example.cc

@@ -58,6 +58,9 @@ using namespace corsika::units::si;
 // The example main program for a particle cascade
 //
 int main() {
+
+  std::cout << "cascade_proton_example" << std::endl;
+
   feenableexcept(FE_INVALID);
   // initialize random number sequence(s)
   random::RNGManager::GetInstance().RegisterRandomStream("cascade");

Documentation/Examples/geometry_example.cc

@@ -21,6 +21,9 @@ using namespace corsika::geometry;
 using namespace corsika::units::si;
 
 int main() {
+
+  std::cout << "geometry_example" << std::endl;    
+  
   // define the root coordinate system
   geometry::CoordinateSystem& root =
       geometry::RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();

Documentation/Examples/helix_example.cc

@@ -20,6 +20,9 @@ using namespace corsika::geometry;
 using namespace corsika::units::si;
 
 int main() {
+
+  std::cout << "helix_example" << std::endl;
+  
   geometry::CoordinateSystem& root =
       geometry::RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
 

Documentation/Examples/particle_list_example.cc

@@ -25,6 +25,9 @@ using namespace std;
 // The example main program for a particle list
 //
 int main() {
+
+  std::cout << "particle_list_example" << std::endl;
+  
   cout << "------------------------------------------"
        << "particles in CORSIKA"
        << "------------------------------------------" << endl;

Documentation/Examples/stack_example.cc

@@ -36,7 +36,7 @@ void fill(corsika::stack::super_stupid::SuperStupidStack& s) {
 }
 
 void read(corsika::stack::super_stupid::SuperStupidStack& s) {
-  assert(s.GetSize() == 11); // stack has 11 particles
+  assert(s.getEntries() == 11); // stack has 11 particles
 
   HEPEnergyType total_energy;
   int i = 0;
@@ -49,6 +49,9 @@ void read(corsika::stack::super_stupid::SuperStupidStack& s) {
 }
 
 int main() {
+
+  std::cout << "stack_example" << std::endl;
+  
   corsika::stack::super_stupid::SuperStupidStack s;
   fill(s);
   read(s);

Documentation/Examples/staticsequence_example.cc

@@ -100,6 +100,9 @@ void modular() {
 }
 
 int main() {
+
+  std::cout << "staticsequence_example" << std::endl;
+  
   modular();
   return 0;
 }

Documentation/Examples/stopping_power.cc

@@ -33,6 +33,9 @@ using namespace corsika::units::si;
 // This example demonstrates the energy loss of muons as function of beta*gamma (=p/m)
 //
 int main() {
+
+  std::cout << "stopping_power" << std::endl;
+  
   feenableexcept(FE_INVALID);
 
   // setup environment, geometry

Documentation/Examples/vertical_EAS.cc

@@ -55,7 +55,7 @@ using namespace corsika::environment;
 using namespace std;
 using namespace corsika::units::si;
 
-void registerRandomStreams() {
+void registerRandomStreams(const int seed) {
   random::RNGManager::GetInstance().RegisterRandomStream("cascade");
   random::RNGManager::GetInstance().RegisterRandomStream("qgsjet");
   random::RNGManager::GetInstance().RegisterRandomStream("sibyll");
@@ -63,17 +63,28 @@ void registerRandomStreams() {
   random::RNGManager::GetInstance().RegisterRandomStream("urqmd");
   random::RNGManager::GetInstance().RegisterRandomStream("proposal");
 
-  random::RNGManager::GetInstance().SeedAll();
+  if (seed==0)
+    random::RNGManager::GetInstance().SeedAll();
+  else 
+    random::RNGManager::GetInstance().SeedAll(seed);
 }
 
 int main(int argc, char** argv) {
-  if (argc != 4) {
-    std::cerr << "usage: vertical_EAS <A> <Z> <energy/GeV>" << std::endl;
+
+  std::cout << "vertical_EAS" << std::endl;
+
+  if (argc < 4) {
+    std::cerr << "usage: vertical_EAS <A> <Z> <energy/GeV> [seed]" << std::endl;
+    std::cerr << "       if no seed is given, a random seed is chosen" << std::endl;
     return 1;
   }
   feenableexcept(FE_INVALID);
+
+  int seed = 0;
+  if (argc>4)
+    seed = std::stoi(std::string(argv[4]));
   // initialize random number sequence(s)
-  registerRandomStreams();
+  registerRandomStreams(seed);
 
   // setup environment, geometry
   using EnvType = Environment<setup::IEnvironmentModel>;

Framework/Cascade/Cascade.h

@@ -111,7 +111,9 @@ namespace corsika::cascade {
       while (!fStack.IsEmpty()) {
         while (!fStack.IsEmpty()) {
           auto pNext = fStack.GetNextParticle();
-          std::cout << "========= next: " << pNext.GetPID() << std::endl;
+          std::cout << "========= next: pid=" << pNext.GetPID()
+                    << ", stack entries=" << fStack.getEntries()
+                    << ", stack deleted=" << fStack.getDeleted() << std::endl;
           Step(pNext);
           std::cout << "========= stack ============" << std::endl;
           fProcessSequence.DoStack(fStack);
@@ -253,7 +255,7 @@ namespace corsika::cascade {
             assert(min_distance == distance_decay);
             decay(vParticle, projectile);
             // make sure particle actually did decay if it should have done so
-            if (secondaries.GetSize() == 1 &&
+            if (secondaries.getSize() == 1 &&
                 projectile.GetPID() == secondaries.GetNextParticle().GetPID())
               throw std::runtime_error(
                   fmt::format("Cascade: {} decayed into itself!",

Framework/Cascade/testCascade.cc

@@ -106,11 +106,10 @@ public:
       if (E < fEcrit) {
         p.Delete();
         fCount++;
-      } else {
-        ++p; // next particle
       }
+      ++p; // next particle
     }
-    cout << "ProcessCut::DoSecondaries size=" << vS.GetSize() << " count=" << fCount
+    cout << "ProcessCut::DoSecondaries size=" << vS.getEntries() << " count=" << fCount
          << endl;
     return EProcessReturn::eOk;
   }
@@ -162,7 +161,7 @@ TEST_CASE("Cascade", "[Cascade]") {
   SECTION("forced interaction") {
     EAS.SetNodes();
     EAS.forceInteraction();
-    CHECK(stack.GetSize() == 2);
+    CHECK(stack.getEntries() == 2);
     CHECK(split.GetCalls() == 1);
   }
 }

Framework/StackInterface/ParticleBase.h

@@ -63,6 +63,11 @@ namespace corsika::stack {
      */
     void Delete() { GetIterator().GetStack().Delete(GetIterator()); }
 
+    /**
+     * Method to retrieve the status of the Particle. Is it already deleted? Or not.
+     */
+    bool isDeleted() const { return GetIterator().GetStack().isDeleted(GetIterator()); }
+
     /**
      * Add a secondary particle based on *this on the stack @param
      * args is a variadic list of input data that has to match the

Framework/StackInterface/SecondaryView.h

@@ -63,6 +63,7 @@ namespace corsika::stack {
      * Helper type for inside this class
      */
     using InnerStackType = Stack<StackDataType&, ParticleInterface>;
+    using InnerStackType::getDeleted;
 
     /**
      * @name We need this "special" types with non-reference StackData for
@@ -104,6 +105,11 @@ namespace corsika::stack {
     template <typename... Args>
     SecondaryView(Args... args) = delete;
 
+  private:
+    InnerStackTypeValue& innerstack_;
+    unsigned int fProjectileIndex;
+    std::vector<unsigned int> fIndices;
+
   public:
     /**
        SecondaryView can only be constructed passing it a valid
@@ -111,6 +117,7 @@ namespace corsika::stack {
      **/
     SecondaryView(StackIteratorValue& vI)
         : Stack<StackDataType&, ParticleInterface>(vI.GetStackData())
+        , innerstack_(vI.GetStack())
         , fProjectileIndex(vI.GetIndex()) {}
 
     StackIterator GetProjectile() {
@@ -128,8 +135,9 @@ namespace corsika::stack {
     auto AddSecondary(StackIterator& proj, const Args... v) {
       // make space on stack
       InnerStackType::GetStackData().IncrementSize();
+      innerstack_.deleted_.push_back(false);
       // get current number of secondaries on stack
-      const unsigned int idSec = GetSize();
+      const unsigned int idSec = getSize();
       // determine index on (inner) stack where new particle will be located
       const unsigned int index = InnerStackType::GetStackData().GetSize() - 1;
       fIndices.push_back(index);
@@ -141,27 +149,65 @@ namespace corsika::stack {
     /**
      * overwrite Stack::GetSize to return actual number of secondaries
      */
-    unsigned int GetSize() const { return fIndices.size(); }
+    unsigned int getSize() const { return fIndices.size(); }
+    unsigned int getEntries() const { return getSize() - getDeleted(); }
+    bool IsEmpty() const { return getEntries() == 0; }
 
     /**
      * @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
+     * 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); }
+    StackIterator begin() {
+      unsigned int i = 0;
+      for (; i < getSize(); ++i) {
+        if (!isDeleted(i)) break;
+      }
+      return StackIterator(*this, i + 1);
+    }
+    auto end() { return StackIterator(*this, getSize() + 1); }
+    auto last() {
+      unsigned int i = 0;
+      for (; i < getSize(); ++i) {
+        if (!isDeleted(getSize() - 1 - i)) break;
+      }
+      return StackIterator(*this, getSize() - 1 - i + 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 begin() const {
+      unsigned int i = 0;
+      for (; i < getSize(); ++i) {
+        if (!isDeleted(i)) break;
+      }
+      return ConstStackIterator(*this, i + 1);
+    }
+    auto end() const { return ConstStackIterator(*this, getSize() + 1); }
+    auto last() const {
+      unsigned int i = 0;
+      for (; i < getSize(); ++i) {
+        if (!isDeleted(getSize() - 1 - i)) break;
+      }
+      return ConstStackIterator(*this, getSize() - 1 - i + 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); }
+    auto cbegin() const {
+      unsigned int i = 0;
+      for (; i < getSize(); ++i) {
+        if (!isDeleted(i)) break;
+      }
+      return ConstStackIterator(*this, i + 1);
+    }
+    auto cend() const { return ConstStackIterator(*this, getSize()); }
+    auto clast() const {
+      unsigned int i = 0;
+      for (; i < getSize(); ++i) {
+        if (!isDeleted(getSize() - 1 - i)) break;
+      }
+      return ConstStackIterator(*this, getSize() - 1 - i + 1);
+    }
     /// @}
 
     /**
@@ -177,42 +223,91 @@ namespace corsika::stack {
      *
      */
     void Delete(StackIterator p) {
-      if (IsEmpty()) { /* error */
+      if (IsEmpty()) { /*error*/
         throw std::runtime_error("Stack, cannot delete entry since size is zero");
       }
-      const int innerSize = InnerStackType::GetSize();
-      const int innerIndex = GetIndexFromIterator(p.GetIndex());
-      if (innerIndex < innerSize - 1)
-        InnerStackType::GetStackData().Copy(innerSize - 1,
-                                            GetIndexFromIterator(p.GetIndex()));
-      DeleteLast();
+      if (isDeleted(p.GetIndex() - 1)) { /*error*/
+        throw std::runtime_error("Stack, cannot delete entry since already deleted");
+      }
+      innerstack_.Delete(GetIndexFromIterator(p.GetIndex()));
+      InnerStackType::nDeleted_++; // also count in SecondaryView
     }
 
     /**
      * need overwrite Stack::Delete, since we want to call SecondaryView::DeleteLast
      */
-    void Delete(ParticleInterfaceType p) { Delete(p.GetIterator()); }
+    // void Delete(ParticleInterfaceType p) { Delete(p.GetIterator()); }
 
     /**
-     * delete last particle on stack by decrementing stack size
+     * return next particle from stack, need to overwrtie Stack::GetNextParticle to get
+     * right reference
      */
-    void DeleteLast() {
-      fIndices.pop_back();
-      InnerStackType::GetStackData().DecrementSize();
+    StackIterator GetNextParticle() {
+      while (purgeLastIfDeleted()) {}
+      return last();
     }
 
     /**
-     * return next particle from stack, need to overwrtie Stack::GetNextParticle to get
-     * right reference
+     * check if this particle was already deleted
+     *
+     * need to re-implement for SecondaryView since StackIterator types are a bit
+     * different
+     */
+    bool isDeleted(const StackIterator& p) const { return isDeleted(p.GetIndex() - 1); }
+    bool isDeleted(const ConstStackIterator& p) const {
+      return isDeleted(p.GetIndex() - 1);
+    }
+    /**
+     * delete this particle
      */
-    StackIterator GetNextParticle() { return last(); }
+    bool isDeleted(const ParticleInterfaceType& p) const {
+      return isDeleted(p.GetIterator());
+    }
+
+    /**
+     * Function to ultimatively remove the last entry from the stack,
+     * if it was marked as deleted before. If this is not the case,
+     * the function will just return false and do nothing.
+     */
+    bool purgeLastIfDeleted() {
+      if (!isDeleted(getSize() - 1))
+        return false; // the last particle is not marked for deletion. Do nothing.
+      innerstack_.purge(GetIndexFromIterator(getSize()));
+      InnerStackType::nDeleted_--;
+      fIndices.pop_back();
+      return true;
+    }
 
     /**
-     * check if there are no further particles on stack
+     * Function to ultimatively remove all entries from the stack
+     * marked as deleted.
+     *
+     * Careful: this will re-order the entries on the stack, since
+     * "gaps" in the stack are filled with entries from the back
+     * (copied).
      */
-    bool IsEmpty() { return GetSize() == 0; }
+    void purge() {
+      unsigned int iStack = 0;
+      unsigned int size = getSize();
+      while (iStack < size) {
+        if (isDeleted(iStack)) {
+          innerstack_.purge(iStack);
+          fIndices.erase(fIndices.begin() + iStack);
+        }
+        size = getSize();
+        iStack++;
+      }
+      InnerStackType::nDeleted_ = 0;
+    }
 
   protected:
+    // forward to inner stack
+    // this also checks the allowed bounds of 'i'
+    bool isDeleted(unsigned int i) const {
+      if (i >= fIndices.size()) return false;
+      return innerstack_.isDeleted(GetIndexFromIterator(i + 1));
+    }
+
     /**
      * We only want to 'see' secondaries indexed in fIndices. In this
      * function the conversion form iterator-index to stack-index is
@@ -222,10 +317,6 @@ namespace corsika::stack {
       if (vI == 0) return fProjectileIndex;
       return fIndices[vI - 1];
     }
-
-  private:
-    unsigned int fProjectileIndex;
-    std::vector<unsigned int> fIndices;
   };
 
   /*

Framework/StackInterface/Stack.h

@@ -10,11 +10,12 @@
 
 #include <corsika/stack/StackIteratorInterface.h>
 // must be after StackIteratorInterface
-#include <corsika/stack/SecondaryView.h>
+// #include <corsika/stack/SecondaryView.h>
 #include <corsika/utl/MetaProgramming.h>
 
 #include <stdexcept>
 #include <type_traits>
+#include <vector>
 
 /**
    All classes around management of particles on a stack.
@@ -51,11 +52,15 @@ namespace corsika::stack {
      loops, ranges, etc.
    */
 
-  template <typename StackDataType, template <typename> typename ParticleInterface>
+  template <typename TStackData, template <typename> typename TParticleInterface>
   class Stack {
-    using StackDataValueType = std::remove_reference_t<StackDataType>;
+    using StackDataValueType = std::remove_reference_t<TStackData>;
 
-    StackDataType fData; ///< this in general holds all the data and can be quite big
+  private:
+    TStackData data_; ///< this in general holds all the data and can be quite big
+    std::vector<bool> deleted_; ///< bit field to flag deleted entries
+  protected:
+    unsigned int nDeleted_ = 0;
 
   private:
     Stack(Stack&) = delete; ///< since Stack can be very big, we don't want to copy it
@@ -64,48 +69,52 @@ namespace corsika::stack {
 
   public:
     /**
-     * if StackDataType is a reference member we *HAVE* to initialize
+     * if TStackData is a reference member we *HAVE* to initialize
      * it in the constructor, this is typically needed for SecondaryView
      */
-    template <typename _ = StackDataType, typename = utl::enable_if<std::is_reference<_>>>
-    Stack(StackDataType vD)
-        : fData(vD) {}
+    template <typename _ = TStackData, typename = utl::enable_if<std::is_reference<_>>>
+    Stack(TStackData vD)
+        : data_(vD)
+        , deleted_(std::vector<bool>(data_.GetSize(), false))
+        , nDeleted_(0) {}
 
     /**
      * This constructor takes any argument and passes it on to the
-     * StackDataType user class. If the user did not provide a suited
+     * TStackData user class. If the user did not provide a suited
      * constructor this will fail with an error message.
      *
      * Furthermore, this is disabled with enable_if for SecondaryView
      * stacks, where the inner data container is always a reference
      * and cannot be initialized here.
      */
-    template <typename... Args, typename _ = StackDataType,
+    template <typename... Args, typename _ = TStackData,
               typename = utl::disable_if<std::is_reference<_>>>
     Stack(Args... args)
-        : fData(args...) {}
+        : data_(args...)
+        , deleted_(std::vector<bool>(data_.GetSize(), false))
+        , nDeleted_(0) {}
 
   public:
-    typedef StackDataType
+    typedef TStackData
         StackImpl; ///< this is the type of the user-provided data structure
 
     template <typename SI>
-    using PIType = ParticleInterface<SI>;
+    using PIType = TParticleInterface<SI>;
 
     /**
      * Via the StackIteratorInterface and ConstStackIteratorInterface
      * specialization, the type of the StackIterator
      * template class is declared for a particular stack data
      * object. Using CRTP, this also determines the type of
-     * ParticleInterface template class simultaneously.
+     * TParticleInterface template class simultaneously.
      */
     using StackIterator =
-        StackIteratorInterface<StackDataValueType, ParticleInterface, Stack>;
+        StackIteratorInterface<StackDataValueType, TParticleInterface, Stack>;
     using ConstStackIterator =
-        ConstStackIteratorInterface<StackDataValueType, ParticleInterface, Stack>;
+        ConstStackIteratorInterface<StackDataValueType, TParticleInterface, Stack>;
 
     /**
-     * this is the full type of the user-declared ParticleInterface
+     * this is the full type of the user-declared TParticleInterface
      */
     using ParticleInterfaceType = typename StackIterator::ParticleInterfaceType;
     /**
@@ -115,21 +124,25 @@ namespace corsika::stack {
     using ParticleType = StackIterator;
 
     // friends are needed since they need access to protected members
-    friend class StackIteratorInterface<StackDataValueType, ParticleInterface, Stack>;
-    friend class ConstStackIteratorInterface<StackDataValueType, ParticleInterface,
+    friend class StackIteratorInterface<StackDataValueType, TParticleInterface, Stack>;
+    friend class ConstStackIteratorInterface<StackDataValueType, TParticleInterface,
                                              Stack>;
+    friend class SecondaryView<StackDataValueType, TParticleInterface>;
 
   public:
     /**
-     * @name Most generic proxy methods for StackDataType fData
+     * @name Most generic proxy methods for TStackData data_
      * @{
      */
-    unsigned int GetCapacity() const { return fData.GetCapacity(); }
-    unsigned int GetSize() const { return fData.GetSize(); }
+    unsigned int GetCapacity() const { return data_.GetCapacity(); }
+    unsigned int getDeleted() const { return nDeleted_; }
+    unsigned int getEntries() const { return getSize() - getDeleted(); }
 
     template <typename... Args>
-    auto Clear(Args... args) {
-      return fData.Clear(args...);
+    void Clear(Args... args) {
+      data_.Clear(args...);
+      deleted_ = std::vector<bool>(data_.GetSize(), false);
+      nDeleted_ = 0;
     }
     ///@}
 
@@ -138,26 +151,68 @@ namespace corsika::stack {
      * @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); }
+    StackIterator begin() {
+      unsigned int i = 0;
+      for (; i < getSize(); ++i) {
+        if (!deleted_[i]) break;
+      }
+      return StackIterator(*this, i);
+    }
+    StackIterator end() { return StackIterator(*this, getSize()); }
+    StackIterator last() {
+      unsigned int i = 0;
+      for (; i < getSize(); ++i) {
+        if (!deleted_[getSize() - 1 - i]) break;
+      }
+      return StackIterator(*this, getSize() - 1 - i);
+    }
 
-    ConstStackIterator begin() const { return ConstStackIterator(*this, 0); }
-    ConstStackIterator end() const { return ConstStackIterator(*this, GetSize()); }
-    ConstStackIterator last() const { return ConstStackIterator(*this, GetSize() - 1); }
+    ConstStackIterator begin() const {
+      unsigned int i = 0;
+      for (; i < getSize(); ++i) {
+        if (!deleted_[i]) break;
+      }
+      return ConstStackIterator(*this, i);
+    }
+    ConstStackIterator end() const { return ConstStackIterator(*this, getSize()); }
+    ConstStackIterator last() const {
+      unsigned int i = 0;
+      for (; i < getSize(); ++i) {
+        if (!deleted_[getSize() - 1 - i]) break;
+      }
+      return ConstStackIterator(*this, getSize() - 1 - i);
+    }
 
-    ConstStackIterator cbegin() const { return ConstStackIterator(*this, 0); }
-    ConstStackIterator cend() const { return ConstStackIterator(*this, GetSize()); }
-    ConstStackIterator clast() const { return ConstStackIterator(*this, GetSize() - 1); }
+    ConstStackIterator cbegin() const {
+      unsigned int i = 0;
+      for (; i < getSize(); ++i) {
+        if (!deleted_[i]) break;
+      }
+      return ConstStackIterator(*this, 0);
+    }
+    ConstStackIterator cend() const { return ConstStackIterator(*this, getSize()); }
+    ConstStackIterator clast() const {
+      unsigned int i = 0;
+      for (; i < getSize(); ++i) {
+        if (!deleted_[getSize() - 1 - i]) break;
+      }
+      return ConstStackIterator(*this, getSize() - 1 - i);
+    }
     /// @}
 
+    StackIterator GetNextParticle() {
+      while (purgeLastIfDeleted()) {}
+      return last();
+    }
+
     /**
      * increase stack size, create new particle at end of stack
      */
     template <typename... Args>
     StackIterator AddParticle(const Args... v) {
-      fData.IncrementSize();
-      return StackIterator(*this, GetSize() - 1, v...);
+      data_.IncrementSize();
+      deleted_.push_back(false);
+      return StackIterator(*this, getSize() - 1, v...);
     }
 
     /**
@@ -166,32 +221,44 @@ namespace corsika::stack {
      */
     template <typename... Args>
     StackIterator AddSecondary(StackIterator& parent, const Args... v) {
-      fData.IncrementSize();
-      return StackIterator(*this, GetSize() - 1, parent, v...);
+      data_.IncrementSize();
+      deleted_.push_back(false);
+      return StackIterator(*this, getSize() - 1, parent, v...);
     }
 
     void Swap(StackIterator a, StackIterator b) {
-      fData.Swap(a.GetIndex(), b.GetIndex());
+      data_.Swap(a.GetIndex(), b.GetIndex());
+      std::swap(deleted_[a.GetIndex()], deleted_[b.GetIndex()]);
     }
     void Swap(ConstStackIterator a, ConstStackIterator b) {
-      fData.Swap(a.GetIndex(), b.GetIndex());
+      data_.Swap(a.GetIndex(), b.GetIndex());
+      std::swap(deleted_[a.GetIndex()], deleted_[b.GetIndex()]);
     }
     void Copy(StackIterator a, StackIterator b) {
-      fData.Copy(a.GetIndex(), b.GetIndex());
+      data_.Copy(a.GetIndex(), b.GetIndex());
+      if (deleted_[b.GetIndex()] && !deleted_[a.GetIndex()]) nDeleted_--;
+      if (!deleted_[b.GetIndex()] && deleted_[a.GetIndex()]) nDeleted_++;
+      deleted_[b.GetIndex()] = deleted_[a.GetIndex()];
     }
     void Copy(ConstStackIterator a, StackIterator b) {
-      fData.Copy(a.GetIndex(), b.GetIndex());
+      data_.Copy(a.GetIndex(), b.GetIndex());
+      if (deleted_[b.GetIndex()] && !deleted_[a.GetIndex()]) nDeleted_--;
+      if (!deleted_[b.GetIndex()] && deleted_[a.GetIndex()]) nDeleted_++;
+      deleted_[b.GetIndex()] = deleted_[a.GetIndex()];
     }
 
     /**
      * delete this particle
      */
+  public:
     void Delete(StackIterator p) {
-      if (GetSize() == 0) { /*error*/
+      if (IsEmpty()) { /*error*/
         throw std::runtime_error("Stack, cannot delete entry since size is zero");
       }
-      if (p.GetIndex() < GetSize() - 1) fData.Copy(GetSize() - 1, p.GetIndex());
-      DeleteLast();
+      if (deleted_[p.GetIndex()]) { /*error*/
+        throw std::runtime_error("Stack, cannot delete entry since already deleted");
+      }
+      Delete(p.GetIndex());
     }
     /**
      * delete this particle
@@ -199,35 +266,99 @@ namespace corsika::stack {
     void Delete(ParticleInterfaceType p) { Delete(p.GetIterator()); }
 
     /**
-     * delete last particle on stack by decrementing stack size
+     * check if there are no further non-deleted particles on stack
      */
-    void DeleteLast() { fData.DecrementSize(); }
+    bool IsEmpty() { return getEntries() == 0; }
 
     /**
-     * check if there are no further particles on stack
+     * check if this particle was already deleted
      */
-    bool IsEmpty() { return GetSize() == 0; }
+  public:
+    bool isDeleted(const StackIterator& p) { return isDeleted(p.GetIndex()); }
+    bool isDeleted(const ConstStackIterator& p) const { return isDeleted(p.GetIndex()); }
+    bool isDeleted(const ParticleInterfaceType& p) { return isDeleted(p.GetIterator()); }
 
     /**
-     * return next particle from stack
+     * Function to ultimatively remove the last entry from the stack,
+     * if it was marked as deleted before. If this is not the case,
+     * the function will just return false and do nothing.
      */
-    StackIterator GetNextParticle() { return last(); }
+    bool purgeLastIfDeleted() {
+      if (!deleted_.back())
+        return false; // the last particle is not marked for deletion. Do nothing.
+      data_.DecrementSize();
+      nDeleted_--;
+      deleted_.pop_back();
+      return true;
+    }
+
+    /**
+     * Function to ultimatively remove all entries from the stack
+     * marked as deleted.
+     *
+     * Careful: this will re-order the entries on the stack, since
+     * "gaps" in the stack are filled with entries from the back
+     * (copied).
+     */
+    void purge() {
+      unsigned int iStackFront = 0;
+      unsigned int iStackBack = getSize() - 1;
+      for (unsigned int iDeleted = 0; iDeleted < getDeleted(); ++iDeleted) {
+        // search first delete entry on stack
+        while (!deleted_[iStackFront]) { iStackFront++; }
+        // search for last non-deleted particle on stack
+        while (deleted_[iStackBack]) { iStackBack--; }
+        // copy entry from iStackBack to iStackFront
+        data_.Copy(iStackBack, iStackFront);
+        data_.DecrementSize();
+      }
+      deleted_.clear();
+      nDeleted_ = 0;
+    }
 
   protected:
+    unsigned int getSize() const { return data_.GetSize(); }
+
+    bool isDeleted(unsigned int i) const {
+      if (i >= deleted_.size()) return false;
+      return deleted_.at(i);
+    }
+
+    void Delete(unsigned int i) {
+      deleted_[i] = true;
+      nDeleted_++;
+    }
+
+    /**
+     * will remove from storage the element i. This is a helper
+     * function for SecondaryView.
+     */
+    void purge(unsigned int i) {
+      unsigned int iStackBack = getSize() - 1;
+      // search for last non-deleted particle on stack
+      while (deleted_[iStackBack]) { iStackBack--; }
+      // copy entry from iStackBack to iStackFront
+      data_.Copy(iStackBack, i);
+      if (deleted_[i]) nDeleted_--;
+      deleted_[i] = deleted_[iStackBack];
+      data_.DecrementSize();
+      deleted_.pop_back();
+    }
+
     /**
      * Function to perform eventual transformation from
      * StackIterator::GetIndex() to index in data stored in
-     * StackDataType fData. By default (and in almost all cases) this
+     * TStackData data_. By default (and in almost all cases) this
      * should just be identiy. See class SecondaryView for an alternative implementation.
      */
     unsigned int GetIndexFromIterator(const unsigned int vI) const { return vI; }
 
     /**
-     * @name Return reference to StackDataType object fData for data access
+     * @name Return reference to TStackData object data_ for data access
      * @{
      */
-    StackDataValueType& GetStackData() { return fData; }
-    const StackDataValueType& GetStackData() const { return fData; }
+    StackDataValueType& GetStackData() { return data_; }
+    const StackDataValueType& GetStackData() const { return data_; }
     ///@}
   };
 

Framework/StackInterface/StackIteratorInterface.h

@@ -17,10 +17,10 @@ namespace corsika::history {
 
 namespace corsika::stack {
 
-  template <typename StackDataType, template <typename> typename ParticleInterface>
+  template <typename TStackData, template <typename> typename TParticleInterface>
   class Stack; // forward decl
 
-  template <typename StackDataType, template <typename> typename ParticleInterface>
+  template <typename TStackData, template <typename> typename TParticleInterface>
   class SecondaryView; // forward decl
 
   /**
@@ -39,54 +39,50 @@ namespace corsika::stack {
      The template argument Stack determines the type of Stack object
      the data is stored in. A pointer to the Stack object is part of
      the StackIteratorInterface. In addition to Stack the iterator only knows
-     the index fIndex in the Stack data.
+     the index index_ in the Stack data.
 
-     The template argument `ParticleInterface` acts as a policy to provide
-     readout function of Particle data from the stack. The ParticleInterface
+     The template argument `TParticleInterface` acts as a policy to provide
+     readout function of Particle data from the stack. The TParticleInterface
      class must know how to retrieve information from the Stack data
-     for a particle entry at any index fIndex.
+     for a particle entry at any index index_.
 
-     The ParticleInterface class must be written and provided by the
+     The TParticleInterface class must be written and provided by the
      user, it contains methods like <code> auto GetData() const {
      return GetStackData().GetData(GetIndex()); }</code>, where
      StackIteratorInterface::GetStackData() return a reference to the
-     object storing the particle data of type StackDataType. And
+     object storing the particle data of type TStackData. And
      StackIteratorInterface::GetIndex() provides the iterator index to
-     be readout. The StackDataType is another user-provided class to
+     be readout. The TStackData is another user-provided class to
      store data and must implement functions compatible with
-     ParticleInterface, in this example StackDataType::GetData(const unsigned int
+     TParticleInterface, in this example TStackData::GetData(const unsigned int
      vIndex).
 
      For two examples see stack_example.cc, or the
      corsika::processes::sibyll::SibStack class
   */
 
-  template <typename StackDataType, template <typename> typename ParticleInterface,
-            typename StackType = Stack<StackDataType, ParticleInterface>>
+  template <typename TStackData, template <typename> typename TParticleInterface,
+            typename StackType = Stack<TStackData, TParticleInterface>>
   class StackIteratorInterface
-      : public ParticleInterface<
-            StackIteratorInterface<StackDataType, ParticleInterface, StackType>> {
+      : public TParticleInterface<
+            StackIteratorInterface<TStackData, TParticleInterface, StackType>> {
 
   public:
     using ParticleInterfaceType =
-        ParticleInterface<corsika::stack::StackIteratorInterface<
-            StackDataType, ParticleInterface, StackType>>;
+        TParticleInterface<corsika::stack::StackIteratorInterface<
+            TStackData, TParticleInterface, StackType>>;
 
     // friends are needed for access to protected methods
-    friend class Stack<StackDataType,
-                       ParticleInterface>; // for access to GetIndex for Stack
-    friend class Stack<StackDataType&, ParticleInterface>; // for access to GetIndex
-                                                           // SecondaryView : public Stack
-    friend class ParticleBase<StackIteratorInterface>; // for access to GetStackDataType
-    friend class SecondaryView<StackDataType,
-                               ParticleInterface>; // access for SecondaryView
-
-    template <typename T>
-    friend class corsika::history::HSecondaryView;
-
+    friend class Stack<TStackData,
+                       TParticleInterface>; // for access to GetIndex for Stack
+    friend class Stack<TStackData&, TParticleInterface>; // for access to GetIndex
+                                                         // SecondaryView : public Stack
+    friend class ParticleBase<StackIteratorInterface>;   // for access to GetStackData
+    friend class SecondaryView<TStackData,
+                               TParticleInterface>; // access for SecondaryView
   private:
-    unsigned int fIndex = 0;
-    StackType* fData = 0; // info: Particles and StackIterators become invalid when parent
+    unsigned int index_ = 0;
+    StackType* data_ = 0; // info: Particles and StackIterators become invalid when parent
                           // Stack is copied or deleted!
 
     // it is not allowed to create a "dangling" stack iterator
@@ -94,12 +90,12 @@ namespace corsika::stack {
 
   public:
     StackIteratorInterface(StackIteratorInterface const& vR)
-        : fIndex(vR.fIndex)
-        , fData(vR.fData) {}
+        : index_(vR.index_)
+        , data_(vR.data_) {}
 
     StackIteratorInterface& operator=(StackIteratorInterface const& vR) {
-      fIndex = vR.fIndex;
-      fData = vR.fData;
+      index_ = vR.index_;
+      data_ = vR.data_;
       return *this;
     }
 
@@ -108,8 +104,8 @@ namespace corsika::stack {
           @param index index on stack
        */
     StackIteratorInterface(StackType& data, const unsigned int index)
-        : fIndex(index)
-        , fData(&data) {}
+        : index_(index)
+        , data_(&data) {}
 
     /** constructor that also sets new values on particle data object
         @param data reference to the stack [rw]
@@ -120,8 +116,8 @@ namespace corsika::stack {
      */
     template <typename... Args>
     StackIteratorInterface(StackType& data, const unsigned int index, const Args... args)
-        : fIndex(index)
-        , fData(&data) {
+        : index_(index)
+        , data_(&data) {
       (**this).SetParticleData(args...);
     }
 
@@ -138,29 +134,37 @@ namespace corsika::stack {
     template <typename... Args>
     StackIteratorInterface(StackType& data, const unsigned int index,
                            StackIteratorInterface& parent, const Args... args)
-        : fIndex(index)
-        , fData(&data) {
+        : index_(index)
+        , data_(&data) {
       (**this).SetParticleData(*parent, args...);
     }
 
+    bool isDeleted() const { return GetStack().isDeleted(*this); }
+
   public:
     /** @name Iterator interface
         @{
     */
     StackIteratorInterface& operator++() {
-      ++fIndex;
+      do {
+        ++index_;
+      } while (
+          GetStack().isDeleted(*this)); // this also check the allowed bounds of index_
       return *this;
     }
     StackIteratorInterface operator++(int) {
       StackIteratorInterface tmp(*this);
-      ++fIndex;
+      do {
+        ++index_;
+      } while (
+          GetStack().isDeleted(*this)); // this also check the allowed bounds of index_
       return tmp;
     }
     StackIteratorInterface operator+(int delta) {
-      return StackIteratorInterface(*fData, fIndex + delta);
+      return StackIteratorInterface(*data_, index_ + delta);
     }
-    bool operator==(const StackIteratorInterface& rhs) { return fIndex == rhs.fIndex; }
-    bool operator!=(const StackIteratorInterface& rhs) { return fIndex != rhs.fIndex; }
+    bool operator==(const StackIteratorInterface& rhs) { return index_ == rhs.index_; }
+    bool operator!=(const StackIteratorInterface& rhs) { return index_ != rhs.index_; }
 
     /**
      * Convert iterator to value type, where value type is the user-provided particle
@@ -184,18 +188,18 @@ namespace corsika::stack {
      * @{
      */
     /// Get current particle index
-    inline unsigned int GetIndex() const { return fIndex; }
+    inline unsigned int GetIndex() const { return index_; }
     /// Get current particle Stack object
-    inline StackType& GetStack() { return *fData; }
+    inline StackType& GetStack() { return *data_; }
     /// Get current particle const Stack object
-    inline const StackType& GetStack() const { return *fData; }
-    /// Get current user particle StackDataType object
-    inline StackDataType& GetStackData() { return fData->GetStackData(); }
-    /// Get current const user particle StackDataType object
-    inline const StackDataType& GetStackData() const { return fData->GetStackData(); }
+    inline const StackType& GetStack() const { return *data_; }
+    /// Get current user particle TStackData object
+    inline TStackData& GetStackData() { return data_->GetStackData(); }
+    /// Get current const user particle TStackData object
+    inline const TStackData& GetStackData() const { return data_->GetStackData(); }
     /// Get data index as mapped in Stack class
     inline unsigned int GetIndexFromIterator() const {
-      return fData->GetIndexFromIterator(fIndex);
+      return data_->GetIndexFromIterator(index_);
     }
     ///@}
   }; // end class StackIterator
@@ -206,24 +210,29 @@ namespace corsika::stack {
      This is the iterator class for const-access to stack data
    */
 
-  template <typename StackDataType, template <typename> typename ParticleInterface,
-            typename StackType = Stack<StackDataType, ParticleInterface>>
+  template <typename TStackData, template <typename> typename TParticleInterface,
+            typename StackType = Stack<TStackData, TParticleInterface>>
   class ConstStackIteratorInterface
-      : public ParticleInterface<
-            ConstStackIteratorInterface<StackDataType, ParticleInterface, StackType>> {
+      : public TParticleInterface<
+            ConstStackIteratorInterface<TStackData, TParticleInterface, StackType>> {
 
   public:
-    typedef ParticleInterface<
-        ConstStackIteratorInterface<StackDataType, ParticleInterface, StackType>>
+    typedef TParticleInterface<
+        ConstStackIteratorInterface<TStackData, TParticleInterface, StackType>>
         ParticleInterfaceType;
 
-    friend class Stack<StackDataType, ParticleInterface>;   // for access to GetIndex
-    friend class ParticleBase<ConstStackIteratorInterface>; // for access to
-                                                            // GetStackDataType
+    // friends are needed for access to protected methods
+    friend class Stack<TStackData,
+                       TParticleInterface>; // for access to GetIndex for Stack
+    friend class Stack<TStackData&, TParticleInterface>; // for access to GetIndex
+                                                         // SecondaryView : public Stack
+    friend class ParticleBase<ConstStackIteratorInterface>; // for access to GetStackData
+    friend class SecondaryView<TStackData,
+                               TParticleInterface>; // access for SecondaryView
 
   private:
-    unsigned int fIndex = 0;
-    const StackType* fData = 0; // info: Particles and StackIterators become invalid when
+    unsigned int index_ = 0;
+    const StackType* data_ = 0; // info: Particles and StackIterators become invalid when
                                 // parent Stack is copied or deleted!
 
     // we don't want to allow dangling iterators to exist
@@ -231,8 +240,8 @@ namespace corsika::stack {
 
   public:
     ConstStackIteratorInterface(const StackType& data, const unsigned int index)
-        : fIndex(index)
-        , fData(&data) {}
+        : index_(index)
+        , data_(&data) {}
 
     /**
        @class ConstStackIteratorInterface
@@ -247,27 +256,36 @@ namespace corsika::stack {
        See documentation of StackIteratorInterface for more details.
     */
 
+    bool isDeleted() const { return GetStack().isDeleted(*this); }
+
   public:
     /** @name Iterator interface
      */
     ///@{
     ConstStackIteratorInterface& operator++() {
-      ++fIndex;
+      do {
+        ++index_;
+      } while (
+          GetStack().isDeleted(*this)); // this also check the allowed bounds of index_
       return *this;
     }
     ConstStackIteratorInterface operator++(int) {
       ConstStackIteratorInterface tmp(*this);
-      ++fIndex;
+      do {
+        ++index_;
+      } while (
+          GetStack().isDeleted(*this)); // this also check the allowed bounds of index_
       return tmp;
     }
-    ConstStackIteratorInterface operator+(int delta) {
-      return ConstStackIteratorInterface(*fData, fIndex + delta);
+    ConstStackIteratorInterface operator+(const int delta) {
+      for (int i = 0; i < delta; ++i) operator++();
+      return ConstStackIteratorInterface(*data_, index_);
     }
     bool operator==(const ConstStackIteratorInterface& rhs) {
-      return fIndex == rhs.fIndex;
+      return index_ == rhs.index_;
     }
     bool operator!=(const ConstStackIteratorInterface& rhs) {
-      return fIndex != rhs.fIndex;
+      return index_ != rhs.index_;
     }
 
     const ParticleInterfaceType& operator*() const {
@@ -280,12 +298,12 @@ namespace corsika::stack {
         Only the const versions for read-only access
      */
     ///@{
-    inline unsigned int GetIndex() const { return fIndex; }
-    inline const StackType& GetStack() const { return *fData; }
-    inline const StackDataType& GetStackData() const { return fData->GetStackData(); }
+    inline unsigned int GetIndex() const { return index_; }
+    inline const StackType& GetStack() const { return *data_; }
+    inline const TStackData& GetStackData() const { return data_->GetStackData(); }
     /// Get data index as mapped in Stack class
     inline unsigned int GetIndexFromIterator() const {
-      return fData->GetIndexFromIterator(fIndex);
+      return data_->GetIndexFromIterator(index_);
     }
     ///@}
   }; // end class ConstStackIterator

Framework/StackInterface/testCombinedStack.cc

@@ -6,6 +6,8 @@
  * the license.
  */
 
+#define protected public // to also test the internal state of objects
+
 #include <corsika/stack/CombinedStack.h>
 #include <corsika/stack/SecondaryView.h>
 #include <corsika/stack/Stack.h>
@@ -107,7 +109,7 @@ TEST_CASE("Combined Stack", "[stack]") {
     s.AddParticle(std::tuple{0.});
     s.Copy(s.cbegin(), s.begin());
     s.Swap(s.begin(), s.begin());
-    REQUIRE(s.GetSize() == 1);
+    CHECK(s.getSize() == 1);
   }
 
   SECTION("construct") {
@@ -120,72 +122,99 @@ TEST_CASE("Combined Stack", "[stack]") {
 
     StackTest s;
     s.AddParticle(std::tuple{9.9});
-    REQUIRE(sum2(s) == 0.);
-    REQUIRE(sum(s) == 9.9);
+    CHECK(sum2(s) == 0.);
+    CHECK(sum(s) == 9.9);
   }
 
   SECTION("delete from stack") {
 
     StackTest s;
-    REQUIRE(s.GetSize() == 0);
+    CHECK(s.getSize() == 0);
     StackTest::StackIterator p =
         s.AddParticle(std::tuple{0.}); // valid way to access particle data
     p.SetData(8.9);
     p.SetData2(3.);
-    REQUIRE(sum2(s) == 3.);
-    REQUIRE(sum(s) == 8.9);
-    REQUIRE(s.GetSize() == 1);
+    CHECK(sum2(s) == 3.);
+    CHECK(sum(s) == 8.9);
+    CHECK(s.getSize() == 1);
+    CHECK(s.getEntries() == 1);
     s.Delete(p);
-    REQUIRE(s.GetSize() == 0);
+    CHECK(s.getSize() == 1);
+    CHECK(s.getEntries() == 0);
   }
 
   SECTION("delete particle") {
 
     StackTest s;
-    REQUIRE(s.GetSize() == 0);
+    CHECK(s.getSize() == 0);
     auto p = s.AddParticle(
         std::tuple{9.9}); // also valid way to access particle data, identical to above
-    REQUIRE(s.GetSize() == 1);
+    CHECK(s.getSize() == 1);
+    CHECK(s.getEntries() == 1);
     p.Delete();
-    REQUIRE(s.GetSize() == 0);
+    CHECK(s.getSize() == 1);
+    CHECK(s.getEntries() == 0);
   }
 
   SECTION("create secondaries") {
     StackTest s;
-    REQUIRE(s.GetSize() == 0);
+    CHECK(s.getSize() == 0);
     auto iter = s.AddParticle(std::tuple{9.9});
     iter.SetData2(2);
-    REQUIRE(s.GetSize() == 1);
+    CHECK(s.getSize() == 1);
+    CHECK(s.getEntries() == 1);
     iter.AddSecondary(std::tuple{4.4});
-    REQUIRE(s.GetSize() == 2);
+    CHECK(s.getSize() == 2);
+    CHECK(s.getEntries() == 2);
     // p.AddSecondary(3.3, 2.2, 1.);
-    // REQUIRE(s.GetSize() == 3);
+    // CHECK(s.getSize() == 3);
     double v = 0;
     for (const auto& i : s) {
       v += i.GetData();
-      REQUIRE(i.GetData2() == 2);
+      CHECK(i.GetData2() == 2);
     }
-    REQUIRE(v == 9.9 + 4.4);
+    CHECK(v == 9.9 + 4.4);
   }
 
   SECTION("get next particle") {
     StackTest s;
-    REQUIRE(s.GetSize() == 0);
+    CHECK(s.getSize() == 0);
+    CHECK(s.getEntries() == 0);
+    CHECK(s.IsEmpty());
+
     auto p1 = s.AddParticle(std::tuple{9.9});
     auto p2 = s.AddParticle(std::tuple{8.8});
     p1.SetData2(20.2);
     p2.SetData2(20.3);
-    auto particle = s.GetNextParticle(); // first particle
-    REQUIRE(particle.GetData() == 8.8);
-    REQUIRE(particle.GetData2() == 20.3);
+    CHECK(s.getSize() == 2);
+    CHECK(s.getEntries() == 2);
+    CHECK(!s.IsEmpty());
 
-    particle.Delete();
+    auto particle = s.GetNextParticle(); // first particle
+    CHECK(particle.GetData() == 8.8);
+    CHECK(particle.GetData2() == 20.3);
+
+    particle.Delete(); // only marks (last) particle as "deleted"
+    CHECK(s.getSize() == 2);
+    CHECK(s.getEntries() == 1);
+    CHECK(!s.IsEmpty());
+
+    /*
+      This following call to GetNextParticle will realize that the
+      current last particle on the stack was marked "deleted" and will
+      purge it: stack size is reduced by one.
+     */
     auto particle2 = s.GetNextParticle(); // first particle
-    REQUIRE(particle2.GetData() == 9.9);
-    REQUIRE(particle2.GetData2() == 20.2);
-    particle2.Delete();
-
-    REQUIRE(s.GetSize() == 0);
+    CHECK(s.getSize() == 1);
+    CHECK(s.getEntries() == 1);
+    CHECK(!s.IsEmpty());
+    CHECK(particle2.GetData() == 9.9);
+    CHECK(particle2.GetData2() == 20.2);
+
+    particle2.Delete(); // also mark this particle as "deleted"
+    CHECK(s.getSize() == 1);
+    CHECK(s.getEntries() == 0);
+    CHECK(s.IsEmpty());
   }
 }
 
@@ -199,7 +228,7 @@ class TestStackData3 {
 public:
   // these functions are needed for the Stack interface
   void Clear() { fData3.clear(); }
-  unsigned int GetSize() const { return fData3.size(); }
+  unsigned int getSize() const { return fData3.size(); }
   unsigned int GetCapacity() const { return fData3.size(); }
   void Copy(const int i1, const int i2) { fData3[i2] = fData3[i1]; }
   void Swap(const int i1, const int i2) {
@@ -259,11 +288,17 @@ TEST_CASE("Combined Stack - multi", "[stack]") {
   SECTION("create secondaries") {
 
     StackTest2 s;
-    REQUIRE(s.GetSize() == 0);
+    CHECK(s.getSize() == 0);
+    CHECK(s.IsEmpty()); // size = entries = 0
+
     // add new particle, only provide tuple data for StackTest
     auto p1 = s.AddParticle(std::tuple{9.9});
     // add new particle, provide tuple data for both StackTest and TestStackData3
     auto p2 = s.AddParticle(std::tuple{8.8}, std::tuple{0.1});
+
+    CHECK(s.getSize() == 2);
+    CHECK(!s.IsEmpty()); // size = entries = 2
+
     // examples to explicitly change data on stack
     p2.SetData2(0.1); // not clear why this is needed, need to check
                       // SetParticleData workflow for more complicated
@@ -271,32 +306,46 @@ TEST_CASE("Combined Stack - multi", "[stack]") {
     p1.SetData3(20.2);
     p2.SetData3(10.3);
 
-    REQUIRE(p1.GetData() == 9.9);
-    REQUIRE(p1.GetData2() == 0.);
+    CHECK(p1.GetData() == 9.9);
+    CHECK(p1.GetData2() == 0.);
     p1.SetData2(10.2);
-    REQUIRE(p1.GetData2() == 10.2);
-    REQUIRE(p1.GetData3() == 20.2);
+    CHECK(p1.GetData2() == 10.2);
+    CHECK(p1.GetData3() == 20.2);
 
-    REQUIRE(p2.GetData() == 8.8);
-    REQUIRE(p2.GetData2() == 0.1);
-    REQUIRE(p2.GetData3() == 10.3);
+    CHECK(p2.GetData() == 8.8);
+    CHECK(p2.GetData2() == 0.1);
+    CHECK(p2.GetData3() == 10.3);
 
     auto particle = s.GetNextParticle(); // first particle
-    REQUIRE(particle.GetData() == 8.8);
-    REQUIRE(particle.GetData2() == 0.1);
-    REQUIRE(particle.GetData3() == 10.3);
+    CHECK(particle.GetData() == 8.8);
+    CHECK(particle.GetData2() == 0.1);
+    CHECK(particle.GetData3() == 10.3);
 
-    REQUIRE(s.GetSize() == 2);
     auto sec = particle.AddSecondary(std::tuple{4.4});
-    REQUIRE(s.GetSize() == 3);
-    REQUIRE(sec.GetData() == 4.4);
-    REQUIRE(sec.GetData2() == 0.1);
-    REQUIRE(sec.GetData3() == 10.3);
-
-    sec.Delete();
-    s.DeleteLast();
-    s.GetNextParticle().Delete();
-    REQUIRE(s.GetSize() == 0);
+    CHECK(s.getSize() == 3);
+    CHECK(s.getEntries() == 3);
+    CHECK(sec.GetData() == 4.4);
+    CHECK(sec.GetData2() == 0.1);
+    CHECK(sec.GetData3() == 10.3);
+
+    sec.Delete(); // mark for deletion: size=3, entries=2
+    CHECK(s.getSize() == 3);
+    CHECK(s.getEntries() == 2);
+    CHECK(!s.IsEmpty());
+
+    s.last().Delete(); // mark for deletion: size=3, entries=1
+    CHECK(s.getSize() == 3);
+    CHECK(s.getEntries() == 1);
+    CHECK(!s.IsEmpty());
+
+    /*
+       GetNextParticle will find two entries marked as "deleted" and
+       will purge this from the end of the stack: size = 1
+    */
+    s.GetNextParticle().Delete(); // mark for deletion: size=3, entries=0
+    CHECK(s.getSize() == 1);
+    CHECK(s.getEntries() == 0);
+    CHECK(s.IsEmpty());
   }
 }
 
@@ -341,6 +390,6 @@ TEST_CASE("Combined Stack - secondary view") {
     auto projectile = view.GetProjectile();
     projectile.AddSecondary(std::tuple{8.8});
 
-    REQUIRE(stack.GetSize() == 2);
+    CHECK(stack.getSize() == 2);
   }
 }

Framework/StackInterface/testSecondaryView.cc

@@ -6,6 +6,8 @@
  * the license.
  */
 
+#define protected public // to also test the internal state of objects
+
 #include <corsika/stack/SecondaryView.h>
 #include <corsika/stack/Stack.h>
 
@@ -45,73 +47,90 @@ TEST_CASE("SecondaryStack", "[stack]") {
   // helper function for sum over stack data
   auto sum = [](const StackTest& stack) {
     double v = 0;
-    for (const auto& p : stack) v += p.GetData();
+    for (const auto& p : stack) { v += p.GetData(); }
     return v;
   };
 
+  auto sumView = [](const StackTestView& stack) {
+    double value = 0;
+    for (const auto& p : stack) { value += p.GetData(); }
+    return value;
+  };
+
   SECTION("secondary view") {
     StackTest s;
-    REQUIRE(s.GetSize() == 0);
+    CHECK(s.getSize() == 0);
+    CHECK(s.IsEmpty());
+
     s.AddParticle(std::tuple{9.9});
     s.AddParticle(std::tuple{8.8});
-    const double sumS = 9.9 + 8.8;
+    const double sumS = 9.9 + 8.8; // helper, see below
+    CHECK(s.getSize() == 2);
+    CHECK(s.getEntries() == 2);
+    CHECK(!s.IsEmpty());
 
     auto particle = s.GetNextParticle();
 
     StackTestView view(particle);
-    REQUIRE(view.GetSize() == 0);
+    CHECK(view.getSize() == 0);
+    CHECK(view.getEntries() == 0);
+    CHECK(view.IsEmpty());
 
     {
       auto proj = view.GetProjectile();
-      REQUIRE(proj.GetData() == particle.GetData());
+      CHECK(proj.GetData() == particle.GetData());
       proj.AddSecondary(std::tuple{4.4});
     }
+    CHECK(view.getSize() == 1);
+    CHECK(view.getEntries() == 1);
+    CHECK(!view.IsEmpty());
+    CHECK(s.getSize() == 3);
+    CHECK(s.getEntries() == 3);
+    CHECK(!s.IsEmpty());
 
     view.AddSecondary(std::tuple{4.5});
     view.AddSecondary(std::tuple{4.6});
+    CHECK(view.getSize() == 3);
+    CHECK(view.getEntries() == 3);
+    CHECK(!view.IsEmpty());
+    CHECK(s.getSize() == 5);
+    CHECK(s.getEntries() == 5);
+    CHECK(!s.IsEmpty());
 
-    REQUIRE(view.GetSize() == 3);
-    REQUIRE(s.GetSize() == 5);
-    REQUIRE(!view.IsEmpty());
-
-    auto sumView = [](const StackTestView& stack) {
-      double value = 0;
-      for (const auto& p : stack) { value += p.GetData(); }
-      return value;
-    };
-
-    REQUIRE(sum(s) == sumS + 4.4 + 4.5 + 4.6);
-    REQUIRE(sumView(view) == 4.4 + 4.5 + 4.6);
+    CHECK(sum(s) == sumS + 4.4 + 4.5 + 4.6);
+    CHECK(sumView(view) == 4.4 + 4.5 + 4.6);
 
-    view.DeleteLast();
-    REQUIRE(view.GetSize() == 2);
-    REQUIRE(s.GetSize() == 4);
+    view.last().Delete();
+    CHECK(view.getSize() == 3);
+    CHECK(view.getEntries() == 2);
+    CHECK(s.getSize() == 5);
+    CHECK(s.getEntries() == 4);
 
-    REQUIRE(sum(s) == sumS + 4.4 + 4.5);
-    REQUIRE(sumView(view) == 4.4 + 4.5);
+    CHECK(sum(s) == sumS + 4.4 + 4.5);
+    CHECK(sumView(view) == 4.4 + 4.5);
 
     auto pDel = view.GetNextParticle();
     view.Delete(pDel);
-    REQUIRE(view.GetSize() == 1);
-    REQUIRE(s.GetSize() == 3);
+    CHECK(view.getSize() == 2);
+    CHECK(s.getSize() == 4);
 
-    REQUIRE(sum(s) == sumS + 4.4 + 4.5 - pDel.GetData());
-    REQUIRE(sumView(view) == 4.4 + 4.5 - pDel.GetData());
+    CHECK(sum(s) == sumS + 4.4 + 4.5 - pDel.GetData());
+    CHECK(sumView(view) == 4.4 + 4.5 - pDel.GetData());
 
     view.Delete(view.GetNextParticle());
-    REQUIRE(sum(s) == sumS);
-    REQUIRE(sumView(view) == 0);
-    REQUIRE(view.IsEmpty());
+    CHECK(sum(s) == sumS);
+    CHECK(sumView(view) == 0);
+    CHECK(view.IsEmpty());
 
     {
       auto proj = view.GetProjectile();
-      REQUIRE(proj.GetData() == particle.GetData());
+      CHECK(proj.GetData() == particle.GetData());
     }
   }
 
   SECTION("secondary view, construct from ParticleType") {
     StackTest s;
-    REQUIRE(s.GetSize() == 0);
+    CHECK(s.getSize() == 0);
     s.AddParticle(std::tuple{9.9});
     s.AddParticle(std::tuple{8.8});
 
@@ -119,11 +138,11 @@ TEST_CASE("SecondaryStack", "[stack]") {
     typename StackTest::ParticleType& particle = iterator; // as in corsika::Cascade
 
     StackTestView view(particle);
-    REQUIRE(view.GetSize() == 0);
+    CHECK(view.getSize() == 0);
 
     view.AddSecondary(std::tuple{4.4});
 
-    REQUIRE(view.GetSize() == 1);
+    CHECK(view.getSize() == 1);
   }
 
   SECTION("deletion") {
@@ -141,21 +160,20 @@ TEST_CASE("SecondaryStack", "[stack]") {
       proj.AddSecondary(std::tuple{1.});
       proj.AddSecondary(std::tuple{2.});
 
-      CHECK(stack.GetSize() == 6); // -99, 0, -2, -1, 1, 2
-      CHECK(view.GetSize() == 4);  // -2, -1, 1, 2
+      CHECK(stack.getSize() == 6); // -99, 0, -2, -1, 1, 2
+      CHECK(view.getSize() == 4);  // -2, -1, 1, 2
 
       // now delete all negative entries, i.e. -1 and -2
       auto p = view.begin();
       while (p != view.end()) {
         auto data = p.GetData();
-        if (data < 0) {
-          p.Delete();
-        } else {
-          ++p;
-        }
+        if (data < 0) { p.Delete(); }
+        ++p;
       }
-      CHECK(stack.GetSize() == 4); // -99, 0, 2, 1 (order changes during deletion)
-      CHECK(view.GetSize() == 2);  // 2, 1
+      CHECK(stack.getSize() == 6);
+      CHECK(stack.getEntries() == 4); // -99, 0, 2, 1 (order changes during deletion)
+      CHECK(view.getSize() == 4);
+      CHECK(view.getEntries() == 2); // 2, 1
     }
 
     // repeat
@@ -175,17 +193,15 @@ TEST_CASE("SecondaryStack", "[stack]") {
       auto p = view.begin();
       while (p != view.end()) {
         auto data = p.GetData();
-        if (data < 0) {
-          p.Delete();
-        } else {
-          ++p;
-        }
+        if (data < 0) { p.Delete(); }
+        ++p;
       }
 
       // stack should contain -99, 0, 2, 1, [2, 1]
       // view should contain 1, 2
 
-      CHECK(stack.GetSize() == 6);
+      CHECK(stack.getEntries() == 6);
+      CHECK(stack.getSize() == 10);
     }
   }
 }

Framework/StackInterface/testStackInterface.cc

@@ -6,6 +6,8 @@
  * the license.
  */
 
+#define protected public // to also test the internal state of objects
+
 #include <corsika/stack/Stack.h>
 
 #include <testTestStack.h> // simple test-stack for testing. This is
@@ -42,7 +44,7 @@ TEST_CASE("Stack", "[Stack]") {
     s.AddParticle(std::tuple{0.});
     s.Copy(s.cbegin(), s.begin());
     s.Swap(s.begin(), s.begin());
-    REQUIRE(s.GetSize() == 1);
+    CHECK(s.getSize() == 1);
   }
 
   SECTION("construct") {
@@ -56,62 +58,108 @@ TEST_CASE("Stack", "[Stack]") {
     StackTest s;
     s.AddParticle(std::tuple{9.9});
     const double v = sum(s);
-    REQUIRE(v == 9.9);
+    CHECK(v == 9.9);
   }
 
   SECTION("delete from stack") {
 
     StackTest s;
-    REQUIRE(s.GetSize() == 0);
+    CHECK(s.getSize() == 0);
     StackTest::StackIterator p =
         s.AddParticle(std::tuple{0.}); // valid way to access particle data
     p.SetData(9.9);
-    REQUIRE(s.GetSize() == 1);
+    CHECK(s.getSize() == 1);
+    CHECK(s.getEntries() == 1);
     s.Delete(p);
-    REQUIRE(s.GetSize() == 0);
+    CHECK(s.getSize() == 1);
+    CHECK(s.getEntries() == 0);
   }
 
   SECTION("delete particle") {
 
     StackTest s;
-    REQUIRE(s.GetSize() == 0);
+    CHECK(s.getSize() == 0);
+    s.AddParticle(std::tuple{8.9});
+    s.AddParticle(std::tuple{7.9});
     auto p = s.AddParticle(
         std::tuple{9.9}); // also valid way to access particle data, identical to above
-    REQUIRE(s.GetSize() == 1);
-    p.Delete();
-    REQUIRE(s.GetSize() == 0);
+
+    CHECK(s.getSize() == 3);
+    CHECK(s.getEntries() == 3);
+    CHECK(!s.IsEmpty());
+
+    p.Delete(); // mark for deletion: size=3, entries=2
+    CHECK(s.getSize() == 3);
+    CHECK(s.getEntries() == 2);
+    CHECK(!s.IsEmpty());
+
+    s.last().Delete(); // mark for deletion: size=3, entries=1
+    CHECK(s.getSize() == 3);
+    CHECK(s.getEntries() == 1);
+    CHECK(!s.IsEmpty());
+
+    /*
+       GetNextParticle will find two entries marked as "deleted" and
+       will purge this from the end of the stack: size = 1
+    */
+    s.GetNextParticle().Delete(); // mark for deletion: size=3, entries=0
+    CHECK(s.getSize() == 1);
+    CHECK(s.getEntries() == 0);
+    CHECK(s.IsEmpty());
   }
 
   SECTION("create secondaries") {
 
     StackTest s;
-    REQUIRE(s.GetSize() == 0);
+    CHECK(s.getSize() == 0);
     auto iter = s.AddParticle(std::tuple{9.9});
     StackTest::ParticleInterfaceType& p =
         *iter; // also this is valid to access particle data
-    REQUIRE(s.GetSize() == 1);
+    CHECK(s.getSize() == 1);
     p.AddSecondary(std::tuple{4.4});
-    REQUIRE(s.GetSize() == 2);
+    CHECK(s.getSize() == 2);
     /*p.AddSecondary(3.3, 2.2);
-    REQUIRE(s.GetSize() == 3);
+    CHECK(s.getSize() == 3);
     double v = 0;
     for (auto& p : s) { v += p.GetData(); }
-    REQUIRE(v == 9.9 + 4.4 + 3.3 + 2.2);*/
+    CHECK(v == 9.9 + 4.4 + 3.3 + 2.2);*/
   }
 
   SECTION("get next particle") {
     StackTest s;
-    REQUIRE(s.GetSize() == 0);
+    CHECK(s.getSize() == 0);
+    CHECK(s.getEntries() == 0);
+    CHECK(s.IsEmpty());
+
     s.AddParticle(std::tuple{9.9});
     s.AddParticle(std::tuple{8.8});
-    auto particle = s.GetNextParticle(); // first particle
-    REQUIRE(particle.GetData() == 8.8);
+    CHECK(s.getSize() == 2);
+    CHECK(s.getEntries() == 2);
+    CHECK(!s.IsEmpty());
 
-    particle.Delete();
+    auto particle = s.GetNextParticle(); // first particle
+    CHECK(particle.GetData() == 8.8);
+
+    particle.Delete(); // only marks (last) particle as deleted
+    CHECK(s.getSize() == 2);
+    CHECK(s.getEntries() == 1);
+    CHECK(!s.IsEmpty());
+
+    /*
+      This following call to GetNextParticle will realize that the
+      current last particle on the stack was marked "deleted" and will
+      purge it: stack size is reduced by one.
+     */
     auto particle2 = s.GetNextParticle(); // first particle
-    REQUIRE(particle2.GetData() == 9.9);
-    particle2.Delete();
+    CHECK(particle2.GetData() == 9.9);
+    CHECK(s.getSize() == 1);
+    CHECK(s.getEntries() == 1);
+    CHECK(!s.IsEmpty());
+
+    particle2.Delete(); // also mark this particle as deleted
 
-    REQUIRE(s.GetSize() == 0);
+    CHECK(s.getSize() == 1);
+    CHECK(s.getEntries() == 0);
+    CHECK(s.IsEmpty());
   }
 }

Processes/CONEXSourceCut/CONEXSourceCut.cc

@@ -29,17 +29,17 @@ corsika::process::EProcessReturn CONEXSourceCut::DoSecondaries(
 
     auto const it = std::find_if(egs_em_codes_.cbegin(), egs_em_codes_.cend(),
                                  [=](auto const& p) { return pid == p.first; });
-    if (it == egs_em_codes_.cend()) {
-      ++p;
-      continue; // no EM particle
-    }
-
-    auto const egs_pid = it->second;
+    if (it != egs_em_codes_.cend()) {
+      // EM particle
+    
+      auto const egs_pid = it->second;
 
-    addParticle(egs_pid, p.GetEnergy(), p.GetMass(), p.GetPosition(),
-                p.GetMomentum().normalized(), p.GetTime());
+      addParticle(egs_pid, p.GetEnergy(), p.GetMass(), p.GetPosition(),
+		  p.GetMomentum().normalized(), p.GetTime());
 
-    p.Delete();
+      p.Delete();
+    }
+    ++p;
   }
 
   return corsika::process::EProcessReturn::eOk;