From 4abfd23d5b0b45ff18549d4551f37b77db6fb4ce Mon Sep 17 00:00:00 2001
From: Felix Riehn <friehn@lip.pt>
Date: Thu, 19 Sep 2024 13:31:39 +0000
Subject: [PATCH] Resolve "ParticleCut for tau leptons"
---
applications/c8_air_shower.cpp | 22 ++++++++++++-------
corsika/detail/modules/ParticleCut.inl | 13 ++++++++---
corsika/modules/ParticleCut.hpp | 8 ++++---
examples/cascade_examples/em_shower.cpp | 4 ++--
examples/cascade_examples/mars.cpp | 3 ++-
examples/cascade_examples/radio_em_shower.cpp | 4 ++--
examples/cascade_examples/water.cpp | 3 ++-
tests/modules/testParticleCut.cpp | 14 ++++++------
8 files changed, 44 insertions(+), 27 deletions(-)
diff --git a/applications/c8_air_shower.cpp b/applications/c8_air_shower.cpp
index e215f6022..998a5bd0a 100644
--- a/applications/c8_air_shower.cpp
+++ b/applications/c8_air_shower.cpp
@@ -195,6 +195,10 @@ int main(int argc, char** argv) {
->default_val(0.3)
->check(CLI::Range(0.000001, 1.e13))
->group("Config");
+ app.add_option("--taucut", "Min. kin. energy of tau leptons in tracking (GeV)")
+ ->default_val(0.3)
+ ->check(CLI::Range(0.000001, 1.e13))
+ ->group("Config");
bool track_neutrinos = false;
app.add_flag("--track-neutrinos", track_neutrinos, "switch on tracking of neutrinos")
->group("Config");
@@ -459,17 +463,19 @@ int main(int argc, char** argv) {
HEPEnergyType const emcut = 1_GeV * app["--emcut"]->as<double>();
HEPEnergyType const hadcut = 1_GeV * app["--hadcut"]->as<double>();
HEPEnergyType const mucut = 1_GeV * app["--mucut"]->as<double>();
- ParticleCut<SubWriter<decltype(dEdX)>> cut(emcut, emcut, hadcut, mucut,
+ HEPEnergyType const taucut = 1_GeV * app["--taucut"]->as<double>();
+ ParticleCut<SubWriter<decltype(dEdX)>> cut(emcut, emcut, hadcut, mucut, taucut,
!track_neutrinos, dEdX);
// tell proposal that we are interested in all energy losses above the particle cut
- set_energy_production_threshold(Code::Electron, std::min({emcut, hadcut, mucut}));
- set_energy_production_threshold(Code::Positron, std::min({emcut, hadcut, mucut}));
- set_energy_production_threshold(Code::Photon, std::min({emcut, hadcut, mucut}));
- set_energy_production_threshold(Code::MuMinus, std::min({emcut, hadcut, mucut}));
- set_energy_production_threshold(Code::MuPlus, std::min({emcut, hadcut, mucut}));
- set_energy_production_threshold(Code::TauMinus, std::min({emcut, hadcut, mucut}));
- set_energy_production_threshold(Code::TauPlus, std::min({emcut, hadcut, mucut}));
+ auto const prod_threshold = std::min({emcut, hadcut, mucut, taucut});
+ set_energy_production_threshold(Code::Electron, prod_threshold);
+ set_energy_production_threshold(Code::Positron, prod_threshold);
+ set_energy_production_threshold(Code::Photon, prod_threshold);
+ set_energy_production_threshold(Code::MuMinus, prod_threshold);
+ set_energy_production_threshold(Code::MuPlus, prod_threshold);
+ set_energy_production_threshold(Code::TauMinus, prod_threshold);
+ set_energy_production_threshold(Code::TauPlus, prod_threshold);
// energy threshold for high energy hadronic model. Affects LE/HE switch for
// hadron interactions and the hadronic photon model in proposal
diff --git a/corsika/detail/modules/ParticleCut.inl b/corsika/detail/modules/ParticleCut.inl
index 514a93650..a5d6a03f2 100644
--- a/corsika/detail/modules/ParticleCut.inl
+++ b/corsika/detail/modules/ParticleCut.inl
@@ -17,19 +17,23 @@ namespace corsika {
inline ParticleCut<TOutput>::ParticleCut(HEPEnergyType const eEleCut,
HEPEnergyType const ePhoCut,
HEPEnergyType const eHadCut,
- HEPEnergyType const eMuCut, bool const inv,
+ HEPEnergyType const eMuCut,
+ HEPEnergyType const eTauCut, bool const inv,
TArgs&&... outputArgs)
: TOutput(std::forward<TArgs>(outputArgs)...)
, cut_electrons_(eEleCut)
, cut_photons_(ePhoCut)
- , cut_muons_(eMuCut)
, cut_hadrons_(eHadCut)
+ , cut_muons_(eMuCut)
+ , cut_tau_(eTauCut)
, doCutInv_(inv) {
for (auto p : get_all_particles()) {
if (is_hadron(p)) // nuclei are also hadrons
set_kinetic_energy_propagation_threshold(p, eHadCut);
else if (is_muon(p))
set_kinetic_energy_propagation_threshold(p, eMuCut);
+ else if (p == Code::TauMinus || p == Code::TauPlus)
+ set_kinetic_energy_propagation_threshold(p, eTauCut);
else if (p == Code::Electron || p == Code::Positron)
set_kinetic_energy_propagation_threshold(p, eEleCut);
else if (p == Code::Photon)
@@ -40,7 +44,8 @@ namespace corsika {
"setting kinetic energy thresholds: electrons = {} GeV, photons = {} GeV, "
"hadrons = {} GeV, "
"muons = {} GeV",
- eEleCut / 1_GeV, ePhoCut / 1_GeV, eHadCut / 1_GeV, eMuCut / 1_GeV);
+ "tau = {} GeV", eEleCut / 1_GeV, ePhoCut / 1_GeV, eHadCut / 1_GeV, eMuCut / 1_GeV,
+ eTauCut / 1_GeV);
}
template <typename TOutput>
@@ -152,6 +157,7 @@ namespace corsika {
CORSIKA_LOG_DEBUG("kinetic energy threshold for photons is {} GeV",
cut_photons_ / 1_GeV);
CORSIKA_LOG_DEBUG("kinetic energy threshold for muons is {} GeV", cut_muons_ / 1_GeV);
+ CORSIKA_LOG_DEBUG("kinetic energy threshold for tau is {} GeV", cut_tau_ / 1_GeV);
CORSIKA_LOG_DEBUG("kinetic energy threshold for hadrons is {} GeV",
cut_hadrons_ / 1_GeV);
@@ -171,6 +177,7 @@ namespace corsika {
node["cut_photons"] = cut_photons_ / 1_GeV;
node["cut_muons"] = cut_muons_ / 1_GeV;
node["cut_hadrons"] = cut_hadrons_ / 1_GeV;
+ node["cut_tau"] = cut_tau_ / 1_GeV;
node["cut_invisibles"] = doCutInv_;
for (auto const& cut : cuts_) {
node[fmt::format("cut_{}", cut.first)] = cut.second / 1_GeV;
diff --git a/corsika/modules/ParticleCut.hpp b/corsika/modules/ParticleCut.hpp
index 431f1bc04..04fdfb231 100644
--- a/corsika/modules/ParticleCut.hpp
+++ b/corsika/modules/ParticleCut.hpp
@@ -44,8 +44,8 @@ namespace corsika {
*/
template <typename... TArgs>
ParticleCut(HEPEnergyType const eEleCut, HEPEnergyType const ePhoCut,
- HEPEnergyType const eHadCut, HEPEnergyType const eMuCut, bool const inv,
- TArgs&&... args);
+ HEPEnergyType const eHadCut, HEPEnergyType const eMuCut,
+ HEPEnergyType const eTauCut, bool const inv, TArgs&&... args);
/**
* particle cut with kinetic energy thresholds for all particles.
@@ -103,6 +103,7 @@ namespace corsika {
HEPEnergyType getElectronKineticECut() const { return cut_electrons_; }
HEPEnergyType getPhotonKineticECut() const { return cut_photons_; }
HEPEnergyType getMuonKineticECut() const { return cut_muons_; }
+ HEPEnergyType getTauKineticECut() const { return cut_tau_; }
HEPEnergyType getHadronKineticECut() const { return cut_hadrons_; }
//! get configuration of this node, for output
@@ -116,8 +117,9 @@ namespace corsika {
private:
HEPEnergyType cut_electrons_;
HEPEnergyType cut_photons_;
- HEPEnergyType cut_muons_;
HEPEnergyType cut_hadrons_;
+ HEPEnergyType cut_muons_;
+ HEPEnergyType cut_tau_;
bool doCutInv_;
std::unordered_map<Code const, HEPEnergyType const> cuts_;
}; // namespace corsika
diff --git a/examples/cascade_examples/em_shower.cpp b/examples/cascade_examples/em_shower.cpp
index c3083ec13..c982bcf46 100644
--- a/examples/cascade_examples/em_shower.cpp
+++ b/examples/cascade_examples/em_shower.cpp
@@ -146,8 +146,8 @@ int main(int argc, char** argv) {
// setup processes, decays and interactions
EnergyLossWriter energyloss{showerAxis, dX};
- ParticleCut<SubWriter<decltype(energyloss)>> cut(5_MeV, 5_MeV, 100_GeV, 100_GeV, true,
- energyloss);
+ ParticleCut<SubWriter<decltype(energyloss)>> cut(5_MeV, 5_MeV, 100_GeV, 100_GeV,
+ 100_GeV, true, energyloss);
corsika::sibyll::Interaction sibyll(corsika::get_all_elements_in_universe(env),
corsika::setup::C7trackedParticles);
diff --git a/examples/cascade_examples/mars.cpp b/examples/cascade_examples/mars.cpp
index b19d1de9e..584202857 100644
--- a/examples/cascade_examples/mars.cpp
+++ b/examples/cascade_examples/mars.cpp
@@ -313,7 +313,8 @@ int main(int argc, char** argv) {
HEPEnergyType const emcut = 1_GeV;
HEPEnergyType const hadcut = 1_GeV;
- ParticleCut<SubWriter<decltype(dEdX)>> cut(emcut, emcut, hadcut, hadcut, true, dEdX);
+ ParticleCut<SubWriter<decltype(dEdX)>> cut(emcut, emcut, hadcut, hadcut, hadcut, true,
+ dEdX);
// tell proposal that we are interested in all energy losses above the particle cut
set_energy_production_threshold(Code::Electron, std::min({emcut, hadcut}));
diff --git a/examples/cascade_examples/radio_em_shower.cpp b/examples/cascade_examples/radio_em_shower.cpp
index 0b89a1f8a..205a25599 100644
--- a/examples/cascade_examples/radio_em_shower.cpp
+++ b/examples/cascade_examples/radio_em_shower.cpp
@@ -231,8 +231,8 @@ int main(int argc, char** argv) {
// setup processes, decays and interactions
EnergyLossWriter energyloss{showerAxis, dX};
- ParticleCut<SubWriter<decltype(energyloss)>> cut(5_MeV, 5_MeV, 100_GeV, 100_GeV, true,
- energyloss);
+ ParticleCut<SubWriter<decltype(energyloss)>> cut(5_MeV, 5_MeV, 100_GeV, 100_GeV,
+ 100_GeV, true, energyloss);
corsika::sibyll::Interaction sibyll(corsika::get_all_elements_in_universe(env),
corsika::setup::C7trackedParticles);
diff --git a/examples/cascade_examples/water.cpp b/examples/cascade_examples/water.cpp
index 26380f0e7..c9f48c618 100644
--- a/examples/cascade_examples/water.cpp
+++ b/examples/cascade_examples/water.cpp
@@ -211,7 +211,8 @@ int main(int argc, char** argv) {
// particle production threshold
HEPEnergyType const emCut = eCut;
HEPEnergyType const hadCut = eCut;
- ParticleCut<SubWriter<decltype(dEdX)>> cut(emCut, emCut, hadCut, hadCut, true, dEdX);
+ ParticleCut<SubWriter<decltype(dEdX)>> cut(emCut, emCut, hadCut, hadCut, hadCut, true,
+ dEdX);
// tell proposal that we are interested in all energy losses above the particle cut
set_energy_production_threshold(Code::Electron, std::min({emCut, hadCut}));
diff --git a/tests/modules/testParticleCut.cpp b/tests/modules/testParticleCut.cpp
index cf61644fd..f305bf28d 100644
--- a/tests/modules/testParticleCut.cpp
+++ b/tests/modules/testParticleCut.cpp
@@ -56,7 +56,7 @@ TEST_CASE("ParticleCut", "process,continuous,secondary") {
SECTION("cut on particle type: inv") {
// particle cut with 20GeV threshold for all, also cut invisible
- ParticleCut cut(20_GeV, 20_GeV, 20_GeV, 20_GeV, true);
+ ParticleCut cut(20_GeV, 20_GeV, 20_GeV, 20_GeV, 20_GeV, true);
CHECK(cut.getHadronKineticECut() == 20_GeV);
// add primary particle to stack
@@ -83,7 +83,7 @@ TEST_CASE("ParticleCut", "process,continuous,secondary") {
SECTION("cut on particle type: em") {
- ParticleCut cut(1_EeV, 1_EeV, 1_GeV, 1_GeV, false);
+ ParticleCut cut(1_EeV, 1_EeV, 1_GeV, 1_GeV, 1_GeV, false);
// add primary particle to stack
auto particle = stack.addParticle(std::make_tuple(
@@ -105,7 +105,7 @@ TEST_CASE("ParticleCut", "process,continuous,secondary") {
}
SECTION("cut low energy") {
- ParticleCut cut(20_GeV, 20_GeV, 20_GeV, 20_GeV, true);
+ ParticleCut cut(20_GeV, 20_GeV, 20_GeV, 20_GeV, 20_GeV, true);
// add primary particle to stack
auto particle = stack.addParticle(std::make_tuple(
@@ -133,8 +133,8 @@ TEST_CASE("ParticleCut", "process,continuous,secondary") {
CHECK(view.getSize() == 13);
}
- SECTION("cut low energy: electrons, photons, hadrons and muons") {
- ParticleCut cut(5_MeV, 5_MeV, 5_GeV, 5_GeV, true);
+ SECTION("cut low energy: electrons, photons, hadrons, muons and tau") {
+ ParticleCut cut(5_MeV, 5_MeV, 5_GeV, 5_GeV, 5_MeV, true);
// add primary particle to stack
auto particle = stack.addParticle(std::make_tuple(Code::Proton, Eabove - Proton::mass,
@@ -177,7 +177,7 @@ TEST_CASE("ParticleCut", "process,continuous,secondary") {
}
SECTION("cut on time") {
- ParticleCut cut(20_GeV, 20_GeV, 20_GeV, 20_GeV, false);
+ ParticleCut cut(20_GeV, 20_GeV, 20_GeV, 20_GeV, 20_GeV, false);
const TimeType too_late = 1_s;
// add primary particle to stack
@@ -205,7 +205,7 @@ TEST_CASE("ParticleCut", "process,continuous,secondary") {
SECTION("cut on doContinous, just invisibles") {
- ParticleCut cut(20_GeV, 20_GeV, 20_GeV, 20_GeV, true);
+ ParticleCut cut(20_GeV, 20_GeV, 20_GeV, 20_GeV, 20_GeV, true);
// add particles, all with energies above the threshold
// only cut is by species
--
GitLab