From 007626fe052ac15ae9f5c4101e1b36a53faa6dc7 Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Mon, 18 Jan 2021 15:25:44 +0000
Subject: [PATCH] added energy threshold to particle properties

---
 .../detail/framework/core/ParticleProperties.inl    |  8 ++++++++
 corsika/framework/core/ParticleProperties.hpp       | 13 +++++++++++++
 src/framework/core/pdxml_reader.py                  |  8 +++++++-
 tests/framework/testParticles.cpp                   |  9 +++++++++
 4 files changed, 37 insertions(+), 1 deletion(-)

diff --git a/corsika/detail/framework/core/ParticleProperties.inl b/corsika/detail/framework/core/ParticleProperties.inl
index 6fc0bea95..b478b1e17 100644
--- a/corsika/detail/framework/core/ParticleProperties.inl
+++ b/corsika/detail/framework/core/ParticleProperties.inl
@@ -13,6 +13,14 @@
 
 namespace corsika {
 
+  HEPEnergyType constexpr get_energy_threshold(Code const p) {
+    return particle::detail::thresholds[static_cast<CodeIntType>(p)];
+  }
+
+  void constexpr set_energy_threshold(Code const p, HEPEnergyType const val) {
+    particle::detail::thresholds[static_cast<CodeIntType>(p)] = val;
+  }
+
   HEPMassType constexpr get_mass(Code const p) {
     if (p == Code::Nucleus)
       throw std::runtime_error("Cannot GetMass() of particle::Nucleus -> unspecified");
diff --git a/corsika/framework/core/ParticleProperties.hpp b/corsika/framework/core/ParticleProperties.hpp
index 2deb03a86..cfc8c5583 100644
--- a/corsika/framework/core/ParticleProperties.hpp
+++ b/corsika/framework/core/ParticleProperties.hpp
@@ -50,6 +50,19 @@ namespace corsika {
   int16_t constexpr get_charge_number(Code);     //!< electric charge in units of e
   ElectricChargeType constexpr get_charge(Code); //!< electric charge
   HEPMassType constexpr get_mass(Code);          //!< mass
+  HEPEnergyType constexpr get_energy_threshold(
+      Code const); //!< get energy threshold below which the particle is discarded, by
+                   //!< default set to particle mass
+  void constexpr set_energy_threshold(
+      Code const, HEPEnergyType const); //!< set energy threshold below which the particle
+                                        //!< is discarded
+
+  inline void set_energy_threshold(std::pair<Code const, HEPEnergyType const>p){
+    set_energy_threshold(p.first, p.second);
+  }
+  inline void set_energy_thresholds(std::unordered_map<Code const,HEPEnergyType const> const& eCuts){
+    for (auto v : eCuts) set_energy_threshold(v);
+  }
 
   //! Particle code according to PDG, "Monte Carlo Particle Numbering Scheme"
   PDGCode constexpr get_PDG(Code);
diff --git a/src/framework/core/pdxml_reader.py b/src/framework/core/pdxml_reader.py
index fcd0b1341..2520cdfe4 100755
--- a/src/framework/core/pdxml_reader.py
+++ b/src/framework/core/pdxml_reader.py
@@ -330,7 +330,13 @@ def gen_properties(particle_db):
     for p in particle_db.values():
         string += "  {mass:e} * 1e9 * corsika::units::si::electronvolt, // {name:s}\n".format(mass = p['mass'], name = p['name'])
     string += "};\n\n"
-                   
+
+    # particle threshold table, initially set to the particle mass
+    string += "static std::array<corsika::units::si::HEPEnergyType, size> thresholds = {\n"    
+    for p in particle_db.values():
+        string += "  {mass:e} * 1e9 * corsika::units::si::electronvolt, // {name:s}\n".format(mass = p['mass'], name = p['name'])
+    string += "};\n\n"
+
     # PDG code table
     string += "static constexpr std::array<PDGCode, size> pdg_codes = {\n"
     for p in particle_db.keys():
diff --git a/tests/framework/testParticles.cpp b/tests/framework/testParticles.cpp
index dd999ecf8..9c80701bd 100644
--- a/tests/framework/testParticles.cpp
+++ b/tests/framework/testParticles.cpp
@@ -80,6 +80,15 @@ TEST_CASE("ParticleProperties", "[Particles]") {
           (Approx(2.1970332555864364e-06).epsilon(1e-5)));
   }
 
+  SECTION("Energy threshold") {
+    //! by default energy thresholds are set to particle mass
+    CHECK(get_energy_threshold(Electron::code) / Electron::mass == Approx(1));
+
+    set_energy_threshold(Electron::code,10_GeV);
+    CHECK_FALSE(get_energy_threshold(Code::Electron) == 1_GeV);
+    CHECK(get_energy_threshold(Code::Electron) == 10_GeV);
+  }
+
   SECTION("Particle groups: electromagnetic") {
     CHECK(is_em(Code::Gamma));
     CHECK(is_em(Code::Electron));
-- 
GitLab