From 7a6ae779ad04aa4022f22a8fc8e928b9bc5573da Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Wed, 13 Jan 2021 18:47:51 +0000
Subject: [PATCH] added separate energy cut thresholds for electrons, photons,
hadrons and muons
---
corsika/detail/modules/ParticleCut.inl | 48 +++++++++++++++++++++++---
corsika/modules/ParticleCut.hpp | 30 ++++++++++++++--
examples/cascade_example.cpp | 2 +-
examples/cascade_proton_example.cpp | 2 +-
examples/em_shower.cpp | 6 ++--
examples/hybrid_MC.cpp | 2 +-
examples/vertical_EAS.cpp | 6 ++--
tests/modules/testParticleCut.cpp | 43 +++++++++++++++++++++--
8 files changed, 120 insertions(+), 19 deletions(-)
diff --git a/corsika/detail/modules/ParticleCut.inl b/corsika/detail/modules/ParticleCut.inl
index 01d99891b..559102cc4 100644
--- a/corsika/detail/modules/ParticleCut.inl
+++ b/corsika/detail/modules/ParticleCut.inl
@@ -12,8 +12,40 @@
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)
+ , doCutInv_(inv)
+ , energy_(0_GeV)
+ , em_energy_(0_GeV)
+ , em_count_(0)
+ , inv_energy_(0_GeV)
+ , inv_count_(0) {}
+
+ 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)
+ , doCutInv_(inv)
+ , energy_(0_GeV)
+ , em_energy_(0_GeV)
+ , em_count_(0)
+ , inv_energy_(0_GeV)
+ , inv_count_(0) {}
+
ParticleCut::ParticleCut(const HEPEnergyType eCut, bool em, bool inv)
- : energy_cut_(eCut)
+ : electron_energy_cut_(eCut)
+ , photon_energy_cut_(eCut)
+ , had_energy_cut_(eCut)
+ , mu_energy_cut_(eCut)
, doCutEm_(em)
, doCutInv_(inv)
, energy_(0_GeV)
@@ -25,13 +57,21 @@ namespace corsika {
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 < 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_);
} else {
- return (energyLab < energy_cut_);
+ // assuming the rest are hadrons
+ return (energyLab < had_energy_cut_);
}
}
diff --git a/corsika/modules/ParticleCut.hpp b/corsika/modules/ParticleCut.hpp
index 098efc005..11d8508ad 100644
--- a/corsika/modules/ParticleCut.hpp
+++ b/corsika/modules/ParticleCut.hpp
@@ -17,11 +17,29 @@
#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:
+ /**
+ * 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(const HEPEnergyType eEleCut, const HEPEnergyType ePhoCut,
+ const HEPEnergyType eHadCut, const HEPEnergyType eMuCut, bool inv);
+
+ //! simple cut. hadrons and muons are cut by threshold. EM particles are all
+ //! discarded.
+ ParticleCut(const HEPEnergyType eHadCut, const HEPEnergyType eMuCut, bool inv);
+
+ //! simplest cut. all particles have same threshold. EM particles can be set to be
+ //! discarded altogether.
ParticleCut(const HEPEnergyType eCut, bool em, bool inv);
void doSecondaries(corsika::setup::StackView&);
@@ -35,7 +53,10 @@ namespace corsika {
void showResults();
void reset();
- HEPEnergyType getECut() const { return energy_cut_; }
+ 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 getInvEnergy() const { return inv_energy_; }
HEPEnergyType getCutEnergy() const { return energy_; }
HEPEnergyType getEmEnergy() const { return em_energy_; }
@@ -50,7 +71,10 @@ namespace corsika {
bool isBelowEnergyCut(TParticle const&) const;
private:
- HEPEnergyType energy_cut_;
+ 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/examples/cascade_example.cpp b/examples/cascade_example.cpp
index d14843c51..eea24768d 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, cut.getElectronECut()};
// 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 003dc8558..e20864685 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, cut.getElectronECut()};
// 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 edb377953..d362ae14f 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, cut.getElectronECut());
+ corsika::proposal::ContinuousProcess em_continuous(env, cut.getElectronECut());
InteractionCounter proposalCounted(proposal);
TrackWriter trackWriter("tracks.dat");
diff --git a/examples/hybrid_MC.cpp b/examples/hybrid_MC.cpp
index 71d8cf828..72e38e49d 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, cut.getElectronECut()};
CONEXhybrid conex_model(center, showerAxis, t, injectionHeight, E0,
get_PDG(Code::Proton));
diff --git a/examples/vertical_EAS.cpp b/examples/vertical_EAS.cpp
index 28488a1bd..46bd64522 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, cut.getElectronECut());
+ corsika::proposal::ContinuousProcess em_continuous(env, cut.getElectronECut());
InteractionCounter proposalCounted(proposal);
OnShellCheck reset_particle_mass(1.e-3, 1.e-1, false);
diff --git a/tests/modules/testParticleCut.cpp b/tests/modules/testParticleCut.cpp
index a88942301..1ac9beec5 100644
--- a/tests/modules/testParticleCut.cpp
+++ b/tests/modules/testParticleCut.cpp
@@ -24,7 +24,7 @@ using namespace corsika;
TEST_CASE("ParticleCut", "[processes]") {
- logging::set_level(logging::level::info);
+ logging::set_level(logging::level::debug);
corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
feenableexcept(FE_INVALID);
@@ -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,43 @@ 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 on time") {
ParticleCut cut(20_GeV, false, false);
const TimeType too_late = 1_s;
@@ -170,7 +207,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
--
GitLab