From da7bee3a07fdc9884e58e2a3580b6a547d150065 Mon Sep 17 00:00:00 2001
From: ralfulrich <ralf.ulrich@kit.edu>
Date: Tue, 2 Apr 2019 11:33:59 +0200
Subject: [PATCH] use new SecondariesProcess, added ObservationLevel (passive
version) for tests
---
Documentation/Examples/CMakeLists.txt | 25 +++-
Documentation/Examples/cascade_example.cc | 108 ++++++++++--------
.../Examples/staticsequence_example.cc | 7 +-
Documentation/Examples/vertical_EAS.cc | 58 +++++-----
4 files changed, 118 insertions(+), 80 deletions(-)
diff --git a/Documentation/Examples/CMakeLists.txt b/Documentation/Examples/CMakeLists.txt
index 06b4467e4..4d3c5fd75 100644
--- a/Documentation/Examples/CMakeLists.txt
+++ b/Documentation/Examples/CMakeLists.txt
@@ -31,7 +31,7 @@ target_link_libraries (cascade_example SuperStupidStack CORSIKAunits
ProcessPythia
CORSIKAcascade
ProcessEnergyLoss
- ProcessStackInspector
+# ProcessStackInspector
ProcessTrackWriter
ProcessTrackingLine
ProcessHadronicElasticModel
@@ -45,6 +45,29 @@ target_link_libraries (cascade_example SuperStupidStack CORSIKAunits
install (TARGETS cascade_example DESTINATION share/examples)
CORSIKA_ADD_TEST (cascade_example)
+add_executable (vertical_EAS vertical_EAS.cc)
+target_compile_options(vertical_EAS PRIVATE -g) # do not skip asserts
+target_link_libraries (vertical_EAS SuperStupidStack CORSIKAunits
+ CORSIKAlogging
+ CORSIKArandom
+ ProcessSibyll
+ ProcessPythia
+ CORSIKAcascade
+ ProcessEnergyLoss
+# ProcessStackInspector
+ ProcessTrackWriter
+ ProcessTrackingLine
+ ProcessHadronicElasticModel
+ CORSIKAprocesses
+ CORSIKAcascade
+ CORSIKAparticles
+ CORSIKAgeometry
+ CORSIKAenvironment
+ CORSIKAprocesssequence
+ )
+install (TARGETS vertical_EAS DESTINATION share/examples)
+CORSIKA_ADD_TEST (vertical_EAS)
+
add_executable (staticsequence_example staticsequence_example.cc)
target_compile_options(staticsequence_example PRIVATE -g) # do not skip asserts
target_link_libraries (staticsequence_example
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index 764eeeba0..6c4036deb 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -58,7 +58,7 @@ using namespace corsika::environment;
using namespace std;
using namespace corsika::units::si;
-class ProcessCut : public process::ContinuousProcess<ProcessCut> {
+class ProcessCut : public process::SecondariesProcess<ProcessCut> {
HEPEnergyType fECut;
@@ -152,47 +152,38 @@ public:
return is_inv;
}
- template <typename Particle>
- LengthType MaxStepLength(Particle& p, setup::Trajectory&) const {
- cout << "ProcessCut: MinStep: pid: " << p.GetPID() << endl;
- cout << "ProcessCut: MinStep: energy (GeV): " << p.GetEnergy() / 1_GeV << endl;
- const Code pid = p.GetPID();
- if (isEmParticle(pid) || isInvisible(pid) || isBelowEnergyCut(p)) {
- cout << "ProcessCut: MinStep: next cut: " << 0. << endl;
- return 0_m;
- } else {
- LengthType next_step = 1_m * std::numeric_limits<double>::infinity();
- cout << "ProcessCut: MinStep: next cut: " << next_step << endl;
- return next_step;
- }
- }
-
- template <typename Particle, typename Stack>
- EProcessReturn DoContinuous(Particle& p, setup::Trajectory&, Stack&) {
- const Code pid = p.GetPID();
- HEPEnergyType energy = p.GetEnergy();
- cout << "ProcessCut: DoContinuous: " << pid << " E= " << energy
- << ", EcutTot=" << (fEmEnergy + fInvEnergy + fEnergy) / 1_GeV << " GeV" << endl;
- EProcessReturn ret = EProcessReturn::eOk;
- if (isEmParticle(pid)) {
- cout << "removing em. particle..." << endl;
- fEmEnergy += energy;
- fEmCount += 1;
- // p.Delete();
- ret = EProcessReturn::eParticleAbsorbed;
- } else if (isInvisible(pid)) {
- cout << "removing inv. particle..." << endl;
- fInvEnergy += energy;
- fInvCount += 1;
- // p.Delete();
- ret = EProcessReturn::eParticleAbsorbed;
- } else if (isBelowEnergyCut(p)) {
- cout << "removing low en. particle..." << endl;
- fEnergy += energy;
- // p.Delete();
- ret = EProcessReturn::eParticleAbsorbed;
+ 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 ret;
+ return EProcessReturn::eOk;
}
void Init() {
@@ -220,6 +211,31 @@ public:
HEPEnergyType GetEmEnergy() const { return fEmEnergy; }
};
+class ObservationLevel : public process::ContinuousProcess<ObservationLevel> {
+
+ LengthType fHeight;
+
+public:
+ ObservationLevel(const LengthType vHeight)
+ : fHeight(vHeight) {}
+
+ template <typename Particle>
+ LengthType MaxStepLength(Particle&, setup::Trajectory&) const {
+ return 1_m * std::numeric_limits<double>::infinity();
+ }
+
+ template <typename TParticle, typename TTrack>
+ EProcessReturn DoContinuous(TParticle&, TTrack& vT) {
+ if ((vT.GetPosition(0).GetZ() <= fHeight && vT.GetPosition(1).GetZ() > fHeight) ||
+ (vT.GetPosition(0).GetZ() > fHeight && vT.GetPosition(1).GetZ() <= fHeight)) {
+ cout << "OBSERVED " << endl;
+ return EProcessReturn::eParticleAbsorbed;
+ }
+ return EProcessReturn::eOk;
+ }
+ void Init() {}
+};
+
//
// The example main program for a particle cascade
//
@@ -252,7 +268,7 @@ int main() {
// setup processes, decays and interactions
tracking_line::TrackingLine<setup::Stack, setup::Trajectory> tracking(env);
- stack_inspector::StackInspector<setup::Stack> stackInspect(true);
+ // stack_inspector::StackInspector<setup::Stack> stackInspect(true);
const std::vector<particles::Code> trackedHadrons = {
particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus,
@@ -265,6 +281,7 @@ int main() {
// random::RNGManager::GetInstance().RegisterRandomStream("pythia");
// process::pythia::Decay decay(trackedHadrons);
ProcessCut cut(20_GeV);
+ ObservationLevel obsLevel(1400_m);
// random::RNGManager::GetInstance().RegisterRandomStream("HadronicElasticModel");
// process::HadronicElasticModel::HadronicElasticInteraction
@@ -277,8 +294,8 @@ int main() {
// auto sequence = stackInspect << sibyll << decay << hadronicElastic << cut <<
// trackWriter; auto sequence = stackInspect << sibyll << sibyllNuc << decay << eLoss <<
// cut << trackWriter;
- // auto sequence = sibyll << sibyllNuc << decay << eLoss << cut;
- auto sequence = stackInspect << sibyll << sibyllNuc << decay << eLoss << cut;
+ auto sequence = sibyll << sibyllNuc << decay << eLoss << cut << obsLevel;
+ // auto sequence = stackInspect << sibyll << sibyllNuc << decay << eLoss << cut;
// cout << "decltype(sequence)=" << type_id_with_cvr<decltype(sequence)>().pretty_name()
// << "\n";
@@ -309,7 +326,8 @@ int main() {
cout << "input particle: " << beamCode << endl;
cout << "input angles: theta=" << theta << " phi=" << phi << endl;
cout << "input momentum: " << plab.GetComponents() / 1_GeV << endl;
- Point pos(rootCS, 0_m, 0_m, 0_m);
+ Point pos(rootCS, 0_m, 0_m,
+ 112.8_km); // this is the CORSIKA 7 start of atmosphere/universe
stack.AddParticle(std::tuple<particles::Code, units::si::HEPEnergyType,
corsika::stack::MomentumVector, geometry::Point,
units::si::TimeType, unsigned short, unsigned short>{
diff --git a/Documentation/Examples/staticsequence_example.cc b/Documentation/Examples/staticsequence_example.cc
index 5e881ebbe..0cd2901c3 100644
--- a/Documentation/Examples/staticsequence_example.cc
+++ b/Documentation/Examples/staticsequence_example.cc
@@ -88,12 +88,11 @@ void modular() {
auto sequence = m1 << m2 << m3 << m4;
- DummyData p;
- DummyStack s;
- DummyTrajectory t;
+ DummyData particle;
+ DummyTrajectory track;
const int n = 1000;
- for (int i = 0; i < n; ++i) { sequence.DoContinuous(p, t, s); }
+ for (int i = 0; i < n; ++i) { sequence.DoContinuous(particle, track); }
for (int i = 0; i < nData; ++i) {
// cout << p.p[i] << endl;
diff --git a/Documentation/Examples/vertical_EAS.cc b/Documentation/Examples/vertical_EAS.cc
index d3bee8145..6c4036deb 100644
--- a/Documentation/Examples/vertical_EAS.cc
+++ b/Documentation/Examples/vertical_EAS.cc
@@ -159,27 +159,28 @@ public:
const Code pid = p.GetPID();
HEPEnergyType energy = p.GetEnergy();
cout << "ProcessCut: DoSecondaries: " << pid << " E= " << energy
- << ", EcutTot=" << (fEmEnergy + fInvEnergy + fEnergy) / 1_GeV << " GeV" << endl;
+ << ", EcutTot=" << (fEmEnergy + fInvEnergy + fEnergy) / 1_GeV << " GeV"
+ << endl;
if (isEmParticle(pid)) {
- cout << "removing em. particle..." << endl;
- fEmEnergy += energy;
- fEmCount += 1;
- p.Delete();
+ 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();
+ 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();
+ 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
+ ++p; // next entry in SecondaryView
}
}
return EProcessReturn::eOk;
@@ -210,33 +211,29 @@ public:
HEPEnergyType GetEmEnergy() const { return fEmEnergy; }
};
-
class ObservationLevel : public process::ContinuousProcess<ObservationLevel> {
LengthType fHeight;
-
+
public:
ObservationLevel(const LengthType vHeight)
- : fHeight(vHeight) {}
+ : fHeight(vHeight) {}
template <typename Particle>
LengthType MaxStepLength(Particle&, setup::Trajectory&) const {
return 1_m * std::numeric_limits<double>::infinity();
}
-
+
template <typename TParticle, typename TTrack>
EProcessReturn DoContinuous(TParticle&, TTrack& vT) {
- if ((vT.GetPosition(0).GetZ()<=fHeight &&
- vT.GetPosition(1).GetZ()>fHeight) ||
- (vT.GetPosition(0).GetZ()>fHeight &&
- vT.GetPosition(1).GetZ()<=fHeight)) {
- cout << "OBSERVED " << endl;
- return EProcessReturn::eParticleAbsorbed;
+ if ((vT.GetPosition(0).GetZ() <= fHeight && vT.GetPosition(1).GetZ() > fHeight) ||
+ (vT.GetPosition(0).GetZ() > fHeight && vT.GetPosition(1).GetZ() <= fHeight)) {
+ cout << "OBSERVED " << endl;
+ return EProcessReturn::eParticleAbsorbed;
}
return EProcessReturn::eOk;
}
void Init() {}
-
};
//
@@ -285,7 +282,7 @@ int main() {
// process::pythia::Decay decay(trackedHadrons);
ProcessCut cut(20_GeV);
ObservationLevel obsLevel(1400_m);
-
+
// random::RNGManager::GetInstance().RegisterRandomStream("HadronicElasticModel");
// process::HadronicElasticModel::HadronicElasticInteraction
// hadronicElastic(env);
@@ -329,7 +326,8 @@ int main() {
cout << "input particle: " << beamCode << endl;
cout << "input angles: theta=" << theta << " phi=" << phi << endl;
cout << "input momentum: " << plab.GetComponents() / 1_GeV << endl;
- Point pos(rootCS, 0_m, 0_m, 112.8_km); // this is the CORSIKA 7 start of atmosphere/universe
+ Point pos(rootCS, 0_m, 0_m,
+ 112.8_km); // this is the CORSIKA 7 start of atmosphere/universe
stack.AddParticle(std::tuple<particles::Code, units::si::HEPEnergyType,
corsika::stack::MomentumVector, geometry::Point,
units::si::TimeType, unsigned short, unsigned short>{
--
GitLab