From 213d7c3c1b3f8120a07856b8ffb862dc49a4c20c Mon Sep 17 00:00:00 2001
From: Maximilian Sackel <maximilian.sackel@udo.edu>
Date: Mon, 26 Apr 2021 12:22:09 +0200
Subject: [PATCH] :sparkles: Adding a particle resolution. Energies below the
resolution limit are treated only continuously and no longer stochastically.
---
corsika/detail/framework/core/ParticleProperties.inl | 9 +++++++++
corsika/detail/modules/proposal/ContinuousProcess.inl | 5 ++---
corsika/detail/modules/proposal/Interaction.inl | 5 ++---
corsika/framework/core/ParticleProperties.hpp | 11 +++++++++++
examples/em_shower.cpp | 9 +++++++++
src/framework/core/code_generator.py | 6 ++++++
6 files changed, 39 insertions(+), 6 deletions(-)
diff --git a/corsika/detail/framework/core/ParticleProperties.inl b/corsika/detail/framework/core/ParticleProperties.inl
index fc9fbf636..76229398e 100644
--- a/corsika/detail/framework/core/ParticleProperties.inl
+++ b/corsika/detail/framework/core/ParticleProperties.inl
@@ -32,6 +32,15 @@ namespace corsika {
return particle::detail::masses[static_cast<CodeIntType>(code)];
}
+ inline HEPEnergyType constexpr get_energy_resolution(Code const p) {
+ return particle::detail::resolutions[static_cast<CodeIntType>(p)];
+ }
+
+ inline void constexpr set_energy_resolution(Code const p, HEPEnergyType const val) {
+ particle::detail::resolutions[static_cast<CodeIntType>(p)] = val;
+ }
+
+
inline bool constexpr is_charged(Code const c) { return get_charge_number(c) != 0; }
inline bool constexpr is_nucleus(Code const code) { return code >= Code::Nucleus; }
diff --git a/corsika/detail/modules/proposal/ContinuousProcess.inl b/corsika/detail/modules/proposal/ContinuousProcess.inl
index 77f77e779..107f743ff 100644
--- a/corsika/detail/modules/proposal/ContinuousProcess.inl
+++ b/corsika/detail/modules/proposal/ContinuousProcess.inl
@@ -29,9 +29,8 @@ namespace corsika::proposal {
// interpolate the crosssection for given media and energy cut. These may
// take some minutes if you have to build the tables and cannot read the
// from disk
- auto const emCut =
- get_kinetic_energy_threshold(code) +
- get_mass(code); //! energy thresholds globally defined for individual particles
+ auto const emCut = get_energy_resolution(
+ code); //! energy resolutions globally defined for individual particles
auto c = p_cross->second(media.at(comp.getHash()), emCut);
// Use higland multiple scattering and deactivate stochastic deflection by
diff --git a/corsika/detail/modules/proposal/Interaction.inl b/corsika/detail/modules/proposal/Interaction.inl
index d8916cbe9..7cc43e298 100644
--- a/corsika/detail/modules/proposal/Interaction.inl
+++ b/corsika/detail/modules/proposal/Interaction.inl
@@ -32,9 +32,8 @@ namespace corsika::proposal {
// interpolate the crosssection for given media and energy cut. These may
// take some minutes if you have to build the tables and cannot read the
// from disk
- auto const emCut =
- get_kinetic_energy_threshold(code) +
- get_mass(code); //! energy thresholds globally defined for individual particles
+ auto const emCut = get_energy_resolution(
+ code); //! energy resolutions globally defined for individual particles
auto c = p_cross->second(media.at(comp.getHash()), emCut);
diff --git a/corsika/framework/core/ParticleProperties.hpp b/corsika/framework/core/ParticleProperties.hpp
index dd03aa106..083005058 100644
--- a/corsika/framework/core/ParticleProperties.hpp
+++ b/corsika/framework/core/ParticleProperties.hpp
@@ -92,6 +92,7 @@ namespace corsika {
int16_t constexpr get_charge_number(Code const); //!< electric charge in units of e
ElectricChargeType constexpr get_charge(Code const); //!< electric charge
HEPMassType constexpr get_mass(Code const); //!< mass
+
HEPEnergyType constexpr get_kinetic_energy_threshold(
Code const); //!< get kinetic energy threshold below which the particle is
//!< discarded, by default set to zero
@@ -102,11 +103,21 @@ namespace corsika {
inline void set_kinetic_energy_threshold(std::pair<Code const, HEPEnergyType const> p) {
set_kinetic_energy_threshold(p.first, p.second);
}
+
inline void set_kinetic_energy_thresholds(
std::unordered_map<Code const, HEPEnergyType const> const& eCuts) {
for (auto v : eCuts) set_kinetic_energy_threshold(v);
}
+ HEPEnergyType constexpr get_energy_resolution(
+ Code const); //!< get energy resolution below which a particle is only
+ //!< handled stoachastically
+ void constexpr set_energy_resolution(
+ Code const, HEPEnergyType const); //!< get energy resolution below which a
+ //!< particle is only handled
+ //!< stoachastically
+
+
//! Particle code according to PDG, "Monte Carlo Particle Numbering Scheme"
PDGCode constexpr get_PDG(Code const);
PDGCode constexpr get_PDG(unsigned int const A, unsigned int const Z);
diff --git a/examples/em_shower.cpp b/examples/em_shower.cpp
index b4a2c39d4..0630bbf2e 100644
--- a/examples/em_shower.cpp
+++ b/examples/em_shower.cpp
@@ -101,6 +101,15 @@ int main(int argc, char** argv) {
env, AtmosphereId::LinsleyUSStd, center, Medium::AirDry1Atm,
MagneticFieldVector{rootCS, 0_T, 50_uT, 0_T});
+ std::unordered_map<Code, HEPEnergyType> energy_resolution = {
+ { Code::Electron, 10_MeV },
+ { Code::Positron, 10_MeV },
+ { Code::Gamma, 10_MeV },
+ };
+ for (auto [pcode, energy] : energy_resolution)
+ set_energy_resolution(pcode, energy);
+
+
// setup particle stack, and add primary particle
setup::Stack stack;
stack.clear();
diff --git a/src/framework/core/code_generator.py b/src/framework/core/code_generator.py
index e06488c7b..8d5f2788c 100755
--- a/src/framework/core/code_generator.py
+++ b/src/framework/core/code_generator.py
@@ -385,6 +385,12 @@ def gen_properties(particle_db):
string += "};\n\n"
string += "static corsika::units::si::HEPEnergyType threshold_nuclei = 0_eV;\n"
+ # particle resolution table, initially set to 1 MeV
+ string += "static std::array<corsika::units::si::HEPEnergyType, size> resolutions = {\n"
+ for p in particle_db.values():
+ string += " 1e6 * corsika::units::si::electronvolt, // {name:s}\n".format(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():
--
GitLab