From 1141920118c7693156bbfc53d0e62383e820f373 Mon Sep 17 00:00:00 2001
From: "felix@home" <felix.riehn@kit.edu>
Date: Mon, 13 May 2019 22:42:09 +0100
Subject: [PATCH] use ParticleCut in Examples
---
Documentation/Examples/CMakeLists.txt | 4 +
Documentation/Examples/boundary_example.cc | 157 +----------------
Documentation/Examples/cascade_example.cc | 158 +-----------------
.../Examples/cascade_proton_example.cc | 156 +----------------
Documentation/Examples/vertical_EAS.cc | 158 +-----------------
5 files changed, 17 insertions(+), 616 deletions(-)
diff --git a/Documentation/Examples/CMakeLists.txt b/Documentation/Examples/CMakeLists.txt
index 104a1513d..f9c631669 100644
--- a/Documentation/Examples/CMakeLists.txt
+++ b/Documentation/Examples/CMakeLists.txt
@@ -31,6 +31,7 @@ if (Pythia8_FOUND)
ProcessSibyll
ProcessPythia
CORSIKAcascade
+ ProcessParticleCut
ProcessEnergyLoss
ProcessStackInspector
ProcessTrackWriter
@@ -53,6 +54,7 @@ target_link_libraries (boundary_example SuperStupidStack CORSIKAunits CORSIKAlog
ProcessSibyll
CORSIKAcascade
ProcessTrackWriter
+ ProcessParticleCut
ProcessTrackingLine
CORSIKAprocesses
CORSIKAparticles
@@ -76,6 +78,7 @@ if (Pythia8_FOUND)
ProcessEnergyLoss
ProcessTrackWriter
ProcessTrackingLine
+ ProcessParticleCut
ProcessHadronicElasticModel
CORSIKAprocesses
CORSIKAcascade
@@ -101,6 +104,7 @@ if (Pythia8_FOUND)
ProcessEnergyLoss
ProcessTrackWriter
ProcessTrackingLine
+ ProcessParticleCut
ProcessHadronicElasticModel
CORSIKAprocesses
CORSIKAcascade
diff --git a/Documentation/Examples/boundary_example.cc b/Documentation/Examples/boundary_example.cc
index 037c8b352..169cb1c6a 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 cb06341c8..400bb07d8 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -33,6 +33,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>
@@ -54,160 +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
//
@@ -261,7 +109,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 e5b4b3456..4e47bf5ce 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 526498b24..f96544ffb 100644
--- a/Documentation/Examples/vertical_EAS.cc
+++ b/Documentation/Examples/vertical_EAS.cc
@@ -30,6 +30,8 @@
#include <corsika/process/pythia/Decay.h>
+#include <corsika/process/particle_cut/ParticleCut.h>
+
#include <corsika/process/track_writer/TrackWriter.h>
#include <corsika/units/PhysicalUnits.h>
@@ -57,159 +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
@@ -255,8 +104,7 @@ int main() {
process::sibyll::Decay decay(trackedHadrons);
// random::RNGManager::GetInstance().RegisterRandomStream("pythia");
// process::pythia::Decay decay(trackedHadrons);
- ProcessCut cut(20_GeV);
-
+ process::particle_cut::ParticleCut cut(20_GeV);
// random::RNGManager::GetInstance().RegisterRandomStream("HadronicElasticModel");
// process::HadronicElasticModel::HadronicElasticInteraction
// hadronicElastic(env);
--
GitLab