diff --git a/corsika/detail/stack/VectorStack.inl b/corsika/detail/stack/VectorStack.inl index 40fb4e61c759e256af9d6fcac6283c6eee44fdb2..9d2a9434c464df9ff1d7c8b1e1b0cdfc4c5143ff 100644 --- a/corsika/detail/stack/VectorStack.inl +++ b/corsika/detail/stack/VectorStack.inl @@ -91,7 +91,13 @@ namespace corsika { time_.clear(); } - inline void VectorStackImpl::copy(size_t i1, size_t i2) { + inline void VectorStackImpl::copy(size_t const i1, size_t const i2) { + // index range check + if (i1 >= getSize() || i2 >= getSize()) { + std::ostringstream err; + err << "VectorStackImpl: trying to access data beyond size of stack !"; + throw std::runtime_error(err.str()); + } dataPID_[i2] = dataPID_[i1]; dataEkin_[i2] = dataEkin_[i1]; direction_[i2] = direction_[i1]; @@ -99,7 +105,13 @@ namespace corsika { time_[i2] = time_[i1]; } - inline void VectorStackImpl::swap(size_t i1, size_t i2) { + inline void VectorStackImpl::swap(size_t const i1, size_t const i2) { + // index range check + if (i1 >= getSize() || i2 >= getSize()) { + std::ostringstream err; + err << "VectorStackImpl: trying to access data beyond size of stack !"; + throw std::runtime_error(err.str()); + } std::swap(dataPID_[i2], dataPID_[i1]); std::swap(dataEkin_[i2], dataEkin_[i1]); std::swap(direction_[i2], direction_[i1]); diff --git a/corsika/stack/VectorStack.hpp b/corsika/stack/VectorStack.hpp index 154a9850bcb11e8157c8c2e9bc29c862f907184c..2a96ef3df912765c7fc81f67aceceb24cb7c6f07 100644 --- a/corsika/stack/VectorStack.hpp +++ b/corsika/stack/VectorStack.hpp @@ -177,7 +177,7 @@ namespace corsika { }; /** - * Memory implementation of the most simple (stupid) particle stack object. + * Memory implementation of the most simple particle stack object. * * @note if we ever want to have off-shell particles, we need to * add momentum as HEPMomentumType, and a lot of care. @@ -222,14 +222,14 @@ namespace corsika { TimeType getTime(size_t i) const { return time_[i]; } /** - * Function to copy particle at location i2 in stack to i1 + * Function to copy particle at location i2 in stack to i1. */ - void copy(size_t i1, size_t i2); + void copy(size_t const i1, size_t const i2); /** - * Function to copy particle at location i2 in stack to i1 + * Function to copy particle at location i2 in stack to i1. */ - void swap(size_t i1, size_t i2); + void swap(size_t const i1, size_t const i2); void incrementSize(); void decrementSize(); diff --git a/tests/stack/testVectorStack.cpp b/tests/stack/testVectorStack.cpp index 7fa0903b456d6bcd26ebf72a728beeb88746eba3..b44e3d55896e1c5d259642210f9b4dfa8de19a68 100644 --- a/tests/stack/testVectorStack.cpp +++ b/tests/stack/testVectorStack.cpp @@ -38,6 +38,17 @@ TEST_CASE("VectorStack", "stack") { CHECK(pout.getEnergy() == 1.5_GeV + Electron::mass); CHECK(pout.getKineticEnergy() == 1.5_GeV); CHECK(pout.getTime() == 100_s); + CHECK(pout.getChargeNumber() == -1); + + // particle with no momentum, has no direction + auto const p0 = s.addParticle( + std::make_tuple(Code::Proton, MomentumVector(dummyCS, {0_GeV, 0_GeV, 0_GeV}), + Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s)); + CHECK(p0.getDirection().getNorm() == 0); + + s.clear(); + CHECK(s.getEntries() == 0); + CHECK(s.getSize() == 0); } SECTION("write+delete") { @@ -56,4 +67,98 @@ TEST_CASE("VectorStack", "stack") { CHECK(s.getEntries() == 0); CHECK(s.getSize() == 1); } -} + + SECTION("stack operations") { + + VectorStack s; + // add 99 particles, each 10th particle is a positron + // i=9, 19, 29, etc. are positron + for (int i = 0; i < 99; ++i) { + Code const pid = ((i + 1) % 10 == 0) ? Code::Positron : Code::Electron; + s.addParticle(std::make_tuple( + pid, i * 1.5_GeV - Electron::mass, DirectionVector(dummyCS, {1, 0, 0}), + Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s)); + } + + // copy + { + s.copy(s.begin() + 9, s.begin() + 10); // positron to electron + const auto& p9 = s.cbegin() + 9; + const auto& p10 = s.cbegin() + 10; + + CHECK(p9.getPID() == Code::Positron); + CHECK(p9.getEnergy() == 9 * 1.5_GeV); + CHECK(p9.getTime() == 100_s); + + CHECK(p10.getPID() == Code::Positron); + CHECK(p10.getEnergy() == 9 * 1.5_GeV); + CHECK(p10.getTime() == 100_s); + } + + // copy + { + s.copy(s.begin() + 93, s.begin() + 9); // electron to positron + const auto& p93 = s.cbegin() + 93; + const auto& p9 = s.cbegin() + 9; + + CHECK(p9.getPID() == Code::Electron); + CHECK(p9.getEnergy() == 93 * 1.5_GeV); + CHECK(p9.getTime() == 100_s); + + CHECK(p93.getPID() == Code::Electron); + CHECK(p93.getEnergy() == 93 * 1.5_GeV); + CHECK(p93.getTime() == 100_s); + } + + // copy + { + s.copy(s.begin() + 89, s.begin() + 79); // positron to positron + const auto& p89 = s.cbegin() + 89; + const auto& p79 = s.cbegin() + 79; + + CHECK(p89.getPID() == Code::Positron); + CHECK(p89.getEnergy() == 89 * 1.5_GeV); + CHECK(p89.getTime() == 100_s); + + CHECK(p79.getPID() == Code::Positron); + CHECK(p79.getEnergy() == 89 * 1.5_GeV); + CHECK(p79.getTime() == 100_s); + } + + // invalid copy + { CHECK_THROWS(s.copy(s.begin(), s.end() + 1000)); } + + // swap + { + s.swap(s.begin() + 11, s.begin() + 10); + const auto& p11 = s.cbegin() + 11; // now: positron + const auto& p10 = s.cbegin() + 10; // now: electron + + CHECK(p11.getPID() == Code::Positron); + CHECK(p11.getEnergy() == 9 * 1.5_GeV); + CHECK(p11.getTime() == 100_s); + + CHECK(p10.getPID() == Code::Electron); + CHECK(p10.getEnergy() == 11 * 1.5_GeV); + CHECK(p10.getTime() == 100_s); + } + + // swap two positrons + { + s.swap(s.begin() + 29, s.begin() + 59); + const auto& p29 = s.cbegin() + 29; + const auto& p59 = s.cbegin() + 59; + + CHECK(p29.getPID() == Code::Positron); + CHECK(p29.getEnergy() == 59 * 1.5_GeV); + CHECK(p29.getTime() == 100_s); + + CHECK(p59.getPID() == Code::Positron); + CHECK(p59.getEnergy() == 29 * 1.5_GeV); + CHECK(p59.getTime() == 100_s); + } + + for (int i = 0; i < 99; ++i) s.last().erase(); + CHECK(s.getEntries() == 0); + } +} \ No newline at end of file diff --git a/tests/stack/testWeightStackExtension.cpp b/tests/stack/testWeightStackExtension.cpp index d984e1e5cc36588fac2fd36623df996cd924ab44..aa3471a821715a7ec9367a776f94c031c8f1d3d9 100644 --- a/tests/stack/testWeightStackExtension.cpp +++ b/tests/stack/testWeightStackExtension.cpp @@ -34,7 +34,6 @@ using TestStack = CombinedStack<typename dummy_stack::DummyStack::stack_data_typ TEST_CASE("WeightStackExtension", "[stack]") { logging::set_level(logging::level::info); - corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v"); dummy_stack::NoData noData; @@ -43,9 +42,10 @@ TEST_CASE("WeightStackExtension", "[stack]") { double const weight = 5.1; TestStack s; - s.addParticle(std::make_tuple(noData), std::tuple<double>{weight}); + auto const p0 = s.addParticle(std::make_tuple(noData), std::tuple<double>{weight}); CHECK(s.getEntries() == 1); + CHECK(p0.getWeight() == weight); } SECTION("write/read weights") { @@ -65,10 +65,10 @@ TEST_CASE("WeightStackExtension", "[stack]") { double const weight = 16; TestStack s; - // add 99 particles, each 10th particle is a nucleus with A=i and Z=A/2! + // add 99 particles, each 10th particle has weigth 100 for (int i = 0; i < 99; ++i) { auto p = s.addParticle(std::tuple<dummy_stack::NoData>{noData}); - p.setWeight(weight); + p.setWeight(((i + 1) % 10 == 0) ? 100 : weight); } CHECK(s.getEntries() == 99); @@ -78,7 +78,77 @@ TEST_CASE("WeightStackExtension", "[stack]") { v += p.getWeight(); p.erase(); } - CHECK(v == Approx(99 * weight)); + CHECK(v == Approx(90 * weight + 9 * 100)); + CHECK(s.getEntries() == 0); + } + + SECTION("stack operations") { + + TestStack s; + + double const weight = 16; + // add 99 particles, each 10th particle is a positron + // i=9, 19, 29, etc. are positron + // add 99 particles, each 10th particle has weigth 100 + for (int i = 0; i < 99; ++i) { + auto p = s.addParticle(std::tuple<dummy_stack::NoData>{noData}); + p.setWeight(((i + 1) % 10 == 0) ? 100 : weight); + } + + // copy + { + s.copy(s.begin() + 9, s.begin() + 10); // 100 to 16 + const auto& p9 = s.cbegin() + 9; + const auto& p10 = s.cbegin() + 10; + + CHECK(p9.getWeight() == 100); + CHECK(p10.getWeight() == 100); + } + + // copy + { + s.copy(s.begin() + 93, s.begin() + 9); // 16 to 100 + const auto& p93 = s.cbegin() + 93; + const auto& p9 = s.cbegin() + 9; + + CHECK(p9.getWeight() == 16); + CHECK(p93.getWeight() == 16); + } + + // copy + { + s.copy(s.begin() + 89, s.begin() + 79); // 100 to 100 + const auto& p89 = s.cbegin() + 89; + const auto& p79 = s.cbegin() + 79; + + CHECK(p89.getWeight() == 100); + CHECK(p79.getWeight() == 100); + } + + // invalid copy + { CHECK_THROWS(s.copy(s.begin(), s.end() + 1000)); } + + // swap + { + s.swap(s.begin() + 11, s.begin() + 10); + const auto& p11 = s.cbegin() + 11; // now: positron + const auto& p10 = s.cbegin() + 10; // now: electron + + CHECK(p11.getWeight() == 100); + CHECK(p10.getWeight() == 16); + } + + // swap two positrons + { + s.swap(s.begin() + 29, s.begin() + 59); + const auto& p29 = s.cbegin() + 29; + const auto& p59 = s.cbegin() + 59; + + CHECK(p29.getWeight() == 100); + CHECK(p59.getWeight() == 100); + } + + for (int i = 0; i < 99; ++i) s.last().erase(); CHECK(s.getEntries() == 0); } }