Resolve "use epos internal decays for exotic resonances"

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;

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);
 

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

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_;

corsika/setup/SetupC7trackedParticles.hpp

+/*
+ * (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

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});
   }
 }