diff --git a/Documentation/Examples/CMakeLists.txt b/Documentation/Examples/CMakeLists.txt
index a3978ae9e04c23e61e02adc7ed25dd681c1bc796..c677f3ab9a3dccd82235fee8d4690be4de697d17 100644
--- a/Documentation/Examples/CMakeLists.txt
+++ b/Documentation/Examples/CMakeLists.txt
@@ -31,6 +31,7 @@ target_link_libraries (cascade_example SuperStupidStack CORSIKAunits
CORSIKAcascade
ProcessEnergyLoss
ProcessStackInspector
+ ProcessParticleCut
ProcessTrackWriter
ProcessTrackingLine
CORSIKAprocesses
@@ -50,6 +51,7 @@ target_link_libraries (boundary_example SuperStupidStack CORSIKAunits CORSIKAlog
ProcessSibyll
CORSIKAcascade
ProcessTrackWriter
+ ProcessParticleCut
ProcessTrackingLine
CORSIKAprocesses
CORSIKAparticles
@@ -72,6 +74,7 @@ if (Pythia8_FOUND)
ProcessEnergyLoss
ProcessTrackWriter
ProcessTrackingLine
+ ProcessParticleCut
ProcessHadronicElasticModel
ProcessStackInspector
CORSIKAprocesses
@@ -97,6 +100,7 @@ target_link_libraries (vertical_EAS SuperStupidStack CORSIKAunits
ProcessTrackWriter
ProcessTrackingLine
ProcessHadronicElasticModel
+ ProcessParticleCut
ProcessStackInspector
CORSIKAprocesses
CORSIKAcascade
diff --git a/Documentation/Examples/boundary_example.cc b/Documentation/Examples/boundary_example.cc
index 037c8b352d2431f83133e85f909d737a2e653d8e..169cb1c6a30a46a16494ac7e37b5bae2e875c256 100644
--- a/Documentation/Examples/boundary_example.cc
+++ b/Documentation/Examples/boundary_example.cc
@@ -29,6 +29,8 @@
#include <corsika/process/track_writer/TrackWriter.h>
+#include <corsika/process/particle_cut/ParticleCut.h>
+
#include <corsika/units/PhysicalUnits.h>
#include <corsika/random/RNGManager.h>
@@ -51,158 +53,6 @@ using namespace corsika::environment;
using namespace std;
using namespace corsika::units::si;
-class ProcessCut : public process::SecondariesProcess<ProcessCut> {
-
- HEPEnergyType fECut;
-
- HEPEnergyType fEnergy = 0_GeV;
- HEPEnergyType fEmEnergy = 0_GeV;
- int fEmCount = 0;
- HEPEnergyType fInvEnergy = 0_GeV;
- int fInvCount = 0;
-
-public:
- ProcessCut(const HEPEnergyType cut)
- : fECut(cut) {}
-
- template <typename Particle>
- bool isBelowEnergyCut(Particle& p) const {
- // nuclei
- if (p.GetPID() == particles::Code::Nucleus) {
- auto const ElabNuc = p.GetEnergy() / p.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. * p.GetEnergy() * 0.93827_GeV);
- if (p.GetEnergy() < fECut || Ecm < 10_GeV)
- return true;
- else
- return false;
- }
- }
-
- bool isEmParticle(Code pCode) const {
- bool is_em = false;
- // FOR NOW: switch
- switch (pCode) {
- 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;
- }
-
- void defineEmParticles() const {
- // create bool array identifying em particles
- }
-
- bool isInvisible(Code pCode) const {
- bool is_inv = false;
- // FOR NOW: switch
- switch (pCode) {
- 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 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 (isEmParticle(pid)) {
- cout << "removing em. particle..." << endl;
- fEmEnergy += energy;
- fEmCount += 1;
- p.Delete();
- } else if (isInvisible(pid)) {
- cout << "removing inv. particle..." << endl;
- fInvEnergy += energy;
- fInvCount += 1;
- p.Delete();
- } else if (isBelowEnergyCut(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 Init() {
- fEmEnergy = 0. * 1_GeV;
- fEmCount = 0;
- fInvEnergy = 0. * 1_GeV;
- fInvCount = 0;
- fEnergy = 0. * 1_GeV;
- // defineEmParticles();
- }
-
- void 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;
- }
-
- HEPEnergyType GetInvEnergy() const { return fInvEnergy; }
- HEPEnergyType GetCutEnergy() const { return fEnergy; }
- HEPEnergyType GetEmEnergy() const { return fEmEnergy; }
-};
template <bool deleteParticle>
struct MyBoundaryCrossingProcess
@@ -276,7 +126,8 @@ int main() {
process::sibyll::Decay decay{{particles::Code::PiPlus, particles::Code::PiMinus,
particles::Code::KPlus, particles::Code::KMinus,
particles::Code::K0Long, particles::Code::K0Short}};
- ProcessCut cut(20_GeV);
+
+ process::particle_cut::ParticleCut cut(20_GeV);
process::track_writer::TrackWriter trackWriter("tracks.dat");
MyBoundaryCrossingProcess<true> boundaryCrossing("crossings.dat");
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index a26f1af688fefea6d381be70c4f957b4d75215d8..7690d664811f4192d224044ab39ebea7b92d6136 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -31,6 +31,8 @@
#include <corsika/process/track_writer/TrackWriter.h>
+#include <corsika/process/particle_cut/ParticleCut.h>
+
#include <corsika/units/PhysicalUnits.h>
#include <corsika/random/RNGManager.h>
@@ -52,159 +54,6 @@ using namespace corsika::environment;
using namespace std;
using namespace corsika::units::si;
-class ProcessCut : public process::SecondariesProcess<ProcessCut> {
-
- HEPEnergyType fECut;
-
- HEPEnergyType fEnergy = 0_GeV;
- HEPEnergyType fEmEnergy = 0_GeV;
- int fEmCount = 0;
- HEPEnergyType fInvEnergy = 0_GeV;
- int fInvCount = 0;
-
-public:
- ProcessCut(const HEPEnergyType cut)
- : fECut(cut) {}
-
- template <typename Particle>
- bool isBelowEnergyCut(Particle& p) const {
- // nuclei
- if (p.GetPID() == particles::Code::Nucleus) {
- auto const ElabNuc = p.GetEnergy() / p.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. * p.GetEnergy() * 0.93827_GeV);
- if (p.GetEnergy() < fECut || Ecm < 10_GeV)
- return true;
- else
- return false;
- }
- }
-
- bool isEmParticle(Code pCode) const {
- bool is_em = false;
- // FOR NOW: switch
- switch (pCode) {
- 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;
- }
-
- void defineEmParticles() const {
- // create bool array identifying em particles
- }
-
- bool isInvisible(Code pCode) const {
- bool is_inv = false;
- // FOR NOW: switch
- switch (pCode) {
- 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 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 (isEmParticle(pid)) {
- cout << "removing em. particle..." << endl;
- fEmEnergy += energy;
- fEmCount += 1;
- p.Delete();
- } else if (isInvisible(pid)) {
- cout << "removing inv. particle..." << endl;
- fInvEnergy += energy;
- fInvCount += 1;
- p.Delete();
- } else if (isBelowEnergyCut(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 Init() {
- fEmEnergy = 0. * 1_GeV;
- fEmCount = 0;
- fInvEnergy = 0. * 1_GeV;
- fInvCount = 0;
- fEnergy = 0. * 1_GeV;
- // defineEmParticles();
- }
-
- void 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;
- }
-
- HEPEnergyType GetInvEnergy() const { return fInvEnergy; }
- HEPEnergyType GetCutEnergy() const { return fEnergy; }
- HEPEnergyType GetEmEnergy() const { return fEmEnergy; }
-};
-
//
// The example main program for a particle cascade
//
@@ -257,7 +106,7 @@ int main() {
process::sibyll::Interaction sibyll;
process::sibyll::NuclearInteraction sibyllNuc(sibyll, env);
process::sibyll::Decay decay(trackedHadrons);
- ProcessCut cut(20_GeV);
+ process::particle_cut::ParticleCut cut(20_GeV);
process::track_writer::TrackWriter trackWriter("tracks.dat");
process::energy_loss::EnergyLoss eLoss;
diff --git a/Documentation/Examples/cascade_proton_example.cc b/Documentation/Examples/cascade_proton_example.cc
index e5b4b3456c42c8e03f28847e41d0babdee9ab464..4e47bf5cecc4c12e04d052655cda951086dea8c3 100644
--- a/Documentation/Examples/cascade_proton_example.cc
+++ b/Documentation/Examples/cascade_proton_example.cc
@@ -32,6 +32,8 @@
#include <corsika/process/track_writer/TrackWriter.h>
+#include <corsika/process/particle_cut/ParticleCut.h>
+
#include <corsika/units/PhysicalUnits.h>
#include <corsika/random/RNGManager.h>
@@ -57,158 +59,6 @@ using namespace corsika::environment;
using namespace std;
using namespace corsika::units::si;
-class ProcessCut : public process::SecondariesProcess<ProcessCut> {
-
- HEPEnergyType fECut;
-
- HEPEnergyType fEnergy = 0_GeV;
- HEPEnergyType fEmEnergy = 0_GeV;
- int fEmCount = 0;
- HEPEnergyType fInvEnergy = 0_GeV;
- int fInvCount = 0;
-
-public:
- ProcessCut(const HEPEnergyType cut)
- : fECut(cut) {}
-
- template <typename Particle>
- bool isBelowEnergyCut(Particle& p) const {
- // nuclei
- if (p.GetPID() == particles::Code::Nucleus) {
- auto const ElabNuc = p.GetEnergy() / p.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. * p.GetEnergy() * 0.93827_GeV);
- if (p.GetEnergy() < fECut || Ecm < 10_GeV)
- return true;
- else
- return false;
- }
- }
-
- bool isEmParticle(Code pCode) const {
- bool is_em = false;
- // FOR NOW: switch
- switch (pCode) {
- 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;
- }
-
- void defineEmParticles() const {
- // create bool array identifying em particles
- }
-
- bool isInvisible(Code pCode) const {
- bool is_inv = false;
- // FOR NOW: switch
- switch (pCode) {
- 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 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 (isEmParticle(pid)) {
- cout << "removing em. particle..." << endl;
- fEmEnergy += energy;
- fEmCount += 1;
- p.Delete();
- } else if (isInvisible(pid)) {
- cout << "removing inv. particle..." << endl;
- fInvEnergy += energy;
- fInvCount += 1;
- p.Delete();
- } else if (isBelowEnergyCut(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 Init() {
- fEmEnergy = 0. * 1_GeV;
- fEmCount = 0;
- fInvEnergy = 0. * 1_GeV;
- fInvCount = 0;
- fEnergy = 0. * 1_GeV;
- // defineEmParticles();
- }
-
- void 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;
- }
-
- HEPEnergyType GetInvEnergy() const { return fInvEnergy; }
- HEPEnergyType GetCutEnergy() const { return fEnergy; }
- HEPEnergyType GetEmEnergy() const { return fEmEnergy; }
-};
//
// The example main program for a particle cascade
@@ -251,7 +101,7 @@ int main() {
// process::sibyll::NuclearInteraction sibyllNuc(env, sibyll);
// process::sibyll::Decay decay(trackedHadrons);
process::pythia::Decay decay(trackedHadrons);
- ProcessCut cut(20_GeV);
+ process::particle_cut::ParticleCut cut(20_GeV);
// random::RNGManager::GetInstance().RegisterRandomStream("HadronicElasticModel");
// process::HadronicElasticModel::HadronicElasticInteraction
diff --git a/Documentation/Examples/vertical_EAS.cc b/Documentation/Examples/vertical_EAS.cc
index 00f38d8b56f4082669b6ba6b995d14d2f5e135cb..872ea2d943bb808edeb295b767cc931dc3e49032 100644
--- a/Documentation/Examples/vertical_EAS.cc
+++ b/Documentation/Examples/vertical_EAS.cc
@@ -28,6 +28,7 @@
#include <corsika/process/sibyll/Interaction.h>
#include <corsika/process/sibyll/NuclearInteraction.h>
+#include <corsika/process/particle_cut/ParticleCut.h>
#include <corsika/process/track_writer/TrackWriter.h>
#include <corsika/units/PhysicalUnits.h>
@@ -55,159 +56,6 @@ using namespace corsika::environment;
using namespace std;
using namespace corsika::units::si;
-class ProcessCut : public process::SecondariesProcess<ProcessCut> {
-
- HEPEnergyType fECut;
-
- HEPEnergyType fEnergy = 0_GeV;
- HEPEnergyType fEmEnergy = 0_GeV;
- int fEmCount = 0;
- HEPEnergyType fInvEnergy = 0_GeV;
- int fInvCount = 0;
-
-public:
- ProcessCut(const HEPEnergyType cut)
- : fECut(cut) {}
-
- template <typename Particle>
- bool isBelowEnergyCut(Particle& p) const {
- // nuclei
- if (p.GetPID() == particles::Code::Nucleus) {
- auto const ElabNuc = p.GetEnergy() / p.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. * p.GetEnergy() * 0.93827_GeV);
- if (p.GetEnergy() < fECut || Ecm < 10_GeV)
- return true;
- else
- return false;
- }
- }
-
- bool isEmParticle(Code pCode) const {
- bool is_em = false;
- // FOR NOW: switch
- switch (pCode) {
- 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;
- }
-
- void defineEmParticles() const {
- // create bool array identifying em particles
- }
-
- bool isInvisible(Code pCode) const {
- bool is_inv = false;
- // FOR NOW: switch
- switch (pCode) {
- 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 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 (isEmParticle(pid)) {
- cout << "removing em. particle..." << endl;
- fEmEnergy += energy;
- fEmCount += 1;
- p.Delete();
- } else if (isInvisible(pid)) {
- cout << "removing inv. particle..." << endl;
- fInvEnergy += energy;
- fInvCount += 1;
- p.Delete();
- } else if (isBelowEnergyCut(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 Init() {
- fEmEnergy = 0. * 1_GeV;
- fEmCount = 0;
- fInvEnergy = 0. * 1_GeV;
- fInvCount = 0;
- fEnergy = 0. * 1_GeV;
- // defineEmParticles();
- }
-
- void 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;
- }
-
- HEPEnergyType GetInvEnergy() const { return fInvEnergy; }
- HEPEnergyType GetCutEnergy() const { return fEnergy; }
- HEPEnergyType GetEmEnergy() const { return fEmEnergy; }
-};
-
//
// The example main program for a particle cascade
//
@@ -250,8 +98,7 @@ int main() {
process::sibyll::Interaction sibyll;
process::sibyll::NuclearInteraction sibyllNuc(sibyll, env);
process::sibyll::Decay decay(trackedHadrons);
- ProcessCut cut(20_GeV);
-
+ process::particle_cut::ParticleCut cut(20_GeV);
// random::RNGManager::GetInstance().RegisterRandomStream("HadronicElasticModel");
// process::HadronicElasticModel::HadronicElasticInteraction
// hadronicElastic(env);
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..28e7013d08d864471de08c742250b3ecb4357760
--- /dev/null
+++ b/Processes/ParticleCut/ParticleCut.cc
@@ -0,0 +1,140 @@
+
+/*
+ * (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 TParticle>
+ bool ParticleCut::ParticleIsBelowEnergyCut(TParticle const& vP) const {
+ auto const energyLab = vP.GetEnergy();
+ // nuclei
+ if (vP.GetPID() == particles::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 {
+ 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;
+ }
+
+ EProcessReturn ParticleCut::DoSecondaries(corsika::setup::StackView& 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._GeV;
+ fEmCount = 0;
+ fInvEnergy = 0._GeV;
+ fInvCount = 0;
+ fEnergy = 0._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..6d5f7a685cccbce7f2e079a3c291f90dbf20039b
--- /dev/null
+++ b/Processes/ParticleCut/ParticleCut.h
@@ -0,0 +1,53 @@
+/*
+ * (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/setup/SetupStack.h>
+#include <corsika/units/PhysicalUnits.h>
+
+namespace corsika::process {
+ namespace particle_cut {
+ class ParticleCut : public process::SecondariesProcess<ParticleCut> {
+
+ units::si::HEPEnergyType const fECut;
+
+ units::si::HEPEnergyType fEnergy = 0 * units::si::electronvolt;
+ units::si::HEPEnergyType fEmEnergy = 0 * units::si::electronvolt;
+ int fEmCount = 0;
+ units::si::HEPEnergyType fInvEnergy = 0 * units::si::electronvolt;
+ int fInvCount = 0;
+
+ public:
+ ParticleCut(const units::si::HEPEnergyType vCut)
+ : fECut(vCut) {}
+
+ bool ParticleIsInvisible(particles::Code) const;
+ EProcessReturn DoSecondaries(corsika::setup::StackView&);
+
+ template <typename TParticle>
+ bool ParticleIsBelowEnergyCut(TParticle const&) const;
+
+ bool ParticleIsEmParticle(particles::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; }
+ };
+ } // namespace particle_cut
+} // namespace corsika::process
+
+#endif
diff --git a/Processes/ParticleCut/testParticleCut.cc b/Processes/ParticleCut/testParticleCut.cc
new file mode 100644
index 0000000000000000000000000000000000000000..c1a69696f06a331289ce20ae901fc3329d46fb34
--- /dev/null
+++ b/Processes/ParticleCut/testParticleCut.cc
@@ -0,0 +1,109 @@
+
+/*
+ * (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();
+
+ // setup empty particle stack
+ setup::Stack stack;
+ stack.Clear();
+ // two energies
+ 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};
+
+ SECTION("cut invisible") {
+
+ ParticleCut cut(20_GeV);
+
+ // 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});
+ // view on secondary particles
+ corsika::stack::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,
+ 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});
+
+ cut.DoSecondaries(view);
+
+ REQUIRE(view.GetSize() == 6);
+ }
+
+ SECTION("cut low energy") {
+ ParticleCut cut(20_GeV);
+
+ // 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});
+ // view on secondary particles
+ corsika::stack::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,
+ 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});
+
+ cut.DoSecondaries(view);
+
+ REQUIRE(view.GetSize() == 0);
+ }
+}