diff --git a/corsika/detail/framework/core/ParticleProperties.inl b/corsika/detail/framework/core/ParticleProperties.inl
index fc9fbf636241f4c83d364cb9a6da80c458a3ecc8..af41af0fc4e8e55ca126777cb541d8595b7a9dd4 100644
--- a/corsika/detail/framework/core/ParticleProperties.inl
+++ b/corsika/detail/framework/core/ParticleProperties.inl
@@ -14,17 +14,18 @@
namespace corsika {
- inline HEPEnergyType constexpr get_kinetic_energy_threshold(Code const code) {
+ inline HEPEnergyType constexpr get_kinetic_energy_propagation_threshold(
+ Code const code) {
if (is_nucleus(code)) return particle::detail::threshold_nuclei;
- return particle::detail::thresholds[static_cast<CodeIntType>(code)];
+ return particle::detail::propagation_thresholds[static_cast<CodeIntType>(code)];
}
- inline void constexpr set_kinetic_energy_threshold(Code const code,
- HEPEnergyType const val) {
+ inline void constexpr set_kinetic_energy_propagation_threshold(
+ Code const code, HEPEnergyType const val) {
if (is_nucleus(code))
particle::detail::threshold_nuclei = val;
else
- particle::detail::thresholds[static_cast<CodeIntType>(code)] = val;
+ particle::detail::propagation_thresholds[static_cast<CodeIntType>(code)] = val;
}
inline HEPMassType constexpr get_mass(Code const code) {
@@ -32,6 +33,15 @@ namespace corsika {
return particle::detail::masses[static_cast<CodeIntType>(code)];
}
+ inline HEPEnergyType constexpr get_energy_production_threshold(Code const p) {
+ return particle::detail::production_thresholds[static_cast<CodeIntType>(p)];
+ }
+
+ inline void constexpr set_energy_production_threshold(Code const p,
+ HEPEnergyType const val) {
+ particle::detail::production_thresholds[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/ParticleCut.inl b/corsika/detail/modules/ParticleCut.inl
index ac72089c166c3a183424b9bb0e6931e565e0634a..80f0244770a5c4b64372e123f7cb0bd0be6747f2 100644
--- a/corsika/detail/modules/ParticleCut.inl
+++ b/corsika/detail/modules/ParticleCut.inl
@@ -27,15 +27,15 @@ namespace corsika {
, doCutInv_(inv) {
for (auto p : get_all_particles()) {
if (is_hadron(p)) // nuclei are also hadrons
- set_kinetic_energy_threshold(p, eHadCut);
+ set_kinetic_energy_propagation_threshold(p, eHadCut);
else if (is_muon(p))
- set_kinetic_energy_threshold(p, eMuCut);
+ set_kinetic_energy_propagation_threshold(p, eMuCut);
else if (p == Code::Electron || p == Code::Positron)
- set_kinetic_energy_threshold(p, eEleCut);
+ set_kinetic_energy_propagation_threshold(p, eEleCut);
else if (p == Code::Photon)
- set_kinetic_energy_threshold(p, ePhoCut);
+ set_kinetic_energy_propagation_threshold(p, ePhoCut);
}
- set_kinetic_energy_threshold(Code::Nucleus, eHadCut);
+ set_kinetic_energy_propagation_threshold(Code::Nucleus, eHadCut);
CORSIKA_LOG_DEBUG(
"setting kinetic energy thresholds: electrons = {} GeV, photons = {} GeV, "
"hadrons = {} GeV, "
@@ -49,8 +49,10 @@ namespace corsika {
TArgs&&... outputArgs)
: TOutput(std::forward<TArgs>(outputArgs)...)
, doCutInv_(inv) {
- for (auto p : get_all_particles()) { set_kinetic_energy_threshold(p, eCut); }
- set_kinetic_energy_threshold(Code::Nucleus, eCut);
+ for (auto p : get_all_particles()) {
+ set_kinetic_energy_propagation_threshold(p, eCut);
+ }
+ set_kinetic_energy_propagation_threshold(Code::Nucleus, eCut);
CORSIKA_LOG_DEBUG("setting kinetic energy threshold {} GeV", eCut / 1_GeV);
}
@@ -61,7 +63,9 @@ namespace corsika {
TArgs&&... args)
: TOutput(std::forward<TArgs>(args)...)
, doCutInv_(inv) {
- set_kinetic_energy_thresholds(eCuts);
+ for (auto const& cut : eCuts) {
+ set_kinetic_energy_propagation_threshold(cut.first, cut.second);
+ }
CORSIKA_LOG_DEBUG("setting threshold particles individually");
}
@@ -74,9 +78,9 @@ namespace corsika {
if (is_nucleus(pid)) {
// calculate energy per nucleon
auto const ElabNuc = energyLab / get_nucleus_A(pid);
- return (ElabNuc < get_kinetic_energy_threshold(pid));
+ return (ElabNuc < get_kinetic_energy_propagation_threshold(pid));
} else {
- return (energyLab < get_kinetic_energy_threshold(pid));
+ return (energyLab < get_kinetic_energy_propagation_threshold(pid));
}
}
diff --git a/corsika/detail/modules/energy_loss/BetheBlochPDG.inl b/corsika/detail/modules/energy_loss/BetheBlochPDG.inl
index 570f830d69c3d257185f830f765190153ba4b518..eb31e373ee2b4fbb2c043fed4d899a43b82de4d2 100644
--- a/corsika/detail/modules/energy_loss/BetheBlochPDG.inl
+++ b/corsika/detail/modules/energy_loss/BetheBlochPDG.inl
@@ -183,7 +183,7 @@ namespace corsika {
auto const energy = vParticle.getKineticEnergy();
auto const energy_lim =
std::max(energy * 0.9, // either 10% relative loss max., or
- get_kinetic_energy_threshold(
+ get_kinetic_energy_propagation_threshold(
vParticle.getPID()) // energy thresholds globally defined
// for individual particles
* 0.99999 // need to go slightly below global e-cut to assure
diff --git a/corsika/detail/modules/proposal/ContinuousProcess.inl b/corsika/detail/modules/proposal/ContinuousProcess.inl
index 77f77e7795a9268e58da251e1b851aad2d0e3d2e..6054a0d5fd0d823779d8bbaeb6dc912be8078a32 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_production_threshold(
+ 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
@@ -134,7 +133,7 @@ namespace corsika::proposal {
auto const energy = vP.getEnergy();
auto const energy_lim =
std::max(energy * 0.9, // either 10% relative loss max., or
- get_kinetic_energy_threshold(
+ get_kinetic_energy_propagation_threshold(
code) // energy thresholds globally defined for individual particles
* 0.9999 // need to go slightly below global e-cut to assure removal
// in ParticleCut. This does not matter since at cut-time
diff --git a/corsika/detail/modules/proposal/Interaction.inl b/corsika/detail/modules/proposal/Interaction.inl
index d8916cbe976e77b8865d7d0edaacc6bfc0135e0b..9b2e23138a060e0d9e4105e9adb6b07b5df8020e 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_production_threshold(
+ 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 dd03aa106c309b6abe3e7f58b69922ebbb8ae875..be6772f7ae6eedf99cb3b8b76de21de0f1f3de31 100644
--- a/corsika/framework/core/ParticleProperties.hpp
+++ b/corsika/framework/core/ParticleProperties.hpp
@@ -31,7 +31,7 @@ namespace corsika {
*
* The properties of all particles are saved in static and flat
* arrays. There is a enum corsika::Code to identify each
- * particles, and each individual particles has its own static class,
+ * particle, and each individual particle has its own static class,
* which can be used to retrieve its physical properties.
*
* The properties of all elementary particles are accessible here. The data
@@ -54,6 +54,25 @@ namespace corsika {
* The names, relations and properties of all particles known to CORSIKA 8 are listed
* below.
*
+ * **Note** on energy threshold on particle production as well as particle propagation.
+ * The functions:
+ * @code {.cpp}
+ * HEPEnergyType constexpr get_energy_production_threshold(Code const);
+ * void constexpr set_energy_production_threshold(Code const, HEPEnergyType const);
+ * @endcode
+ * can be used to tune the transition where explicit production of new particles, e.g.
+ * in Bremsstrahlung, is simulated versus a continuous handling of low-energy particles
+ * as generic energy losses. The default value for all particle types is 1 MeV.
+ *
+ * Furthermore, the functions:
+ * @code {.cpp}
+ * HEPEnergyType constexpr get_kinetic_energy_propagation_threshold(Code const);
+ * void constexpr set_kinetic_energy_propagation_threshold(Code const, HEPEnergyType
+ * const);
+ * @endcode
+ * are used to discard low energy particle during tracking. The default value for all
+ * particle types is 1 GeV.
+ *
* @addtogroup Particles
* @{
*/
@@ -89,23 +108,39 @@ namespace corsika {
namespace corsika {
// forward declarations to be used in GeneratedParticleProperties
+
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
- void constexpr set_kinetic_energy_threshold(
- Code const, HEPEnergyType const); //!< set kinetic energy threshold below which the
- //!< particle is discarded
-
- 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);
- }
+
+ /**
+ * Get the kinetic energy propagation threshold.
+ *
+ * Particles are tracked only above the kinetic energy propagation threshold. Below
+ * this, they are discarded and removed. Sensible default values must be configured for
+ * a simulation.
+ */
+ HEPEnergyType constexpr get_kinetic_energy_propagation_threshold(Code const);
+
+ /**
+ * Set the kinetic energy propagation threshold object.
+ */
+ void constexpr set_kinetic_energy_propagation_threshold(Code const,
+ HEPEnergyType const);
+
+ /**
+ * Get the particle production energy threshold.
+ *
+ * The (total) energy below which a particle is only handled stoachastically (no
+ * production below this energy). This is for example important for stachastic discrete
+ * Bremsstrahlung versus low-enregy Bremsstrahlung as part of continuous energy losses.
+ */
+ HEPEnergyType constexpr get_energy_production_threshold(Code const); //!<
+
+ /**
+ * Set the particle production energy threshold.
+ */
+ void constexpr set_energy_production_threshold(Code const, HEPEnergyType const);
//! Particle code according to PDG, "Monte Carlo Particle Numbering Scheme"
PDGCode constexpr get_PDG(Code const);
diff --git a/examples/em_shower.cpp b/examples/em_shower.cpp
index b4a2c39d414290318a9ceed010a96cf25ddfbf02..615b840c18d53d1d9d6bfc4ec92c0203e9c3619b 100644
--- a/examples/em_shower.cpp
+++ b/examples/em_shower.cpp
@@ -101,6 +101,14 @@ 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::Photon, 10_MeV},
+ };
+ for (auto [pcode, energy] : energy_resolution)
+ set_energy_production_threshold(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 e06488c7b4da21840a26986e4514741d17a6c197..7fa7499166edec35f4aaf6c1e8f13bb19388e8db 100755
--- a/src/framework/core/code_generator.py
+++ b/src/framework/core/code_generator.py
@@ -20,8 +20,9 @@ mproton = 0.9382720813 # GeV
namespace = "corsika"
# IDs of Nuclei are 10LZZZAAAI
-nucleusIdStr = "10{L:01d}{Z:03d}{A:03d}{I:01d}" # used with .format(L=,Z=,A=,I=)
-nucleusIdOffset = int(nucleusIdStr.format(L=0,A=0,Z=0,I=0))
+# used with .format(L=,Z=,A=,I=)
+nucleusIdStr = "10{L:01d}{Z:03d}{A:03d}{I:01d}"
+nucleusIdOffset = int(nucleusIdStr.format(L=0, A=0, Z=0, I=0))
##############################################################
@@ -243,7 +244,7 @@ def read_nuclei_db(filename, particle_db, classnames):
"mass": mass, # in GeV
"electric_charge": electric_charge, # in e/3
"lifetime": lifetime,
- "ngc_code": int(nucleusIdStr.format(L=0,A=A,Z=Z,I=0)),
+ "ngc_code": int(nucleusIdStr.format(L=0, A=A, Z=Z, I=0)),
"A": A,
"Z": Z,
"isNucleus": True,
@@ -308,26 +309,28 @@ def gen_internal_enum(particle_db, nuclei_db):
" FirstParticle = 1, // if you want to loop over particles, you want to start with \"1\" \n") # identifier for eventual loops...
# non-nuclei loop
- for k in filter(lambda k: "ngc_code" in particle_db[k] and not particle_db[k]["isNucleus"],
+ for k in filter(lambda k: "ngc_code" in particle_db[k] and not particle_db[k]["isNucleus"],
particle_db):
last_ngc_id = particle_db[k]['ngc_code']
string += " {key:s} = {code:d},\n".format(key=k, code=last_ngc_id)
- string += " LastParticle = {:d},\n".format(last_ngc_id + 1) # identifier for eventual loops...
+ # identifier for eventual loops...
+ string += " LastParticle = {:d},\n".format(last_ngc_id + 1)
if last_ngc_id > 0x7fffffff: # does not fit into int32_t, this will never happen....
raise Exception(
"Integer overflow in internal particle code definition prevented!")
- if last_ngc_id + 1 >= nucleusIdOffset:
+ if last_ngc_id + 1 >= nucleusIdOffset:
raise Exception(
"Too many particles. Fix conflict with Code::Nuclear == {:d} (just increase id...) !").format(nucleusIdOffset)
# marker to mark the beginning of generic nuclear IDs: 1aaazzz
- string += " Nucleus = {:d},\n".format(nucleusIdOffset) # identifier for Nuclei
-
+ # identifier for Nuclei
+ string += " Nucleus = {:d},\n".format(nucleusIdOffset)
+
# nuclei loop
- for k in filter(lambda k: "ngc_code" in nuclei_db[k] and nuclei_db[k]["isNucleus"],
+ for k in filter(lambda k: "ngc_code" in nuclei_db[k] and nuclei_db[k]["isNucleus"],
nuclei_db):
last_ngc_id = nuclei_db[k]['ngc_code']
string += " {key:s} = {code:d},\n".format(key=k, code=last_ngc_id)
@@ -372,18 +375,26 @@ def gen_properties(particle_db):
string += "\n"
# particle masses table
- string += "static constexpr std::array<corsika::units::si::HEPMassType const, size> masses = {\n"
+ string += "static constexpr std::array<HEPMassType const, size> masses = {\n"
for p in particle_db.values():
- string += " {mass:e} * 1e9 * corsika::units::si::electronvolt, // {name:s}\n".format(
+ string += " {mass:e} * 1e9 * electronvolt, // {name:s}\n".format(
mass=p['mass'], name=p['name'])
string += "};\n\n"
# particle threshold table, initially set to 0
- string += "static std::array<corsika::units::si::HEPEnergyType, size> thresholds = {\n"
+ string += "static std::array<HEPEnergyType, size> propagation_thresholds = {\n"
for k in particle_db:
- string += " 0 * corsika::units::si::electronvolt, // {name:s}\n".format( name = k)
+ string += " 1e9 * electronvolt, // {name:s}\n".format(
+ name=k)
+ string += "};\n\n"
+ string += "static HEPEnergyType threshold_nuclei = 0_eV;\n"
+
+ # particle production_threshold table, initially set to 1 MeV
+ string += "static std::array<HEPEnergyType, size> production_thresholds = {\n"
+ for p in particle_db.values():
+ string += " 1e6 * electronvolt, // {name:s}\n".format(
+ name=p['name'])
string += "};\n\n"
- string += "static corsika::units::si::HEPEnergyType threshold_nuclei = 0_eV;\n"
# PDG code table
string += "static constexpr std::array<PDGCode, size> pdg_codes = {\n"
@@ -410,15 +421,15 @@ def gen_properties(particle_db):
# string += "};\n"
# lifetime
- #string += "static constexpr std::array<corsika::units::si::TimeType const, size> lifetime = {\n"
+ #string += "static constexpr std::array<TimeType const, size> lifetime = {\n"
string += "static constexpr std::array<double const, size> lifetime = {\n"
for p in particle_db.values():
if p['lifetime'] == float("Inf"):
- # * corsika::units::si::second, \n"
+ # * second, \n"
string += " std::numeric_limits<double>::infinity(), \n"
else:
string += " {tau:e}, \n".format(tau=p['lifetime'])
- #string += " {tau:e} * corsika::units::si::second, \n".format(tau = p['lifetime'])
+ #string += " {tau:e} * second, \n".format(tau = p['lifetime'])
string += "};\n"
# is Hadron flag
@@ -443,7 +454,8 @@ def gen_classes(particle_db, nuclei_db):
"/** @defgroup ParticleClasses \n"
" @{ */\n")
- common_db = OrderedDict(list(particle_db.items()) + list(nuclei_db.items()))
+ common_db = OrderedDict(
+ list(particle_db.items()) + list(nuclei_db.items()))
for cname in common_db:
if cname == "Nucleus":
string += "// skipping Nucleus"
@@ -538,7 +550,8 @@ def inc_end():
# Serialize particle_db into file
#
def serialize_particle_db(particle_db, nuclei_db, file):
- common_db = OrderedDict(list(particle_db.items()) + list(nuclei_db.items()))
+ common_db = OrderedDict(
+ list(particle_db.items()) + list(nuclei_db.items()))
pickle.dump(common_db, file)
@@ -555,11 +568,11 @@ if __name__ == "__main__":
print("\n code_generator.py: automatically produce particle properties from input files\n")
- names = read_class_names(sys.argv[3]) # re-names and conventions
- particle_db = OrderedDict() # the DB for pythia8 pdg particles
- read_pythia_db(sys.argv[1], particle_db, names) # pythia8 pdg DB
- nuclei_db = OrderedDict() # the DB for specific nuclei
- read_nuclei_db(sys.argv[2], nuclei_db, names) # list of nuclei
+ names = read_class_names(sys.argv[3]) # re-names and conventions
+ particle_db = OrderedDict() # the DB for pythia8 pdg particles
+ read_pythia_db(sys.argv[1], particle_db, names) # pythia8 pdg DB
+ nuclei_db = OrderedDict() # the DB for specific nuclei
+ read_nuclei_db(sys.argv[2], nuclei_db, names) # list of nuclei
with open("GeneratedParticleProperties.inc", "w") as f:
print(inc_start(), file=f)
diff --git a/tests/framework/testParticles.cpp b/tests/framework/testParticles.cpp
index eaa22922212be3badea612ecb262e13249025eb1..52648ba7120aabb8705bc679323f8322b96ac64c 100644
--- a/tests/framework/testParticles.cpp
+++ b/tests/framework/testParticles.cpp
@@ -88,13 +88,18 @@ TEST_CASE("ParticleProperties", "[Particles]") {
(Approx(2.1970332555864364e-06).epsilon(1e-5)));
}
- SECTION("Energy threshold") {
+ SECTION("Energy thresholds") {
//! by default energy thresholds are set to zero
- CHECK(get_kinetic_energy_threshold(Electron::code) == 0_GeV);
+ CHECK(get_kinetic_energy_propagation_threshold(Electron::code) == 1_GeV);
+ set_kinetic_energy_propagation_threshold(Electron::code, 10_GeV);
+ CHECK_FALSE(get_kinetic_energy_propagation_threshold(Code::Electron) == 1_GeV);
+ CHECK(get_kinetic_energy_propagation_threshold(Code::Electron) == 10_GeV);
- set_kinetic_energy_threshold(Electron::code, 10_GeV);
- CHECK_FALSE(get_kinetic_energy_threshold(Code::Electron) == 1_GeV);
- CHECK(get_kinetic_energy_threshold(Code::Electron) == 10_GeV);
+ //! by default energy thresholds are set to zero
+ CHECK(get_energy_production_threshold(Neutron::code) == 1_MeV);
+ set_energy_production_threshold(Neutron::code, 1_GeV);
+ CHECK_FALSE(get_energy_production_threshold(Code::Neutron) == 1_MeV);
+ CHECK(get_energy_production_threshold(Code::Neutron) == 1_GeV);
}
SECTION("Particle groups: electromagnetic") {
diff --git a/tests/modules/testParticleCut.cpp b/tests/modules/testParticleCut.cpp
index e12ac460144b571c70951498179b464c40917819..3d3c7e666e6c29928a6880119ef267ec8cbc26c5 100644
--- a/tests/modules/testParticleCut.cpp
+++ b/tests/modules/testParticleCut.cpp
@@ -164,11 +164,12 @@ TEST_CASE("ParticleCut", "process,continuous,secondary") {
SECTION("cut low energy: reset thresholds of arbitrary set of particles") {
ParticleCut cut({{Code::Electron, 5_MeV}, {Code::Positron, 50_MeV}}, false);
- CHECK(get_kinetic_energy_threshold(Code::Electron) !=
- get_kinetic_energy_threshold(Code::Positron));
- CHECK_FALSE(get_kinetic_energy_threshold(Code::Electron) == Electron::mass);
+ CHECK(get_kinetic_energy_propagation_threshold(Code::Electron) !=
+ get_kinetic_energy_propagation_threshold(Code::Positron));
+ CHECK_FALSE(get_kinetic_energy_propagation_threshold(Code::Electron) ==
+ Electron::mass);
// test default values still correct
- CHECK(get_kinetic_energy_threshold(Code::Proton) == 5_GeV);
+ CHECK(get_kinetic_energy_propagation_threshold(Code::Proton) == 5_GeV);
}
SECTION("cut on time") {