diff --git a/corsika/detail/framework/core/ParticleProperties.inl b/corsika/detail/framework/core/ParticleProperties.inl
index 47504f4ea5f0dd02672418e260eb61dd6971c5c3..6e56d357a7f3891637c0bf89b1fe9593b342137f 100644
--- a/corsika/detail/framework/core/ParticleProperties.inl
+++ b/corsika/detail/framework/core/ParticleProperties.inl
@@ -13,11 +13,12 @@
namespace corsika {
- inline HEPEnergyType constexpr get_energy_threshold(Code const p) {
+ inline HEPEnergyType constexpr get_kinetic_energy_threshold(Code const p) {
return particle::detail::thresholds[static_cast<CodeIntType>(p)];
}
- inline void constexpr set_energy_threshold(Code const p, HEPEnergyType const val) {
+ inline void constexpr set_kinetic_energy_threshold(Code const p,
+ HEPEnergyType const val) {
particle::detail::thresholds[static_cast<CodeIntType>(p)] = val;
}
diff --git a/corsika/detail/modules/ParticleCut.inl b/corsika/detail/modules/ParticleCut.inl
index ffbbbed8ddd51cb99ffe45f92217ed426fe2ee91..04aff2d69f0e237541c73389a10b5fc7823004b0 100644
--- a/corsika/detail/modules/ParticleCut.inl
+++ b/corsika/detail/modules/ParticleCut.inl
@@ -18,25 +18,27 @@ namespace corsika {
bool const inv)
: doCutEm_(false)
, doCutInv_(inv)
- , energy_(0_GeV)
- , em_energy_(0_GeV)
+ , energy_cut_(0_GeV)
+ , energy_timecut_(0_GeV)
+ , energy_emcut_(0_GeV)
+ , energy_invcut_(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);
+ set_kinetic_energy_threshold(p, eHadCut);
else if (is_muon(p))
- set_energy_threshold(p, eMuCut);
+ set_kinetic_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);
+ set_kinetic_energy_threshold(p, eEleCut);
+ else if (p == Code::Photon)
+ set_kinetic_energy_threshold(p, ePhoCut);
else if (p == Code::Nucleus)
// nuclei have same threshold as hadrons on the nucleon level.
- set_energy_threshold(p, eHadCut);
+ set_kinetic_energy_threshold(p, eHadCut);
CORSIKA_LOG_DEBUG(
- "setting thresholds: electrons = {} GeV, photons = {} GeV, hadrons = {} GeV, "
+ "setting kinetic energy thresholds: electrons = {} GeV, photons = {} GeV, "
+ "hadrons = {} GeV, "
"muons = {} GeV",
eEleCut / 1_GeV, ePhoCut / 1_GeV, eHadCut / 1_GeV, eMuCut / 1_GeV);
printThresholds();
@@ -46,17 +48,18 @@ namespace corsika {
bool const inv)
: doCutEm_(true)
, doCutInv_(inv)
- , energy_(0_GeV)
- , em_energy_(0_GeV)
+ , energy_cut_(0_GeV)
+ , energy_timecut_(0_GeV)
+ , energy_emcut_(0_GeV)
+ , energy_invcut_(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);
+ set_kinetic_energy_threshold(p, eHadCut);
else if (is_muon(p))
- set_energy_threshold(p, eMuCut);
+ set_kinetic_energy_threshold(p, eMuCut);
CORSIKA_LOG_DEBUG(
"setting thresholds: hadrons = {} GeV, "
"muons = {} GeV",
@@ -67,13 +70,15 @@ namespace corsika {
inline ParticleCut::ParticleCut(HEPEnergyType const eCut, bool const em, bool const inv)
: doCutEm_(em)
, doCutInv_(inv)
- , energy_(0_GeV)
- , em_energy_(0_GeV)
+ , energy_cut_(0_GeV)
+ , energy_timecut_(0_GeV)
+ , energy_emcut_(0_GeV)
+ , energy_invcut_(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);
+ for (auto p : get_all_particles()) set_kinetic_energy_threshold(p, eCut);
+ CORSIKA_LOG_DEBUG("setting kinetic energy threshold for all particles to {} GeV",
+ eCut / 1_GeV);
printThresholds();
}
@@ -82,27 +87,28 @@ namespace corsika {
bool const inv)
: doCutEm_(em)
, doCutInv_(inv)
- , energy_(0_GeV)
- , em_energy_(0_GeV)
+ , energy_cut_(0_GeV)
+ , energy_timecut_(0_GeV)
+ , energy_emcut_(0_GeV)
+ , energy_invcut_(0_GeV)
, em_count_(0)
- , inv_energy_(0_GeV)
, inv_count_(0) {
- set_energy_thresholds(eCuts);
+ set_kinetic_energy_thresholds(eCuts);
CORSIKA_LOG_DEBUG("setting threshold particles individually");
printThresholds();
}
template <typename TParticle>
inline bool ParticleCut::isBelowEnergyCut(TParticle const& vP) const {
- auto const energyLab = vP.getEnergy();
+ auto const energyLab = vP.getKineticEnergy();
auto const pid = vP.getPID();
// nuclei
if (pid == Code::Nucleus) {
// calculate energy per nucleon
auto const ElabNuc = energyLab / vP.getNuclearA();
- return (ElabNuc < get_energy_threshold(pid));
+ return (ElabNuc < get_kinetic_energy_threshold(pid));
} else {
- return (energyLab < get_energy_threshold(pid));
+ return (energyLab < get_kinetic_energy_threshold(pid));
}
}
@@ -114,26 +120,29 @@ namespace corsika {
inline bool ParticleCut::checkCutParticle(TParticle const& particle) {
Code const pid = particle.getPID();
- HEPEnergyType energy = particle.getEnergy();
- CORSIKA_LOG_DEBUG("ParticleCut: checking {}, E= {} GeV, EcutTot={} GeV", pid,
- energy / 1_GeV, (em_energy_ + inv_energy_ + energy_) / 1_GeV);
+ HEPEnergyType const kine_energy = particle.getKineticEnergy();
+ HEPEnergyType const energy = particle.getEnergy();
+ CORSIKA_LOG_DEBUG(
+ "ParticleCut: checking {}, E_kin= {} GeV, EcutTot={} GeV", pid,
+ kine_energy / 1_GeV,
+ (energy_emcut_ + energy_invcut_ + energy_cut_ + energy_timecut_) / 1_GeV);
if (doCutEm_ && is_em(pid)) {
CORSIKA_LOG_DEBUG("removing em. particle...");
- em_energy_ += energy;
+ energy_emcut_ += energy;
em_count_ += 1;
return true;
} else if (doCutInv_ && is_neutrino(pid)) {
CORSIKA_LOG_DEBUG("removing inv. particle...");
- inv_energy_ += energy;
+ energy_invcut_ += energy;
inv_count_ += 1;
return true;
} else if (isBelowEnergyCut(particle)) {
CORSIKA_LOG_DEBUG("removing low en. particle...");
- energy_ += energy;
+ energy_cut_ += energy;
return true;
} else if (particle.getTime() > 10_ms) {
CORSIKA_LOG_DEBUG("removing OLD particle...");
- energy_ += energy;
+ energy_timecut_ += energy;
return true;
}
return false; // this particle will not be removed/cut
@@ -161,29 +170,33 @@ namespace corsika {
inline 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);
+ auto const Eth = get_kinetic_energy_threshold(p);
+ CORSIKA_LOG_INFO("kinetic energy threshold for particle {} is {} GeV", p,
+ Eth / 1_GeV);
}
}
inline void ParticleCut::showResults() {
CORSIKA_LOG_INFO(
" ******************************\n"
- " energy in em. component (GeV): {} \n "
- " no. of em. particles injected: {} \n "
- " energy in inv. component (GeV): {} \n "
- " no. of inv. particles injected: {} \n "
- " energy below particle cut (GeV): {} \n"
+ " energy removed by cut of electromagnetic (GeV): {} \n "
+ " no. of em. particles removed : {} \n "
+ " energy removed by cut of invisible (GeV): {} \n "
+ " no. of invisible particles removed : {} \n "
+ " energy removed by kinetic energy cut (GeV): {} \n"
+ " energy removed by time cut (GeV): {} \n"
" ******************************",
- em_energy_ / 1_GeV, em_count_, inv_energy_ / 1_GeV, inv_count_, energy_ / 1_GeV);
+ energy_emcut_ / 1_GeV, em_count_, energy_invcut_ / 1_GeV, inv_count_,
+ energy_cut_ / 1_GeV, energy_timecut_ / 1_GeV);
}
inline void ParticleCut::reset() {
- em_energy_ = 0_GeV;
+ energy_emcut_ = 0_GeV;
em_count_ = 0;
- inv_energy_ = 0_GeV;
+ energy_invcut_ = 0_GeV;
inv_count_ = 0;
- energy_ = 0_GeV;
+ energy_cut_ = 0_GeV;
+ energy_timecut_ = 0_GeV;
}
} // namespace corsika
diff --git a/corsika/framework/core/ParticleProperties.hpp b/corsika/framework/core/ParticleProperties.hpp
index d0d17d264d0dbe8e7a7338ecce743482d0e700c1..bb87d66e3ba7159793c61a4bdec9713c4515d8c2 100644
--- a/corsika/framework/core/ParticleProperties.hpp
+++ b/corsika/framework/core/ParticleProperties.hpp
@@ -71,19 +71,19 @@ namespace corsika {
int16_t constexpr get_charge_number(Code const); //!< electric charge in units of e
ElectricChargeType constexpr get_charge(Code const); //!< electric charge
HEPMassType constexpr get_mass(Code const); //!< 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);
+ HEPEnergyType constexpr get_kinetic_energy_threshold(
+ Code const); //!< get kinetic energy threshold below which the particle is
+ //!< discarded, by default set to zero
+ void constexpr set_kinetic_energy_threshold(
+ Code const, HEPEnergyType const); //!< set kinetic energy threshold below which the
+ //!< particle is discarded
+
+ inline void set_kinetic_energy_threshold(std::pair<Code const, HEPEnergyType const> p) {
+ set_kinetic_energy_threshold(p.first, p.second);
}
- inline void set_energy_thresholds(
+ inline void set_kinetic_energy_thresholds(
std::unordered_map<Code const, HEPEnergyType const> const& eCuts) {
- for (auto v : eCuts) set_energy_threshold(v);
+ for (auto v : eCuts) set_kinetic_energy_threshold(v);
}
//! Particle code according to PDG, "Monte Carlo Particle Numbering Scheme"
diff --git a/corsika/modules/ParticleCut.hpp b/corsika/modules/ParticleCut.hpp
index 9842dd034530d7528e1359e63983d8adc85b25f0..7857c691f53024263ab2dd8e284cab8777a9e317 100644
--- a/corsika/modules/ParticleCut.hpp
+++ b/corsika/modules/ParticleCut.hpp
@@ -21,15 +21,18 @@
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.
+ removes particles according to their kinetic energy. Particles with a time delay of
+ more than 10ms are removed as well. Invisible particles (neutrinos) can be removed if
+ selected. The threshold value is set to 0 by default but in principle can be configured
+ for each particle. Special constructors for cuts by the following groups are
+ implemented: (electrons,positrons), photons, hadrons and muons.
**/
class ParticleCut : public SecondariesProcess<ParticleCut>,
public ContinuousProcess<ParticleCut> {
public:
/**
- * particle cut with energy thresholds for electrons, photons,
+ * particle cut with kinetic energy thresholds for electrons, photons,
* hadrons (including nuclei with energy per nucleon) and muons
* invisible particles (neutrinos) can be cut or not
**/
@@ -63,14 +66,30 @@ namespace corsika {
void showResults(); // LCOV_EXCL_LINE
void reset();
- 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_; }
+ HEPEnergyType getElectronKineticECut() const {
+ return get_kinetic_energy_threshold(Code::Electron);
+ }
+ HEPEnergyType getPhotonKineticECut() const {
+ return get_kinetic_energy_threshold(Code::Photon);
+ }
+ HEPEnergyType getMuonKineticECut() const {
+ return get_kinetic_energy_threshold(Code::MuPlus);
+ }
+ HEPEnergyType getHadronKineticECut() const {
+ return get_kinetic_energy_threshold(Code::Proton);
+ }
+ //! returns total energy of particles that were removed by cut for invisible particles
+ HEPEnergyType getInvEnergy() const { return energy_invcut_; }
+ //! returns total energy of particles that were removed by cut in time
+ HEPEnergyType getTimeCutEnergy() const { return energy_timecut_; }
+ //! returns total energy of particles that were removed by cut in kinetic energy
+ HEPEnergyType getCutEnergy() const { return energy_cut_; }
+ //! returns total energy of particles that were removed by cut for electromagnetic
+ //! particles
+ HEPEnergyType getEmEnergy() const { return energy_emcut_; }
+ //! returns number of electromagnetic particles
unsigned int getNumberEmParticles() const { return em_count_; }
+ //! returns number of invisible particles
unsigned int getNumberInvParticles() const { return inv_count_; }
private:
@@ -86,10 +105,11 @@ namespace corsika {
private:
bool doCutEm_;
bool doCutInv_;
- HEPEnergyType energy_ = 0 * electronvolt;
- HEPEnergyType em_energy_ = 0 * electronvolt;
+ HEPEnergyType energy_cut_ = 0 * electronvolt;
+ HEPEnergyType energy_timecut_ = 0 * electronvolt;
+ HEPEnergyType energy_emcut_ = 0 * electronvolt;
+ HEPEnergyType energy_invcut_ = 0 * electronvolt;
unsigned int em_count_ = 0;
- HEPEnergyType inv_energy_ = 0 * electronvolt;
unsigned int inv_count_ = 0;
};
diff --git a/corsika/stack/NuclearStackExtension.hpp b/corsika/stack/NuclearStackExtension.hpp
index 411c353810d9c84c0207a6445e83acd6b0297c87..350a343ca34f25632691dc4829d67b1e358405fa 100644
--- a/corsika/stack/NuclearStackExtension.hpp
+++ b/corsika/stack/NuclearStackExtension.hpp
@@ -62,9 +62,8 @@ namespace corsika::nuclear_stack {
void setParticleData(super_type& p, altenative_particle_data_type const& v);
-
typedef std::tuple<Code, HEPEnergyType, MomentumVector, Point, TimeType>
- particle_data_momentum_type;
+ particle_data_momentum_type;
typedef std::tuple<Code, HEPEnergyType, MomentumVector, Point, TimeType,
unsigned short, unsigned short>
@@ -111,7 +110,7 @@ namespace corsika::nuclear_stack {
/**
* Overwrite normal getMomentum function with nuclear version
- */
+ */
MomentumVector getMomentum() const;
/**
diff --git a/src/framework/core/pdxml_reader.py b/src/framework/core/pdxml_reader.py
index a099dc427ff835ee6bad75bdd5169513ea100d05..e9f5808c3e90b49d76affc2e220197e4b71209b4 100755
--- a/src/framework/core/pdxml_reader.py
+++ b/src/framework/core/pdxml_reader.py
@@ -333,10 +333,10 @@ def gen_properties(particle_db):
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
+ # particle threshold table, initially set to 0
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'])
+ for k in particle_db:
+ string += " 0 * corsika::units::si::electronvolt, // {name:s}\n".format( name = k)
string += "};\n\n"
# PDG code table
diff --git a/tests/framework/testParticles.cpp b/tests/framework/testParticles.cpp
index d8c71e6e787ce465d367c29dfa5d06b0b2a51faf..461dae438e9670a0727b0fa0eb039214f343f37d 100644
--- a/tests/framework/testParticles.cpp
+++ b/tests/framework/testParticles.cpp
@@ -80,12 +80,12 @@ TEST_CASE("ParticleProperties", "[Particles]") {
}
SECTION("Energy threshold") {
- //! by default energy thresholds are set to particle mass
- CHECK(get_energy_threshold(Electron::code) / Electron::mass == Approx(1));
+ //! by default energy thresholds are set to zero
+ CHECK(get_kinetic_energy_threshold(Electron::code) == 0_GeV);
- set_energy_threshold(Electron::code, 10_GeV);
- CHECK_FALSE(get_energy_threshold(Code::Electron) == 1_GeV);
- CHECK(get_energy_threshold(Code::Electron) == 10_GeV);
+ set_kinetic_energy_threshold(Electron::code, 10_GeV);
+ CHECK_FALSE(get_kinetic_energy_threshold(Code::Electron) == 1_GeV);
+ CHECK(get_kinetic_energy_threshold(Code::Electron) == 10_GeV);
}
SECTION("Particle groups: electromagnetic") {
diff --git a/tests/modules/testParticleCut.cpp b/tests/modules/testParticleCut.cpp
index 0a2388772f960667cf99f7431089e5938ae1f66b..cacf4840c572417e9c9c431ad25afa1a9c456c73 100644
--- a/tests/modules/testParticleCut.cpp
+++ b/tests/modules/testParticleCut.cpp
@@ -50,8 +50,9 @@ TEST_CASE("ParticleCut", "processes") {
SECTION("cut on particle type: inv") {
+ // particle cut with 20GeV threshold for all, also cut invisible
ParticleCut cut(20_GeV, false, true);
- CHECK(cut.getHadronECut() == 20_GeV);
+ CHECK(cut.getHadronKineticECut() == 20_GeV);
// add primary particle to stack
auto particle = stack.addParticle(std::make_tuple(
@@ -173,10 +174,11 @@ TEST_CASE("ParticleCut", "processes") {
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);
+ CHECK(get_kinetic_energy_threshold(Code::Electron) !=
+ get_kinetic_energy_threshold(Code::Positron));
+ CHECK_FALSE(get_kinetic_energy_threshold(Code::Electron) == Electron::mass);
// test default values still correct
- CHECK(get_energy_threshold(Code::Proton) == 5_GeV);
+ CHECK(get_kinetic_energy_threshold(Code::Proton) == 5_GeV);
}
SECTION("cut on time") {
@@ -193,7 +195,7 @@ TEST_CASE("ParticleCut", "processes") {
// only this way the secondary view is populated
auto projectile = view.getProjectile();
// add secondaries, all with energies above the threshold
- // only cut is by species
+ // only cut is by time
for (auto proType : particleList) {
projectile.addSecondary(
std::make_tuple(proType, Eabove, MomentumVector(rootCS, {0_GeV, 0_GeV, 0_GeV}),
@@ -202,9 +204,10 @@ TEST_CASE("ParticleCut", "processes") {
cut.doSecondaries(view);
CHECK(view.getEntries() == 0);
- CHECK(cut.getCutEnergy() / 1_GeV == 11000);
- cut.reset();
+ CHECK(cut.getTimeCutEnergy() == 11 * Eabove);
CHECK(cut.getCutEnergy() == 0_GeV);
+ cut.reset();
+ CHECK(cut.getTimeCutEnergy() == 0_GeV);
}
setup::Trajectory const track = setup::testing::make_track<setup::Trajectory>(
@@ -214,7 +217,7 @@ TEST_CASE("ParticleCut", "processes") {
SECTION("cut on DoContinous, just invisibles") {
ParticleCut cut(20_GeV, false, true);
- CHECK(cut.getHadronECut() == 20_GeV);
+ CHECK(cut.getHadronKineticECut() == 20_GeV);
// add particles, all with energies above the threshold
// only cut is by species