From c157f6ac497a748e595bf8b29d401c36117aed9d Mon Sep 17 00:00:00 2001
From: Felix Riehn <friehn@lip.pt>
Date: Wed, 7 Aug 2024 11:07:28 +0000
Subject: [PATCH] Resolve "use epos internal decays for exotic resonances"
---
applications/c8_air_shower.cpp | 3 +-
.../detail/modules/epos/InteractionModel.inl | 46 ++++++++++++-
corsika/modules/Epos.hpp | 6 +-
corsika/modules/epos/InteractionModel.hpp | 5 +-
corsika/setup/SetupC7trackedParticles.hpp | 25 +++++++
tests/modules/testEpos.cpp | 68 ++++---------------
6 files changed, 92 insertions(+), 61 deletions(-)
create mode 100644 corsika/setup/SetupC7trackedParticles.hpp
diff --git a/applications/c8_air_shower.cpp b/applications/c8_air_shower.cpp
index 10a9f4867..40a762726 100644
--- a/applications/c8_air_shower.cpp
+++ b/applications/c8_air_shower.cpp
@@ -76,6 +76,7 @@
#include <corsika/setup/SetupStack.hpp>
#include <corsika/setup/SetupTrajectory.hpp>
+#include <corsika/setup/SetupC7trackedParticles.hpp>
#include <boost/filesystem.hpp>
@@ -425,7 +426,7 @@ int main(int argc, char** argv) {
std::make_shared<corsika::qgsjetII::Interaction>()};
} else if (modelStr == "EPOS-LHC") {
heModel = DynamicInteractionProcess<StackType>{
- std::make_shared<corsika::epos::Interaction>()};
+ std::make_shared<corsika::epos::Interaction>(corsika::setup::C7trackedParticles)};
} else {
CORSIKA_LOG_CRITICAL("invalid choice \"{}\"; also check argument parser", modelStr);
return EXIT_FAILURE;
diff --git a/corsika/detail/modules/epos/InteractionModel.inl b/corsika/detail/modules/epos/InteractionModel.inl
index 6d2b94c61..cf35f1aca 100644
--- a/corsika/detail/modules/epos/InteractionModel.inl
+++ b/corsika/detail/modules/epos/InteractionModel.inl
@@ -24,7 +24,8 @@
namespace corsika::epos {
- inline InteractionModel::InteractionModel(std::string const& dataPath,
+ inline InteractionModel::InteractionModel(std::set<Code> vList,
+ std::string const& dataPath,
bool const epos_printout_on)
: data_path_(dataPath)
, epos_listing_(epos_printout_on) {
@@ -36,7 +37,45 @@ namespace corsika::epos {
data_path_ = (std::string(corsika_data("EPOS").c_str()) + "/").c_str();
}
initialize();
+ if (vList.empty()) {
+ CORSIKA_LOGGER_DEBUG(logger_,
+ "set all particles known to CORSIKA stable inside EPOS..");
+ setParticleListStable(get_all_particles());
+ } else {
+ CORSIKA_LOGGER_DEBUG(logger_, "set specific particles stable inside EPOS..");
+ setParticleListStable(vList);
+ }
+ }
+ }
+
+ inline void InteractionModel::setParticleListStable(std::set<Code> vPartList) const {
+ for (auto& p : vPartList) {
+ int const eid = convertToEposRaw(p);
+ if (eid != 0) {
+ // LCOV_EXCL_START
+ // this is only a safeguard against messing up the epos internals by initializing
+ // more than once.
+ unsigned int const n_particles_stable_epos =
+ ::epos::nodcy_.nrnody; // avoid waring -Wsign-compare
+ if (n_particles_stable_epos < ::epos::mxnody) {
+ CORSIKA_LOGGER_TRACE(logger_, "setting {} with EposId={} stable inside EPOS.",
+ p, eid);
+ ::epos::nodcy_.nrnody = ::epos::nodcy_.nrnody + 1;
+ ::epos::nodcy_.nody[::epos::nodcy_.nrnody - 1] = eid;
+ } else {
+ CORSIKA_LOGGER_ERROR(logger_, "List of stable particles too long for Epos!");
+ throw std::runtime_error("Epos initialization error!");
+ }
+ // LCOV_EXCL_STOP
+ } else {
+ CORSIKA_LOG_TRACE(
+ "particle conversion Corsika-->Epos not known for {}. Using {}. Setting "
+ "unstable in Epos!",
+ p, eid);
+ }
}
+ CORSIKA_LOGGER_DEBUG(logger_, "set {} particles stable inside Epos",
+ ::epos::nodcy_.nrnody);
}
inline bool InteractionModel::isValid(Code const projectileId, Code const targetId,
@@ -103,7 +142,8 @@ namespace corsika::epos {
::epos::othe2_.iframe = 11; // cms frame
// decay settings
- ::epos::othe2_.idecay = 0; // no decays in epos
+ // activate decays in epos for particles defined by set_stable/set_unstable
+ // ::epos::othe2_.idecay = 0; // no decays in epos
// set paths to tables in corsika data
::epos::datadir BASE(data_path_);
@@ -417,7 +457,7 @@ namespace corsika::epos {
if (is_nucleus(targetId)) {
targetA = get_nucleus_A(targetId);
targetZ = get_nucleus_Z(targetId);
- CORSIKA_LOGGER_DEBUG(logger_, "target: A={}, Z={} ", beamA, beamZ);
+ CORSIKA_LOGGER_DEBUG(logger_, "target: A={}, Z={} ", targetA, targetZ);
}
initializeEventCoM(projectileId, beamA, beamZ, targetId, targetA, targetZ, sqrtSNN);
diff --git a/corsika/modules/Epos.hpp b/corsika/modules/Epos.hpp
index 0fb3db624..9610edb38 100644
--- a/corsika/modules/Epos.hpp
+++ b/corsika/modules/Epos.hpp
@@ -26,5 +26,9 @@ namespace corsika::epos {
* The epos::InteractionModel is wrapped as an InteractionProcess here in order
* to provide all the functions for ProcessSequence.
*/
- class Interaction : public InteractionModel, public InteractionProcess<Interaction> {};
+ class Interaction : public InteractionModel, public InteractionProcess<Interaction> {
+ public:
+ Interaction(std::set<Code> const& stableList)
+ : InteractionModel{stableList} {};
+ };
} // namespace corsika::epos
diff --git a/corsika/modules/epos/InteractionModel.hpp b/corsika/modules/epos/InteractionModel.hpp
index dcb841895..388298a4b 100644
--- a/corsika/modules/epos/InteractionModel.hpp
+++ b/corsika/modules/epos/InteractionModel.hpp
@@ -9,6 +9,7 @@
#pragma once
#include <tuple>
+#include <set>
#include <corsika/framework/core/ParticleProperties.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp>
@@ -20,7 +21,7 @@ namespace corsika::epos {
class InteractionModel {
public:
- InteractionModel(std::string const& dataPath = "",
+ InteractionModel(std::set<Code> = {}, std::string const& dataPath = "",
bool const epos_printout_on = false);
~InteractionModel();
@@ -101,7 +102,7 @@ namespace corsika::epos {
// initialize and setParticlesStable are private since they can only be called once at
// the beginning and are already called in the constructor!
void initialize() const;
- void setParticlesStable() const;
+ void setParticleListStable(std::set<Code>) const;
inline static bool isInitialized_ = false;
std::string data_path_;
diff --git a/corsika/setup/SetupC7trackedParticles.hpp b/corsika/setup/SetupC7trackedParticles.hpp
new file mode 100644
index 000000000..82c284431
--- /dev/null
+++ b/corsika/setup/SetupC7trackedParticles.hpp
@@ -0,0 +1,25 @@
+/*
+ * (c) Copyright 2024 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+namespace corsika::setup {
+
+ /**
+ \file SetupC7TrackedParticles.hpp
+
+ set of particles that should be tracked by corsika. all others should decay.
+ this is the corsika 7 configuration
+ */
+
+ std::set<Code> C7trackedParticles{
+ Code::Proton, Code::Neutron, Code::AntiProton, Code::AntiNeutron,
+ Code::PiPlus, Code::PiMinus, Code::Pi0, Code::KPlus,
+ Code::KMinus, Code::K0Long, Code::K0Short, Code::Lambda0,
+ Code::Lambda0Bar, Code::SigmaPlus, Code::SigmaPlusBar, Code::SigmaMinus,
+ Code::SigmaMinusBar, Code::Xi0, Code::Xi0Bar, Code::OmegaMinus,
+ Code::OmegaPlusBar, Code::MuPlus, Code::MuMinus};
+} // namespace corsika::setup
\ No newline at end of file
diff --git a/tests/modules/testEpos.cpp b/tests/modules/testEpos.cpp
index 4b984bbca..3235d8012 100644
--- a/tests/modules/testEpos.cpp
+++ b/tests/modules/testEpos.cpp
@@ -225,71 +225,31 @@ TEST_CASE("Epos", "modules") {
{Oxygen::mass, {cs, 0_GeV, 0_GeV, 0_GeV}}));
}
- SECTION("InteractionInterface - nuclear projectile") {
+ SECTION("InteractionInterface - valid projectile target combinations") {
+
+ HEPMomentumType const P0 = 10_TeV;
+ Code const projectileId =
+ GENERATE(Code::Proton, Code::PiPlus, Code::KPlus, Code::Iron);
+ Code const targetId = GENERATE(Code::Proton, Code::Neutron, Code::Nitrogen);
- HEPEnergyType const P0 = 10_TeV;
- Code const pid = get_nucleus_code(40, 20);
auto [stack, viewPtr] = setup::testing::setup_stack(
- pid, P0, (DummyEnvironment::BaseNodeType* const)nodePtr, cs);
- MomentumVector plab =
- MomentumVector(cs, {P0, 0_eV, 0_eV}); // this is secret knowledge about
+ projectileId, P0, (DummyEnvironment::BaseNodeType* const)nodePtr, cs);
test::StackView& view = *viewPtr;
// @todo This is very obscure since it fails for -O2, but for both clang and gcc ???
- model.doInteraction(view, pid, Code::Oxygen,
- {sqrt(static_pow<2>(P0) + static_pow<2>(get_mass(pid))), plab},
- {Oxygen::mass, {cs, 0_GeV, 0_GeV, 0_GeV}});
+ model.doInteraction(
+ view, projectileId, targetId,
+ {calculate_total_energy(P0, get_mass(projectileId)), {cs, {P0, 0_eV, 0_eV}}},
+ {get_mass(targetId), {cs, 0_GeV, 0_GeV, 0_GeV}});
// simply check if stack is not empty after the event. Energy and momentum
// conservation will be tested elsewhere
CHECK(view.getSize() > 0);
-
- // auto const pSum = sumMomentum(view, cs);
-
- // CHECK(pSum.getComponents(cs).getX() / P0 == Approx(1).margin(0.05));
- // CHECK(pSum.getComponents(cs).getY() / 1_GeV ==
- // Approx(0).margin(0.5)); // this is not physics validation
- // CHECK(pSum.getComponents(cs).getZ() / 1_GeV ==
- // Approx(0).margin(0.5)); // this is not physics validation
-
- // CHECK((pSum - plab).getNorm() / 1_GeV ==
- // Approx(0).margin(plab.getNorm() * 0.05 / 1_GeV));
- // CHECK(pSum.getNorm() / P0 == Approx(1).margin(0.05));
- // [[maybe_unused]] const GrammageType length =
- // model.getInteractionLength(particle);
- // CHECK(length / 1_g * 1_cm * 1_cm ==
- // Approx(30).margin(20)); // this is no physics validation
}
- // SECTION("InteractionInterface")
- {
- HEPEnergyType const P0 = 10_TeV;
- Code const pid = Code::Proton;
- auto [stack, viewPtr] = setup::testing::setup_stack(
- pid, P0, (DummyEnvironment::BaseNodeType* const)nodePtr, cs);
- MomentumVector plab =
- MomentumVector(cs, {P0, 0_eV, 0_eV}); // this is secret knowledge about
- test::StackView& view = *viewPtr;
+ SECTION("Decay config") {
+ logging::set_level(logging::level::debug);
- // @todo This is very obscure since it fails for -O2, but for both clang and gcc ???
- model.doInteraction(view, pid, Code::Oxygen,
- {sqrt(static_pow<2>(P0) + static_pow<2>(get_mass(pid))), plab},
- {Oxygen::mass, {cs, 0_GeV, 0_GeV, 0_GeV}});
-
- auto const pSum = sumMomentum(view, cs);
-
- CHECK(pSum.getComponents(cs).getX() / P0 == Approx(1).margin(0.05));
- CHECK(pSum.getComponents(cs).getY() / 1_GeV ==
- Approx(0).margin(0.5)); // this is not physics validation
- CHECK(pSum.getComponents(cs).getZ() / 1_GeV ==
- Approx(0).margin(0.5)); // this is not physics validation
-
- CHECK((pSum - plab).getNorm() / 1_GeV ==
- Approx(0).margin(plab.getNorm() * 0.05 / 1_GeV));
- CHECK(pSum.getNorm() / P0 == Approx(1).margin(0.05));
- // [[maybe_unused]] const GrammageType length =
- // model.getInteractionLength(particle);
- // CHECK(length / 1_g * 1_cm * 1_cm ==
- // Approx(30).margin(20)); // this is no physics validation
+ InteractionModel model(std::set<Code>{Code::Proton, Code::PiPlus, Code::KPlus});
}
}
--
GitLab