Merge branch '183-something-wrong-with-secondaryview-delete' into 'master'

Resolve "something wrong with SecondaryView::Delete?"

Closes #183

See merge request AirShowerPhysics/corsika!117

Documentation/Examples/cascade_example.cc

@@ -207,7 +207,6 @@ public:
   HEPEnergyType GetEmEnergy() const { return fEmEnergy; }
 };
 
-
 //
 // The example main program for a particle cascade
 //

Documentation/Examples/vertical_EAS.cc

@@ -210,7 +210,6 @@ public:
   HEPEnergyType GetEmEnergy() const { return fEmEnergy; }
 };
 
-
 //
 // The example main program for a particle cascade
 //

Framework/StackInterface/SecondaryView.h

@@ -43,14 +43,18 @@ namespace corsika::stack {
      All operations possible on a Stack object are also possible on a
      SecondaryView object. This means you can add, delete, copy, swap,
      iterate, etc.
-   */
 
-  /* INTERNAL NOTE FOR DEVELOPERS The secondary particle indices are
-     stored in a std::vector fIndices, the index of the primary projectle
-     particle is explicitly stored in fProjectileIndex. StackIterator
-     indices are refering to those numbers, where
-     StackIterator::GetIndex()==0 refers to the fProjectileIndex and
-     StackIterator::GetIndex()>0 to fIndices[i+1], see function GetIndexFromIterator.
+     *Further information about implementation (for developers):* All
+     data is stored in the original stack privided at construction
+     time. The secondary particle (view) indices are stored in an
+     extra std::vector of SecondaryView class 'fIndices' referring to
+     the original stack slot indices. The index of the primary
+     projectle particle is also explicitly stored in
+     'fProjectileIndex'. StackIterator indices
+     'i = StackIterator::GetIndex()' are referring to those numbers,
+     where 'i==0' refers to the 'fProjectileIndex', and
+     'StackIterator::GetIndex()>0' to 'fIndices[i-1]', see function
+     GetIndexFromIterator.
    */
 
   template <typename StackDataType, template <typename> typename ParticleInterface>
@@ -125,8 +129,11 @@ namespace corsika::stack {
 
     template <typename... Args>
     auto AddSecondary(StackIterator& proj, const Args... v) {
+      // make space on stack
       InnerStackType::GetStackData().IncrementSize();
+      // get current number of secondaries on stack
       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);
       // NOTE: "+1" is since "0" is special marker here for PROJECTILE, see
@@ -161,14 +168,26 @@ namespace corsika::stack {
     /// @}
 
     /**
-     * need overwrite Stack::Delete, since we want to call SecondaryView::DeleteLast
+     * need overwrite Stack::Delete, since we want to call
+     * SecondaryView::DeleteLast
+     *
+     * The particle is deleted on the underlying (internal) stack. The
+     * local references in SecondaryView in fIndices must be fixed,
+     * too.  The approach is to a) check if the particle 'p' is at the
+     * very end of the internal stack, b) if not: move it there by
+     * copying the last particle to the current particle location, c)
+     * remove the last particle.
+     *
      */
     void Delete(StackIterator p) {
       if (IsEmpty()) { /* error */
         throw std::runtime_error("Stack, cannot delete entry since size is zero");
       }
-      if (p.GetIndex() < GetSize() - 1)
-        InnerStackType::GetStackData().Copy(GetSize() - 1, p.GetIndex());
+      const int innerSize = InnerStackType::GetSize();
+      const int innerIndex = GetIndexFromIterator(p.GetIndex());
+      if (innerIndex < innerSize - 1)
+        InnerStackType::GetStackData().Copy(innerSize - 1,
+                                            GetIndexFromIterator(p.GetIndex()));
       DeleteLast();
     }
 

Framework/StackInterface/testCombinedStack.cc

@@ -333,7 +333,8 @@ using StackTest2 = CombinedStack<typename StackTest::StackImpl, TestStackData3,
                                  CombinedTestInterfaceType2>;
 
 #if defined(__clang__)
-using StackTestView = SecondaryView<typename StackTest2::StackImpl, CombinedTestInterfaceType2>;
+using StackTestView =
+    SecondaryView<typename StackTest2::StackImpl, CombinedTestInterfaceType2>;
 #elif defined(__GNUC__) || defined(__GNUG__)
 using StackTestView = corsika::stack::MakeView<StackTest2>::type;
 #endif

Framework/StackInterface/testSecondaryView.cc

@@ -134,4 +134,67 @@ TEST_CASE("SecondaryStack", "[stack]") {
 
     REQUIRE(view.GetSize() == 1);
   }
+
+  SECTION("deletion") {
+    StackTest stack;
+    stack.AddParticle(std::tuple{-99.});
+    stack.AddParticle(std::tuple{0.});
+
+    {
+      auto particle = stack.GetNextParticle();
+      StackTestView view(particle);
+
+      auto proj = view.GetProjectile();
+      proj.AddSecondary(std::tuple{-2.});
+      proj.AddSecondary(std::tuple{-1.});
+      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
+
+      // 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;
+        }
+      }
+      CHECK(stack.GetSize() == 4); // -99, 0, 2, 1 (order changes during deletion)
+      CHECK(view.GetSize() == 2);  // 2, 1
+    }
+
+    // repeat
+
+    {
+      auto particle = stack.GetNextParticle();
+      StackTestView view(particle);
+
+      // put -2,...,+2 on stack
+      auto proj = view.GetProjectile();
+      proj.AddSecondary(std::tuple{-2.});
+      proj.AddSecondary(std::tuple{-1.});
+      proj.AddSecondary(std::tuple{1.});
+      proj.AddSecondary(std::tuple{2.});
+      // stack should contain -99, 0, 2, 1, [-2, -1, 1, 2]
+
+      auto p = view.begin();
+      while (p != view.end()) {
+        auto data = p.GetData();
+        if (data < 0) {
+          p.Delete();
+        } else {
+          ++p;
+        }
+      }
+
+      // stack should contain -99, 0, 2, 1, [2, 1]
+      // view should contain 1, 2
+
+      CHECK(stack.GetSize() == 6);
+    }
+  }
 }

Framework/Units/PhysicalConstants.h

@@ -50,7 +50,7 @@ namespace corsika::units::constants {
   constexpr quantity<mass_d> u{Rep(1.6605402e-27L) * kilogram};
 
   auto constexpr nucleonMass = 0.5 * (0.93827 + 0.93957) * 1e9 * electronvolt;
-  
+
   // molar gas constant
   auto constexpr R = Rep(8.314'459'8) * joule / (mole * kelvin);
 

Framework/Utilities/sgn.h

@@ -3,12 +3,12 @@
 
 namespace corsika::utl {
 
-//! sign function without branches
-template <typename T>
-static int sgn(T val) {
-  return (T(0) < val) - (val < T(0));
-}
+  //! sign function without branches
+  template <typename T>
+  static int sgn(T val) {
+    return (T(0) < val) - (val < T(0));
+  }
 
-}
+} // namespace corsika::utl
 
 #endif

Processes/EnergyLoss/EnergyLoss.h

@@ -44,7 +44,8 @@ namespace corsika::process::energy_loss {
     static units::si::HEPEnergyType BetheBloch(setup::Stack::ParticleType const& p,
                                                const units::si::GrammageType dX);
 
-    int GetXbin(setup::Stack::ParticleType const&, setup::Trajectory const&, units::si::HEPEnergyType);
+    int GetXbin(setup::Stack::ParticleType const&, setup::Trajectory const&,
+                units::si::HEPEnergyType);
 
     units::si::HEPEnergyType fEnergyLossTot;
     units::si::GrammageType fdX;    // profile binning

Processes/Pythia/Decay.cc

@@ -147,7 +147,6 @@ namespace corsika::process::pythia {
 
     // set particle stable
     Decay::SetStable(vP.GetPID());
-
   }
 
 } // namespace corsika::process::pythia