diff --git a/Processes/CMakeLists.txt b/Processes/CMakeLists.txt
index 6afe1cc51acdc689b97a90dec334d727b61c7f60..b54f76d8fd3d16d37de019ffb202015ab15bd663 100644
--- a/Processes/CMakeLists.txt
+++ b/Processes/CMakeLists.txt
@@ -11,6 +11,7 @@ add_subdirectory (HadronicElasticModel)
add_subdirectory (EnergyLoss)
add_subdirectory (UrQMD)
add_subdirectory (SwitchProcess)
+add_subdirectory (ParticleCut)
#add_custom_target(CORSIKAprocesses)
add_library (CORSIKAprocesses INTERFACE)
@@ -23,3 +24,4 @@ add_dependencies(CORSIKAprocesses ProcessStackInspector)
add_dependencies(CORSIKAprocesses ProcessTrackingLine)
add_dependencies(CORSIKAprocesses ProcessEnergyLoss)
add_dependencies(CORSIKAprocesses ProcessUrQMD)
+add_dependencies(CORSIKAprocesses ProcessParticleCut)
diff --git a/Processes/ParticleCut/CMakeLists.txt b/Processes/ParticleCut/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..48b539b47d273f0c1e359cb3f2d4f68cc07f1ea9
--- /dev/null
+++ b/Processes/ParticleCut/CMakeLists.txt
@@ -0,0 +1,67 @@
+set (
+ MODEL_SOURCES
+ ParticleCut.cc
+)
+
+set (
+ MODEL_HEADERS
+ ParticleCut.h
+ )
+
+set (
+ MODEL_NAMESPACE
+ corsika/process/particle_cut
+ )
+
+add_library (ProcessParticleCut STATIC ${MODEL_SOURCES})
+CORSIKA_COPY_HEADERS_TO_NAMESPACE (ProcessParticleCut ${MODEL_NAMESPACE} ${MODEL_HEADERS})
+
+set_target_properties (
+ ProcessParticleCut
+ PROPERTIES
+ VERSION ${PROJECT_VERSION}
+ SOVERSION 1
+# PUBLIC_HEADER "${MODEL_HEADERS}"
+ )
+
+# target dependencies on other libraries (also the header onlys)
+target_link_libraries (
+ ProcessParticleCut
+ CORSIKAunits
+ CORSIKAparticles
+ CORSIKAsetup
+ )
+
+target_include_directories (
+ ProcessParticleCut
+ INTERFACE
+ $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
+ $<INSTALL_INTERFACE:include/include>
+ )
+
+install (
+ TARGETS ProcessParticleCut
+ LIBRARY DESTINATION lib
+ ARCHIVE DESTINATION lib
+# PUBLIC_HEADER DESTINATION include/${MODEL_NAMESPACE}
+ )
+
+# --------------------
+# code unit testing
+add_executable (testParticleCut
+ testParticleCut.cc
+ ${MODEL_HEADERS}
+ )
+
+target_link_libraries (
+ testParticleCut
+ ProcessParticleCut
+ CORSIKAunits
+ CORSIKAstackinterface
+ CORSIKAprocesssequence
+ CORSIKAsetup
+ CORSIKAgeometry
+ CORSIKAenvironment
+ CORSIKAthirdparty # for catch2
+ )
+CORSIKA_ADD_TEST(testParticleCut)
diff --git a/Processes/ParticleCut/ParticleCut.cc b/Processes/ParticleCut/ParticleCut.cc
new file mode 100644
index 0000000000000000000000000000000000000000..455a9b8e66ea53e455b797e914ef931d2ff50f47
--- /dev/null
+++ b/Processes/ParticleCut/ParticleCut.cc
@@ -0,0 +1,156 @@
+
+/*
+ * (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.
+ */
+
+#include <corsika/process/particle_cut/ParticleCut.h>
+
+using namespace std;
+
+using namespace corsika;
+using namespace corsika::process;
+using namespace corsika::units::si;
+using namespace corsika::particles;
+using namespace corsika::setup;
+
+namespace corsika::process {
+ namespace particle_cut {
+
+ template <typename Particle>
+ bool ParticleCut::ParticleIsBelowEnergyCut(Particle& vP) const {
+ auto const energyLab = vP.GetEnergy();
+ // nuclei
+ if (vP.GetPID() == particles::Code::Nucleus) {
+ auto const ElabNuc = energyLab / vP.GetNuclearA();
+ auto const EcmNN = sqrt(2. * ElabNuc * 0.93827_GeV);
+ if (ElabNuc < fECut || EcmNN < 10_GeV)
+ return true;
+ else
+ return false;
+ } else {
+ // TODO: center-of-mass energy hard coded
+ const HEPEnergyType Ecm = sqrt(2. * energyLab * 0.93827_GeV);
+ if (energyLab < fECut || Ecm < 10_GeV)
+ return true;
+ else
+ return false;
+ }
+ }
+
+ bool ParticleCut::ParticleIsEmParticle(Code vCode) const {
+ bool is_em = false;
+ // FOR NOW: switch
+ switch (vCode) {
+ case Code::Electron:
+ is_em = true;
+ break;
+ case Code::Positron:
+ is_em = true;
+ break;
+ case Code::Gamma:
+ is_em = true;
+ break;
+ default:
+ break;
+ }
+ return is_em;
+ }
+
+ bool ParticleCut::ParticleIsInvisible(Code vCode) const {
+ bool is_inv = false;
+ // FOR NOW: switch
+ switch (vCode) {
+ case Code::NuE:
+ is_inv = true;
+ break;
+ case Code::NuEBar:
+ is_inv = true;
+ break;
+ case Code::NuMu:
+ is_inv = true;
+ break;
+ case Code::NuMuBar:
+ is_inv = true;
+ break;
+ case Code::MuPlus:
+ is_inv = true;
+ break;
+ case Code::MuMinus:
+ is_inv = true;
+ break;
+
+ case Code::Neutron:
+ is_inv = true;
+ break;
+
+ case Code::AntiNeutron:
+ is_inv = true;
+ break;
+
+ default:
+ break;
+ }
+ return is_inv;
+ }
+
+ template <typename TSecondaries>
+ EProcessReturn ParticleCut::DoSecondaries(TSecondaries& vS) {
+ auto p = vS.begin();
+ while (p != vS.end()) {
+ const Code pid = p.GetPID();
+ HEPEnergyType energy = p.GetEnergy();
+ cout << "ProcessCut: DoSecondaries: " << pid << " E= " << energy
+ << ", EcutTot=" << (fEmEnergy + fInvEnergy + fEnergy) / 1_GeV << " GeV"
+ << endl;
+ if (ParticleIsEmParticle(pid)) {
+ cout << "removing em. particle..." << endl;
+ fEmEnergy += energy;
+ fEmCount += 1;
+ p.Delete();
+ } else if (ParticleIsInvisible(pid)) {
+ cout << "removing inv. particle..." << endl;
+ fInvEnergy += energy;
+ fInvCount += 1;
+ p.Delete();
+ } else if (ParticleIsBelowEnergyCut(p)) {
+ cout << "removing low en. particle..." << endl;
+ fEnergy += energy;
+ p.Delete();
+ } else if (p.GetTime() > 10_ms) {
+ cout << "removing OLD particle..." << endl;
+ fEnergy += energy;
+ p.Delete();
+ } else {
+ ++p; // next entry in SecondaryView
+ }
+ }
+ return EProcessReturn::eOk;
+ }
+
+ void ParticleCut::Init() {
+ fEmEnergy = 0. * 1_GeV;
+ fEmCount = 0;
+ fInvEnergy = 0. * 1_GeV;
+ fInvCount = 0;
+ fEnergy = 0. * 1_GeV;
+ // defineEmParticles();
+ }
+
+ void ParticleCut::ShowResults() {
+ cout << " ******************************" << endl
+ << " ParticleCut: " << endl
+ << " energy in em. component (GeV): " << fEmEnergy / 1_GeV << endl
+ << " no. of em. particles injected: " << fEmCount << endl
+ << " energy in inv. component (GeV): " << fInvEnergy / 1_GeV << endl
+ << " no. of inv. particles injected: " << fInvCount << endl
+ << " energy below particle cut (GeV): " << fEnergy / 1_GeV << endl
+ << " ******************************" << endl;
+ }
+ } // namespace particle_cut
+} // namespace corsika::process
diff --git a/Processes/ParticleCut/ParticleCut.h b/Processes/ParticleCut/ParticleCut.h
new file mode 100644
index 0000000000000000000000000000000000000000..e715d85c453969da44b3bdd8a72c95933537267d
--- /dev/null
+++ b/Processes/ParticleCut/ParticleCut.h
@@ -0,0 +1,57 @@
+/*
+ * (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.
+ */
+
+#ifndef _corsika_process_particle_cut_ParticleCut_h_
+#define _corsika_process_particle_cut_ParticleCut_h_
+
+#include <corsika/particles/ParticleProperties.h>
+#include <corsika/process/SecondariesProcess.h>
+#include <corsika/units/PhysicalUnits.h>
+
+using namespace corsika;
+using namespace corsika::units;
+using namespace corsika::units::si;
+
+namespace corsika::process {
+ namespace particle_cut {
+ class ParticleCut : public process::SecondariesProcess<ParticleCut> {
+
+ HEPEnergyType fECut;
+
+ HEPEnergyType fEnergy = 0_GeV;
+ HEPEnergyType fEmEnergy = 0_GeV;
+ int fEmCount = 0;
+ HEPEnergyType fInvEnergy = 0_GeV;
+ int fInvCount = 0;
+
+ public:
+ ParticleCut(const HEPEnergyType vCut)
+ : fECut(vCut) {}
+
+ bool ParticleIsInvisible(particles::Code) const;
+ template <typename TSecondaries>
+ EProcessReturn DoSecondaries(TSecondaries&);
+
+ template <typename Particle>
+ bool ParticleIsBelowEnergyCut(Particle&) const;
+
+ bool ParticleIsEmParticle(particles::Code) const;
+
+ void Init();
+ void ShowResults();
+
+ HEPEnergyType GetInvEnergy() const { return fInvEnergy; }
+ HEPEnergyType GetCutEnergy() const { return fEnergy; }
+ HEPEnergyType GetEmEnergy() const { return fEmEnergy; }
+ };
+ } // namespace particlecut
+} // namespace corsika::process
+
+#endif
diff --git a/Processes/ParticleCut/testParticleCut.cc b/Processes/ParticleCut/testParticleCut.cc
new file mode 100644
index 0000000000000000000000000000000000000000..f315cb17e98d6a51f9a653ab4a39d976ab685f75
--- /dev/null
+++ b/Processes/ParticleCut/testParticleCut.cc
@@ -0,0 +1,72 @@
+
+/*
+ * (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.
+ */
+
+#include <corsika/process/particle_cut/ParticleCut.h>
+
+#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/setup/SetupStack.h>
+
+#define CATCH_CONFIG_MAIN // This tells Catch to provide a main() - only do this in one
+ // cpp file
+#include <catch2/catch.hpp>
+
+using namespace corsika;
+using namespace corsika::process::particle_cut;
+using namespace corsika::units;
+using namespace corsika::units::si;
+
+TEST_CASE("ParticleCut", "[processes]") {
+ feenableexcept(FE_INVALID);
+ using EnvType = environment::Environment<setup::IEnvironmentModel>;
+ EnvType env;
+ const geometry::CoordinateSystem& rootCS = env.GetCoordinateSystem();
+ SECTION("cut") {
+
+ ParticleCut cut(20_GeV);
+
+ // setup particle stack, and add primary particle
+ setup::Stack stack;
+ stack.Clear();
+ const HEPEnergyType Eabove = 1_TeV;
+ const HEPEnergyType Ebelow = 10_GeV;
+ 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};
+ 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});
+
+ for (auto proType : particleList)
+
+ particle.AddSecondary(std::tuple<particles::Code, units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, geometry::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});
+
+ // auto secViewPtr = corsika::stack::SecondaryView(particle);
+ corsika::stack::SecondaryView view(particle);
+ cut.DoSecondaries(view);
+ REQUIRE(stack.GetSize() == 1);
+
+ }
+}