added separate energy cut thresholds for electrons, photons, hadrons and muons

corsika/detail/modules/ParticleCut.inl

@@ -12,8 +12,40 @@
 
 namespace corsika {
 
+  ParticleCut::ParticleCut(const HEPEnergyType eEleCut, const HEPEnergyType ePhoCut,
+                           const HEPEnergyType eHadCut, const HEPEnergyType eMuCut,
+                           bool inv)
+      : electron_energy_cut_(eEleCut)
+      , photon_energy_cut_(ePhoCut)
+      , had_energy_cut_(eHadCut)
+      , mu_energy_cut_(eMuCut)
+      , doCutEm_(false)
+      , doCutInv_(inv)
+      , energy_(0_GeV)
+      , em_energy_(0_GeV)
+      , em_count_(0)
+      , inv_energy_(0_GeV)
+      , inv_count_(0) {}
+
+  ParticleCut::ParticleCut(const HEPEnergyType eHadCut, const HEPEnergyType eMuCut,
+                           bool inv)
+      : electron_energy_cut_(0_eV)
+      , photon_energy_cut_(0_eV)
+      , had_energy_cut_(eHadCut)
+      , mu_energy_cut_(eMuCut)
+      , doCutEm_(true)
+      , doCutInv_(inv)
+      , energy_(0_GeV)
+      , em_energy_(0_GeV)
+      , em_count_(0)
+      , inv_energy_(0_GeV)
+      , inv_count_(0) {}
+
   ParticleCut::ParticleCut(const HEPEnergyType eCut, bool em, bool inv)
-      : energy_cut_(eCut)
+      : electron_energy_cut_(eCut)
+      , photon_energy_cut_(eCut)
+      , had_energy_cut_(eCut)
+      , mu_energy_cut_(eCut)
       , doCutEm_(em)
       , doCutInv_(inv)
       , energy_(0_GeV)
@@ -25,13 +57,21 @@ namespace corsika {
   template <typename TParticle>
   bool ParticleCut::isBelowEnergyCut(TParticle const& vP) const {
     auto const energyLab = vP.getEnergy();
+    auto const pid = vP.getPID();
     // nuclei
-    if (vP.getPID() == Code::Nucleus) {
+    if (pid == Code::Nucleus) {
       // calculate energy per nucleon
       auto const ElabNuc = energyLab / vP.getNuclearA();
-      return (ElabNuc < energy_cut_);
+      return (ElabNuc < had_energy_cut_);
+    } else if (pid == Code::Gamma) {
+      return (energyLab < photon_energy_cut_);
+    } else if (pid == Code::Electron || pid == Code::Positron) {
+      return (energyLab < electron_energy_cut_);
+    } else if (is_muon(pid)) {
+      return (energyLab < mu_energy_cut_);
     } else {
-      return (energyLab < energy_cut_);
+      // assuming the rest are hadrons
+      return (energyLab < had_energy_cut_);
     }
   }
 

corsika/modules/ParticleCut.hpp

@@ -17,11 +17,29 @@
 #include <corsika/setup/SetupTrajectory.hpp>
 
 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.
+   **/
   class ParticleCut : public SecondariesProcess<ParticleCut>,
                       public ContinuousProcess<ParticleCut> {
 
   public:
+    /**
+     * particle cut with energy thresholds for electrons, photons,
+     *    hadrons (including nuclei with energy per nucleon) and muons
+     *    invisible particles (neutrinos) can be cut or not
+     **/
+    ParticleCut(const HEPEnergyType eEleCut, const HEPEnergyType ePhoCut,
+                const HEPEnergyType eHadCut, const HEPEnergyType eMuCut, bool inv);
+
+    //! simple cut. hadrons and muons are cut by threshold. EM particles are all
+    //! discarded.
+    ParticleCut(const HEPEnergyType eHadCut, const HEPEnergyType eMuCut, bool inv);
+
+    //! simplest cut. all particles have same threshold. EM particles can be set to be
+    //! discarded altogether.
     ParticleCut(const HEPEnergyType eCut, bool em, bool inv);
 
     void doSecondaries(corsika::setup::StackView&);
@@ -35,7 +53,10 @@ namespace corsika {
     void showResults();
     void reset();
 
-    HEPEnergyType getECut() const { return energy_cut_; }
+    HEPEnergyType getElectronECut() const { return electron_energy_cut_; }
+    HEPEnergyType getPhotonECut() const { return photon_energy_cut_; }
+    HEPEnergyType getMuonECut() const { return mu_energy_cut_; }
+    HEPEnergyType getHadronECut() const { return had_energy_cut_; }
     HEPEnergyType getInvEnergy() const { return inv_energy_; }
     HEPEnergyType getCutEnergy() const { return energy_; }
     HEPEnergyType getEmEnergy() const { return em_energy_; }
@@ -50,7 +71,10 @@ namespace corsika {
     bool isBelowEnergyCut(TParticle const&) const;
 
   private:
-    HEPEnergyType energy_cut_;
+    HEPEnergyType electron_energy_cut_;
+    HEPEnergyType photon_energy_cut_;
+    HEPEnergyType had_energy_cut_;
+    HEPEnergyType mu_energy_cut_;
     bool doCutEm_;
     bool doCutInv_;
     HEPEnergyType energy_ = 0 * electronvolt;

examples/cascade_example.cpp

@@ -142,7 +142,7 @@ int main() {
   ParticleCut cut(80_GeV, true, true);
 
   TrackWriter trackWriter("tracks.dat");
-  BetheBlochPDG eLoss{showerAxis, cut.getECut()};
+  BetheBlochPDG eLoss{showerAxis, cut.getElectronECut()};
 
   // assemble all processes into an ordered process list
   auto sequence =

examples/cascade_proton_example.cpp

@@ -132,7 +132,7 @@ int main() {
 
   TrackWriter trackWriter("tracks.dat");
   ShowerAxis const showerAxis{injectionPos, Vector{rootCS, 0_m, 0_m, -100_km}, env};
-  BetheBlochPDG eLoss{showerAxis, cut.getECut()};
+  BetheBlochPDG eLoss{showerAxis, cut.getElectronECut()};
 
   // assemble all processes into an ordered process list
   // auto sequence = sibyll << decay << hadronicElastic << cut << trackWriter;

examples/em_shower.cpp

@@ -142,9 +142,9 @@ int main(int argc, char** argv) {
   // setup processes, decays and interactions
 
   // PROPOSAL processs proposal{...};
-  ParticleCut cut(10_GeV, false, true);
-  corsika::proposal::Interaction proposal(env, cut.getECut());
-  corsika::proposal::ContinuousProcess em_continuous(env, cut.getECut());
+  ParticleCut cut(10_GeV, 10_GeV, 100_PeV, 100_PeV, true);
+  corsika::proposal::Interaction proposal(env, cut.getElectronECut());
+  corsika::proposal::ContinuousProcess em_continuous(env, cut.getElectronECut());
   InteractionCounter proposalCounted(proposal);
 
   TrackWriter trackWriter("tracks.dat");

examples/hybrid_MC.cpp

@@ -207,7 +207,7 @@ int main(int argc, char** argv) {
   decaySibyll.printDecayConfig();
 
   ParticleCut cut{60_GeV, false, true};
-  BetheBlochPDG eLoss{showerAxis, cut.getECut()};
+  BetheBlochPDG eLoss{showerAxis, cut.getElectronECut()};
 
   CONEXhybrid conex_model(center, showerAxis, t, injectionHeight, E0,
                           get_PDG(Code::Proton));

examples/vertical_EAS.cpp

@@ -221,9 +221,9 @@ int main(int argc, char** argv) {
 
   decaySibyll.printDecayConfig();
 
-  ParticleCut cut{60_GeV, false, true};
-  corsika::proposal::Interaction proposal(env, cut.getECut());
-  corsika::proposal::ContinuousProcess em_continuous(env, cut.getECut());
+  ParticleCut cut{60_GeV, 60_GeV, 60_GeV, 60_GeV, true};
+  corsika::proposal::Interaction proposal(env, cut.getElectronECut());
+  corsika::proposal::ContinuousProcess em_continuous(env, cut.getElectronECut());
   InteractionCounter proposalCounted(proposal);
 
   OnShellCheck reset_particle_mass(1.e-3, 1.e-1, false);

tests/modules/testParticleCut.cpp

@@ -24,7 +24,7 @@ using namespace corsika;
 
 TEST_CASE("ParticleCut", "[processes]") {
 
-  logging::set_level(logging::level::info);
+  logging::set_level(logging::level::debug);
   corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
 
   feenableexcept(FE_INVALID);
@@ -51,7 +51,7 @@ TEST_CASE("ParticleCut", "[processes]") {
   SECTION("cut on particle type: inv") {
 
     ParticleCut cut(20_GeV, false, true);
-    CHECK(cut.getECut() == 20_GeV);
+    CHECK(cut.getHadronECut() == 20_GeV);
 
     // add primary particle to stack
     auto particle = stack.addParticle(std::make_tuple(
@@ -135,6 +135,43 @@ TEST_CASE("ParticleCut", "[processes]") {
     CHECK(view.getSize() == 13);
   }
 
+  SECTION("cut low energy: electrons, photons, hadrons and muons") {
+    ParticleCut cut(5_MeV,5_MeV,5_GeV,5_GeV, true);
+
+    // add primary particle to stack
+    auto particle = stack.addParticle(
+        std::make_tuple(Code::Proton, Eabove,
+                        MomentumVector(rootCS, {0_GeV, 0_GeV, 0_GeV}), point0, 0_ns));
+    // view on secondary particles
+    setup::StackView view(particle);
+    // ref. to primary particle through the secondary view.
+    // only this way the secondary view is populated
+    auto projectile = view.getProjectile();
+    // add secondaries 
+    projectile.addSecondary(std::make_tuple(Code::Gamma, 3_MeV,
+                                            MomentumVector(rootCS, {0_GeV, 0_GeV, 0_GeV}),
+                                            point0, 0_ns));
+    projectile.addSecondary(std::make_tuple(Code::Electron, 3_MeV,
+                                            MomentumVector(rootCS, {0_GeV, 0_GeV, 0_GeV}),
+                                            point0, 0_ns));
+    projectile.addSecondary(std::make_tuple(Code::PiPlus, 4_GeV,
+                                            MomentumVector(rootCS, {0_GeV, 0_GeV, 0_GeV}),
+                                            point0, 0_ns));
+    unsigned short A = 18;
+    unsigned short Z = 8;
+    projectile.addSecondary(std::make_tuple(Code::Nucleus, 4_GeV * A,
+                                            MomentumVector(rootCS, {0_GeV, 0_GeV, 0_GeV}),
+                                            point0, 0_ns, A, Z));
+    projectile.addSecondary(std::make_tuple(Code::Nucleus, 6_GeV * A,
+                                            MomentumVector(rootCS, {0_GeV, 0_GeV, 0_GeV}),
+                                            point0, 0_ns, A, Z));
+
+    cut.doSecondaries(view);
+
+    CHECK(view.getEntries() == 1);
+    CHECK(view.getSize() == 5);
+  }
+
   SECTION("cut on time") {
     ParticleCut cut(20_GeV, false, false);
     const TimeType too_late = 1_s;
@@ -170,7 +207,7 @@ TEST_CASE("ParticleCut", "[processes]") {
   SECTION("cut on DoContinous, just invisibles") {
 
     ParticleCut cut(20_GeV, false, true);
-    CHECK(cut.getECut() == 20_GeV);
+    CHECK(cut.getHadronECut() == 20_GeV);
 
     // add particles, all with energies above the threshold
     // only cut is by species