added ParticleCut proess

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)

Processes/ParticleCut/CMakeLists.txt

+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)

Processes/ParticleCut/ParticleCut.cc

+
+/*
+ * (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

Processes/ParticleCut/ParticleCut.h

+/*
+ * (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

Processes/ParticleCut/testParticleCut.cc

+
+/*
+ * (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);
+    
+  }
+}