diff --git a/corsika/detail/framework/core/ParticleProperties.inl b/corsika/detail/framework/core/ParticleProperties.inl index fc9fbf636241f4c83d364cb9a6da80c458a3ecc8..76229398e7172b911ce12c0fa87644138bb0ba3c 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 77f77e7795a9268e58da251e1b851aad2d0e3d2e..107f743ff5d0360d4167a927ff2ba9797a87695a 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 d8916cbe976e77b8865d7d0edaacc6bfc0135e0b..7cc43e298e985091232910ac5383c170f5996d04 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 dd03aa106c309b6abe3e7f58b69922ebbb8ae875..083005058ea8fc67d7069626e2db752fcb9ae899 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 b4a2c39d414290318a9ceed010a96cf25ddfbf02..0630bbf2e7dd4a91fd65a113e2530391c72836b3 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 e06488c7b4da21840a26986e4514741d17a6c197..8d5f2788c628bbea8ddb1200293eba273ea771b0 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():