From cd4188fa716b8dedf2e303eee81fb1222d6e8f2c Mon Sep 17 00:00:00 2001 From: ralfulrich <ralf.ulrich@kit.edu> Date: Mon, 20 Dec 2021 16:11:57 +0100 Subject: [PATCH] rename quantities --- .../framework/core/ParticleProperties.inl | 21 +++--- corsika/detail/modules/ParticleCut.inl | 24 ++++--- .../modules/energy_loss/BetheBlochPDG.inl | 2 +- .../modules/proposal/ContinuousProcess.inl | 4 +- .../detail/modules/proposal/Interaction.inl | 2 +- corsika/framework/core/ParticleProperties.hpp | 71 ++++++++++++------- examples/em_shower.cpp | 9 ++- src/framework/core/code_generator.py | 61 +++++++++------- tests/framework/testParticles.cpp | 8 +-- tests/modules/testParticleCut.cpp | 9 +-- 10 files changed, 123 insertions(+), 88 deletions(-) diff --git a/corsika/detail/framework/core/ParticleProperties.inl b/corsika/detail/framework/core/ParticleProperties.inl index 76229398e..af41af0fc 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,15 +33,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 HEPEnergyType constexpr get_energy_production_threshold(Code const p) { + return particle::detail::production_thresholds[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 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 ac72089c1..80f024477 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 570f830d6..eb31e373e 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 107f743ff..6054a0d5f 100644 --- a/corsika/detail/modules/proposal/ContinuousProcess.inl +++ b/corsika/detail/modules/proposal/ContinuousProcess.inl @@ -29,7 +29,7 @@ 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_energy_resolution( + 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); @@ -133,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 7cc43e298..9b2e23138 100644 --- a/corsika/detail/modules/proposal/Interaction.inl +++ b/corsika/detail/modules/proposal/Interaction.inl @@ -32,7 +32,7 @@ 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_energy_resolution( + 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 083005058..855d261f9 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,24 @@ 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. + * + * 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. + * * @addtogroup Particles * @{ */ @@ -89,34 +107,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); - } - - 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 + /** + * 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 0630bbf2e..615b840c1 100644 --- a/examples/em_shower.cpp +++ b/examples/em_shower.cpp @@ -102,13 +102,12 @@ int main(int argc, char** argv) { 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 }, + {Code::Electron, 10_MeV}, + {Code::Positron, 10_MeV}, + {Code::Photon, 10_MeV}, }; for (auto [pcode, energy] : energy_resolution) - set_energy_resolution(pcode, energy); - + set_energy_production_threshold(pcode, energy); // setup particle stack, and add primary particle setup::Stack stack; diff --git a/src/framework/core/code_generator.py b/src/framework/core/code_generator.py index 8d5f2788c..9e8c41a37 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,23 +375,25 @@ 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 += " 0 * electronvolt, // {name:s}\n".format( + name=k) string += "};\n\n" - string += "static corsika::units::si::HEPEnergyType threshold_nuclei = 0_eV;\n" + string += "static 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" + # 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 * corsika::units::si::electronvolt, // {name:s}\n".format(name = p['name']) + string += " 1e6 * electronvolt, // {name:s}\n".format( + name=p['name']) string += "};\n\n" # PDG code table @@ -416,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 @@ -449,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" @@ -544,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) @@ -561,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 eaa229222..9c25eef2e 100644 --- a/tests/framework/testParticles.cpp +++ b/tests/framework/testParticles.cpp @@ -90,11 +90,11 @@ TEST_CASE("ParticleProperties", "[Particles]") { SECTION("Energy threshold") { //! 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) == 0_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); + 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); } SECTION("Particle groups: electromagnetic") { diff --git a/tests/modules/testParticleCut.cpp b/tests/modules/testParticleCut.cpp index e12ac4601..3d3c7e666 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") { -- GitLab