From f7582114494d39c4e28714bcbce5ae350effafa0 Mon Sep 17 00:00:00 2001
From: Fan <fan_hu@pku.edu.cn>
Date: Sat, 30 Oct 2021 15:24:06 +0800
Subject: [PATCH] book kinetic energy in particle cut
---
corsika/detail/modules/ParticleCut.inl | 24 ++++++++++++------------
tests/modules/testParticleCut.cpp | 13 +++++++------
2 files changed, 19 insertions(+), 18 deletions(-)
diff --git a/corsika/detail/modules/ParticleCut.inl b/corsika/detail/modules/ParticleCut.inl
index 0e6163af5..44bcdf424 100644
--- a/corsika/detail/modules/ParticleCut.inl
+++ b/corsika/detail/modules/ParticleCut.inl
@@ -133,26 +133,26 @@ namespace corsika {
CORSIKA_LOG_DEBUG("p={}", particle.asString());
if (doCutEm_ && is_em(pid)) {
CORSIKA_LOG_DEBUG("removing em. particle...");
- energy_emcut_ += energy;
+ energy_emcut_ += kine_energy;
em_count_ += 1;
- energy_event_ += energy;
+ energy_event_ += kine_energy;
return true;
} else if (doCutInv_ && is_neutrino(pid)) {
CORSIKA_LOG_DEBUG("removing inv. particle...");
- energy_invcut_ += energy;
+ energy_invcut_ += kine_energy;
inv_count_ += 1;
- energy_event_ += energy;
+ energy_event_ += kine_energy;
return true;
} else if (isBelowEnergyCut(particle)) {
CORSIKA_LOG_DEBUG("removing low en. particle...");
- energy_cut_ += energy;
+ energy_cut_ += kine_energy;
energy_count_ += 1;
- energy_event_ += energy;
+ energy_event_ += kine_energy;
return true;
} else if (particle.getTime() > 10_ms) {
CORSIKA_LOG_DEBUG("removing OLD particle...");
- energy_timecut_ += energy;
- energy_event_ += energy;
+ energy_timecut_ += kine_energy;
+ energy_event_ += kine_energy;
return true;
}
return false; // this particle will not be removed/cut
@@ -192,10 +192,10 @@ namespace corsika {
inline void ParticleCut::showResults() {
CORSIKA_LOG_INFO(
"\n ******************************\n "
- " energy removed by cut of electromagnetic (GeV): {} (number: {})\n "
- " energy removed by cut of invisible (GeV): {} (number: {})\n "
- " energy removed by kinetic energy cut (GeV): {} (number: {}) \n "
- " energy removed by time cut (GeV): {} \n"
+ " kinetic energy removed by cut of electromagnetic (GeV): {} (number: {})\n "
+ " kinetic energy removed by cut of invisible (GeV): {} (number: {})\n "
+ " kinetic energy removed by kinetic energy cut (GeV): {} (number: {}) \n "
+ " kinetic energy removed by time cut (GeV): {} \n"
" ******************************",
energy_emcut_ / 1_GeV, em_count_, energy_invcut_ / 1_GeV, inv_count_,
energy_cut_ / 1_GeV, energy_count_, energy_timecut_ / 1_GeV);
diff --git a/tests/modules/testParticleCut.cpp b/tests/modules/testParticleCut.cpp
index 453764066..661e961e9 100644
--- a/tests/modules/testParticleCut.cpp
+++ b/tests/modules/testParticleCut.cpp
@@ -76,7 +76,7 @@ TEST_CASE("ParticleCut", "process,continuous,secondary") {
CHECK(view.getEntries() == 9);
CHECK(cut.getNumberInvParticles() == 2);
- CHECK(cut.getInvEnergy() / 1_GeV == 2000);
+ CHECK(cut.getInvEnergy() / 1_GeV == Approx(2000.));
}
SECTION("cut on particle type: em") {
@@ -101,14 +101,14 @@ TEST_CASE("ParticleCut", "process,continuous,secondary") {
CHECK(view.getEntries() == 10);
CHECK(cut.getNumberEmParticles() == 1);
- CHECK(cut.getEmEnergy() / 1_GeV == Approx(1000. + 0.000511));
+ CHECK(cut.getEmEnergy() / 1_GeV == 1000.);
}
SECTION("cut low energy") {
ParticleCut cut(20_GeV, true, true);
// add primary particle to stack
- auto particle = stack.addParticle(std::make_tuple(Code::Proton, Eabove - Proton::mass,
+ auto particle = stack.addParticle(std::make_tuple(Code::Proton, Eabove,
DirectionVector(rootCS, {1, 0, 0}),
point0, 0_ns));
// view on secondary particles
@@ -134,6 +134,7 @@ TEST_CASE("ParticleCut", "process,continuous,secondary") {
CHECK(view.getEntries() == 1);
CHECK(view.getSize() == 13);
+ CHECK(cut.getEmEnergy() == Ebelow * A);
}
SECTION("cut low energy: electrons, photons, hadrons and muons") {
@@ -151,10 +152,10 @@ TEST_CASE("ParticleCut", "process,continuous,secondary") {
// add secondaries
projectile.addSecondary(std::make_tuple(
Code::Photon, 3_MeV, DirectionVector(rootCS, {1, 0, 0}), point0, 0_ns));
- projectile.addSecondary(std::make_tuple(Code::Electron, 3_MeV - Electron::mass,
+ projectile.addSecondary(std::make_tuple(Code::Electron, 3_MeV,
DirectionVector(rootCS, {1, 0, 0}), point0,
0_ns));
- projectile.addSecondary(std::make_tuple(Code::PiPlus, 4_GeV - PiPlus::mass,
+ projectile.addSecondary(std::make_tuple(Code::PiPlus, 4_GeV,
DirectionVector(rootCS, {1, 0, 0}), point0,
0_ns));
@@ -197,7 +198,7 @@ TEST_CASE("ParticleCut", "process,continuous,secondary") {
// add secondaries, all with energies above the threshold
// only cut is by time
for (auto proType : particleList) {
- projectile.addSecondary(std::make_tuple(proType, Eabove - get_mass(proType),
+ projectile.addSecondary(std::make_tuple(proType, Eabove,
DirectionVector(rootCS, {1, 0, 0}), point0,
too_late));
}
--
GitLab