IAP GITLAB

Commit b8946967 authored by Maximilian Sackel's avatar Maximilian Sackel

Adding a particle resolution.

Energies below the resolution limit are treated only continuously and no longer stochastically.
parent c24339ce
Pipeline #4072 passed with stages
in 73 minutes and 11 seconds
...@@ -21,6 +21,14 @@ namespace corsika { ...@@ -21,6 +21,14 @@ namespace corsika {
particle::detail::thresholds[static_cast<CodeIntType>(p)] = val; particle::detail::thresholds[static_cast<CodeIntType>(p)] = val;
} }
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 HEPMassType constexpr get_mass(Code const p) { inline HEPMassType constexpr get_mass(Code const p) {
if (p == Code::Nucleus) if (p == Code::Nucleus)
throw std::runtime_error("Cannot GetMass() of particle::Nucleus -> unspecified"); throw std::runtime_error("Cannot GetMass() of particle::Nucleus -> unspecified");
......
...@@ -33,8 +33,8 @@ namespace corsika::proposal { ...@@ -33,8 +33,8 @@ namespace corsika::proposal {
// interpolate the crosssection for given media and energy cut. These may // 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 // take some minutes if you have to build the tables and cannot read the
// from disk // from disk
auto const emCut = get_energy_threshold( auto const emCut = get_energy_resolution(
code); //! energy thresholds globally defined for individual particles code); //! energy resolutions globally defined for individual particles
auto c = p_cross->second(media.at(comp.getHash()), emCut); auto c = p_cross->second(media.at(comp.getHash()), emCut);
// Build displacement integral and scattering object and interpolate them too and // Build displacement integral and scattering object and interpolate them too and
......
...@@ -36,8 +36,8 @@ namespace corsika::proposal { ...@@ -36,8 +36,8 @@ namespace corsika::proposal {
// interpolate the crosssection for given media and energy cut. These may // 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 // take some minutes if you have to build the tables and cannot read the
// from disk // from disk
auto const emCut = get_energy_threshold( auto const emCut = get_energy_resolution(
code); //! energy thresholds globally defined for individual particles code); //! energy resolutions globally defined for individual particles
auto c = p_cross->second(media.at(comp.getHash()), emCut); auto c = p_cross->second(media.at(comp.getHash()), emCut);
......
...@@ -77,6 +77,13 @@ namespace corsika { ...@@ -77,6 +77,13 @@ namespace corsika {
void constexpr set_energy_threshold( void constexpr set_energy_threshold(
Code const, HEPEnergyType const); //!< set energy threshold below which the particle Code const, HEPEnergyType const); //!< set energy threshold below which the particle
//!< is discarded //!< is discarded
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
inline void set_energy_threshold(std::pair<Code const, HEPEnergyType const> p) { inline void set_energy_threshold(std::pair<Code const, HEPEnergyType const> p) {
set_energy_threshold(p.first, p.second); set_energy_threshold(p.first, p.second);
......
...@@ -97,6 +97,15 @@ int main(int argc, char** argv) { ...@@ -97,6 +97,15 @@ int main(int argc, char** argv) {
builder.addLinearLayer(1e9_cm, 112.8_km); builder.addLinearLayer(1e9_cm, 112.8_km);
builder.assemble(env); builder.assemble(env);
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 particle stack, and add primary particle
setup::Stack stack; setup::Stack stack;
stack.clear(); stack.clear();
......
...@@ -339,6 +339,12 @@ def gen_properties(particle_db): ...@@ -339,6 +339,12 @@ def gen_properties(particle_db):
string += " {mass:e} * 1e9 * corsika::units::si::electronvolt, // {name:s}\n".format(mass = p['mass'], name = p['name']) string += " {mass:e} * 1e9 * corsika::units::si::electronvolt, // {name:s}\n".format(mass = p['mass'], name = p['name'])
string += "};\n\n" string += "};\n\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 # PDG code table
string += "static constexpr std::array<PDGCode, size> pdg_codes = {\n" string += "static constexpr std::array<PDGCode, size> pdg_codes = {\n"
for p in particle_db.keys(): for p in particle_db.keys():
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment