added ParticleCut

corsika/detail/modules/particle_cut/ParticleCut.inl

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

corsika/modules/particle_cut/ParticleCut.hpp

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

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

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