conanfile.txt

@@ -6,7 +6,7 @@ boost/1.76.0
 zlib/1.2.11
 proposal/7.4.2
 yaml-cpp/0.7.0
-arrow/8.0.1
+arrow/10.0.0
 cli11/1.9.1
 
 [build_requires]

corsika/detail/modules/epos/InteractionModel.inl

@@ -10,9 +10,8 @@
 
 #include <corsika/modules/epos/InteractionModel.hpp>
 #include <corsika/modules/epos/EposStack.hpp>
-
 #include <corsika/framework/geometry/Point.hpp>
-
+#include <corsika/framework/core/ParticleProperties.hpp>
 #include <corsika/framework/utility/COMBoost.hpp>
 #include <corsika/framework/utility/CorsikaData.hpp>
 
@@ -35,8 +34,8 @@ namespace corsika::epos {
         data_path_ = (std::string(corsika_data("EPOS").c_str()) + "/").c_str();
       }
       initialize();
+      setParticlesStable();
     }
-    setParticlesStable();
   }
 
   inline void InteractionModel::setParticlesStable() const {
@@ -46,10 +45,30 @@ namespace corsika::epos {
       if (!is_hadron(p)) continue;
       int const eid = convertToEposRaw(p);
       if (eid != 0) {
-        ::epos::nodcy_.nrnody = ::epos::nodcy_.nrnody + 1;
-        ::epos::nodcy_.nody[::epos::nodcy_.nrnody - 1] = eid;
+        // 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,
@@ -59,7 +78,7 @@ namespace corsika::epos {
     if (!is_nucleus(targetId) && targetId != Code::Neutron && targetId != Code::Proton) {
       return false;
     }
-    if (is_nucleus(targetId) && (get_nucleus_A(targetId) >= maxTargetMassNumber_)) {
+    if (is_nucleus(targetId) && (get_nucleus_A(targetId) >= get_nucleus_A(maxNucleus_))) {
       return false;
     }
     if ((minEnergyCoM_ > sqrtS) || (sqrtS > maxEnergyCoM_)) { return false; }
@@ -136,9 +155,10 @@ namespace corsika::epos {
     strcpy(::epos::fname_.fncs, CS.data);
     ::epos::nfname_.nfncs = CS.length;
 
-    // initialiazes maximum energy and mass
-    initializeEventCoM(Code::Lead, Lead::nucleus_A, Lead::nucleus_Z, Code::Lead,
-                       Lead::nucleus_A, Lead::nucleus_Z, 1_PeV);
+    // initializes maximum energy and mass
+    initializeEventCoM(
+        maxNucleus_, get_nucleus_A(maxNucleus_), get_nucleus_Z(maxNucleus_), maxNucleus_,
+        get_nucleus_A(maxNucleus_), get_nucleus_Z(maxNucleus_), maxEnergyCoM_);
   }
 
   inline void InteractionModel::initializeEventCoM(Code const idBeam, int const iBeamA,
@@ -284,7 +304,7 @@ namespace corsika::epos {
                          "beamA={}, "
                          "beamZ={} "
                          "targetId={}, "
-                         "ELab={:4.3f} GeV,",
+                         "ELab={:12.2f} GeV,",
                          BeamId, BeamA, BeamZ, TargetId, EnergyLab / 1_GeV);
 
     // read cross section from epos internal tables
@@ -300,7 +320,7 @@ namespace corsika::epos {
       int const iBeam = epos::getEposXSCode(
           BeamId); // 0 (can not interact, 1: pion-like, 2: proton-like, 3:kaon-like)
       CORSIKA_LOGGER_TRACE(logger_,
-                           "projectile cross section type={} "
+                           "readCrossSectionTableLab: projectile cross section type={} "
                            "(0: cannot interact, 1:pion, 2:baryon, 3:kaon)",
                            iBeam);
 
@@ -314,16 +334,25 @@ namespace corsika::epos {
 
     int iMode = 3; // 0: air, >0 not air
 
-    CORSIKA_LOGGER_DEBUG(logger_,
-                         "inside Epos "
+    CORSIKA_LOGGER_TRACE(logger_,
+                         "readCrossSectionTableLab: inside Epos "
                          "beamId={}, beamXS={}",
                          ::epos::hadr2_.idproj, ::epos::had10_.iclpro);
 
+    CORSIKA_LOGGER_TRACE(logger_,
+                         "readCrossSectionTableLab: calling Epos cross section with"
+                         "Ekin = {}, Abeam = {}, Atarget = {}, iMode = {}",
+                         Ekin, Abeam, Atarget, iMode);
+
     // cross section from table, FAST
     float sigProdEpos = ::epos::eposcrse_(Ekin, Abeam, Atarget, iMode);
     // sig-el from analytic calculation, no fast
     float sigElaEpos = ::epos::eposelacrse_(Ekin, Abeam, Atarget, iMode);
 
+    CORSIKA_LOGGER_TRACE(logger_,
+                         "readCrossSectionTableLab: result: sigProd = {}, sigEla = {}",
+                         sigProdEpos, sigElaEpos);
+
     return std::make_tuple(sigProdEpos * 1_mb, sigElaEpos * 1_mb);
   }
 

corsika/modules/epos/InteractionModel.hpp

@@ -89,16 +89,18 @@ namespace corsika::epos {
     void doInteraction(TSecondaries&, Code const projectileId, Code const targetId,
                        FourMomentum const& projectileP4, FourMomentum const& targetP4);
 
-    void initialize() const;
     void initializeEventCoM(Code const, int const, int const, Code const, int const,
                             int const, HEPEnergyType const) const;
     void initializeEventLab(Code const, int const, int const, Code const, int const,
                             int const, HEPEnergyType const) const;
     void configureParticles(Code const, int const, int const, Code const, int const,
                             int const) const;
-    void setParticlesStable() const;
 
   private:
+    // 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;
     inline static bool isInitialized_ = false;
 
     std::string data_path_;
@@ -109,8 +111,7 @@ namespace corsika::epos {
     std::shared_ptr<spdlog::logger> logger_ = get_logger("corsika_epos_Interaction");
     HEPEnergyType const minEnergyCoM_ = 6 * 1e9 * electronvolt;
     HEPEnergyType const maxEnergyCoM_ = 2.e6 * 1e9 * electronvolt;
-    static unsigned int constexpr maxTargetMassNumber_ = 20;
-    static unsigned int constexpr minNuclearTargetA_ = 4;
+    static Code constexpr maxNucleus_ = Code::Lead;
   };
 
 } // namespace corsika::epos

modules/epos/epos-sem-lhc.f

@@ -2430,6 +2430,7 @@ c      p=sqrt(max(0.,egy**2-ampro**2))
         maproj=mapr
         matarg=matg
         engy=egy
+        call Class("pscrse-1")
         if(matg.eq.1.and.mapr.eq.1)then
           if(iqq.eq.1)then !sig ela
             call psfz(1,gz2,0.)
@@ -2479,6 +2480,7 @@ c      p=sqrt(max(0.,egy**2-ampro**2))
         maproj=maprojsave
         matarg=matargsave
         engy=engysave
+        call Class("pscrse-2")
       else
         ye=log10(max(1.,egy/1.5))+1.
         je=min(5,int(ye))
@@ -2597,7 +2599,7 @@ c------------------------------------------------------------------------------
       else
         eposelacrse=pscrse(ek,mapro,matar,1)
       endif
-
+      
       return
       end
 

src/modules/epos/epos_codes.dat

@@ -19,20 +19,30 @@ PiPlus                   120     1       Pion
 PiMinus                 -120     1       Pion
 KPlus                    130     1       Kaon
 KMinus                  -130     1       Kaon
+K0                      230      0       CannotInteract
+K0Bar                  -230      0       CannotInteract
+K0Long                  -20      1       Kaon
+K0Short                  20      1       Kaon
 Neutron                1220      1       Baryon
 AntiNeutron            -1220     1       Baryon
 Proton                  1120     1       Baryon
 AntiProton             -1120     1       Baryon
-Eta                      220	 0       CannotInteract
-Lambda0                 2130	 0       CannotInteract
-Lambda0Bar             -2130	 0       CannotInteract
+Lambda0                 2130	 1       Baryon
+Lambda0Bar             -2130	 1       Baryon
+
+SigmaPlus               1130     1       Baryon
+SigmaMinusBar          -1130     1       Baryon
+Sigma0                  1230     1       Baryon
+Sigma0Bar              -1230     1       Baryon
+SigmaMinus              2230     1       Baryon
+SigmaPlusBar           -2230     1       Baryon
 
-SigmaPlus               1130     0       CannotInteract
-SigmaMinusBar          -1130     0       CannotInteract
-Sigma0                  1230     0       CannotInteract
-Sigma0Bar              -1230     0       CannotInteract
-SigmaMinus              2230     0       CannotInteract
-SigmaPlusBar           -2230     0       CannotInteract
+SigmaStarPlus           1131     0       CannotInteract
+SigmaStarMinusBar      -1131     0       CannotInteract
+SigmaStar0              1231     0       CannotInteract
+SigmaStar0Bar          -1231     0       CannotInteract
+SigmaStarMinus          2231     0       CannotInteract
+SigmaStarPlusBar       -2231     0       CannotInteract
 
 DeltaPlus               1121     0       CannotInteract
 DeltaMinusBar          -1121     0       CannotInteract
@@ -41,14 +51,24 @@ Delta0Bar              -1221     0       CannotInteract
 DeltaMinus              2221     0       CannotInteract
 DeltaPlusBar           -2221     0       CannotInteract
 
+DeltaPlusPlus           1111     0       CannotInteract
+DeltaMinusMinusBar     -1111     0       CannotInteract
+
 Xi0                     1330     0       CannotInteract
 XiMinus	                2330 	 0       CannotInteract
 Xi0Bar	               -1330	 0       CannotInteract
 XiPlusBar              -2330 	 0       CannotInteract
 
+XiStar0                 1331     0       CannotInteract
+XiStarMinus	            2331 	 0       CannotInteract
+XiStar0Bar	           -1331	 0       CannotInteract
+XiStarPlusBar          -2331 	 0       CannotInteract
+
 OmegaMinus              3331     0       CannotInteract
 OmegaPlusBar           -3331     0       CannotInteract
 
+
+Eta                      220	 0       CannotInteract
 EtaPrime                 330     0       CannotInteract
 Phi                      331     0       CannotInteract
 Omega                    221     0       CannotInteract
@@ -56,9 +76,6 @@ Rho0                     111     1       Pion
 RhoPlus                  121     0       CannotInteract
 RhoMinus                -121     0       CannotInteract
 
-DeltaPlusPlus           1111     0       CannotInteract
-DeltaMinusMinusBar     -1111     0       CannotInteract
-
 KStarPlus                131     0       CannotInteract
 KStarMinus              -131     0       CannotInteract
 KStar0                   231     0       CannotInteract

tests/modules/testEpos.cpp

@@ -66,8 +66,8 @@ TEST_CASE("EposBasics", "module,process") {
 
   SECTION("cross-section type") {
     CHECK(corsika::epos::getEposXSCode(Code::Electron) == 0);
-    CHECK(corsika::epos::getEposXSCode(Code::K0Long) == 0);
-    CHECK(corsika::epos::getEposXSCode(Code::SigmaPlus) == 0);
+    CHECK(corsika::epos::getEposXSCode(Code::K0Long) == 3);
+    CHECK(corsika::epos::getEposXSCode(Code::SigmaPlus) == 2);
     CHECK(corsika::epos::getEposXSCode(Code::KMinus) == 3);
     CHECK(corsika::epos::getEposXSCode(Code::PiMinus) == 1);
     CHECK(corsika::epos::getEposXSCode(Code::Proton) == 2);
@@ -148,7 +148,8 @@ TEST_CASE("Epos", "modules") {
     CHECK_FALSE(model.isValid(Code::Proton, Code::Electron, 100_GeV));
     CHECK(model.isValid(Code::Proton, Code::Hydrogen, 100_GeV));
     CHECK(model.isValid(Code::Proton, Code::Helium, 100_GeV));
-    CHECK_FALSE(model.isValid(Code::Proton, Code::Iron, 100_GeV));
+    CHECK_FALSE(model.isValid(Code::Proton, Code::Iron, 10_EeV));
+    CHECK_FALSE(model.isValid(Code::Proton, get_nucleus_code(240, 120), 10_EeV));
     CHECK(model.isValid(Code::Proton, Code::Oxygen, 100_GeV));
   }
 
@@ -239,51 +240,53 @@ TEST_CASE("Epos", "modules") {
     CHECK(xs_prod2 / 1_mb == Approx(1076.7).margin(3.1));
   }
 
-  /*
-    SECTION("InteractionInterface - invalid") {
-      Code const pid = Code::Electron;
-      HEPEnergyType const P0 = 10_TeV;
-      auto [stack, viewPtr] = setup::testing::setup_stack(
-          pid, P0, (setup::Environment::BaseNodeType* const)nodePtr, cs);
-      setup::StackView& view = *viewPtr;
-      CHECK_THROWS(model.doInteraction(
-          view, pid, Code::Oxygen,
-          {sqrt(static_pow<2>(P0) + static_pow<2>(get_mass(pid))), {cs, P0, 0_GeV,
-    0_GeV}}, {Oxygen::mass, {cs, 0_GeV, 0_GeV, 0_GeV}}));
-    }
-  */
-  /*
-    SECTION("InteractionInterface - nuclear projectile") {
-
-      HEPEnergyType const P0 = 10_TeV;
-      Code const pid = get_nucleus_code(40, 20);
-      auto [stack, viewPtr] = setup::testing::setup_stack(
-          pid, P0, (setup::Environment::BaseNodeType* const)nodePtr, cs);
-      MomentumVector plab =
-          MomentumVector(cs, {P0, 0_eV, 0_eV}); // this is secret knowledge about
-    setupStack setup::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}});
-
-      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 - invalid") {
+    Code const pid = Code::Electron;
+    HEPEnergyType const P0 = 10_TeV;
+    auto [stack, viewPtr] = setup::testing::setup_stack(
+        pid, P0, (DummyEnvironment::BaseNodeType* const)nodePtr, cs);
+    test::StackView& view = *viewPtr;
+    CHECK_THROWS(model.doInteraction(
+        view, pid, Code::Oxygen,
+        {sqrt(static_pow<2>(P0) + static_pow<2>(get_mass(pid))), {cs, P0, 0_GeV, 0_GeV}},
+        {Oxygen::mass, {cs, 0_GeV, 0_GeV, 0_GeV}}));
+  }
+
+  SECTION("InteractionInterface - nuclear projectile") {
+
+    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
+    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}});
+
+    //  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")
   {
@@ -316,4 +319,4 @@ TEST_CASE("Epos", "modules") {
     //  CHECK(length / 1_g * 1_cm * 1_cm ==
     //      Approx(30).margin(20)); // this is no physics validation
   }
-}
\ No newline at end of file
+}