diff --git a/corsika/detail/modules/ParticleCut.inl b/corsika/detail/modules/ParticleCut.inl
index 559102cc43278178050323eb84f3fba2eb7b7b30..4afa7b37f339b919c663754db1511a7eec69d7f0 100644
--- a/corsika/detail/modules/ParticleCut.inl
+++ b/corsika/detail/modules/ParticleCut.inl
@@ -15,44 +15,80 @@ namespace corsika {
ParticleCut::ParticleCut(const HEPEnergyType eEleCut, const HEPEnergyType ePhoCut,
const HEPEnergyType eHadCut, const HEPEnergyType eMuCut,
bool inv)
- : electron_energy_cut_(eEleCut)
- , photon_energy_cut_(ePhoCut)
- , had_energy_cut_(eHadCut)
- , mu_energy_cut_(eMuCut)
- , doCutEm_(false)
+ : 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(const HEPEnergyType eHadCut, const HEPEnergyType eMuCut,
bool inv)
- : electron_energy_cut_(0_eV)
- , photon_energy_cut_(0_eV)
- , had_energy_cut_(eHadCut)
- , mu_energy_cut_(eMuCut)
- , doCutEm_(true)
+ : doCutEm_(true)
, 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);
+ CORSIKA_LOG_DEBUG(
+ "setting thresholds: hadrons = {} GeV, "
+ "muons = {} GeV",
+ eHadCut / 1_GeV, eMuCut / 1_GeV);
+ printThresholds();
+ }
ParticleCut::ParticleCut(const HEPEnergyType eCut, bool em, bool inv)
- : electron_energy_cut_(eCut)
- , photon_energy_cut_(eCut)
- , had_energy_cut_(eCut)
- , mu_energy_cut_(eCut)
- , doCutEm_(em)
+ : 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 em, bool inv)
+ : doCutEm_(em)
, doCutInv_(inv)
, energy_(0_GeV)
, em_energy_(0_GeV)
, em_count_(0)
, inv_energy_(0_GeV)
- , inv_count_(0) {}
+ , 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 {
@@ -62,16 +98,9 @@ namespace corsika {
if (pid == Code::Nucleus) {
// calculate energy per nucleon
auto const ElabNuc = energyLab / vP.getNuclearA();
- return (ElabNuc < had_energy_cut_);
- } else if (pid == Code::Gamma) {
- return (energyLab < photon_energy_cut_);
- } else if (pid == Code::Electron || pid == Code::Positron) {
- return (energyLab < electron_energy_cut_);
- } else if (is_muon(pid)) {
- return (energyLab < mu_energy_cut_);
+ return (ElabNuc < get_energy_threshold(pid));
} else {
- // assuming the rest are hadrons
- return (energyLab < had_energy_cut_);
+ return (energyLab < get_energy_threshold(pid));
}
}
@@ -125,6 +154,11 @@ namespace corsika {
return ProcessReturn::Ok;
}
+ void ParticleCut::printThresholds() {
+ for(auto p : get_all_particles())
+ CORSIKA_LOG_DEBUG("energy threshold for particle {} is {} GeV", p, get_energy_threshold(p) / 1_GeV);
+ }
+
void ParticleCut::showResults() {
CORSIKA_LOG_INFO(
" ******************************\n"
diff --git a/corsika/modules/ParticleCut.hpp b/corsika/modules/ParticleCut.hpp
index 11d8508adca3637e4b3d77b061e50acb44181ef1..fe709ee27d470ad385919856c102c09ff6c9e00c 100644
--- a/corsika/modules/ParticleCut.hpp
+++ b/corsika/modules/ParticleCut.hpp
@@ -42,6 +42,10 @@ namespace corsika {
//! discarded altogether.
ParticleCut(const HEPEnergyType eCut, bool em, bool 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 em, bool inv);
+
void doSecondaries(corsika::setup::StackView&);
ProcessReturn doContinuous(corsika::setup::Stack::particle_type& vParticle,
corsika::setup::Trajectory const& vTrajectory);
@@ -50,13 +54,14 @@ namespace corsika {
return meter * std::numeric_limits<double>::infinity();
}
+ void printThresholds();
void showResults();
void reset();
- HEPEnergyType getElectronECut() const { return electron_energy_cut_; }
- HEPEnergyType getPhotonECut() const { return photon_energy_cut_; }
- HEPEnergyType getMuonECut() const { return mu_energy_cut_; }
- HEPEnergyType getHadronECut() const { return had_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_; }
@@ -71,10 +76,6 @@ namespace corsika {
bool isBelowEnergyCut(TParticle const&) const;
private:
- HEPEnergyType electron_energy_cut_;
- HEPEnergyType photon_energy_cut_;
- HEPEnergyType had_energy_cut_;
- HEPEnergyType mu_energy_cut_;
bool doCutEm_;
bool doCutInv_;
HEPEnergyType energy_ = 0 * electronvolt;
diff --git a/tests/modules/testParticleCut.cpp b/tests/modules/testParticleCut.cpp
index 97eebdc8e8177577d2f4ab5814a0205de0ee627d..5d05f04aff38cd17922e9caca0854d57c704a714 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;
@@ -171,6 +171,14 @@ TEST_CASE("ParticleCut", "[processes]") {
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;