diff --git a/corsika/detail/framework/core/ParticleProperties.inl b/corsika/detail/framework/core/ParticleProperties.inl
index 6fc0bea950f383c6e12e537eabcd95cf66cbe4ca..b478b1e17f3f12f850482c6c652c8c884f104656 100644
--- a/corsika/detail/framework/core/ParticleProperties.inl
+++ b/corsika/detail/framework/core/ParticleProperties.inl
@@ -13,6 +13,14 @@
namespace corsika {
+ HEPEnergyType constexpr get_energy_threshold(Code const p) {
+ return particle::detail::thresholds[static_cast<CodeIntType>(p)];
+ }
+
+ void constexpr set_energy_threshold(Code const p, HEPEnergyType const val) {
+ particle::detail::thresholds[static_cast<CodeIntType>(p)] = val;
+ }
+
HEPMassType constexpr get_mass(Code const p) {
if (p == Code::Nucleus)
throw std::runtime_error("Cannot GetMass() of particle::Nucleus -> unspecified");
diff --git a/corsika/detail/framework/geometry/FourVector.inl b/corsika/detail/framework/geometry/FourVector.inl
index 6c0a039e5466ef46d2d3014b7e91a6f0621c922c..135790bfcfa83efaf1a78f218228d0c2557f5853 100644
--- a/corsika/detail/framework/geometry/FourVector.inl
+++ b/corsika/detail/framework/geometry/FourVector.inl
@@ -95,7 +95,7 @@ namespace corsika {
template <typename TTimeType, typename TSpaceVecType>
typename FourVector<TTimeType, TSpaceVecType>::norm_type
- FourVector<TTimeType, TSpaceVecType>::operator*(FourVector const& b) {
+ FourVector<TTimeType, TSpaceVecType>::operator*(FourVector const& b) {
if constexpr (std::is_same<time_type, decltype(std::declval<space_type>() / meter *
second)>::value)
return timeLike_ * b.timeLike_ * constants::cSquared - spaceLike_.norm();
diff --git a/corsika/detail/framework/geometry/QuantityVector.inl b/corsika/detail/framework/geometry/QuantityVector.inl
index 151b35daec64d4909ef3492e80ee3c0c57a64429..994f70689a97d2cc6f183f7d6e03e9558f08b7ba 100644
--- a/corsika/detail/framework/geometry/QuantityVector.inl
+++ b/corsika/detail/framework/geometry/QuantityVector.inl
@@ -17,8 +17,8 @@
namespace corsika {
template <typename TDimension>
- inline typename QuantityVector<TDimension>::quantity_type QuantityVector<TDimension>::
- operator[](size_t const index) const {
+ inline typename QuantityVector<TDimension>::quantity_type
+ QuantityVector<TDimension>::operator[](size_t const index) const {
return quantity_type(phys::units::detail::magnitude_tag, eigenVector_[index]);
}
diff --git a/corsika/detail/media/WeightProvider.inl b/corsika/detail/media/WeightProvider.inl
index af9849aa3ca75b9bdaeefc661d56232f1ab895f4..c075b236350127c9978fb7c08d87ed035292e5a8 100644
--- a/corsika/detail/media/WeightProvider.inl
+++ b/corsika/detail/media/WeightProvider.inl
@@ -20,7 +20,7 @@ namespace corsika {
template <class AConstIterator, class BConstIterator>
typename WeightProviderIterator<AConstIterator, BConstIterator>::value_type
- WeightProviderIterator<AConstIterator, BConstIterator>::operator*() const {
+ WeightProviderIterator<AConstIterator, BConstIterator>::operator*() const {
return ((*aIter_) * (*bIter_)).magnitude();
}
diff --git a/corsika/detail/modules/ParticleCut.inl b/corsika/detail/modules/ParticleCut.inl
index da19bce2727550a09a3ccede6b08a8ead7ca51ec..7a22b4d719beddef778acab9861b65ba7578217f 100644
--- a/corsika/detail/modules/ParticleCut.inl
+++ b/corsika/detail/modules/ParticleCut.inl
@@ -12,68 +12,114 @@
namespace corsika {
- ParticleCut::ParticleCut(const HEPEnergyType eCut, bool em, bool inv)
- : energy_cut_(eCut)
- , doCutEm_(em)
+ ParticleCut::ParticleCut(HEPEnergyType const eEleCut, HEPEnergyType const ePhoCut,
+ HEPEnergyType const eHadCut, HEPEnergyType const eMuCut,
+ bool const inv)
+ : doCutEm_(false)
, doCutInv_(inv)
, energy_(0_GeV)
, em_energy_(0_GeV)
, em_count_(0)
, inv_energy_(0_GeV)
- , inv_count_(0) {}
+ , inv_count_(0) {
+ for (auto p : get_all_particles())
+ if (is_hadron(p))
+ set_energy_threshold(p, eHadCut);
+ else if (is_muon(p))
+ set_energy_threshold(p, eMuCut);
+ else if (p == Code::Electron || p == Code::Positron)
+ set_energy_threshold(p, eEleCut);
+ else if (p == Code::Gamma)
+ set_energy_threshold(p, ePhoCut);
+ else if (p == Code::Nucleus)
+ // nuclei have same threshold as hadrons on the nucleon level.
+ set_energy_threshold(p, eHadCut);
+ CORSIKA_LOG_DEBUG(
+ "setting thresholds: electrons = {} GeV, photons = {} GeV, hadrons = {} GeV, "
+ "muons = {} GeV",
+ eEleCut / 1_GeV, ePhoCut / 1_GeV, eHadCut / 1_GeV, eMuCut / 1_GeV);
+ printThresholds();
+ }
+
+ ParticleCut::ParticleCut(HEPEnergyType const eHadCut, HEPEnergyType const eMuCut,
+ bool const inv)
+ : doCutEm_(true)
+ , doCutInv_(inv)
+ , energy_(0_GeV)
+ , em_energy_(0_GeV)
+ , em_count_(0)
+ , inv_energy_(0_GeV)
+ , inv_count_(0) {
+
+ for (auto p : get_all_particles())
+ if (is_hadron(p))
+ set_energy_threshold(p, eHadCut);
+ else if (is_muon(p))
+ set_energy_threshold(p, eMuCut);
+ CORSIKA_LOG_DEBUG(
+ "setting thresholds: hadrons = {} GeV, "
+ "muons = {} GeV",
+ eHadCut / 1_GeV, eMuCut / 1_GeV);
+ printThresholds();
+ }
+
+ ParticleCut::ParticleCut(HEPEnergyType const eCut, bool const em, bool const inv)
+ : doCutEm_(em)
+ , doCutInv_(inv)
+ , energy_(0_GeV)
+ , em_energy_(0_GeV)
+ , em_count_(0)
+ , inv_energy_(0_GeV)
+ , inv_count_(0) {
+ for (auto p : get_all_particles()) set_energy_threshold(p, eCut);
+ CORSIKA_LOG_DEBUG("setting threshold for all particles to {} GeV", eCut / 1_GeV);
+ printThresholds();
+ }
+
+ ParticleCut::ParticleCut(
+ std::unordered_map<Code const, HEPEnergyType const> const& eCuts, bool const em,
+ bool const inv)
+ : doCutEm_(em)
+ , doCutInv_(inv)
+ , energy_(0_GeV)
+ , em_energy_(0_GeV)
+ , em_count_(0)
+ , inv_energy_(0_GeV)
+ , inv_count_(0) {
+ set_energy_thresholds(eCuts);
+ CORSIKA_LOG_DEBUG("setting threshold particles individually");
+ printThresholds();
+ }
template <typename TParticle>
bool ParticleCut::isBelowEnergyCut(TParticle const& vP) const {
auto const energyLab = vP.getEnergy();
+ auto const pid = vP.getPID();
// nuclei
- if (vP.getPID() == Code::Nucleus) {
+ if (pid == Code::Nucleus) {
// calculate energy per nucleon
auto const ElabNuc = energyLab / vP.getNuclearA();
- return (ElabNuc < energy_cut_);
+ return (ElabNuc < get_energy_threshold(pid));
} else {
- return (energyLab < energy_cut_);
- }
- }
-
- bool ParticleCut::isEmParticle(Code vCode) const {
- // FOR NOW: switch
- switch (vCode) {
- case Code::Gamma:
- case Code::Electron:
- case Code::Positron:
- return true;
- default:
- return false;
+ return (energyLab < get_energy_threshold(pid));
}
}
- bool ParticleCut::isInvisible(Code vCode) const {
- switch (vCode) {
- case Code::NuE:
- case Code::NuEBar:
- case Code::NuMu:
- case Code::NuMuBar:
- return true;
-
- default:
- return false;
- }
- }
+ bool ParticleCut::isInvisible(Code const& vCode) const { return is_neutrino(vCode); }
template <typename TParticle>
- bool ParticleCut::checkCutParticle(const TParticle& particle) {
+ bool ParticleCut::checkCutParticle(TParticle const& particle) {
- const Code pid = particle.getPID();
+ Code const pid = particle.getPID();
HEPEnergyType energy = particle.getEnergy();
- CORSIKA_LOG_DEBUG(fmt::format("ParticleCut: checking {}, E= {} GeV, EcutTot={} GeV",
- pid, energy / 1_GeV,
- (em_energy_ + inv_energy_ + energy_) / 1_GeV));
- if (doCutEm_ && isEmParticle(pid)) {
+ CORSIKA_LOG_DEBUG("ParticleCut: checking {}, E= {} GeV, EcutTot={} GeV", pid,
+ energy / 1_GeV, (em_energy_ + inv_energy_ + energy_) / 1_GeV);
+ if (doCutEm_ && is_em(pid)) {
CORSIKA_LOG_DEBUG("removing em. particle...");
em_energy_ += energy;
em_count_ += 1;
return true;
- } else if (doCutInv_ && isInvisible(pid)) {
+ } else if (doCutInv_ && is_neutrino(pid)) {
CORSIKA_LOG_DEBUG("removing inv. particle...");
inv_energy_ += energy;
inv_count_ += 1;
@@ -110,6 +156,13 @@ namespace corsika {
return ProcessReturn::Ok;
}
+ void ParticleCut::printThresholds() {
+ for (auto p : get_all_particles()) {
+ auto const Eth = get_energy_threshold(p);
+ CORSIKA_LOG_INFO("energy threshold for particle {} is {} GeV", p, Eth / 1_GeV);
+ }
+ }
+
void ParticleCut::showResults() {
CORSIKA_LOG_INFO(
" ******************************\n"
diff --git a/corsika/detail/modules/energy_loss/BetheBlochPDG.inl b/corsika/detail/modules/energy_loss/BetheBlochPDG.inl
index 776531540cce808be3fb017c40865d6d1bddc373..1f3d0e9303cb60f40fb588ecad70034c634dd1c8 100644
--- a/corsika/detail/modules/energy_loss/BetheBlochPDG.inl
+++ b/corsika/detail/modules/energy_loss/BetheBlochPDG.inl
@@ -26,11 +26,10 @@ namespace corsika {
return sqrt((Elab - m) * (Elab + m));
};
- BetheBlochPDG::BetheBlochPDG(ShowerAxis const& shower_axis, HEPEnergyType emCut)
+ BetheBlochPDG::BetheBlochPDG(ShowerAxis const& shower_axis)
: dX_(10_g / square(1_cm)) // profile binning
, dX_threshold_(0.0001_g / square(1_cm))
, shower_axis_(shower_axis)
- , emCut_(emCut)
, profile_(int(shower_axis.getMaximumX() / dX_) + 1) {}
HEPEnergyType BetheBlochPDG::getBetheBloch(setup::Stack::particle_type const& p,
@@ -170,8 +169,9 @@ namespace corsika {
auto const dEdX = -getTotalEnergyLoss(vParticle, dX) / dX; // dE > 0
//~ auto const Ekin = vParticle.getEnergy() - vParticle.getMass();
- // in any case: never go below 0.99*emCut_ This needs to be
- // slightly smaller than emCut_ since, either this Step is limited
+ auto const emCut = get_energy_threshold(vParticle.getPID());
+ // in any case: never go below 0.99*emCut This needs to be
+ // slightly smaller than emCut since, either this Step is limited
// by energy_lim, then the particle is stopped in a very short
// range (before doing anythin else) and is then removed
// instantly. The exact position where it reaches emCut is not
@@ -179,7 +179,7 @@ namespace corsika {
// afterwards.
//
const auto energy = vParticle.getEnergy();
- auto energy_lim = std::max(0.9 * energy, 0.99 * emCut_);
+ auto energy_lim = std::max(0.9 * energy, 0.99 * emCut);
auto const maxGrammage = (energy - energy_lim) / dEdX;
diff --git a/corsika/detail/modules/proposal/ContinuousProcess.inl b/corsika/detail/modules/proposal/ContinuousProcess.inl
index 16d419b0be200bb0d65b8748770818491d5e53b2..8bfbd705b8bb0dfc19d1e7efe79252aa5a6a69be 100644
--- a/corsika/detail/modules/proposal/ContinuousProcess.inl
+++ b/corsika/detail/modules/proposal/ContinuousProcess.inl
@@ -32,7 +32,9 @@ 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 c = p_cross->second(media.at(comp.getHash()), emCut_);
+ auto const emCut = get_energy_threshold(
+ code); //! energy thresholds globally defined for individual particles
+ auto c = p_cross->second(media.at(comp.getHash()), emCut);
// Build displacement integral and scattering object and interpolate them too and
// saved in the calc map by a key build out of a hash of composed of the component and
@@ -45,9 +47,8 @@ namespace corsika::proposal {
}
template <>
- ContinuousProcess::ContinuousProcess(setup::Environment const& _env,
- HEPEnergyType _emCut)
- : ProposalProcessBase(_env, _emCut) {}
+ ContinuousProcess::ContinuousProcess(setup::Environment const& _env)
+ : ProposalProcessBase(_env) {}
template <>
void ContinuousProcess::scatter(setup::Stack::particle_type& vP,
@@ -115,21 +116,24 @@ namespace corsika::proposal {
template <>
LengthType ContinuousProcess::getMaxStepLength(setup::Stack::particle_type const& vP,
setup::Trajectory const& vT) {
-
- if (!canInteract(vP.getPID())) return meter * std::numeric_limits<double>::infinity();
+ auto const code = vP.getPID();
+ if (!canInteract(code)) return meter * std::numeric_limits<double>::infinity();
// Limit the step size of a conitnuous loss. The maximal continuous loss seems to be a
// hyper parameter which must be adjusted.
//
- // in any case: never go below 0.99*emCut_ This needs to be
- // slightly smaller than emCut_ since, either this Step is limited
+ auto const emCut = get_energy_threshold(
+ code); //! energy thresholds globally defined for individual particles
+
+ // in any case: never go below 0.99*emCut This needs to be
+ // slightly smaller than emCut since, either this Step is limited
// by energy_lim, then the particle is stopped in a very short
// range (before doing anythin else) and is then removed
// instantly. The exact position where it reaches emCut is not
// important, the important fact is that its E_kin is zero
// afterwards.
//
- auto energy_lim = std::max(0.9 * vP.getEnergy(), 0.99 * emCut_);
+ auto energy_lim = std::max(0.9 * vP.getEnergy(), 0.99 * emCut);
// solving the track integral for giving energy lim
auto c = getCalculator(vP, calc);
diff --git a/corsika/detail/modules/proposal/Interaction.inl b/corsika/detail/modules/proposal/Interaction.inl
index 853302967fcc0bd460f6181f593254c9335b10bb..4874bda1f7fb7f6ca89d898cf7f1e01e8a51b791 100644
--- a/corsika/detail/modules/proposal/Interaction.inl
+++ b/corsika/detail/modules/proposal/Interaction.inl
@@ -24,8 +24,8 @@
namespace corsika::proposal {
template <>
- Interaction::Interaction(setup::Environment const& _env, HEPEnergyType _emCut)
- : ProposalProcessBase(_env, _emCut) {}
+ Interaction::Interaction(setup::Environment const& _env)
+ : ProposalProcessBase(_env) {}
void Interaction::buildCalculator(Code code, NuclearComposition const& comp) {
// search crosssection builder for given particle
@@ -36,7 +36,10 @@ 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 c = p_cross->second(media.at(comp.getHash()), emCut_);
+ auto const emCut = get_energy_threshold(
+ code); //! energy thresholds globally defined for individual particles
+
+ auto c = p_cross->second(media.at(comp.getHash()), emCut);
// Look which interactions take place and build the corresponding
// interaction and secondarie builder. The interaction integral will
diff --git a/corsika/detail/modules/proposal/ProposalProcessBase.inl b/corsika/detail/modules/proposal/ProposalProcessBase.inl
index b71b0fa9b1f9aff2784796588fb4ea7017a17d1f..21d0cce576f306055d745a8d50a6a2baa3f21466 100644
--- a/corsika/detail/modules/proposal/ProposalProcessBase.inl
+++ b/corsika/detail/modules/proposal/ProposalProcessBase.inl
@@ -30,10 +30,8 @@ namespace corsika::proposal {
return false;
}
- ProposalProcessBase::ProposalProcessBase(setup::Environment const& _env,
- HEPEnergyType _emCut)
- : emCut_(_emCut)
- , RNG_(RNGManager::getInstance().getRandomStream("proposal")) {
+ ProposalProcessBase::ProposalProcessBase(setup::Environment const& _env)
+ : RNG_(RNGManager::getInstance().getRandomStream("proposal")) {
_env.getUniverse()->walk([&](auto& vtn) {
if (vtn.hasModelProperties()) {
const auto& prop = vtn.getModelProperties();
diff --git a/corsika/framework/core/ParticleProperties.hpp b/corsika/framework/core/ParticleProperties.hpp
index 2deb03a865c6d32e37f4c274c83a604f72fecb04..85227425c42bf584ffcde2cea20c75d11570709e 100644
--- a/corsika/framework/core/ParticleProperties.hpp
+++ b/corsika/framework/core/ParticleProperties.hpp
@@ -25,6 +25,7 @@
#include <iosfwd>
#include <string_view>
#include <type_traits>
+#include <unordered_map>
#include <corsika/framework/core/PhysicalUnits.hpp>
@@ -50,6 +51,20 @@ namespace corsika {
int16_t constexpr get_charge_number(Code); //!< electric charge in units of e
ElectricChargeType constexpr get_charge(Code); //!< electric charge
HEPMassType constexpr get_mass(Code); //!< mass
+ HEPEnergyType constexpr get_energy_threshold(
+ Code const); //!< get energy threshold below which the particle is discarded, by
+ //!< default set to particle mass
+ void constexpr set_energy_threshold(
+ Code const, HEPEnergyType const); //!< set energy threshold below which the particle
+ //!< is discarded
+
+ inline void set_energy_threshold(std::pair<Code const, HEPEnergyType const> p) {
+ set_energy_threshold(p.first, p.second);
+ }
+ inline void set_energy_thresholds(
+ std::unordered_map<Code const, HEPEnergyType const> const& eCuts) {
+ for (auto v : eCuts) set_energy_threshold(v);
+ }
//! Particle code according to PDG, "Monte Carlo Particle Numbering Scheme"
PDGCode constexpr get_PDG(Code);
diff --git a/corsika/modules/ParticleCut.hpp b/corsika/modules/ParticleCut.hpp
index 7abc7d6be02517cf78c8bc19749387cbf7055a23..ca7f87013d59cf2cbab4194c7e0fa02b31101c42 100644
--- a/corsika/modules/ParticleCut.hpp
+++ b/corsika/modules/ParticleCut.hpp
@@ -8,6 +8,8 @@
#pragma once
+#include <unordered_map>
+
#include <corsika/framework/core/ParticleProperties.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/process/SecondariesProcess.hpp>
@@ -17,12 +19,35 @@
#include <corsika/setup/SetupTrajectory.hpp>
namespace corsika {
-
+ /**
+ simple ParticleCut process. Goes through the secondaries of an interaction and
+ removes particles according to their energy. Particles with a time delay of more than
+ 10ms are removed as well. Invisible particles (neutrinos) can be removed if selected.
+ **/
class ParticleCut : public SecondariesProcess<ParticleCut>,
public ContinuousProcess<ParticleCut> {
public:
- ParticleCut(const HEPEnergyType eCut, bool em, bool inv);
+ /**
+ * particle cut with energy thresholds for electrons, photons,
+ * hadrons (including nuclei with energy per nucleon) and muons
+ * invisible particles (neutrinos) can be cut or not
+ **/
+ ParticleCut(HEPEnergyType const eEleCut, HEPEnergyType const ePhoCut,
+ HEPEnergyType const eHadCut, HEPEnergyType const eMuCut, bool const inv);
+
+ //! simple cut. hadrons and muons are cut by threshold. EM particles are all
+ //! discarded.
+ ParticleCut(HEPEnergyType const eHadCut, HEPEnergyType const euCut, bool const inv);
+
+ //! simplest cut. all particles have same threshold. EM particles can be set to be
+ //! discarded altogether.
+ ParticleCut(HEPEnergyType const eCut, bool const em, bool const inv);
+
+ //! threshold for specific particles redefined. EM and invisible particles can be set
+ //! to be discarded altogether.
+ ParticleCut(std::unordered_map<Code const, HEPEnergyType const> const& eCuts,
+ bool const em, bool const inv);
void doSecondaries(corsika::setup::StackView&);
ProcessReturn doContinuous(corsika::setup::Stack::particle_type& vParticle,
@@ -32,10 +57,14 @@ namespace corsika {
return meter * std::numeric_limits<double>::infinity();
}
+ void printThresholds();
void showResults();
void reset();
- HEPEnergyType getECut() const { return energy_cut_; }
+ HEPEnergyType getElectronECut() const { return get_energy_threshold(Code::Electron); }
+ HEPEnergyType getPhotonECut() const { return get_energy_threshold(Code::Gamma); }
+ HEPEnergyType getMuonECut() const { return get_energy_threshold(Code::MuPlus); }
+ HEPEnergyType getHadronECut() const { return get_energy_threshold(Code::Proton); }
HEPEnergyType getInvEnergy() const { return inv_energy_; }
HEPEnergyType getCutEnergy() const { return energy_; }
HEPEnergyType getEmEnergy() const { return em_energy_; }
@@ -44,15 +73,15 @@ namespace corsika {
private:
template <typename TParticle>
- bool checkCutParticle(const TParticle& p);
+ bool checkCutParticle(TParticle const& p);
template <typename TParticle>
bool isBelowEnergyCut(TParticle const&) const;
- bool isEmParticle(Code) const;
- bool isInvisible(Code) const;
+
+ //! defines which particles are invisible, by default only neutrinos
+ bool isInvisible(Code const&) const;
private:
- HEPEnergyType energy_cut_;
bool doCutEm_;
bool doCutInv_;
HEPEnergyType energy_ = 0 * electronvolt;
diff --git a/corsika/modules/energy_loss/BetheBlochPDG.hpp b/corsika/modules/energy_loss/BetheBlochPDG.hpp
index 3e64567b974c9dfdb71b6f67c99e7877dbd6d39d..f8e8759376a3833bd3f7df438b2dbe4132c5e3c8 100644
--- a/corsika/modules/energy_loss/BetheBlochPDG.hpp
+++ b/corsika/modules/energy_loss/BetheBlochPDG.hpp
@@ -41,11 +41,13 @@ namespace corsika {
using MeVgcm2 = decltype(1e6 * electronvolt / gram * square(1e-2 * meter));
public:
- BetheBlochPDG(ShowerAxis const& showerAxis, HEPEnergyType emCut);
+ BetheBlochPDG(ShowerAxis const& showerAxis);
ProcessReturn doContinuous(setup::Stack::particle_type&, setup::Trajectory const&);
LengthType getMaxStepLength(setup::Stack::particle_type const&,
- setup::Trajectory const&) const;
+ setup::Trajectory const&)
+ const; //! limited by the energy threshold! By default the limit is the particle
+ //! rest mass, i.e. kinetic energy is zero
static HEPEnergyType getBetheBloch(setup::Stack::particle_type const&,
const GrammageType);
static HEPEnergyType getRadiationLosses(setup::Stack::particle_type const&,
@@ -66,7 +68,6 @@ namespace corsika {
GrammageType const dX_ = 10_g / square(1_cm); // profile binning
GrammageType const dX_threshold_ = 0.0001_g / square(1_cm);
ShowerAxis const& shower_axis_;
- HEPEnergyType emCut_;
HEPEnergyType energy_lost_ = HEPEnergyType::zero();
std::vector<HEPEnergyType> profile_; // longitudinal profile
};
diff --git a/corsika/modules/proposal/ContinuousProcess.hpp b/corsika/modules/proposal/ContinuousProcess.hpp
index 3bf53b145246392848f6ba09cbd293ac722c0a80..3b8976f3ac6bfe7a6d5251740782e4ca8fec89bc 100644
--- a/corsika/modules/proposal/ContinuousProcess.hpp
+++ b/corsika/modules/proposal/ContinuousProcess.hpp
@@ -52,7 +52,7 @@ namespace corsika::proposal {
//! compositions and stochastic description limited by the particle cut.
//!
template <typename TEnvironment>
- ContinuousProcess(TEnvironment const&, HEPEnergyType _emCut);
+ ContinuousProcess(TEnvironment const&);
//!
//! Multiple Scattering of the lepton. Stochastic deflection is not yet taken into
diff --git a/corsika/modules/proposal/Interaction.hpp b/corsika/modules/proposal/Interaction.hpp
index 69f27cca42b728cd52a77883cd377c6c3ccafc36..c6a389b1e2dcfcbb5308dcd53b9dd7b33589ef00 100644
--- a/corsika/modules/proposal/Interaction.hpp
+++ b/corsika/modules/proposal/Interaction.hpp
@@ -48,7 +48,7 @@ namespace corsika::proposal {
//! compositions and stochastic description limited by the particle cut.
//!
template <typename TEnvironment>
- Interaction(TEnvironment const& env, HEPEnergyType emCut);
+ Interaction(TEnvironment const& env);
//!
//! Calculate the rates for the different targets and interactions. Sample a
diff --git a/corsika/modules/proposal/ProposalProcessBase.hpp b/corsika/modules/proposal/ProposalProcessBase.hpp
index 89037f9fead31d255b76e40180a265aaa00d3cb7..cedad07efa87d7266764ebd61d9b808983e68de5 100644
--- a/corsika/modules/proposal/ProposalProcessBase.hpp
+++ b/corsika/modules/proposal/ProposalProcessBase.hpp
@@ -44,7 +44,10 @@ namespace corsika::proposal {
//!
template <typename T>
static auto cross_builder =
- [](PROPOSAL::Medium& m, corsika::units::si::HEPEnergyType emCut) {
+ [](PROPOSAL::Medium& m,
+ corsika::units::si::HEPEnergyType
+ emCut) { //!< Stochastic losses smaller than the given cut
+ //!< will be handeled continuously.
using namespace corsika::units::si;
auto p_cut =
std::make_shared<const PROPOSAL::EnergyCutSettings>(emCut / 1_MeV, 1, true);
@@ -73,8 +76,6 @@ namespace corsika::proposal {
//!
class ProposalProcessBase {
protected:
- HEPEnergyType emCut_; //!< Stochastic losses smaller than the given cut
- //!< will be handeled continuously.
RNGManager::prng_type RNG_; //!< random number generator used by proposal
std::unordered_map<std::size_t, PROPOSAL::Medium>
@@ -85,7 +86,7 @@ namespace corsika::proposal {
//! Store cut and nuclear composition of the whole universe in media which are
//! required for creating crosssections by proposal.
//!
- ProposalProcessBase(corsika::setup::Environment const& _env, HEPEnergyType _emCut);
+ ProposalProcessBase(corsika::setup::Environment const& _env);
//!
//! Checks if a particle can be processed by proposal
diff --git a/examples/cascade_example.cpp b/examples/cascade_example.cpp
index d14843c51730d8c05dd2e8d46c7374ba5e4224b7..91bc3f7cc99640a10e067da3360eee6e50911020 100644
--- a/examples/cascade_example.cpp
+++ b/examples/cascade_example.cpp
@@ -142,7 +142,7 @@ int main() {
ParticleCut cut(80_GeV, true, true);
TrackWriter trackWriter("tracks.dat");
- BetheBlochPDG eLoss{showerAxis, cut.getECut()};
+ BetheBlochPDG eLoss{showerAxis};
// assemble all processes into an ordered process list
auto sequence =
diff --git a/examples/cascade_proton_example.cpp b/examples/cascade_proton_example.cpp
index 003dc85584ebbe6911348b74109471fda5b3c096..24ed9e77d3be401229568a818511c9b2f2735563 100644
--- a/examples/cascade_proton_example.cpp
+++ b/examples/cascade_proton_example.cpp
@@ -132,7 +132,7 @@ int main() {
TrackWriter trackWriter("tracks.dat");
ShowerAxis const showerAxis{injectionPos, Vector{rootCS, 0_m, 0_m, -100_km}, env};
- BetheBlochPDG eLoss{showerAxis, cut.getECut()};
+ BetheBlochPDG eLoss{showerAxis};
// assemble all processes into an ordered process list
// auto sequence = sibyll << decay << hadronicElastic << cut << trackWriter;
diff --git a/examples/em_shower.cpp b/examples/em_shower.cpp
index edb377953427cb89120828193319b41ce152d8b4..2afddc84ee9744909d4378eb3cf2456267655308 100644
--- a/examples/em_shower.cpp
+++ b/examples/em_shower.cpp
@@ -142,9 +142,9 @@ int main(int argc, char** argv) {
// setup processes, decays and interactions
// PROPOSAL processs proposal{...};
- ParticleCut cut(10_GeV, false, true);
- corsika::proposal::Interaction proposal(env, cut.getECut());
- corsika::proposal::ContinuousProcess em_continuous(env, cut.getECut());
+ ParticleCut cut(10_GeV, 10_GeV, 100_PeV, 100_PeV, true);
+ corsika::proposal::Interaction proposal(env);
+ corsika::proposal::ContinuousProcess em_continuous(env);
InteractionCounter proposalCounted(proposal);
TrackWriter trackWriter("tracks.dat");
diff --git a/examples/hybrid_MC.cpp b/examples/hybrid_MC.cpp
index 71d8cf82875f39c85a7f6b11dd39f7ce2faac5bf..83ca7bd57098a6af751f430077dffc7414fa31e1 100644
--- a/examples/hybrid_MC.cpp
+++ b/examples/hybrid_MC.cpp
@@ -207,7 +207,7 @@ int main(int argc, char** argv) {
decaySibyll.printDecayConfig();
ParticleCut cut{60_GeV, false, true};
- BetheBlochPDG eLoss{showerAxis, cut.getECut()};
+ BetheBlochPDG eLoss{showerAxis};
CONEXhybrid conex_model(center, showerAxis, t, injectionHeight, E0,
get_PDG(Code::Proton));
diff --git a/examples/stopping_power.cpp b/examples/stopping_power.cpp
index ef8252501703205734ed361e4667df9fabdcbbed..d419efbb86b9cbfb3bc5a2b8db3cc8adb3a6e2bb 100644
--- a/examples/stopping_power.cpp
+++ b/examples/stopping_power.cpp
@@ -49,7 +49,7 @@ int main() {
112.8_km); // this is the CORSIKA 7 start of atmosphere/universe
ShowerAxis showerAxis{injectionPos, Vector<length_d>{rootCS, 0_m, 0_m, 1_m}, env};
- BetheBlochPDG eLoss{showerAxis, 300_MeV};
+ BetheBlochPDG eLoss{showerAxis};
setup::Stack stack;
diff --git a/examples/vertical_EAS.cpp b/examples/vertical_EAS.cpp
index 28488a1bdb53a8cc9a6b8f234152acf289919dfa..f012440e82deb19c3ab2b5220622c0255da9dbc5 100644
--- a/examples/vertical_EAS.cpp
+++ b/examples/vertical_EAS.cpp
@@ -221,9 +221,9 @@ int main(int argc, char** argv) {
decaySibyll.printDecayConfig();
- ParticleCut cut{60_GeV, false, true};
- corsika::proposal::Interaction proposal(env, cut.getECut());
- corsika::proposal::ContinuousProcess em_continuous(env, cut.getECut());
+ ParticleCut cut{60_GeV, 60_GeV, 60_GeV, 60_GeV, true};
+ corsika::proposal::Interaction proposal(env);
+ corsika::proposal::ContinuousProcess em_continuous(env);
InteractionCounter proposalCounted(proposal);
OnShellCheck reset_particle_mass(1.e-3, 1.e-1, false);
diff --git a/modules/conex b/modules/conex
index 73e3d1fa850d2d3602a15b3b948ac1789fbc968d..820f042b6a055276d465437c74160ef7c199b646 160000
--- a/modules/conex
+++ b/modules/conex
@@ -1 +1 @@
-Subproject commit 73e3d1fa850d2d3602a15b3b948ac1789fbc968d
+Subproject commit 820f042b6a055276d465437c74160ef7c199b646
diff --git a/modules/data b/modules/data
index afe83dc19c1464deb176f38c3ddab80a64e081f4..8b76a9ca2599cd0ce1f204b17362eb06bbcf5277 160000
--- a/modules/data
+++ b/modules/data
@@ -1 +1 @@
-Subproject commit afe83dc19c1464deb176f38c3ddab80a64e081f4
+Subproject commit 8b76a9ca2599cd0ce1f204b17362eb06bbcf5277
diff --git a/src/framework/core/pdxml_reader.py b/src/framework/core/pdxml_reader.py
index fcd0b134187647d81bf747c3634e4dce2eb76a57..2520cdfe4adb3f75bf316043425ca47e3e3742d7 100755
--- a/src/framework/core/pdxml_reader.py
+++ b/src/framework/core/pdxml_reader.py
@@ -330,7 +330,13 @@ def gen_properties(particle_db):
for p in particle_db.values():
string += " {mass:e} * 1e9 * corsika::units::si::electronvolt, // {name:s}\n".format(mass = p['mass'], name = p['name'])
string += "};\n\n"
-
+
+ # particle threshold table, initially set to the particle mass
+ string += "static std::array<corsika::units::si::HEPEnergyType, size> thresholds = {\n"
+ for p in particle_db.values():
+ string += " {mass:e} * 1e9 * corsika::units::si::electronvolt, // {name:s}\n".format(mass = p['mass'], 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():
diff --git a/tests/framework/testParticles.cpp b/tests/framework/testParticles.cpp
index dd999ecf8c6be288546fd7b32ec5770b414c39b4..b74d836bc47d5715ec4f2504b30fb056af545af4 100644
--- a/tests/framework/testParticles.cpp
+++ b/tests/framework/testParticles.cpp
@@ -80,6 +80,15 @@ TEST_CASE("ParticleProperties", "[Particles]") {
(Approx(2.1970332555864364e-06).epsilon(1e-5)));
}
+ SECTION("Energy threshold") {
+ //! by default energy thresholds are set to particle mass
+ CHECK(get_energy_threshold(Electron::code) / Electron::mass == Approx(1));
+
+ set_energy_threshold(Electron::code, 10_GeV);
+ CHECK_FALSE(get_energy_threshold(Code::Electron) == 1_GeV);
+ CHECK(get_energy_threshold(Code::Electron) == 10_GeV);
+ }
+
SECTION("Particle groups: electromagnetic") {
CHECK(is_em(Code::Gamma));
CHECK(is_em(Code::Electron));
diff --git a/tests/modules/testParticleCut.cpp b/tests/modules/testParticleCut.cpp
index a8894230196e5041482cd53d267d1786a5a05a21..0dede160e563c6abab6a9a3fb2b1f0318f3dc1c7 100644
--- a/tests/modules/testParticleCut.cpp
+++ b/tests/modules/testParticleCut.cpp
@@ -25,7 +25,7 @@ using namespace corsika;
TEST_CASE("ParticleCut", "[processes]") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
+ corsika_logger->set_pattern("[%n:%^%-8l%$] %v");
feenableexcept(FE_INVALID);
using EnvType = setup::Environment;
@@ -51,7 +51,7 @@ TEST_CASE("ParticleCut", "[processes]") {
SECTION("cut on particle type: inv") {
ParticleCut cut(20_GeV, false, true);
- CHECK(cut.getECut() == 20_GeV);
+ CHECK(cut.getHadronECut() == 20_GeV);
// add primary particle to stack
auto particle = stack.addParticle(std::make_tuple(
@@ -135,6 +135,50 @@ TEST_CASE("ParticleCut", "[processes]") {
CHECK(view.getSize() == 13);
}
+ SECTION("cut low energy: electrons, photons, hadrons and muons") {
+ ParticleCut cut(5_MeV, 5_MeV, 5_GeV, 5_GeV, true);
+
+ // add primary particle to stack
+ auto particle = stack.addParticle(
+ std::make_tuple(Code::Proton, Eabove,
+ MomentumVector(rootCS, {0_GeV, 0_GeV, 0_GeV}), point0, 0_ns));
+ // view on secondary particles
+ setup::StackView view(particle);
+ // ref. to primary particle through the secondary view.
+ // only this way the secondary view is populated
+ auto projectile = view.getProjectile();
+ // add secondaries
+ projectile.addSecondary(std::make_tuple(
+ Code::Gamma, 3_MeV, MomentumVector(rootCS, {0_GeV, 0_GeV, 0_GeV}), point0, 0_ns));
+ projectile.addSecondary(std::make_tuple(Code::Electron, 3_MeV,
+ MomentumVector(rootCS, {0_GeV, 0_GeV, 0_GeV}),
+ point0, 0_ns));
+ projectile.addSecondary(std::make_tuple(Code::PiPlus, 4_GeV,
+ MomentumVector(rootCS, {0_GeV, 0_GeV, 0_GeV}),
+ point0, 0_ns));
+ unsigned short A = 18;
+ unsigned short Z = 8;
+ projectile.addSecondary(std::make_tuple(Code::Nucleus, 4_GeV * A,
+ MomentumVector(rootCS, {0_GeV, 0_GeV, 0_GeV}),
+ point0, 0_ns, A, Z));
+ projectile.addSecondary(std::make_tuple(Code::Nucleus, 6_GeV * A,
+ MomentumVector(rootCS, {0_GeV, 0_GeV, 0_GeV}),
+ point0, 0_ns, A, Z));
+
+ cut.doSecondaries(view);
+
+ CHECK(view.getEntries() == 1);
+ CHECK(view.getSize() == 5);
+ }
+
+ SECTION("cut low energy: reset thresholds of arbitrary set of particles") {
+ ParticleCut cut({{Code::Electron, 5_MeV}, {Code::Positron, 50_MeV}}, false, true);
+ CHECK(get_energy_threshold(Code::Electron) != get_energy_threshold(Code::Positron));
+ CHECK_FALSE(get_energy_threshold(Code::Electron) == Electron::mass);
+ // test default values still correct
+ CHECK(get_energy_threshold(Code::Proton) == 5_GeV);
+ }
+
SECTION("cut on time") {
ParticleCut cut(20_GeV, false, false);
const TimeType too_late = 1_s;
@@ -170,7 +214,7 @@ TEST_CASE("ParticleCut", "[processes]") {
SECTION("cut on DoContinous, just invisibles") {
ParticleCut cut(20_GeV, false, true);
- CHECK(cut.getECut() == 20_GeV);
+ CHECK(cut.getHadronECut() == 20_GeV);
// add particles, all with energies above the threshold
// only cut is by species