diff --git a/corsika/detail/stack/NuclearStackExtension.inl b/corsika/detail/stack/NuclearStackExtension.inl index 5dab28ab2506181580fd510e28164c507b70b878..e8b72b1288cd529eca69bbe03577fcb513f81386 100644 --- a/corsika/detail/stack/NuclearStackExtension.inl +++ b/corsika/detail/stack/NuclearStackExtension.inl @@ -121,6 +121,13 @@ namespace corsika::nuclear_stack { return super_type::getCharge(); } + template <template <typename> class InnerParticleInterface, + typename StackIteratorInterface> + inline HEPEnergyType NuclearParticleInterface< + InnerParticleInterface, StackIteratorInterface>::getKineticEnergy() const { + return this->getEnergy() - this->getMass(); + } + template <template <typename> class InnerParticleInterface, typename StackIteratorInterface> inline int16_t NuclearParticleInterface< diff --git a/corsika/stack/NuclearStackExtension.hpp b/corsika/stack/NuclearStackExtension.hpp index 19e9a44c32aa401147d57585c89a7c1ab3ec6f4e..445350506811ece92db3697c72d96b4b52c2c5f5 100644 --- a/corsika/stack/NuclearStackExtension.hpp +++ b/corsika/stack/NuclearStackExtension.hpp @@ -92,10 +92,17 @@ namespace corsika::nuclear_stack { * Overwrite normal getParticleMass function with nuclear version */ HEPMassType getMass() const; + /** * Overwrite normal getParticleCharge function with nuclear version */ ElectricChargeType getCharge() const; + + /** + * Overwrite normal getKineticEnergy function with nuclear version + */ + HEPMassType getKineticEnergy() const; + /** * Overwirte normal getChargeNumber function with nuclear version **/ diff --git a/corsika/stack/VectorStack.hpp b/corsika/stack/VectorStack.hpp index e5a765781ae44f38b38da7d4135a6363c697f5d3..23211e4cb95f762cf6ae9c328b05a7a54a004473 100644 --- a/corsika/stack/VectorStack.hpp +++ b/corsika/stack/VectorStack.hpp @@ -95,6 +95,8 @@ namespace corsika { ElectricChargeType getCharge() const { return get_charge(this->getPID()); } + HEPEnergyType getKineticEnergy() const { return this->getEnergy() - this->getMass(); } + int16_t getChargeNumber() const { return get_charge_number(this->getPID()); } ///@} }; diff --git a/tests/stack/testNuclearStackExtension.cpp b/tests/stack/testNuclearStackExtension.cpp index 1e33f6c0570230bed7fc31baa13f0f2de1cb0c91..1b310ba4dee4f08f9b189f1e41da1e276181ffc6 100644 --- a/tests/stack/testNuclearStackExtension.cpp +++ b/tests/stack/testNuclearStackExtension.cpp @@ -17,6 +17,14 @@ using namespace corsika; #include <iostream> using namespace std; +template <typename TParticle> +HEPEnergyType kineticEnergy(TParticle const p) { + if (p.getPID() == Code::Nucleus) + return p.getEnergy() - get_nucleus_mass(p.getNuclearA(), p.getNuclearZ()); + else + return p.getEnergy() - get_mass(p.getPID()); +} + TEST_CASE("NuclearStackExtension", "[stack]") { logging::set_level(logging::level::info); @@ -64,13 +72,18 @@ TEST_CASE("NuclearStackExtension", "[stack]") { } SECTION("read nucleus") { + auto const A = 10; + auto const Z = 9; nuclear_stack::ParticleDataStack s; s.addParticle(std::make_tuple( Code::Nucleus, 1.5_GeV, MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}), - Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s, 10, 9)); + Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s, A, Z)); const auto pout = s.getNextParticle(); CHECK(pout.getPID() == Code::Nucleus); CHECK(pout.getEnergy() == 1.5_GeV); + CHECK(pout.getMass() == get_nucleus_mass(A, Z)); + CHECK(pout.getKineticEnergy() == kineticEnergy(pout)); + CHECK(pout.getKineticEnergy() > 0_GeV); CHECK(pout.getTime() == 100_s); CHECK(pout.getNuclearA() == 10); CHECK(pout.getNuclearZ() == 9); diff --git a/tests/stack/testVectorStack.cpp b/tests/stack/testVectorStack.cpp index 9767b6c847beff9097e873aed8a390ccb62a45fd..074c34dead76cf4988fb28b6dbc279f5a64ca333 100644 --- a/tests/stack/testVectorStack.cpp +++ b/tests/stack/testVectorStack.cpp @@ -17,6 +17,11 @@ using namespace corsika; using namespace std; +template <typename TParticle> +HEPEnergyType kineticEnergy(TParticle const p) { + return p.getEnergy() - get_mass(p.getPID()); +} + TEST_CASE("VectorStack", "[stack]") { logging::set_level(logging::level::info); @@ -37,6 +42,7 @@ TEST_CASE("VectorStack", "[stack]") { auto pout = s.getNextParticle(); CHECK(pout.getPID() == Code::Electron); CHECK(pout.getEnergy() == 1.5_GeV); + CHECK(pout.getKineticEnergy() == kineticEnergy(pout)); CHECK(pout.getTime() == 100_s); }