diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc index cb06341c84be104542fd64451d1b72e6f5d8539c..2ef7cdcf0d03a973e5e57099d0fa11f48eabd58f 100644 --- a/Documentation/Examples/cascade_example.cc +++ b/Documentation/Examples/cascade_example.cc @@ -207,7 +207,6 @@ public: HEPEnergyType GetEmEnergy() const { return fEmEnergy; } }; - // // The example main program for a particle cascade // diff --git a/Documentation/Examples/vertical_EAS.cc b/Documentation/Examples/vertical_EAS.cc index 526498b2421055dc21596c76db08adf23842c38a..d029a35cf6c9eb0d1af938c5522d2ee68a4bf75a 100644 --- a/Documentation/Examples/vertical_EAS.cc +++ b/Documentation/Examples/vertical_EAS.cc @@ -210,7 +210,6 @@ public: HEPEnergyType GetEmEnergy() const { return fEmEnergy; } }; - // // The example main program for a particle cascade // diff --git a/Framework/StackInterface/SecondaryView.h b/Framework/StackInterface/SecondaryView.h index 1bbbe02ccaad07fcc533d5782c753427c0dddbf1..fae9d463f9c6808b69d346092fe6277a7c5fec7a 100644 --- a/Framework/StackInterface/SecondaryView.h +++ b/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(); } diff --git a/Framework/StackInterface/testCombinedStack.cc b/Framework/StackInterface/testCombinedStack.cc index 8fe244b667bf736cc418368ee2e06116ef845608..aceccd1b472a99c57168eef4d78851b3641c092c 100644 --- a/Framework/StackInterface/testCombinedStack.cc +++ b/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 diff --git a/Framework/StackInterface/testSecondaryView.cc b/Framework/StackInterface/testSecondaryView.cc index 2d2064b64e2e09c96ec5f77ebbf25d50f5d60e9e..95c196d5ecef139339b21bb41623f8a52e846d76 100644 --- a/Framework/StackInterface/testSecondaryView.cc +++ b/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); + } + } } diff --git a/Framework/Units/PhysicalConstants.h b/Framework/Units/PhysicalConstants.h index 78f9ecd35a5872ccc05ddebab3dd5f6873a65d1f..23494f036ed2f692de71984c8fc1b81a6d919112 100644 --- a/Framework/Units/PhysicalConstants.h +++ b/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); diff --git a/Framework/Utilities/sgn.h b/Framework/Utilities/sgn.h index 71d5bec5bbb59fe8c71ef35d0c3f9da1a40b27b7..43d9a9fad728785ec851424b2113f5e7dac60af1 100644 --- a/Framework/Utilities/sgn.h +++ b/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 diff --git a/Processes/EnergyLoss/EnergyLoss.h b/Processes/EnergyLoss/EnergyLoss.h index 711e135f716099fc3204a9be14239239a97a3016..d0726d9d44220f27c356f5128a43a0092cecc558 100644 --- a/Processes/EnergyLoss/EnergyLoss.h +++ b/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 diff --git a/Processes/Pythia/Decay.cc b/Processes/Pythia/Decay.cc index 670d24196423606bbce79329ce95272630412633..20d2aaa5044e484b2727e5cdc03f011fe36be83b 100644 --- a/Processes/Pythia/Decay.cc +++ b/Processes/Pythia/Decay.cc @@ -147,7 +147,6 @@ namespace corsika::process::pythia { // set particle stable Decay::SetStable(vP.GetPID()); - } } // namespace corsika::process::pythia