From a1736558a7b3ffd48033d468764fb8db633c8bf9 Mon Sep 17 00:00:00 2001
From: ralfulrich <ralf.ulrich@kit.edu>
Date: Tue, 20 Oct 2020 12:13:42 +0200
Subject: [PATCH] added ParticleCut
---
.../modules/particle_cut/ParticleCut.inl | 113 ++++++++++++++++++
corsika/modules/particle_cut/ParticleCut.hpp | 54 +++++++++
tests/modules/CMakeLists.txt | 2 +-
tests/modules/testParticleCut.cpp | 73 ++++++-----
4 files changed, 204 insertions(+), 38 deletions(-)
create mode 100644 corsika/detail/modules/particle_cut/ParticleCut.inl
create mode 100644 corsika/modules/particle_cut/ParticleCut.hpp
diff --git a/corsika/detail/modules/particle_cut/ParticleCut.inl b/corsika/detail/modules/particle_cut/ParticleCut.inl
new file mode 100644
index 000000000..bd739a5f5
--- /dev/null
+++ b/corsika/detail/modules/particle_cut/ParticleCut.inl
@@ -0,0 +1,113 @@
+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#pragma once
+
+#include <corsika/modules/particle_cut/ParticleCut.hpp>
+
+namespace corsika::particle_cut {
+
+ template <typename TParticle>
+ bool ParticleCut::ParticleIsBelowEnergyCut(TParticle const& vP) const {
+ auto const energyLab = vP.GetEnergy();
+ // nuclei
+ if (vP.GetPID() == corsika::Code::Nucleus) {
+ // calculate energy per nucleon
+ auto const ElabNuc = energyLab / vP.GetNuclearA();
+ return (ElabNuc < fECut);
+ } else {
+ return (energyLab < fECut);
+ }
+ }
+
+ bool ParticleCut::ParticleIsEmParticle(Code vCode) const {
+ // FOR NOW: switch
+ switch (vCode) {
+ case Code::Gamma:
+ case Code::Electron:
+ case Code::Positron:
+ return true;
+ default:
+ return false;
+ }
+ }
+
+ bool ParticleCut::ParticleIsInvisible(Code vCode) const {
+ switch (vCode) {
+ case Code::NuE:
+ case Code::NuEBar:
+ case Code::NuMu:
+ case Code::NuMuBar:
+ return true;
+
+ default:
+ return false;
+ }
+ }
+
+ EProcessReturn ParticleCut::DoSecondaries(corsika::setup::StackView& vS) {
+
+ using namespace units::si;
+
+ auto p = vS.begin();
+ while (p != vS.end()) {
+ const Code pid = p.GetPID();
+ units::si::HEPEnergyType energy = p.GetEnergy();
+ std::cout << "ProcessCut: DoSecondaries: " << pid << " E= " << energy
+ << ", EcutTot=" << (fEmEnergy + fInvEnergy + fEnergy) / 1_GeV << " GeV"
+ << std::endl;
+ if (ParticleIsEmParticle(pid)) {
+ std::cout << "removing em. particle..." << std::endl;
+ fEmEnergy += energy;
+ fEmCount += 1;
+ p.Delete();
+ } else if (ParticleIsInvisible(pid)) {
+ std::cout << "removing inv. particle..." << std::endl;
+ fInvEnergy += energy;
+ fInvCount += 1;
+ p.Delete();
+ } else if (ParticleIsBelowEnergyCut(p)) {
+ std::cout << "removing low en. particle..." << std::endl;
+ fEnergy += energy;
+ p.Delete();
+ } else if (p.GetTime() > 10_ms) {
+ std::cout << "removing OLD particle..." << std::endl;
+ fEnergy += energy;
+ p.Delete();
+ } else {
+ ++p; // next entry in SecondaryView
+ }
+ }
+ return EProcessReturn::eOk;
+ }
+
+ void ParticleCut::Init() {
+ using namespace units::si;
+ fEmEnergy = 0_GeV;
+ fEmCount = 0;
+ fInvEnergy = 0_GeV;
+ fInvCount = 0;
+ fEnergy = 0_GeV;
+ // defineEmParticles();
+ }
+
+ void ParticleCut::ShowResults() {
+ using namespace units::si;
+ std::cout << " ******************************" << std::endl
+ << " ParticleCut: " << std::endl
+ << " energy in em. component (GeV): " << fEmEnergy / 1_GeV << std::endl
+ << " no. of em. particles injected: " << fEmCount << std::endl
+ << " energy in inv. component (GeV): " << fInvEnergy / 1_GeV << std::endl
+ << " no. of inv. particles injected: " << fInvCount << std::endl
+ << " energy below particle cut (GeV): " << fEnergy / 1_GeV << std::endl
+ << " ******************************" << std::endl;
+ }
+
+} // namespace corsika::particle_cut
diff --git a/corsika/modules/particle_cut/ParticleCut.hpp b/corsika/modules/particle_cut/ParticleCut.hpp
new file mode 100644
index 000000000..57ef280bd
--- /dev/null
+++ b/corsika/modules/particle_cut/ParticleCut.hpp
@@ -0,0 +1,54 @@
+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/sequence/SecondariesProcess.hpp>
+#include <corsika/setup/SetupStack.hpp>
+
+namespace corsika::particle_cut {
+
+ class ParticleCut : public corsika::SecondariesProcess<ParticleCut> {
+
+ units::si::HEPEnergyType const fECut;
+
+ units::si::HEPEnergyType fEnergy = 0 * units::si::electronvolt;
+ units::si::HEPEnergyType fEmEnergy = 0 * units::si::electronvolt;
+ unsigned int fEmCount = 0;
+ units::si::HEPEnergyType fInvEnergy = 0 * units::si::electronvolt;
+ unsigned int fInvCount = 0;
+
+ public:
+ ParticleCut(const units::si::HEPEnergyType vCut)
+ : fECut(vCut) {}
+
+ bool ParticleIsInvisible(corsika::Code) const;
+ EProcessReturn DoSecondaries(corsika::setup::StackView&);
+
+ template <typename TParticle>
+ bool ParticleIsBelowEnergyCut(TParticle const&) const;
+
+ bool ParticleIsEmParticle(corsika::Code) const;
+
+ void Init();
+ void ShowResults();
+
+ units::si::HEPEnergyType GetInvEnergy() const { return fInvEnergy; }
+ units::si::HEPEnergyType GetCutEnergy() const { return fEnergy; }
+ units::si::HEPEnergyType GetEmEnergy() const { return fEmEnergy; }
+ unsigned int GetNumberEmParticles() const { return fEmCount; }
+ unsigned int GetNumberInvParticles() const { return fInvCount; }
+ };
+
+} // namespace corsika::particle_cut
+
+#include <corsika/detail/modules/particle_cut/ParticleCut.inl>
diff --git a/tests/modules/CMakeLists.txt b/tests/modules/CMakeLists.txt
index e4449b878..7aa9031a4 100644
--- a/tests/modules/CMakeLists.txt
+++ b/tests/modules/CMakeLists.txt
@@ -4,7 +4,7 @@ set (test_modules_sources
testSibyll.cpp
testNullModel.cpp
testObservationPlane.cpp
- #testParticleCut.cpp
+ testParticleCut.cpp
#testPythia.cpp
#testQGSJetII.cpp
#testSibyll.cpp
diff --git a/tests/modules/testParticleCut.cpp b/tests/modules/testParticleCut.cpp
index 2375600c8..bc003a779 100644
--- a/tests/modules/testParticleCut.cpp
+++ b/tests/modules/testParticleCut.cpp
@@ -9,29 +9,28 @@
* the license.
*/
-#include <corsika/process/particle_cut/ParticleCut.h>
+#include <corsika/modules/particle_cut/ParticleCut.hpp>
-#include <corsika/environment/Environment.h>
-#include <corsika/geometry/Point.h>
-#include <corsika/geometry/RootCoordinateSystem.h>
-#include <corsika/geometry/Vector.h>
-#include <corsika/units/PhysicalUnits.h>
-#include <corsika/utl/CorsikaFenv.h>
+#include <corsika/media/Environment.hpp>
+#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/framework/geometry/RootCoordinateSystem.hpp>
+#include <corsika/framework/geometry/Vector.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/utility/CorsikaFenv.hpp>
-#include <corsika/setup/SetupStack.h>
+#include <corsika/setup/SetupStack.hpp>
#include <catch2/catch.hpp>
using namespace corsika;
-using namespace corsika::process::particle_cut;
-using namespace corsika::units;
+using namespace corsika::particle_cut;
using namespace corsika::units::si;
TEST_CASE("ParticleCut", "[processes]") {
feenableexcept(FE_INVALID);
- using EnvType = environment::Environment<setup::IEnvironmentModel>;
+ using EnvType = corsika::Environment<setup::IEnvironmentModel>;
EnvType env;
- const geometry::CoordinateSystem& rootCS = env.GetCoordinateSystem();
+ const corsika::CoordinateSystem& rootCS = env.GetCoordinateSystem();
// setup empty particle stack
setup::Stack stack;
@@ -40,11 +39,11 @@ TEST_CASE("ParticleCut", "[processes]") {
const HEPEnergyType Eabove = 1_TeV;
const HEPEnergyType Ebelow = 10_GeV;
// list of arbitrary particles
- std::vector<particles::Code> particleList = {
- particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus,
- particles::Code::KMinus, particles::Code::K0Long, particles::Code::K0Short,
- particles::Code::Electron, particles::Code::MuPlus, particles::Code::NuE,
- particles::Code::Neutron};
+ std::vector<corsika::Code> particleList = {
+ corsika::Code::PiPlus, corsika::Code::PiMinus, corsika::Code::KPlus,
+ corsika::Code::KMinus, corsika::Code::K0Long, corsika::Code::K0Short,
+ corsika::Code::Electron, corsika::Code::MuPlus, corsika::Code::NuE,
+ corsika::Code::Neutron};
SECTION("cut on particle type") {
@@ -52,24 +51,24 @@ TEST_CASE("ParticleCut", "[processes]") {
// add primary particle to stack
auto particle = stack.AddParticle(
- std::tuple<particles::Code, units::si::HEPEnergyType,
- corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
- particles::Code::Proton, Eabove,
- corsika::stack::MomentumVector(rootCS, {0_GeV, 0_GeV, 0_GeV}),
- geometry::Point(rootCS, 0_m, 0_m, 0_m), 0_ns});
+ std::tuple<corsika::Code, units::si::HEPEnergyType,
+ corsika::MomentumVector, corsika::Point, units::si::TimeType>{
+ corsika::Code::Proton, Eabove,
+ corsika::MomentumVector(rootCS, {0_GeV, 0_GeV, 0_GeV}),
+ corsika::Point(rootCS, 0_m, 0_m, 0_m), 0_ns});
// view on secondary particles
- corsika::stack::SecondaryView view(particle);
+ corsika::SecondaryView view(particle);
// ref. to primary particle through the secondary view.
// 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
for (auto proType : particleList)
- projectile.AddSecondary(std::tuple<particles::Code, units::si::HEPEnergyType,
- corsika::stack::MomentumVector, geometry::Point,
+ projectile.AddSecondary(std::tuple<corsika::Code, units::si::HEPEnergyType,
+ corsika::MomentumVector, corsika::Point,
units::si::TimeType>{
- proType, Eabove, corsika::stack::MomentumVector(rootCS, {0_GeV, 0_GeV, 0_GeV}),
- geometry::Point(rootCS, 0_m, 0_m, 0_m), 0_ns});
+ proType, Eabove, corsika::MomentumVector(rootCS, {0_GeV, 0_GeV, 0_GeV}),
+ corsika::Point(rootCS, 0_m, 0_m, 0_m), 0_ns});
cut.DoSecondaries(view);
@@ -81,24 +80,24 @@ TEST_CASE("ParticleCut", "[processes]") {
// add primary particle to stack
auto particle = stack.AddParticle(
- std::tuple<particles::Code, units::si::HEPEnergyType,
- corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
- particles::Code::Proton, Eabove,
- corsika::stack::MomentumVector(rootCS, {0_GeV, 0_GeV, 0_GeV}),
- geometry::Point(rootCS, 0_m, 0_m, 0_m), 0_ns});
+ std::tuple<corsika::Code, units::si::HEPEnergyType,
+ corsika::MomentumVector, corsika::Point, units::si::TimeType>{
+ corsika::Code::Proton, Eabove,
+ corsika::MomentumVector(rootCS, {0_GeV, 0_GeV, 0_GeV}),
+ corsika::Point(rootCS, 0_m, 0_m, 0_m), 0_ns});
// view on secondary particles
- corsika::stack::SecondaryView view(particle);
+ corsika::SecondaryView view(particle);
// ref. to primary particle through the secondary view.
// only this way the secondary view is populated
auto projectile = view.GetProjectile();
// add secondaries, all with energies below the threshold
// only cut is by species
for (auto proType : particleList)
- projectile.AddSecondary(std::tuple<particles::Code, units::si::HEPEnergyType,
- corsika::stack::MomentumVector, geometry::Point,
+ projectile.AddSecondary(std::tuple<corsika::Code, units::si::HEPEnergyType,
+ corsika::MomentumVector, corsika::Point,
units::si::TimeType>{
- proType, Ebelow, corsika::stack::MomentumVector(rootCS, {0_GeV, 0_GeV, 0_GeV}),
- geometry::Point(rootCS, 0_m, 0_m, 0_m), 0_ns});
+ proType, Ebelow, corsika::MomentumVector(rootCS, {0_GeV, 0_GeV, 0_GeV}),
+ corsika::Point(rootCS, 0_m, 0_m, 0_m), 0_ns});
cut.DoSecondaries(view);
--
GitLab