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