diff --git a/Documentation/Examples/CMakeLists.txt b/Documentation/Examples/CMakeLists.txt
index 6529b5153565a74ddc39fc16b3e7c90f882f7aac..ffc14426f790c8c760d77fdc056d782836bc1a94 100644
--- a/Documentation/Examples/CMakeLists.txt
+++ b/Documentation/Examples/CMakeLists.txt
@@ -47,8 +47,8 @@ CORSIKA_ADD_TEST (cascade_example)
add_executable (boundary_example boundary_example.cc)
target_compile_options(boundary_example PRIVATE -g) # do not skip asserts
target_link_libraries (boundary_example SuperStupidStack CORSIKAunits CORSIKAlogging
- CORSIKArandom
- ProcessSibyll
+ CORSIKArandom
+ ProcessSibyll
CORSIKAcascade
ProcessTrackWriter
ProcessTrackingLine
@@ -70,11 +70,12 @@ target_link_libraries (cascade_proton_example SuperStupidStack CORSIKAunits
ProcessSibyll
ProcessPythia
CORSIKAcascade
- ProcessStackInspector
+ ProcessEnergyLoss
ProcessTrackWriter
ProcessTrackingLine
ProcessHadronicElasticModel
CORSIKAprocesses
+ CORSIKAcascade
CORSIKAparticles
CORSIKAgeometry
CORSIKAenvironment
@@ -84,6 +85,28 @@ target_link_libraries (cascade_proton_example SuperStupidStack CORSIKAunits
install (TARGETS cascade_proton_example DESTINATION share/examples)
CORSIKA_ADD_TEST (cascade_proton_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
+ ProcessTrackWriter
+ ProcessTrackingLine
+ ProcessHadronicElasticModel
+ CORSIKAprocesses
+ CORSIKAcascade
+ CORSIKAparticles
+ CORSIKAgeometry
+ CORSIKAenvironment
+ CORSIKAprocesssequence
+ )
+install (TARGETS vertical_EAS DESTINATION share/examples)
+# CORSIKA_ADD_TEST (vertical_EAS) not yet...
+
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/boundary_example.cc b/Documentation/Examples/boundary_example.cc
index 98e01b89f0aaeaa556e8d0364b2be3fff4498ecc..b1d540eebb540c10063ef79648086f749d779ba9 100644
--- a/Documentation/Examples/boundary_example.cc
+++ b/Documentation/Examples/boundary_example.cc
@@ -51,7 +51,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;
@@ -145,47 +145,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() {
@@ -251,7 +242,7 @@ int main() {
random::RNGManager::GetInstance().RegisterRandomStream("cascade");
// setup environment, geometry
- using EnvType = environment::Environment<setup::IEnvironmentModel>;
+ using EnvType = Environment<setup::IEnvironmentModel>;
EnvType env;
auto& universe = *(env.GetUniverse());
@@ -266,7 +257,7 @@ int main() {
->SetModelProperties<environment::HomogeneousMedium<setup::IEnvironmentModel>>(
1_kg / (1_m * 1_m * 1_m),
environment::NuclearComposition(
- std::vector<particles::Code>{particles::Code::Oxygen},
+ std::vector<particles::Code>{particles::Code::Proton},
std::vector<float>{1.f}));
auto innerMedium = EnvType::CreateNode<Sphere>(Point{rootCS, 0_m, 0_m, 0_m}, 5_km);
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index 9855dfcb8bfa227c1b5da8487d73fe0d64b4941f..03feb8f6eb88b8d5f023bbcaec7ede771e9fadfd 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -41,7 +41,6 @@
#include <iostream>
#include <limits>
-#include <typeinfo>
using namespace corsika;
using namespace corsika::process;
@@ -55,7 +54,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;
@@ -149,47 +148,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() {
@@ -217,10 +207,14 @@ public:
HEPEnergyType GetEmEnergy() const { return fEmEnergy; }
};
+
//
// The example main program for a particle cascade
//
int main() {
+
+ const LengthType height_atmosphere = 112.8_km;
+
feenableexcept(FE_INVALID);
// initialize random number sequence(s)
random::RNGManager::GetInstance().RegisterRandomStream("cascade");
@@ -256,7 +250,7 @@ int main() {
// setup processes, decays and interactions
tracking_line::TrackingLine tracking;
- stack_inspector::StackInspector<setup::Stack> stackInspect(true);
+ stack_inspector::StackInspector<setup::Stack> stackInspect(1, true);
const std::vector<particles::Code> trackedHadrons = {
particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus,
@@ -273,7 +267,8 @@ int main() {
process::EnergyLoss::EnergyLoss eLoss;
// assemble all processes into an ordered process list
- auto sequence = sibyll << sibyllNuc << decay << cut << trackWriter;
+ auto sequence = stackInspect << sibyll << sibyllNuc << decay << eLoss << cut
+ << trackWriter;
// setup particle stack, and add primary particle
setup::Stack stack;
@@ -282,7 +277,7 @@ int main() {
const int nuclA = 4;
const int nuclZ = int(nuclA / 2.15 + 0.7);
const HEPMassType mass = GetNucleusMass(nuclA, nuclZ);
- const HEPEnergyType E0 = nuclA * 10_TeV;
+ const HEPEnergyType E0 = nuclA * 1_TeV;
double theta = 0.;
double phi = 0.;
@@ -301,7 +296,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,
+ height_atmosphere); // 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/cascade_proton_example.cc b/Documentation/Examples/cascade_proton_example.cc
index 97e6c40ac5f4b80fd2a64430eec8b19fe81867ed..743f8f007b0fbfd594cb209d1d580444a25aba21 100644
--- a/Documentation/Examples/cascade_proton_example.cc
+++ b/Documentation/Examples/cascade_proton_example.cc
@@ -12,7 +12,6 @@
#include <corsika/cascade/Cascade.h>
#include <corsika/process/ProcessSequence.h>
#include <corsika/process/hadronic_elastic_model/HadronicElasticModel.h>
-#include <corsika/process/stack_inspector/StackInspector.h>
#include <corsika/process/tracking_line/TrackingLine.h>
#include <corsika/setup/SetupStack.h>
@@ -58,7 +57,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 +151,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() {
@@ -229,7 +219,7 @@ int main() {
random::RNGManager::GetInstance().RegisterRandomStream("cascade");
// setup environment, geometry
- using EnvType = environment::Environment<setup::IEnvironmentModel>;
+ using EnvType = Environment<setup::IEnvironmentModel>;
EnvType env;
auto& universe = *(env.GetUniverse());
@@ -237,12 +227,11 @@ int main() {
EnvType::CreateNode<Sphere>(Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m},
1_km * std::numeric_limits<double>::infinity());
- using MyHomogeneousModel = environment::HomogeneousMedium<environment::IMediumModel>;
+ using MyHomogeneousModel = HomogeneousMedium<IMediumModel>;
theMedium->SetModelProperties<MyHomogeneousModel>(
1_kg / (1_m * 1_m * 1_m),
- environment::NuclearComposition(
- std::vector<particles::Code>{particles::Code::Hydrogen},
- std::vector<float>{(float)1.}));
+ NuclearComposition(std::vector<particles::Code>{particles::Code::Hydrogen},
+ std::vector<float>{(float)1.}));
universe.AddChild(std::move(theMedium));
@@ -250,7 +239,6 @@ int main() {
// setup processes, decays and interactions
tracking_line::TrackingLine tracking;
- stack_inspector::StackInspector<setup::Stack> stackInspect(true);
const std::vector<particles::Code> trackedHadrons = {
particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus,
@@ -272,10 +260,8 @@ int main() {
process::TrackWriter::TrackWriter trackWriter("tracks.dat");
// assemble all processes into an ordered process list
- // auto sequence = p0 << sibyll << decay << hadronicElastic << cut << trackWriter;
- // auto sequence = p0 << sibyll << sibyllNuc << decay << cut << trackWriter;
-
- auto sequence = stackInspect << pythia << decay << cut << trackWriter;
+ // auto sequence = sibyll << decay << hadronicElastic << cut << trackWriter;
+ auto sequence = pythia << decay << cut << trackWriter;
// cout << "decltype(sequence)=" << type_id_with_cvr<decltype(sequence)>().pretty_name()
// << "\n";
diff --git a/Documentation/Examples/staticsequence_example.cc b/Documentation/Examples/staticsequence_example.cc
index 5e881ebbe3e38ae028d2c263ff27bcc41b0490b5..0cd2901c3520d203371a7103f3dff96243328f9f 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
new file mode 100644
index 0000000000000000000000000000000000000000..e2d381914668b83431cccac1c2a51caa117ed22f
--- /dev/null
+++ b/Documentation/Examples/vertical_EAS.cc
@@ -0,0 +1,322 @@
+
+/*
+ * (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/cascade/Cascade.h>
+#include <corsika/process/ProcessSequence.h>
+#include <corsika/process/energy_loss/EnergyLoss.h>
+#include <corsika/process/hadronic_elastic_model/HadronicElasticModel.h>
+#include <corsika/process/tracking_line/TrackingLine.h>
+
+#include <corsika/setup/SetupStack.h>
+#include <corsika/setup/SetupTrajectory.h>
+
+#include <corsika/environment/Environment.h>
+#include <corsika/environment/HomogeneousMedium.h>
+#include <corsika/environment/NuclearComposition.h>
+
+#include <corsika/geometry/Sphere.h>
+
+#include <corsika/process/sibyll/Decay.h>
+#include <corsika/process/sibyll/Interaction.h>
+#include <corsika/process/sibyll/NuclearInteraction.h>
+
+#include <corsika/process/pythia/Decay.h>
+
+#include <corsika/process/track_writer/TrackWriter.h>
+
+#include <corsika/units/PhysicalUnits.h>
+
+#include <corsika/random/RNGManager.h>
+
+#include <corsika/utl/CorsikaFenv.h>
+
+#include <boost/type_index.hpp>
+using boost::typeindex::type_id_with_cvr;
+
+#include <iostream>
+#include <limits>
+#include <typeinfo>
+
+using namespace corsika;
+using namespace corsika::process;
+using namespace corsika::units;
+using namespace corsika::particles;
+using namespace corsika::random;
+using namespace corsika::setup;
+using namespace corsika::geometry;
+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
+//
+int main() {
+ feenableexcept(FE_INVALID);
+ // initialize random number sequence(s)
+ random::RNGManager::GetInstance().RegisterRandomStream("cascade");
+
+ // setup environment, geometry
+ using EnvType = Environment<setup::IEnvironmentModel>;
+ EnvType env;
+ auto& universe = *(env.GetUniverse());
+
+ auto theMedium =
+ EnvType::CreateNode<Sphere>(Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m},
+ 1_km * std::numeric_limits<double>::infinity());
+
+ // fraction of oxygen
+ const float fox = 0.20946;
+ using MyHomogeneousModel = HomogeneousMedium<environment::IMediumModel>;
+ theMedium->SetModelProperties<MyHomogeneousModel>(
+ 1_kg / (1_m * 1_m * 1_m),
+ environment::NuclearComposition(
+ std::vector<particles::Code>{particles::Code::Nitrogen,
+ particles::Code::Oxygen},
+ std::vector<float>{(float)1. - fox, fox}));
+
+ universe.AddChild(std::move(theMedium));
+
+ const CoordinateSystem& rootCS = env.GetCoordinateSystem();
+
+ // setup processes, decays and interactions
+ tracking_line::TrackingLine tracking;
+
+ const std::vector<particles::Code> trackedHadrons = {
+ particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus,
+ particles::Code::KMinus, particles::Code::K0Long, particles::Code::K0Short};
+
+ random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
+ process::sibyll::Interaction sibyll;
+ process::sibyll::NuclearInteraction sibyllNuc(sibyll, env);
+ process::sibyll::Decay decay(trackedHadrons);
+ // random::RNGManager::GetInstance().RegisterRandomStream("pythia");
+ // process::pythia::Decay decay(trackedHadrons);
+ ProcessCut cut(20_GeV);
+
+ // random::RNGManager::GetInstance().RegisterRandomStream("HadronicElasticModel");
+ // process::HadronicElasticModel::HadronicElasticInteraction
+ // hadronicElastic(env);
+
+ process::TrackWriter::TrackWriter trackWriter("tracks.dat");
+ process::EnergyLoss::EnergyLoss eLoss;
+
+ // assemble all processes into an ordered process list
+ // cut << trackWriter;
+ auto sequence = sibyll << sibyllNuc << decay << eLoss << cut;
+
+ // cout << "decltype(sequence)=" << type_id_with_cvr<decltype(sequence)>().pretty_name()
+ // << "\n";
+
+ // setup particle stack, and add primary particle
+ setup::Stack stack;
+ stack.Clear();
+ const Code beamCode = Code::Nucleus;
+ const int nuclA = 4;
+ const int nuclZ = int(nuclA / 2.15 + 0.7);
+ const HEPMassType mass = GetNucleusMass(nuclA, nuclZ);
+ const HEPEnergyType E0 = nuclA * 10_TeV;
+ double theta = 0.;
+ double phi = 0.;
+
+ {
+ auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
+ return sqrt((Elab - m) * (Elab + m));
+ };
+ HEPMomentumType P0 = elab2plab(E0, mass);
+ auto momentumComponents = [](double theta, double phi, HEPMomentumType ptot) {
+ return std::make_tuple(ptot * sin(theta) * cos(phi), ptot * sin(theta) * sin(phi),
+ -ptot * cos(theta));
+ };
+ auto const [px, py, pz] =
+ momentumComponents(theta / 180. * M_PI, phi / 180. * M_PI, P0);
+ auto plab = corsika::stack::MomentumVector(rootCS, {px, py, pz});
+ 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
+ stack.AddParticle(std::tuple<particles::Code, units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, geometry::Point,
+ units::si::TimeType, unsigned short, unsigned short>{
+ beamCode, E0, plab, pos, 0_ns, nuclA, nuclZ});
+ }
+
+ // define air shower object, run simulation
+ cascade::Cascade EAS(env, tracking, sequence, stack);
+ EAS.Init();
+ EAS.Run();
+
+ eLoss.PrintProfile(); // print longitudinal profile
+
+ cut.ShowResults();
+ const HEPEnergyType Efinal =
+ cut.GetCutEnergy() + cut.GetInvEnergy() + cut.GetEmEnergy();
+ cout << "total cut energy (GeV): " << Efinal / 1_GeV << endl
+ << "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl;
+ cout << "total dEdX energy (GeV): " << eLoss.GetTotal() / 1_GeV << endl
+ << "relative difference (%): " << eLoss.GetTotal() / E0 * 100 << endl;
+}
diff --git a/Environment/FlatExponential.h b/Environment/FlatExponential.h
index 19e42c725c4b5204dc9c6fc166f9be00f44bdb7a..61fa2dced76357946bccfc9fc2951752361c9a42 100644
--- a/Environment/FlatExponential.h
+++ b/Environment/FlatExponential.h
@@ -1,5 +1,6 @@
+
/*
- * (c) Copyright 2019 CORSIKA Project, corsika-project@lists.kit.edu
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
diff --git a/Environment/VolumeTreeNode.h b/Environment/VolumeTreeNode.h
index faadc94c6a070876d34653fcf5cdebfb9482b7af..7a8b2e785f134ba59546e5622366a03c8f7554a4 100644
--- a/Environment/VolumeTreeNode.h
+++ b/Environment/VolumeTreeNode.h
@@ -108,7 +108,7 @@ namespace corsika::environment {
auto const& GetModelProperties() const { return *fModelProperties; }
- auto HasModelProperties() const { return fModelProperties != nullptr; }
+ bool HasModelProperties() const { return fModelProperties.get() != nullptr; }
template <typename ModelProperties, typename... Args>
auto SetModelProperties(Args&&... args) {
diff --git a/Framework/Cascade/CMakeLists.txt b/Framework/Cascade/CMakeLists.txt
index 3e294b214dd325ee9c082c32e2e0d46cdd70de3b..cce178637582bb48a8b17d1017bc489802c192aa 100644
--- a/Framework/Cascade/CMakeLists.txt
+++ b/Framework/Cascade/CMakeLists.txt
@@ -45,6 +45,7 @@ target_link_libraries (
CORSIKAcascade
ProcessStackInspector
ProcessTrackingLine
+ ProcessNullModel
CORSIKAstackinterface
CORSIKAprocesses
CORSIKAparticles
diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h
index c7a187447a5d3f6034600a05bffe09e3543cd7f0..75f29e4ed51bc70d7d31d4917ef4c186baee4bef 100644
--- a/Framework/Cascade/Cascade.h
+++ b/Framework/Cascade/Cascade.h
@@ -17,15 +17,28 @@
#include <corsika/random/ExponentialDistribution.h>
#include <corsika/random/RNGManager.h>
#include <corsika/random/UniformRealDistribution.h>
-#include <corsika/setup/SetupTrajectory.h>
+#include <corsika/stack/SecondaryView.h>
#include <corsika/units/PhysicalUnits.h>
+#include <corsika/setup/SetupTrajectory.h>
+
+/* see Issue 161, we need to include SetupStack only because we need
+ to globally define StackView. This is clearly not nice and should
+ be changed, when possible. It might be that StackView needs to be
+ templated in Cascade, but this would be even worse... so we don't
+ do that until it is really needed.
+ */
+#include <corsika/setup/SetupStack.h>
+
#include <cassert>
#include <cmath>
#include <iostream>
#include <limits>
#include <type_traits>
+#include <boost/type_index.hpp>
+using boost::typeindex::type_id_with_cvr;
+
/**
* The cascade namespace assembles all objects needed to simulate full particles cascades.
*/
@@ -39,23 +52,27 @@ namespace corsika::cascade {
* it very versatile. Via the template arguments physics models are
* plugged into the cascade simulation.
*
- * <b>Tracking</b> must be a class according to the
+ * <b>TTracking</b> must be a class according to the
* TrackingInterface providing the functions:
* <code>auto GetTrack(Particle const& p)</auto>,
* with the return type <code>geometry::Trajectory<corsika::geometry::Line>
* </code>
*
- * <b>ProcessList</b> must be a ProcessSequence. *
+ * <b>TProcessList</b> must be a ProcessSequence. *
* <b>Stack</b> is the storage object for particle data, i.e. with
* Particle class type <code>Stack::ParticleType</code>
*
*
-
*/
- template <typename Tracking, typename ProcessList, typename Stack>
+ template <typename TTracking, typename TProcessList, typename TStack,
+ /*
+ TStackView is needed as template parameter because of issue 161 and the
+ inability of clang to understand "MakeView" so far.
+ */
+ typename TStackView = corsika::setup::StackView>
class Cascade {
- using Particle = typename Stack::ParticleType;
+ using Particle = typename TStack::ParticleType;
using VolumeTreeNode =
std::remove_pointer_t<decltype(((Particle*)nullptr)->GetNode())>;
using MediumInterface = typename VolumeTreeNode::IModelProperties;
@@ -68,8 +85,8 @@ namespace corsika::cascade {
* Cascade class cannot be default constructed, but needs a valid
* list of physics processes for configuration at construct time.
*/
- Cascade(corsika::environment::Environment<MediumInterface> const& env, Tracking& tr,
- ProcessList& pl, Stack& stack)
+ Cascade(corsika::environment::Environment<MediumInterface> const& env, TTracking& tr,
+ TProcessList& pl, TStack& stack)
: fEnvironment(env)
, fTracking(tr)
, fProcessSequence(pl)
@@ -107,6 +124,7 @@ namespace corsika::cascade {
while (!fStack.IsEmpty()) {
auto pNext = fStack.GetNextParticle();
Step(pNext);
+ fProcessSequence.DoStack(fStack);
}
// do cascade equations, which can put new particles on Stack,
// thus, the double loop
@@ -125,15 +143,17 @@ namespace corsika::cascade {
* New particles produced in one step are subject to further
* processing, e.g. thinning, etc.
*/
- void Step(Particle& particle) {
+ void Step(Particle& vParticle) {
+ using namespace corsika;
using namespace corsika::units::si;
// determine geometric tracking
- auto [step, geomMaxLength, nextVol] = fTracking.GetTrack(particle);
+ auto [step, geomMaxLength, nextVol] = fTracking.GetTrack(vParticle);
+ [[maybe_unused]] auto const& dummy_nextVol = nextVol;
// determine combined total interaction length (inverse)
InverseGrammageType const total_inv_lambda =
- fProcessSequence.GetTotalInverseInteractionLength(particle, step);
+ fProcessSequence.GetTotalInverseInteractionLength(vParticle, step);
// sample random exponential step length in grammage
corsika::random::ExponentialDistribution expDist(1 / total_inv_lambda);
@@ -142,7 +162,7 @@ namespace corsika::cascade {
std::cout << "total_inv_lambda=" << total_inv_lambda
<< ", next_interact=" << next_interact << std::endl;
- auto const* currentLogicalNode = particle.GetNode();
+ auto const* currentLogicalNode = vParticle.GetNode();
// assert that particle stays outside void Universe if it has no
// model properties set
@@ -155,12 +175,12 @@ namespace corsika::cascade {
next_interact);
// determine the maximum geometric step length
- LengthType const distance_max = fProcessSequence.MaxStepLength(particle, step);
+ LengthType const distance_max = fProcessSequence.MaxStepLength(vParticle, step);
std::cout << "distance_max=" << distance_max << std::endl;
// determine combined total inverse decay time
InverseTimeType const total_inv_lifetime =
- fProcessSequence.GetTotalInverseLifetime(particle);
+ fProcessSequence.GetTotalInverseLifetime(vParticle);
// sample random exponential decay time
corsika::random::ExponentialDistribution expDistDecay(1 / total_inv_lifetime);
@@ -169,9 +189,8 @@ namespace corsika::cascade {
<< ", next_decay=" << next_decay << std::endl;
// convert next_decay from time to length [m]
- LengthType const distance_decay = next_decay * particle.GetMomentum().norm() /
- particle.GetEnergy() *
- corsika::units::constants::c;
+ LengthType const distance_decay = next_decay * vParticle.GetMomentum().norm() /
+ vParticle.GetEnergy() * units::constants::c;
// take minimum of geometry, interaction, decay for next step
auto const min_distance =
@@ -180,73 +199,98 @@ namespace corsika::cascade {
std::cout << " move particle by : " << min_distance << std::endl;
// here the particle is actually moved along the trajectory to new position:
- // std::visit(corsika::setup::ParticleUpdate<Particle>{particle}, step);
- particle.SetPosition(step.PositionFromArclength(min_distance));
+ // std::visit(setup::ParticleUpdate<Particle>{vParticle}, step);
+ vParticle.SetPosition(step.PositionFromArclength(min_distance));
// .... also update time, momentum, direction, ...
+ vParticle.SetTime(vParticle.GetTime() + min_distance / units::constants::c);
step.LimitEndTo(min_distance);
// apply all continuous processes on particle + track
- corsika::process::EProcessReturn status =
- fProcessSequence.DoContinuous(particle, step, fStack);
+ process::EProcessReturn status = fProcessSequence.DoContinuous(vParticle, step);
- if (status == corsika::process::EProcessReturn::eParticleAbsorbed) {
- std::cout << "Cascade: delete absorbed particle " << particle.GetPID() << " "
- << particle.GetEnergy() / 1_GeV << "GeV" << std::endl;
- particle.Delete();
+ if (status == process::EProcessReturn::eParticleAbsorbed) {
+ std::cout << "Cascade: delete absorbed particle " << vParticle.GetPID() << " "
+ << vParticle.GetEnergy() / 1_GeV << "GeV" << std::endl;
+ vParticle.Delete();
return;
}
std::cout << "sth. happening before geometric limit ? "
<< ((min_distance < geomMaxLength) ? "yes" : "no") << std::endl;
+ /*
+ Create SecondaryView object on Stack. The data container
+ remains untouched and identical, and 'projectil' is identical
+ to 'vParticle' above this line. However,
+ projectil.AddSecondaries populate the SecondaryView, which can
+ then be used afterwards for further processing. Thus: it is
+ important to use projectle (and not vParticle) for Interaction,
+ and Decay!
+ */
+
+ TStackView secondaries(vParticle);
+ [[maybe_unused]] auto projectile = secondaries.GetProjectile();
+
if (min_distance < geomMaxLength) { // interaction to happen within geometric limit
+
+ // check whether decay or interaction limits this step the
+ // outcome of decay or interaction MAY be a) new particles in
+ // secondaries, b) the projectile particle deleted (or
+ // changed)
+
if (min_distance == distance_interact) {
std::cout << "collide" << std::endl;
InverseGrammageType const actual_inv_length =
- fProcessSequence.GetTotalInverseInteractionLength(particle, step);
+ fProcessSequence.GetTotalInverseInteractionLength(vParticle, step);
- corsika::random::UniformRealDistribution<InverseGrammageType> uniDist(
- actual_inv_length);
+ random::UniformRealDistribution<InverseGrammageType> uniDist(actual_inv_length);
const auto sample_process = uniDist(fRNG);
InverseGrammageType inv_lambda_count = 0. * meter * meter / gram;
- fProcessSequence.SelectInteraction(particle, step, fStack, sample_process,
+ fProcessSequence.SelectInteraction(vParticle, projectile, step, sample_process,
inv_lambda_count);
} else if (min_distance == distance_decay) {
std::cout << "decay" << std::endl;
InverseTimeType const actual_decay_time =
- fProcessSequence.GetTotalInverseLifetime(particle);
+ fProcessSequence.GetTotalInverseLifetime(vParticle);
- corsika::random::UniformRealDistribution<InverseTimeType> uniDist(
- actual_decay_time);
+ random::UniformRealDistribution<InverseTimeType> uniDist(actual_decay_time);
const auto sample_process = uniDist(fRNG);
InverseTimeType inv_decay_count = 0 / second;
- fProcessSequence.SelectDecay(particle, fStack, sample_process, inv_decay_count);
+ fProcessSequence.SelectDecay(vParticle, projectile, sample_process,
+ inv_decay_count);
} else { // step-length limitation within volume
std::cout << "step-length limitation" << std::endl;
}
+ fProcessSequence.DoSecondaries(secondaries);
+ vParticle.Delete(); // todo: this should be reviewed. Where
+ // exactly are particles best deleted, and
+ // where they should NOT be
+ // deleted... maybe Delete function should
+ // be "protected" and not accessible to physics
+
auto const assertion = [&] {
auto const* numericalNodeAfterStep =
- fEnvironment.GetUniverse()->GetContainingNode(particle.GetPosition());
+ fEnvironment.GetUniverse()->GetContainingNode(vParticle.GetPosition());
return numericalNodeAfterStep == currentLogicalNode;
};
assert(assertion()); // numerical and logical nodes don't match
} else { // boundary crossing, step is limited by volume boundary
std::cout << "boundary crossing! next node = " << nextVol << std::endl;
- particle.SetNode(nextVol);
+ vParticle.SetNode(nextVol);
// DoBoundary may delete the particle (or not)
- fProcessSequence.DoBoundaryCrossing(particle, *currentLogicalNode, *nextVol);
+ fProcessSequence.DoBoundaryCrossing(vParticle, *currentLogicalNode, *nextVol);
}
}
private:
corsika::environment::Environment<MediumInterface> const& fEnvironment;
- Tracking& fTracking;
- ProcessList& fProcessSequence;
- Stack& fStack;
+ TTracking& fTracking;
+ TProcessList& fProcessSequence;
+ TStack& fStack;
corsika::random::RNG& fRNG =
corsika::random::RNGManager::GetInstance().GetRandomStream("cascade");
}; // namespace corsika::cascade
diff --git a/Framework/Cascade/testCascade.cc b/Framework/Cascade/testCascade.cc
index 5ec522a670de7ea048d33397e868590ea9e6faed..89fd65ef42150149a9c0ffbd7183fed865ba659d 100644
--- a/Framework/Cascade/testCascade.cc
+++ b/Framework/Cascade/testCascade.cc
@@ -9,13 +9,12 @@
* the license.
*/
-#include <limits>
-
#include <corsika/cascade/testCascade.h>
#include <corsika/cascade/Cascade.h>
#include <corsika/process/ProcessSequence.h>
+#include <corsika/process/null_model/NullModel.h>
#include <corsika/process/stack_inspector/StackInspector.h>
#include <corsika/process/tracking_line/TrackingLine.h>
@@ -28,10 +27,6 @@
#include <corsika/environment/HomogeneousMedium.h>
#include <corsika/environment/NuclearComposition.h>
-#include <corsika/setup/SetupStack.h>
-#include <corsika/setup/SetupTrajectory.h>
-using corsika::setup::Trajectory;
-
#define CATCH_CONFIG_MAIN // This tells Catch to provide a main() - only do this in one
// cpp file
#include <catch2/catch.hpp>
@@ -43,6 +38,7 @@ using namespace corsika::units::si;
using namespace corsika::geometry;
#include <iostream>
+#include <limits>
using namespace std;
auto MakeDummyEnv() {
@@ -51,11 +47,11 @@ auto MakeDummyEnv() {
auto theMedium = TestEnvironmentType::CreateNode<Sphere>(
Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m},
- 1_km * std::numeric_limits<double>::infinity());
+ 100_km * std::numeric_limits<double>::infinity());
using MyHomogeneousModel = environment::HomogeneousMedium<environment::IMediumModel>;
theMedium->SetModelProperties<MyHomogeneousModel>(
- 1_g / (1_m * 1_m * 1_m),
+ 1_g / (1_cm * 1_cm * 1_cm),
environment::NuclearComposition(
std::vector<particles::Code>{particles::Code::Proton}, std::vector<float>{1.}));
@@ -64,47 +60,75 @@ auto MakeDummyEnv() {
return env;
}
-class ProcessSplit : public process::ContinuousProcess<ProcessSplit> {
+class ProcessSplit : public process::InteractionProcess<ProcessSplit> {
+
+ int fCalls = 0;
+ GrammageType fX0;
+
+public:
+ ProcessSplit(GrammageType const X0)
+ : fX0(X0) {}
+
+ template <typename Particle, typename Track>
+ corsika::units::si::GrammageType GetInteractionLength(Particle&, Track&) const {
+ return fX0;
+ }
+
+ template <typename TProjectile>
+ corsika::process::EProcessReturn DoInteraction(TProjectile& vP) {
+ fCalls++;
+ const HEPEnergyType E = vP.GetEnergy();
+ vP.AddSecondary(
+ std::tuple<particles::Code, units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
+ vP.GetPID(), E / 2, vP.GetMomentum(), vP.GetPosition(), vP.GetTime()});
+ vP.AddSecondary(
+ std::tuple<particles::Code, units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
+ vP.GetPID(), E / 2, vP.GetMomentum(), vP.GetPosition(), vP.GetTime()});
+ return EProcessReturn::eInteracted;
+ }
+
+ void Init() { fCalls = 0; }
+
+ int GetCalls() const { return fCalls; }
+};
+
+class ProcessCut : public process::SecondariesProcess<ProcessCut> {
int fCount = 0;
int fCalls = 0;
HEPEnergyType fEcrit;
public:
- ProcessSplit(HEPEnergyType e)
+ ProcessCut(HEPEnergyType e)
: fEcrit(e) {}
- template <typename Particle, typename T>
- LengthType MaxStepLength(Particle&, T&) const {
- return 1_m;
- }
-
- template <typename Particle, typename T, typename Stack>
- EProcessReturn DoContinuous(Particle& p, T&, Stack&) {
+ template <typename TStack>
+ EProcessReturn DoSecondaries(TStack& vS) {
fCalls++;
- HEPEnergyType E = p.GetEnergy();
- if (E < fEcrit) {
- p.Delete();
- fCount++;
- } else {
- p.SetEnergy(E / 2);
- p.AddSecondary(std::tuple<particles::Code, units::si::HEPEnergyType,
- corsika::stack::MomentumVector, geometry::Point,
- units::si::TimeType>{p.GetPID(), E / 2, p.GetMomentum(),
- p.GetPosition(), p.GetTime()});
+ auto p = vS.begin();
+ while (p != vS.end()) {
+ HEPEnergyType E = p.GetEnergy();
+ if (E < fEcrit) {
+ p.Delete();
+ fCount++;
+ } else {
+ ++p; // next particle
+ }
}
+ cout << "ProcessCut::DoSecondaries size=" << vS.GetSize() << " count=" << fCount
+ << endl;
return EProcessReturn::eOk;
}
void Init() {
- fCount = 0;
fCalls = 0;
+ fCount = 0;
}
int GetCount() const { return fCount; }
int GetCalls() const { return fCalls; }
-
-private:
};
TEST_CASE("Cascade", "[Cascade]") {
@@ -114,14 +138,19 @@ TEST_CASE("Cascade", "[Cascade]") {
auto env = MakeDummyEnv();
tracking_line::TrackingLine tracking;
- stack_inspector::StackInspector<TestCascadeStack> p0(true);
+ stack_inspector::StackInspector<TestCascadeStack> stackInspect(1, true);
+ null_model::NullModel nullModel;
+ const GrammageType X0 = 20_g / square(1_cm);
const HEPEnergyType Ecrit = 85_MeV;
- ProcessSplit p1(Ecrit);
- auto sequence = p0 << p1;
+ ProcessSplit split(X0);
+ ProcessCut cut(Ecrit);
+ auto sequence = nullModel << stackInspect << split << cut;
TestCascadeStack stack;
- cascade::Cascade EAS(env, tracking, sequence, stack);
+ cascade::Cascade<tracking_line::TrackingLine, decltype(sequence), TestCascadeStack,
+ TestCascadeStackView>
+ EAS(env, tracking, sequence, stack);
CoordinateSystem const& rootCS =
RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
@@ -136,6 +165,7 @@ TEST_CASE("Cascade", "[Cascade]") {
EAS.Init();
EAS.Run();
- CHECK(p1.GetCount() == 2048);
- CHECK(p1.GetCalls() == 4095);
+ CHECK(cut.GetCount() == 2048);
+ CHECK(cut.GetCalls() == 2047);
+ CHECK(split.GetCalls() == 2047);
}
diff --git a/Framework/Cascade/testCascade.h b/Framework/Cascade/testCascade.h
index c84d1e14bff303274e648a4d018e711c4ff9dfd5..386a7281f429ad1b9a8ce83de4ed1e5e830ef040 100644
--- a/Framework/Cascade/testCascade.h
+++ b/Framework/Cascade/testCascade.h
@@ -26,8 +26,20 @@ template <typename StackIter>
using StackWithGeometryInterface = corsika::stack::CombinedParticleInterface<
corsika::setup::detail::ParticleDataStack::PIType, SetupGeometryDataInterface,
StackIter>;
+
using TestCascadeStack = corsika::stack::CombinedStack<
typename corsika::setup::detail::ParticleDataStack::StackImpl,
GeometryData<TestEnvironmentType>, StackWithGeometryInterface>;
+/*
+ See also Issue 161
+*/
+#if defined(__clang__)
+using TestCascadeStackView =
+ corsika::stack::SecondaryView<typename TestCascadeStack::StackImpl,
+ StackWithGeometryInterface>;
+#elif defined(__GNUC__) || defined(__GNUG__)
+using TestCascadeStackView = corsika::stack::MakeView<TestCascadeStack>::type;
+#endif
+
#endif
diff --git a/Framework/Geometry/Point.h b/Framework/Geometry/Point.h
index f05bd31d1ca628f1f32834a9e6d6300e6ef2578c..55dfa2950a1d5452bb9065ca4a814d4ee7562dee 100644
--- a/Framework/Geometry/Point.h
+++ b/Framework/Geometry/Point.h
@@ -37,6 +37,9 @@ namespace corsika::geometry {
// TODO: this should be private or protected, we don NOT want to expose numbers
// without reference to outside:
auto GetCoordinates() const { return BaseVector<length_d>::qVector; }
+ auto GetX() const { return BaseVector<length_d>::qVector.GetX(); }
+ auto GetY() const { return BaseVector<length_d>::qVector.GetY(); }
+ auto GetZ() const { return BaseVector<length_d>::qVector.GetZ(); }
/// this always returns a QuantityVector as triple
auto GetCoordinates(CoordinateSystem const& pCS) const {
diff --git a/Framework/Geometry/QuantityVector.h b/Framework/Geometry/QuantityVector.h
index 2712c188d8daf8c65f354bf416711c7d4d71e32c..c606a2934d6415aecefc541c71d2e0335be02948 100644
--- a/Framework/Geometry/QuantityVector.h
+++ b/Framework/Geometry/QuantityVector.h
@@ -55,6 +55,12 @@ namespace corsika::geometry {
return Quantity(phys::units::detail::magnitude_tag, eVector[index]);
}
+ auto GetX() const { return Quantity(phys::units::detail::magnitude_tag, eVector[0]); }
+
+ auto GetY() const { return Quantity(phys::units::detail::magnitude_tag, eVector[1]); }
+
+ auto GetZ() const { return Quantity(phys::units::detail::magnitude_tag, eVector[2]); }
+
auto norm() const {
return Quantity(phys::units::detail::magnitude_tag, eVector.norm());
}
diff --git a/Framework/ProcessSequence/BaseProcess.h b/Framework/ProcessSequence/BaseProcess.h
index 8b4014478f5ebef33e4acd76f6b7315f2c1e5041..6733fcf9a9f66a53c3a4d55326cdc37eabbcc314 100644
--- a/Framework/ProcessSequence/BaseProcess.h
+++ b/Framework/ProcessSequence/BaseProcess.h
@@ -13,6 +13,7 @@
#define _include_corsika_baseprocess_h_
#include <corsika/process/ProcessReturn.h> // for convenience
+#include <type_traits>
namespace corsika::process {
diff --git a/Framework/ProcessSequence/CMakeLists.txt b/Framework/ProcessSequence/CMakeLists.txt
index 9ec9072dabb1dab6ecb17213ddc84a45bb73cd5f..33aad8236f6b2943d083fd95addcea1698700d43 100644
--- a/Framework/ProcessSequence/CMakeLists.txt
+++ b/Framework/ProcessSequence/CMakeLists.txt
@@ -13,7 +13,9 @@ set (
BaseProcess.h
BoundaryCrossingProcess.h
ContinuousProcess.h
+ SecondariesProcess.h
InteractionProcess.h
+ StackProcess.h
DecayProcess.h
ProcessSequence.h
ProcessReturn.h
diff --git a/Framework/ProcessSequence/ContinuousProcess.h b/Framework/ProcessSequence/ContinuousProcess.h
index 9a3424a35d4637a9e21212a0909f2456f71b3597..58d378aa3eb5bdad5f8029555ea8825afa450007 100644
--- a/Framework/ProcessSequence/ContinuousProcess.h
+++ b/Framework/ProcessSequence/ContinuousProcess.h
@@ -33,8 +33,8 @@ namespace corsika::process {
// here starts the interface part
// -> enforce derived to implement DoContinuous...
- template <typename Particle, typename Track, typename Stack>
- EProcessReturn DoContinuous(Particle&, Track&, Stack&) const;
+ template <typename Particle, typename Track>
+ EProcessReturn DoContinuous(Particle&, Track&) const;
// -> enforce derived to implement MaxStepLength...
template <typename Particle, typename Track>
diff --git a/Framework/ProcessSequence/DecayProcess.h b/Framework/ProcessSequence/DecayProcess.h
index 82be572d6fc8e19260fa50a5267f1171c20e042b..16d63ae667775848db49076d242fa451f38060e6 100644
--- a/Framework/ProcessSequence/DecayProcess.h
+++ b/Framework/ProcessSequence/DecayProcess.h
@@ -34,15 +34,15 @@ namespace corsika::process {
/// here starts the interface-definition part
// -> enforce derived to implement DoDecay...
- template <typename Particle, typename Stack>
- EProcessReturn DoDecay(Particle&, Stack&);
+ template <typename Particle>
+ EProcessReturn DoDecay(Particle&);
template <typename Particle>
corsika::units::si::TimeType GetLifetime(Particle& p);
template <typename Particle>
- corsika::units::si::InverseTimeType GetInverseLifetime(Particle& p) {
- return 1. / GetRef().GetLifetime(p);
+ corsika::units::si::InverseTimeType GetInverseLifetime(Particle& vP) {
+ return 1. / GetRef().GetLifetime(vP);
}
};
diff --git a/Framework/ProcessSequence/InteractionProcess.h b/Framework/ProcessSequence/InteractionProcess.h
index 723cbf02fcc7ce443d968300293ccc44d1b40811..27ad9c222d5ce644d6b344a469edaf8e717fcf87 100644
--- a/Framework/ProcessSequence/InteractionProcess.h
+++ b/Framework/ProcessSequence/InteractionProcess.h
@@ -35,8 +35,8 @@ namespace corsika::process {
/// here starts the interface-definition part
// -> enforce derived to implement DoInteraction...
- template <typename P, typename S>
- EProcessReturn DoInteraction(P&, S&);
+ template <typename Particle>
+ EProcessReturn DoInteraction(Particle&);
template <typename Particle, typename Track>
corsika::units::si::GrammageType GetInteractionLength(Particle& p, Track& t);
diff --git a/Framework/ProcessSequence/ProcessSequence.h b/Framework/ProcessSequence/ProcessSequence.h
index c697d09f416d49025a0d4763f00851c047cedb62..3bd759c04a97e55360c87473cc4457159bd967fa 100644
--- a/Framework/ProcessSequence/ProcessSequence.h
+++ b/Framework/ProcessSequence/ProcessSequence.h
@@ -18,6 +18,8 @@
#include <corsika/process/DecayProcess.h>
#include <corsika/process/InteractionProcess.h>
#include <corsika/process/ProcessReturn.h>
+#include <corsika/process/SecondariesProcess.h>
+#include <corsika/process/StackProcess.h>
#include <corsika/units/PhysicalUnits.h>
#include <cmath>
@@ -89,78 +91,120 @@ namespace corsika::process {
return ret;
}
- template <typename Particle, typename Track, typename Stack>
- EProcessReturn DoContinuous(Particle& p, Track& t, Stack& s) {
+ template <typename TParticle, typename TTrack>
+ EProcessReturn DoContinuous(TParticle& vP, TTrack& vT) {
EProcessReturn ret = EProcessReturn::eOk;
if constexpr (std::is_base_of<ContinuousProcess<T1type>, T1type>::value ||
is_process_sequence<T1>::value) {
- ret |= A.DoContinuous(p, t, s);
+ ret |= A.DoContinuous(vP, vT);
}
if constexpr (std::is_base_of<ContinuousProcess<T2type>, T2type>::value ||
is_process_sequence<T2>::value) {
- ret |= B.DoContinuous(p, t, s);
+ ret |= B.DoContinuous(vP, vT);
}
return ret;
}
- template <typename Particle, typename Track>
- corsika::units::si::LengthType MaxStepLength(Particle& p, Track& track) {
+ template <typename TSecondaries>
+ EProcessReturn DoSecondaries(TSecondaries& vS) {
+ EProcessReturn ret = EProcessReturn::eOk;
+ if constexpr (std::is_base_of<SecondariesProcess<T1type>, T1type>::value ||
+ is_process_sequence<T1>::value) {
+ ret |= A.DoSecondaries(vS);
+ }
+ if constexpr (std::is_base_of<SecondariesProcess<T2type>, T2type>::value ||
+ is_process_sequence<T2>::value) {
+ ret |= B.DoSecondaries(vS);
+ }
+ return ret;
+ }
+
+ bool CheckStep() {
+ bool ret = false;
+ if constexpr (std::is_base_of<StackProcess<T1type>, T1type>::value ||
+ is_process_sequence<T1>::value) {
+ ret |= A.CheckStep();
+ }
+ if constexpr (std::is_base_of<StackProcess<T2type>, T2type>::value ||
+ is_process_sequence<T2>::value) {
+ ret |= B.CheckStep();
+ }
+ return ret;
+ }
+
+ template <typename TStack>
+ EProcessReturn DoStack(TStack& vS) {
+ EProcessReturn ret = EProcessReturn::eOk;
+ if constexpr (std::is_base_of<StackProcess<T1type>, T1type>::value ||
+ is_process_sequence<T1>::value) {
+ if (A.CheckStep()) { ret |= A.DoStack(vS); }
+ }
+ if constexpr (std::is_base_of<StackProcess<T2type>, T2type>::value ||
+ is_process_sequence<T2>::value) {
+ if (B.CheckStep()) { ret |= B.DoStack(vS); }
+ }
+ return ret;
+ }
+
+ template <typename TParticle, typename TTrack>
+ corsika::units::si::LengthType MaxStepLength(TParticle& vP, TTrack& vTrack) {
corsika::units::si::LengthType
max_length = // if no other process in the sequence implements it
std::numeric_limits<double>::infinity() * corsika::units::si::meter;
if constexpr (std::is_base_of<ContinuousProcess<T1type>, T1type>::value ||
is_process_sequence<T1>::value) {
- corsika::units::si::LengthType const len = A.MaxStepLength(p, track);
+ corsika::units::si::LengthType const len = A.MaxStepLength(vP, vTrack);
max_length = std::min(max_length, len);
}
if constexpr (std::is_base_of<ContinuousProcess<T2type>, T2type>::value ||
is_process_sequence<T2>::value) {
- corsika::units::si::LengthType const len = B.MaxStepLength(p, track);
+ corsika::units::si::LengthType const len = B.MaxStepLength(vP, vTrack);
max_length = std::min(max_length, len);
}
return max_length;
}
- template <typename Particle, typename Track>
- corsika::units::si::GrammageType GetTotalInteractionLength(Particle& p, Track& t) {
- return 1. / GetInverseInteractionLength(p, t);
+ template <typename TParticle, typename TTrack>
+ corsika::units::si::GrammageType GetTotalInteractionLength(TParticle& vP,
+ TTrack& vT) {
+ return 1. / GetInverseInteractionLength(vP, vT);
}
- template <typename Particle, typename Track>
- corsika::units::si::InverseGrammageType GetTotalInverseInteractionLength(Particle& p,
- Track& t) {
- return GetInverseInteractionLength(p, t);
+ template <typename TParticle, typename TTrack>
+ corsika::units::si::InverseGrammageType GetTotalInverseInteractionLength(
+ TParticle& vP, TTrack& vT) {
+ return GetInverseInteractionLength(vP, vT);
}
- template <typename Particle, typename Track>
- corsika::units::si::InverseGrammageType GetInverseInteractionLength(Particle& p,
- Track& t) {
+ template <typename TParticle, typename TTrack>
+ corsika::units::si::InverseGrammageType GetInverseInteractionLength(TParticle& vP,
+ TTrack& vT) {
using namespace corsika::units::si;
InverseGrammageType tot = 0 * meter * meter / gram;
if constexpr (std::is_base_of<InteractionProcess<T1type>, T1type>::value ||
is_process_sequence<T1>::value) {
- tot += A.GetInverseInteractionLength(p, t);
+ tot += A.GetInverseInteractionLength(vP, vT);
}
if constexpr (std::is_base_of<InteractionProcess<T2type>, T2type>::value ||
is_process_sequence<T2>::value) {
- tot += B.GetInverseInteractionLength(p, t);
+ tot += B.GetInverseInteractionLength(vP, vT);
}
return tot;
}
- template <typename Particle, typename Track, typename Stack>
+ template <typename TParticle, typename TSecondaries, typename TTrack>
EProcessReturn SelectInteraction(
- Particle& vP, Track& vT, Stack& vS,
+ TParticle& vP, TSecondaries& vS, TTrack& vT,
[[maybe_unused]] corsika::units::si::InverseGrammageType lambda_select,
corsika::units::si::InverseGrammageType& lambda_inv_count) {
if constexpr (is_process_sequence<T1type>::value) {
// if A is a process sequence --> check inside
const EProcessReturn ret =
- A.SelectInteraction(vP, vT, vS, lambda_select, lambda_inv_count);
+ A.SelectInteraction(vP, vS, vT, lambda_select, lambda_inv_count);
// if A did succeed, stop routine
if (ret != EProcessReturn::eOk) { return ret; }
} else if constexpr (std::is_base_of<InteractionProcess<T1type>, T1type>::value) {
@@ -168,7 +212,7 @@ namespace corsika::process {
lambda_inv_count += A.GetInverseInteractionLength(vP, vT);
// check if we should execute THIS process and then EXIT
if (lambda_select < lambda_inv_count) {
- A.DoInteraction(vP, vS);
+ A.DoInteraction(vS);
return EProcessReturn::eInteracted;
}
} // end branch A
@@ -176,7 +220,7 @@ namespace corsika::process {
if constexpr (is_process_sequence<T2>::value) {
// if A is a process sequence --> check inside
const EProcessReturn ret =
- B.SelectInteraction(vP, vT, vS, lambda_select, lambda_inv_count);
+ B.SelectInteraction(vP, vS, vT, lambda_select, lambda_inv_count);
// if A did succeed, stop routine
if (ret != EProcessReturn::eOk) { return ret; }
} else if constexpr (std::is_base_of<InteractionProcess<T2type>, T2type>::value) {
@@ -184,25 +228,25 @@ namespace corsika::process {
lambda_inv_count += B.GetInverseInteractionLength(vP, vT);
// check if we should execute THIS process and then EXIT
if (lambda_select < lambda_inv_count) {
- B.DoInteraction(vP, vS);
+ B.DoInteraction(vS);
return EProcessReturn::eInteracted;
}
} // end branch A
return EProcessReturn::eOk;
}
- template <typename Particle>
- corsika::units::si::TimeType GetTotalLifetime(Particle& p) {
+ template <typename TParticle>
+ corsika::units::si::TimeType GetTotalLifetime(TParticle& p) {
return 1. / GetInverseLifetime(p);
}
- template <typename Particle>
- corsika::units::si::InverseTimeType GetTotalInverseLifetime(Particle& p) {
+ template <typename TParticle>
+ corsika::units::si::InverseTimeType GetTotalInverseLifetime(TParticle& p) {
return GetInverseLifetime(p);
}
- template <typename Particle>
- corsika::units::si::InverseTimeType GetInverseLifetime(Particle& p) {
+ template <typename TParticle>
+ corsika::units::si::InverseTimeType GetInverseLifetime(TParticle& p) {
using namespace corsika::units::si;
corsika::units::si::InverseTimeType tot = 0 / second;
@@ -219,38 +263,38 @@ namespace corsika::process {
}
// select decay process
- template <typename Particle, typename Stack>
+ template <typename TParticle, typename TSecondaries>
EProcessReturn SelectDecay(
- Particle& p, Stack& s,
+ TParticle& vP, TSecondaries& vS,
[[maybe_unused]] corsika::units::si::InverseTimeType decay_select,
corsika::units::si::InverseTimeType& decay_inv_count) {
if constexpr (is_process_sequence<T1>::value) {
// if A is a process sequence --> check inside
- const EProcessReturn ret = A.SelectDecay(p, s, decay_select, decay_inv_count);
+ const EProcessReturn ret = A.SelectDecay(vP, vS, decay_select, decay_inv_count);
// if A did succeed, stop routine
if (ret != EProcessReturn::eOk) { return ret; }
} else if constexpr (std::is_base_of<DecayProcess<T1type>, T1type>::value) {
// if this is not a ContinuousProcess --> evaluate probability
- decay_inv_count += A.GetInverseLifetime(p);
+ decay_inv_count += A.GetInverseLifetime(vP);
// check if we should execute THIS process and then EXIT
if (decay_select < decay_inv_count) { // more pedagogical: rndm_select <
// decay_inv_count / decay_inv_tot
- A.DoDecay(p, s);
+ A.DoDecay(vS);
return EProcessReturn::eDecayed;
}
} // end branch A
if constexpr (is_process_sequence<T2>::value) {
// if A is a process sequence --> check inside
- const EProcessReturn ret = B.SelectDecay(p, s, decay_select, decay_inv_count);
+ const EProcessReturn ret = B.SelectDecay(vP, vS, decay_select, decay_inv_count);
// if A did succeed, stop routine
if (ret != EProcessReturn::eOk) { return ret; }
} else if constexpr (std::is_base_of<DecayProcess<T2type>, T2type>::value) {
// if this is not a ContinuousProcess --> evaluate probability
- decay_inv_count += B.GetInverseLifetime(p);
+ decay_inv_count += B.GetInverseLifetime(vP);
// check if we should execute THIS process and then EXIT
if (decay_select < decay_inv_count) {
- B.DoDecay(p, s);
+ B.DoDecay(vS);
return EProcessReturn::eDecayed;
}
} // end branch B
@@ -272,8 +316,8 @@ namespace corsika::process {
typename P1, typename P2,
typename std::enable_if<is_process<typename std::decay<P1>::type>::value &&
is_process<typename std::decay<P2>::type>::value>::type...>
- inline auto operator<<(P1&& A, P2&& B) -> ProcessSequence<P1, P2> {
- return ProcessSequence<P1, P2>(A.GetRef(), B.GetRef());
+ inline auto operator<<(P1&& vA, P2&& vB) -> ProcessSequence<P1, P2> {
+ return ProcessSequence<P1, P2>(vA.GetRef(), vB.GetRef());
}
/// marker to identify objectas ProcessSequence
diff --git a/Framework/ProcessSequence/SecondariesProcess.h b/Framework/ProcessSequence/SecondariesProcess.h
new file mode 100644
index 0000000000000000000000000000000000000000..7e4f11e442e39ecfdde376332ceb276095fc8484
--- /dev/null
+++ b/Framework/ProcessSequence/SecondariesProcess.h
@@ -0,0 +1,48 @@
+
+/*
+ * (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 _include_corsika_secondariesprocess_h_
+#define _include_corsika_secondariesprocess_h_
+
+#include <corsika/process/ProcessReturn.h> // for convenience
+#include <corsika/setup/SetupTrajectory.h>
+#include <corsika/units/PhysicalUnits.h>
+
+namespace corsika::process {
+
+ /**
+ \class SecondariesProcess
+
+ The structural base type of a process object in a
+ ProcessSequence. Both, the ProcessSequence and all its elements
+ are of type SecondariesProcess<T>
+
+ */
+
+ template <typename derived>
+ struct SecondariesProcess {
+
+ derived& GetRef() { return static_cast<derived&>(*this); }
+ const derived& GetRef() const { return static_cast<const derived&>(*this); }
+
+ /// here starts the interface-definition part
+ // -> enforce derived to implement DoSecondaries...
+ template <typename TSecondaries>
+ inline EProcessReturn DoSecondaries(TSecondaries&);
+ };
+
+ // overwrite the default trait class, to mark BaseProcess<T> as useful process
+ template <class T>
+ std::true_type is_process_impl(const SecondariesProcess<T>* impl);
+
+} // namespace corsika::process
+
+#endif
diff --git a/Framework/ProcessSequence/SecondaryProcess.h b/Framework/ProcessSequence/SecondaryProcess.h
new file mode 100644
index 0000000000000000000000000000000000000000..523ce4d88ff4b2a3684958f6015540ab415682d3
--- /dev/null
+++ b/Framework/ProcessSequence/SecondaryProcess.h
@@ -0,0 +1,48 @@
+
+/*
+ * (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 _include_corsika_secondariesprocess_h_
+#define _include_corsika_secondariesprocess_h_
+
+#include <corsika/process/ProcessReturn.h> // for convenience
+#include <corsika/setup/SetupTrajectory.h>
+#include <corsika/units/PhysicalUnits.h>
+
+namespace corsika::process {
+
+ /**
+ \class SecondariesProcess
+
+ The structural base type of a process object in a
+ ProcessSequence. Both, the ProcessSequence and all its elements
+ are of type SecondariesProcess<T>
+
+ */
+
+ template <typename derived>
+ struct SecondariesProcess {
+
+ derived& GetRef() { return static_cast<derived&>(*this); }
+ const derived& GetRef() const { return static_cast<const derived&>(*this); }
+
+ /// here starts the interface-definition part
+ // -> enforce derived to implement DoSecondaries...
+ template <typename Particle>
+ inline EProcessReturn DoSecondaries(Particle&);
+ };
+
+ // overwrite the default trait class, to mark BaseProcess<T> as useful process
+ template <class T>
+ std::true_type is_process_impl(const SecondariesProcess<T>* impl);
+
+} // namespace corsika::process
+
+#endif
diff --git a/Framework/ProcessSequence/StackProcess.h b/Framework/ProcessSequence/StackProcess.h
new file mode 100644
index 0000000000000000000000000000000000000000..1458fa75ca5938cdd42d50d0b4d22af46f641302
--- /dev/null
+++ b/Framework/ProcessSequence/StackProcess.h
@@ -0,0 +1,66 @@
+
+/*
+ * (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 _include_corsika_stackprocess_h_
+#define _include_corsika_stackprocess_h_
+
+#include <corsika/process/ProcessReturn.h> // for convenience
+#include <corsika/setup/SetupTrajectory.h>
+#include <corsika/units/PhysicalUnits.h>
+
+namespace corsika::process {
+
+ /**
+ \class StackProcess
+
+ The structural base type of a process object in a
+ ProcessSequence. Both, the ProcessSequence and all its elements
+ are of type StackProcess<T>
+
+ */
+
+ template <typename derived>
+ class StackProcess {
+
+ public:
+ StackProcess() = delete;
+ StackProcess(const unsigned int nStep)
+ : fNStep(nStep) {}
+
+ derived& GetRef() { return static_cast<derived&>(*this); }
+ const derived& GetRef() const { return static_cast<const derived&>(*this); }
+
+ /// here starts the interface-definition part
+ // -> enforce derived to implement DoStack...
+ template <typename TStack>
+ inline EProcessReturn DoStack(TStack&);
+
+ bool CheckStep() { return !((++fIStep) % fNStep); }
+
+ private:
+ /**
+ @name The number of "steps" during the cascade processing after
+ which this StackProcess is going to be executed. The logic is
+ "fIStep modulo fNStep"
+ @{
+ */
+ unsigned int fNStep = 0;
+ unsigned long int fIStep = 0;
+ //! @}
+ };
+
+ // overwrite the default trait class, to mark BaseProcess<T> as useful process
+ template <class T>
+ std::true_type is_process_impl(const StackProcess<T>* impl);
+
+} // namespace corsika::process
+
+#endif
diff --git a/Framework/ProcessSequence/testProcessSequence.cc b/Framework/ProcessSequence/testProcessSequence.cc
index 8fd46418ac838ca1af0f2beb3691bf599a40eb05..aa78ee89c4c012eb2df3d322b6b1ddf6fa44ba1e 100644
--- a/Framework/ProcessSequence/testProcessSequence.cc
+++ b/Framework/ProcessSequence/testProcessSequence.cc
@@ -39,8 +39,8 @@ public:
assert(globalCount == fV);
globalCount++;
}
- template <typename D, typename T, typename S>
- inline EProcessReturn DoContinuous(D& d, T&, S&) const {
+ template <typename D, typename T>
+ inline EProcessReturn DoContinuous(D& d, T&) const {
cout << "ContinuousProcess1::DoContinuous" << endl;
for (int i = 0; i < nData; ++i) d.p[i] += 0.933;
return EProcessReturn::eOk;
@@ -58,8 +58,8 @@ public:
assert(globalCount == fV);
globalCount++;
}
- template <typename D, typename T, typename S>
- inline EProcessReturn DoContinuous(D& d, T&, S&) const {
+ template <typename D, typename T>
+ inline EProcessReturn DoContinuous(D& d, T&) const {
cout << "ContinuousProcess2::DoContinuous" << endl;
for (int i = 0; i < nData; ++i) d.p[i] += 0.933;
return EProcessReturn::eOk;
@@ -95,8 +95,8 @@ public:
assert(globalCount == fV);
globalCount++;
}
- template <typename Particle, typename Stack>
- inline EProcessReturn DoInteraction(Particle&, Stack&) const {
+ template <typename Particle>
+ inline EProcessReturn DoInteraction(Particle&) const {
cout << "Process2::DoInteraction" << endl;
return EProcessReturn::eOk;
}
@@ -118,8 +118,8 @@ public:
assert(globalCount == fV);
globalCount++;
}
- template <typename Particle, typename Stack>
- inline EProcessReturn DoInteraction(Particle&, Stack&) const {
+ template <typename Particle>
+ inline EProcessReturn DoInteraction(Particle&) const {
cout << "Process3::DoInteraction" << endl;
return EProcessReturn::eOk;
}
@@ -141,14 +141,14 @@ public:
assert(globalCount == fV);
globalCount++;
}
- template <typename D, typename T, typename S>
- inline EProcessReturn DoContinuous(D& d, T&, S&) const {
+ template <typename D, typename T>
+ inline EProcessReturn DoContinuous(D& d, T&) const {
for (int i = 0; i < nData; ++i) { d.p[i] /= 1.2; }
return EProcessReturn::eOk;
}
// inline double MinStepLength(D& d) {
- template <typename Particle, typename Stack>
- EProcessReturn DoInteraction(Particle&, Stack&) const {
+ template <typename Particle>
+ EProcessReturn DoInteraction(Particle&) const {
return EProcessReturn::eOk;
}
};
@@ -168,16 +168,30 @@ public:
TimeType GetLifetime(Particle&) const {
return 1_s;
}
- template <typename Particle, typename Stack>
- EProcessReturn DoDecay(Particle&, Stack&) const {
+ template <typename Particle>
+ EProcessReturn DoDecay(Particle&) const {
+ return EProcessReturn::eOk;
+ }
+};
+
+class Stack1 : public StackProcess<Stack1> {
+ int fCount = 0;
+
+public:
+ Stack1(const int n)
+ : StackProcess(n) {}
+ template <typename TStack>
+ EProcessReturn DoStack(TStack&) {
+ fCount++;
return EProcessReturn::eOk;
}
+ int GetCount() const { return fCount; }
};
+struct DummyStack {};
struct DummyData {
double p[nData] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
};
-struct DummyStack {};
struct DummyTrajectory {};
TEST_CASE("Process Sequence", "[Process Sequence]") {
@@ -205,12 +219,13 @@ TEST_CASE("Process Sequence", "[Process Sequence]") {
Process2 m2(1);
Process3 m3(2);
- DummyStack s;
- DummyTrajectory t;
+ DummyData particle;
+ DummyTrajectory track;
auto sequence2 = cp1 << m2 << m3;
- GrammageType const tot = sequence2.GetTotalInteractionLength(s, t);
- InverseGrammageType const tot_inv = sequence2.GetTotalInverseInteractionLength(s, t);
+ GrammageType const tot = sequence2.GetTotalInteractionLength(particle, track);
+ InverseGrammageType const tot_inv =
+ sequence2.GetTotalInverseInteractionLength(particle, track);
cout << "lambda_tot=" << tot << "; lambda_tot_inv=" << tot_inv << endl;
}
@@ -220,11 +235,11 @@ TEST_CASE("Process Sequence", "[Process Sequence]") {
Process3 m3(2);
Decay1 d3(2);
- DummyStack s;
+ DummyData particle;
auto sequence2 = cp1 << m2 << m3 << d3;
- TimeType const tot = sequence2.GetTotalLifetime(s);
- InverseTimeType const tot_inv = sequence2.GetTotalInverseLifetime(s);
+ TimeType const tot = sequence2.GetTotalLifetime(particle);
+ InverseTimeType const tot_inv = sequence2.GetTotalInverseLifetime(particle);
cout << "lambda_tot=" << tot << "; lambda_tot_inv=" << tot_inv << endl;
}
@@ -237,27 +252,44 @@ TEST_CASE("Process Sequence", "[Process Sequence]") {
auto sequence2 = cp1 << m2 << m3 << cp2;
- DummyData p;
- DummyStack s;
- DummyTrajectory t;
+ DummyData particle;
+ DummyTrajectory track;
cout << "-->init sequence2" << endl;
globalCount = 0;
sequence2.Init();
cout << "-->docont" << endl;
- sequence2.DoContinuous(p, t, s);
+ sequence2.DoContinuous(particle, track);
cout << "-->dodisc" << endl;
- // sequence2.DoInteraction(p, s);
cout << "-->done" << endl;
const int nLoop = 5;
cout << "Running loop with n=" << nLoop << endl;
- for (int i = 0; i < nLoop; ++i) {
- sequence2.DoContinuous(p, t, s);
- // sequence2.DoInteraction(p, s);
+ for (int i = 0; i < nLoop; ++i) { sequence2.DoContinuous(particle, track); }
+ for (int i = 0; i < nData; i++) {
+ cout << "data[" << i << "]=" << particle.p[i] << endl;
}
- for (int i = 0; i < nData; i++) { cout << "data[" << i << "]=" << p.p[i] << endl; }
cout << "done" << endl;
}
+
+ SECTION("StackProcess") {
+
+ ContinuousProcess1 cp1(0);
+ ContinuousProcess2 cp2(3);
+ Process2 m2(1);
+ Process3 m3(2);
+ Stack1 s1(1);
+ Stack1 s2(2);
+
+ auto sequence = s1 << s2;
+
+ DummyStack stack;
+
+ const int nLoop = 20;
+ for (int i = 0; i < nLoop; ++i) { sequence.DoStack(stack); }
+
+ CHECK(s1.GetCount() == 20);
+ CHECK(s2.GetCount() == 10);
+ }
}
diff --git a/Framework/StackInterface/CombinedStack.h b/Framework/StackInterface/CombinedStack.h
index 44a29a6280bbf7aa61e05301d3bd9fd9adc8dd5b..8b987d5100732a7a339a23d6fadd3530870bb2ef 100644
--- a/Framework/StackInterface/CombinedStack.h
+++ b/Framework/StackInterface/CombinedStack.h
@@ -19,51 +19,92 @@
namespace corsika::stack {
/**
+ * @class CombinedParticleInterface
*
+ * You may combine two StackData object, see class CombinedStackImpl
+ * below, into one Stack, using a combined StackIterator (aka
+ * CombinedParticleInterface) interface class.
+ *
+ * This allows to add specific information to a given Stack, could
+ * be special information on a subset of entries
+ * (e.g. NuclearStackExtension) or also (multi) thinning weights for
+ * all particles.
+ *
+ * Many Stacks can be combined into more complex object.
+ *
+ * The two sub-stacks must both provide their independent
+ * ParticleInterface classes.
*
*/
- template <template <typename> typename ParticleInterface,
- template <typename> typename ParticleInterfaceAdd, typename StackIterator>
+ template <template <typename> typename ParticleInterfaceA,
+ template <typename> typename ParticleInterfaceB, typename StackIterator>
class CombinedParticleInterface
- : public ParticleInterfaceAdd<ParticleInterface<StackIterator>> {
+ : public ParticleInterfaceB<ParticleInterfaceA<StackIterator>> {
+
+ // template<template <typename> typename _PI>
+ // template <typename StackDataType, template <typename> typename ParticleInterface>
+ // template<typename T1, template <typename> typename T2> friend class Stack<T1, T2>;
- using C =
- CombinedParticleInterface<ParticleInterface, ParticleInterfaceAdd, StackIterator>;
- using T = ParticleInterfaceAdd<ParticleInterface<StackIterator>>;
- using I = ParticleInterface<StackIterator>;
+ using PI_C =
+ CombinedParticleInterface<ParticleInterfaceA, ParticleInterfaceB, StackIterator>;
+ using PI_A = ParticleInterfaceA<StackIterator>;
+ using PI_B = ParticleInterfaceB<ParticleInterfaceA<StackIterator>>;
protected:
- using T::GetIndex;
- using T::GetStackData;
+ using PI_B::GetIndex; // choose B, A would also work
+ using PI_B::GetStackData; // choose B, A would also work
public:
+ /**
+ * @name wrapper for user functions
+ * @{
+ *
+ * In this set of functions we call the user-provide
+ * ParticleInterface SetParticleData(...) methods, either with
+ * parent particle reference, or w/o.
+ *
+ * There is one implicit assumption here: if only one data tuple
+ * is provided for SetParticleData, the data is passed on to
+ * ParticleInterfaceA and the ParticleInterfaceB is
+ * default-initialized. There are many occasions where this is the
+ * desired behaviour, e.g. for thinning etc.
+ *
+ */
+
template <typename... Args1>
void SetParticleData(const std::tuple<Args1...> vA) {
- I::SetParticleData(vA);
- T::SetParticleData();
+ PI_A::SetParticleData(vA);
+ PI_B::SetParticleData();
}
template <typename... Args1, typename... Args2>
void SetParticleData(const std::tuple<Args1...> vA, const std::tuple<Args2...> vB) {
- I::SetParticleData(vA);
- T::SetParticleData(vB);
+ PI_A::SetParticleData(vA);
+ PI_B::SetParticleData(vB);
}
template <typename... Args1>
- void SetParticleData(C& p, const std::tuple<Args1...> vA) {
+ void SetParticleData(PI_C& p, const std::tuple<Args1...> vA) {
// static_assert(MT<I>::has_not, "error");
- I::SetParticleData(static_cast<I&>(p), vA); // original stack
- T::SetParticleData(static_cast<T&>(p)); // addon stack
+ PI_A::SetParticleData(static_cast<PI_A&>(p), vA); // original stack
+ PI_B::SetParticleData(static_cast<PI_B&>(p)); // addon stack
}
template <typename... Args1, typename... Args2>
- void SetParticleData(C& p, const std::tuple<Args1...> vA,
+ void SetParticleData(PI_C& p, const std::tuple<Args1...> vA,
const std::tuple<Args2...> vB) {
- I::SetParticleData(static_cast<I&>(p), vA);
- T::SetParticleData(static_cast<T&>(p), vB);
+ PI_A::SetParticleData(static_cast<PI_A&>(p), vA);
+ PI_B::SetParticleData(static_cast<PI_B&>(p), vB);
}
+ ///@}
};
/**
- * Memory implementation of the most simple (stupid) particle stack object.
+ * @class CombinedStackImpl
+ *
+ * Memory implementation of a combined data stack.
+ *
+ * The two stack data user objects Stack1Impl and Stack2Impl are
+ * merged into one consistent Stack container object providing
+ * access to the combined number of data entries.
*/
template <typename Stack1Impl, typename Stack2Impl>
class CombinedStackImpl : public Stack1Impl, public Stack2Impl {
@@ -118,11 +159,17 @@ namespace corsika::stack {
Stack2Impl::DecrementSize();
}
- private:
- /// the actual memory to store particle data
-
}; // end class CombinedStackImpl
+ /**
+ * Helper template alias `CombinedStack` to construct new combined
+ * stack from two stack data objects and a particle readout interface.
+ *
+ * Note that the Stack2Impl provides only /additional/ data to
+ * Stack1Impl. This is important (see above) since tuple data for
+ * initialization are forwarded to Stack1Impl (first).
+ */
+
template <typename Stack1Impl, typename Stack2Impl, template <typename> typename _PI>
using CombinedStack = Stack<CombinedStackImpl<Stack1Impl, Stack2Impl>, _PI>;
diff --git a/Framework/StackInterface/ParticleBase.h b/Framework/StackInterface/ParticleBase.h
index b5704e762a667b8ff7d4c7442a068f561d2dee8a..04c9184ba5187111e44ffae0902a7352950498da 100644
--- a/Framework/StackInterface/ParticleBase.h
+++ b/Framework/StackInterface/ParticleBase.h
@@ -52,7 +52,6 @@ namespace corsika::stack {
ParticleBase() = default;
private:
- /*
// those copy constructors and assigments should never be implemented
ParticleBase(ParticleBase&) = delete;
ParticleBase operator=(ParticleBase&) = delete;
@@ -60,7 +59,7 @@ namespace corsika::stack {
ParticleBase operator=(ParticleBase&&) = delete;
ParticleBase(const ParticleBase&) = delete;
ParticleBase operator=(const ParticleBase&) = delete;
- */
+
public:
/**
* Delete this particle on the stack. The corresponding iterator
@@ -73,8 +72,8 @@ namespace corsika::stack {
* args is a variadic list of input data that has to match the
* function description in the user defined ParticleInterface::AddSecondary(...)
*/
- template <typename... Args>
- StackIterator AddSecondary(const Args... args) {
+ template <typename... TArgs>
+ StackIterator AddSecondary(const TArgs... args) {
return GetStack().AddSecondary(GetIterator(), args...);
}
@@ -88,7 +87,7 @@ namespace corsika::stack {
return static_cast<const StackIterator&>(*this);
}
- // protected:
+ protected:
/**
@name Access to underlying stack fData, these are service
function for user classes. User code can only rely on GetIndex
@@ -107,26 +106,6 @@ namespace corsika::stack {
///@}
};
- template <typename T>
- class ParticleBaseAdd {
-
- public:
- ParticleBaseAdd() = default;
-
- using T::GetIndex;
- using T::GetStackData;
-
- public:
- /*
- template <typename... Args1, typename... Args2>
- void SetParticleData(Args1... args1, Args2... args2) {
- T::SetParticleData(args1...);
- }
- template <typename... Args1, typename... Args2>
- void SetParticleData(T& p, Args1... args1, Args2... args2) {}
- */
- };
-
} // namespace corsika::stack
#endif
diff --git a/Framework/StackInterface/SecondaryView.h b/Framework/StackInterface/SecondaryView.h
index ac2e3ed19e462ecaca7e7def0bfd2abcb58e31ac..1bbbe02ccaad07fcc533d5782c753427c0dddbf1 100644
--- a/Framework/StackInterface/SecondaryView.h
+++ b/Framework/StackInterface/SecondaryView.h
@@ -14,6 +14,7 @@
#include <corsika/stack/Stack.h>
+#include <stdexcept>
#include <vector>
namespace corsika::stack {
@@ -68,10 +69,10 @@ namespace corsika::stack {
* the constructor of the SecondaryView class
* @{
*/
- using InnerStackTypeV = Stack<StackDataType, ParticleInterface>;
- typedef StackIteratorInterface<typename std::remove_reference<StackDataType>::type,
- ParticleInterface, InnerStackTypeV>
- StackIteratorV;
+ using InnerStackTypeValue = Stack<StackDataType, ParticleInterface>;
+ using StackIteratorValue =
+ StackIteratorInterface<typename std::remove_reference<StackDataType>::type,
+ ParticleInterface, InnerStackTypeValue>;
/// @}
public:
@@ -85,7 +86,8 @@ namespace corsika::stack {
/**
* this is the full type of the declared ParticleInterface: typedef typename
*/
- using ParticleType = typename StackIterator::ParticleInterfaceType;
+ using ParticleType = StackIterator;
+ using ParticleInterfaceType = typename StackIterator::ParticleInterfaceType;
friend class StackIteratorInterface<
typename std::remove_reference<StackDataType>::type, ParticleInterface, ViewType>;
@@ -99,25 +101,34 @@ namespace corsika::stack {
* new stack.
*/
template <typename... Args>
- SecondaryView(Args... args);
+ SecondaryView(Args... args) = delete;
public:
- SecondaryView(StackIteratorV& vP)
- : Stack<StackDataType&, ParticleInterface>(vP.GetStackData())
- , fProjectileIndex(vP.GetIndex()) {}
-
- auto GetProjectile() {
+ /**
+ SecondaryView can only be constructed passing it a valid
+ StackIterator to another Stack object
+ **/
+ SecondaryView(StackIteratorValue& vI)
+ : Stack<StackDataType&, ParticleInterface>(vI.GetStackData())
+ , fProjectileIndex(vI.GetIndex()) {}
+
+ StackIterator GetProjectile() {
// NOTE: 0 is special marker here for PROJECTILE, see GetIndexFromIterator
return StackIterator(*this, 0);
}
template <typename... Args>
auto AddSecondary(const Args... v) {
+ StackIterator proj = GetProjectile();
+ return AddSecondary(proj, v...);
+ }
+
+ template <typename... Args>
+ auto AddSecondary(StackIterator& proj, const Args... v) {
InnerStackType::GetStackData().IncrementSize();
const unsigned int idSec = GetSize();
const unsigned int index = InnerStackType::GetStackData().GetSize() - 1;
fIndices.push_back(index);
- StackIterator proj = GetProjectile();
// NOTE: "+1" is since "0" is special marker here for PROJECTILE, see
// GetIndexFromIterator
return StackIterator(*this, idSec + 1, proj, v...);
@@ -164,7 +175,7 @@ namespace corsika::stack {
/**
* need overwrite Stack::Delete, since we want to call SecondaryView::DeleteLast
*/
- void Delete(ParticleType p) { Delete(p.GetIterator()); }
+ void Delete(ParticleInterfaceType p) { Delete(p.GetIterator()); }
/**
* delete last particle on stack by decrementing stack size
@@ -201,6 +212,21 @@ namespace corsika::stack {
std::vector<unsigned int> fIndices;
};
+ /*
+ See Issue 161
+
+ unfortunately clang does not support this in the same way (yet) as
+ gcc, so we have to distinguish here. If clang cataches up, we
+ could remove the #if here and elsewhere. The gcc code is much more
+ generic and universal.
+ */
+#if not defined(__clang__) && defined(__GNUC__) || defined(__GNUG__)
+ template <typename S, template <typename> typename _PIType = S::template PIType>
+ struct MakeView {
+ using type = corsika::stack::SecondaryView<typename S::StackImpl, _PIType>;
+ };
+#endif
+
} // namespace corsika::stack
#endif
diff --git a/Framework/StackInterface/Stack.h b/Framework/StackInterface/Stack.h
index 72f3a72fe2d9c5078f893b3f815ed880fd515118..2f183af2a9d7b64def393ea698966cfdabe0d413 100644
--- a/Framework/StackInterface/Stack.h
+++ b/Framework/StackInterface/Stack.h
@@ -13,12 +13,12 @@
#define _include_Stack_h__
#include <corsika/stack/StackIteratorInterface.h>
+// must be after StackIteratorInterface
+#include <corsika/stack/SecondaryView.h>
#include <stdexcept>
#include <type_traits>
-#include <corsika/stack/SecondaryView.h>
-
/**
All classes around management of particles on a stack.
*/
@@ -85,6 +85,10 @@ namespace corsika::stack {
* This constructor takes any argument and passes it on to the
* StackDataType user class. If the user did not provide a suited
* constructor this will fail with an error message.
+ *
+ * Furthermore, this is disabled with enable_if for SecondaryView
+ * stacks, where the inner data container is always a reference
+ * and cannot be initialized here.
*/
template <
typename... Args, typename _StackDataType = StackDataType,
diff --git a/Framework/StackInterface/StackIteratorInterface.h b/Framework/StackInterface/StackIteratorInterface.h
index cf1b6d7b53efb14f3fbdd4283831e1f34fe751de..03c59cb44c1d9b7cc74f18fce20cd0ceaaac687f 100644
--- a/Framework/StackInterface/StackIteratorInterface.h
+++ b/Framework/StackInterface/StackIteratorInterface.h
@@ -14,8 +14,6 @@
#include <corsika/stack/ParticleBase.h>
-#include <type_traits>
-
namespace corsika::stack {
template <typename StackDataType, template <typename> typename ParticleInterface>
@@ -58,20 +56,20 @@ namespace corsika::stack {
ParticleInterface, in this example StackDataType::GetData(const unsigned int
vIndex).
- For an example see stack_example.cc, or the
+ For two examples see stack_example.cc, or the
corsika::processes::sibyll::SibStack class
*/
template <typename StackDataType, template <typename> typename ParticleInterface,
- typename StackType =
- Stack<StackDataType, ParticleInterface>> //, bool IsBase=true >
+ typename StackType = Stack<StackDataType, ParticleInterface>>
class StackIteratorInterface
- : public ParticleInterface<StackIteratorInterface<StackDataType, ParticleInterface,
- StackType>> { //,IsBase> {
+ : public ParticleInterface<
+ StackIteratorInterface<StackDataType, ParticleInterface, StackType>> {
public:
- using ParticleInterfaceType = ParticleInterface<
- StackIteratorInterface<StackDataType, ParticleInterface, StackType>>; //,IsBase>;
+ using ParticleInterfaceType =
+ ParticleInterface<corsika::stack::StackIteratorInterface<
+ StackDataType, ParticleInterface, StackType>>;
// friends are needed for access to protected methods
friend class Stack<StackDataType,
@@ -91,10 +89,20 @@ namespace corsika::stack {
StackIteratorInterface() = delete;
public:
+ StackIteratorInterface(StackIteratorInterface const& vR)
+ : fIndex(vR.fIndex)
+ , fData(vR.fData) {}
+
+ StackIteratorInterface& operator=(StackIteratorInterface const& vR) {
+ fIndex = vR.fIndex;
+ fData = vR.fData;
+ return *this;
+ }
+
/** iterator must always point to data, with an index:
- @param data reference to the stack [rw]
- @param index index on stack
- */
+ @param data reference to the stack [rw]
+ @param index index on stack
+ */
StackIteratorInterface(StackType& data, const unsigned int index)
: fIndex(index)
, fData(&data) {}
@@ -110,8 +118,7 @@ namespace corsika::stack {
StackIteratorInterface(StackType& data, const unsigned int index, const Args... args)
: fIndex(index)
, fData(&data) {
- ParticleInterfaceType& p = **this;
- p.SetParticleData(args...);
+ (**this).SetParticleData(args...);
}
/** constructor that also sets new values on particle data object, including reference
@@ -129,9 +136,7 @@ namespace corsika::stack {
StackIteratorInterface& parent, const Args... args)
: fIndex(index)
, fData(&data) {
- ParticleInterfaceType& p = **this;
- ParticleInterfaceType& pa = *parent;
- p.SetParticleData(pa, args...);
+ (**this).SetParticleData(*parent, args...);
}
public:
@@ -200,12 +205,12 @@ namespace corsika::stack {
template <typename StackDataType, template <typename> typename ParticleInterface,
typename StackType = Stack<StackDataType, ParticleInterface>>
class ConstStackIteratorInterface
- : public ParticleInterface<ConstStackIteratorInterface<
- StackDataType, ParticleInterface, StackType>> { //,IsBase> {
+ : public ParticleInterface<
+ ConstStackIteratorInterface<StackDataType, ParticleInterface, StackType>> {
public:
- typedef ParticleInterface<ConstStackIteratorInterface<
- StackDataType, ParticleInterface, StackType>> //,IsBase>
+ typedef ParticleInterface<
+ ConstStackIteratorInterface<StackDataType, ParticleInterface, StackType>>
ParticleInterfaceType;
friend class Stack<StackDataType, ParticleInterface>; // for access to GetIndex
diff --git a/Framework/StackInterface/testCombinedStack.cc b/Framework/StackInterface/testCombinedStack.cc
index 11a3314a2e60a517b9c1b197f913ffec1e7cc421..fbe4a031c928da4ea5a6c326ad1a98c903933972 100644
--- a/Framework/StackInterface/testCombinedStack.cc
+++ b/Framework/StackInterface/testCombinedStack.cc
@@ -10,6 +10,7 @@
*/
#include <corsika/stack/CombinedStack.h>
+#include <corsika/stack/SecondaryView.h>
#include <corsika/stack/Stack.h>
#include <testTestStack.h> // for testing: simple stack. This is a
@@ -31,6 +32,10 @@ using namespace corsika;
using namespace corsika::stack;
using namespace std;
+////////////////////////////////////////////////////////////
+// first level test: combine two stacks:
+// StackTest = (TestStackData + TestStackData2)
+
// definition of stack-data object
class TestStackData2 {
@@ -82,7 +87,7 @@ public:
double GetData2() const { return GetStackData().GetData2(GetIndex()); }
};
-// combined stack
+// combined stack: StackTest = (TestStackData + TestStackData2)
template <typename StackIter>
using CombinedTestInterfaceType =
corsika::stack::CombinedParticleInterface<TestParticleInterface,
@@ -196,6 +201,8 @@ TEST_CASE("Combined Stack", "[stack]") {
}
////////////////////////////////////////////////////////////
+// next level: combine three stacks:
+// combined stack: StackTest2 = ((TestStackData + TestStackData2) + TestStackData3)
// definition of stack-data object
class TestStackData3 {
@@ -228,6 +235,7 @@ private:
std::vector<double> fData3;
};
+// ---------------------------------------
// defintion of a stack-readout object, the iteractor dereference
// operator will deliver access to these function
template <typename T>
@@ -248,7 +256,7 @@ public:
double GetData3() const { return GetStackData().GetData3(GetIndex()); }
};
-// double combined stack
+// double combined stack:
// combined stack
template <typename StackIter>
using CombinedTestInterfaceType2 =
@@ -264,12 +272,14 @@ TEST_CASE("Combined Stack - multi", "[stack]") {
StackTest2 s;
REQUIRE(s.GetSize() == 0);
+ // add new particle, only provide tuple data for StackTest
auto p1 = s.AddParticle(std::tuple{9.9});
+ // add new particle, provide tuple data for both StackTest and TestStackData3
auto p2 = s.AddParticle(std::tuple{8.8}, std::tuple{0.1});
+ // examples to explicitly change data on stack
p2.SetData2(0.1); // not clear why this is needed, need to check
// SetParticleData workflow for more complicated
// settings
- // auto p3 = s.AddParticle( std::tuple {8.8}, std::tuple{1.}, std::tuple{0.1} );
p1.SetData3(20.2);
p2.SetData3(10.3);
@@ -301,3 +311,47 @@ TEST_CASE("Combined Stack - multi", "[stack]") {
REQUIRE(s.GetSize() == 0);
}
}
+
+////////////////////////////////////////////////////////////
+
+// final level test, create SecondaryView on StackTest2
+
+/*
+ See Issue 161
+
+ unfortunately clang does not support this in the same way (yet) as
+ gcc, so we have to distinguish here. If clang cataches up, we could
+ remove the clang branch here and also in corsika::Cascade. The gcc
+ code is much more generic and universal.
+ */
+template <typename StackIter>
+using CombinedTestInterfaceType2 =
+ corsika::stack::CombinedParticleInterface<StackTest::PIType, TestParticleInterface3,
+ StackIter>;
+
+using StackTest2 = CombinedStack<typename StackTest::StackImpl, TestStackData3,
+ CombinedTestInterfaceType2>;
+
+#if defined(__clang__)
+using StackTestView = SecondaryView<TestStackData, TestParticleInterface>;
+#elif defined(__GNUC__) || defined(__GNUG__)
+using StackTestView = corsika::stack::MakeView<StackTest2>::type;
+#endif
+
+using Particle2 = typename StackTest2::ParticleType;
+
+TEST_CASE("Combined Stack - secondary view") {
+
+ SECTION("create secondaries via secondaryview") {
+
+ StackTest2 stack;
+ auto particle = stack.AddParticle(std::tuple{9.9});
+ // cout << boost::typeindex::type_id_runtime(particle).pretty_name() << endl;
+ StackTestView view(particle);
+
+ auto projectile = view.GetProjectile();
+ projectile.AddSecondary(std::tuple{8.8});
+
+ REQUIRE(stack.GetSize() == 2);
+ }
+}
diff --git a/Framework/StackInterface/testSecondaryView.cc b/Framework/StackInterface/testSecondaryView.cc
index 08d27c169bfa6144abd203603eae74dc1e134c98..2d2064b64e2e09c96ec5f77ebbf25d50f5d60e9e 100644
--- a/Framework/StackInterface/testSecondaryView.cc
+++ b/Framework/StackInterface/testSecondaryView.cc
@@ -33,6 +33,22 @@ using namespace std;
typedef Stack<TestStackData, TestParticleInterface> StackTest;
+/*
+ See Issue 161
+
+ unfortunately clang does not support this in the same way (yet) as
+ gcc, so we have to distinguish here. If clang cataches up, we could
+ remove the clang branch here and also in corsika::Cascade. The gcc
+ code is much more generic and universal.
+ */
+#if defined(__clang__)
+using StackTestView = SecondaryView<TestStackData, TestParticleInterface>;
+#elif defined(__GNUC__) || defined(__GNUG__)
+using StackTestView = MakeView<StackTest>::type;
+#endif
+
+using Particle = typename StackTest::ParticleType;
+
TEST_CASE("SecondaryStack", "[stack]") {
// helper function for sum over stack data
@@ -51,55 +67,71 @@ TEST_CASE("SecondaryStack", "[stack]") {
auto particle = s.GetNextParticle();
- typedef SecondaryView<TestStackData, TestParticleInterface> StackTestView;
- StackTestView v(particle);
- REQUIRE(v.GetSize() == 0);
+ StackTestView view(particle);
+ REQUIRE(view.GetSize() == 0);
{
- auto proj = v.GetProjectile();
+ auto proj = view.GetProjectile();
REQUIRE(proj.GetData() == particle.GetData());
+ proj.AddSecondary(std::tuple{4.4});
}
- v.AddSecondary(std::tuple{4.4});
- v.AddSecondary(std::tuple{4.5});
- v.AddSecondary(std::tuple{4.6});
+ view.AddSecondary(std::tuple{4.5});
+ view.AddSecondary(std::tuple{4.6});
- REQUIRE(v.GetSize() == 3);
+ REQUIRE(view.GetSize() == 3);
REQUIRE(s.GetSize() == 5);
- REQUIRE(!v.IsEmpty());
+ REQUIRE(!view.IsEmpty());
auto sumView = [](const StackTestView& stack) {
- double v = 0;
- for (const auto& p : stack) { v += p.GetData(); }
- return v;
+ double value = 0;
+ for (const auto& p : stack) { value += p.GetData(); }
+ return value;
};
REQUIRE(sum(s) == sumS + 4.4 + 4.5 + 4.6);
- REQUIRE(sumView(v) == 4.4 + 4.5 + 4.6);
+ REQUIRE(sumView(view) == 4.4 + 4.5 + 4.6);
- v.DeleteLast();
- REQUIRE(v.GetSize() == 2);
+ view.DeleteLast();
+ REQUIRE(view.GetSize() == 2);
REQUIRE(s.GetSize() == 4);
REQUIRE(sum(s) == sumS + 4.4 + 4.5);
- REQUIRE(sumView(v) == 4.4 + 4.5);
+ REQUIRE(sumView(view) == 4.4 + 4.5);
- auto pDel = v.GetNextParticle();
- v.Delete(pDel);
- REQUIRE(v.GetSize() == 1);
+ auto pDel = view.GetNextParticle();
+ view.Delete(pDel);
+ REQUIRE(view.GetSize() == 1);
REQUIRE(s.GetSize() == 3);
REQUIRE(sum(s) == sumS + 4.4 + 4.5 - pDel.GetData());
- REQUIRE(sumView(v) == 4.4 + 4.5 - pDel.GetData());
+ REQUIRE(sumView(view) == 4.4 + 4.5 - pDel.GetData());
- v.Delete(v.GetNextParticle());
+ view.Delete(view.GetNextParticle());
REQUIRE(sum(s) == sumS);
- REQUIRE(sumView(v) == 0);
- REQUIRE(v.IsEmpty());
+ REQUIRE(sumView(view) == 0);
+ REQUIRE(view.IsEmpty());
{
- auto proj = v.GetProjectile();
+ auto proj = view.GetProjectile();
REQUIRE(proj.GetData() == particle.GetData());
}
}
+
+ SECTION("secondary view, construct from ParticleType") {
+ StackTest s;
+ REQUIRE(s.GetSize() == 0);
+ s.AddParticle(std::tuple{9.9});
+ s.AddParticle(std::tuple{8.8});
+
+ auto iterator = s.GetNextParticle();
+ typename StackTest::ParticleType& particle = iterator; // as in corsika::Cascade
+
+ StackTestView view(particle);
+ REQUIRE(view.GetSize() == 0);
+
+ view.AddSecondary(std::tuple{4.4});
+
+ REQUIRE(view.GetSize() == 1);
+ }
}
diff --git a/Processes/EnergyLoss/EnergyLoss.cc b/Processes/EnergyLoss/EnergyLoss.cc
index 960bbd3235e0de25b4898d5a68dd84813ca27ea4..dde5cf442c175d167cd14f605a04f4e5fd288736 100644
--- a/Processes/EnergyLoss/EnergyLoss.cc
+++ b/Processes/EnergyLoss/EnergyLoss.cc
@@ -148,7 +148,7 @@ namespace corsika::process::EnergyLoss {
(0.5 * log(aux) - beta2 - Cadj / Z - delta / 2 + barkas + bloch) * dX;
}
- process::EProcessReturn EnergyLoss::DoContinuous(Particle& p, Track& t, Stack&) {
+ process::EProcessReturn EnergyLoss::DoContinuous(Particle& p, Track& t) {
if (p.GetChargeNumber() == 0) return process::EProcessReturn::eOk;
GrammageType const dX =
p.GetNode()->GetModelProperties().IntegratedGrammage(t, t.GetLength());
@@ -178,16 +178,16 @@ namespace corsika::process::EnergyLoss {
return units::si::meter * std::numeric_limits<double>::infinity();
}
- void EnergyLoss::MomentumUpdate(corsika::setup::Stack::ParticleType& p,
+ void EnergyLoss::MomentumUpdate(corsika::setup::Stack::ParticleType& vP,
corsika::units::si::HEPEnergyType Enew) {
- HEPMomentumType Pnew = elab2plab(Enew, p.GetMass());
- auto pnew = p.GetMomentum();
- p.SetMomentum(pnew * Pnew / pnew.GetNorm());
+ HEPMomentumType Pnew = elab2plab(Enew, vP.GetMass());
+ auto pnew = vP.GetMomentum();
+ vP.SetMomentum(pnew * Pnew / pnew.GetNorm());
}
#include <corsika/geometry/CoordinateSystem.h>
- int EnergyLoss::GetXbin(corsika::setup::Stack::ParticleType& p,
+ int EnergyLoss::GetXbin(corsika::setup::Stack::ParticleType& vP,
const HEPEnergyType dE) {
using namespace corsika::geometry;
@@ -195,12 +195,12 @@ namespace corsika::process::EnergyLoss {
CoordinateSystem const& rootCS =
RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
Point pos1(rootCS, 0_m, 0_m, 0_m);
- Point pos2(rootCS, 0_m, 0_m, p.GetPosition().GetCoordinates()[2]);
+ Point pos2(rootCS, 0_m, 0_m, vP.GetPosition().GetCoordinates()[2]);
Vector delta = (pos2 - pos1) / 1_s;
Trajectory t(Line(pos1, delta), 1_s);
GrammageType const grammage =
- p.GetNode()->GetModelProperties().IntegratedGrammage(t, t.GetLength());
+ vP.GetNode()->GetModelProperties().IntegratedGrammage(t, t.GetLength());
const int bin = grammage / fdX;
diff --git a/Processes/EnergyLoss/EnergyLoss.h b/Processes/EnergyLoss/EnergyLoss.h
index 5eb2a6b6782b12bb7b336e67d43e2882d666b80a..323ef85735eb1525ef3514329bbc112d7c44cf97 100644
--- a/Processes/EnergyLoss/EnergyLoss.h
+++ b/Processes/EnergyLoss/EnergyLoss.h
@@ -35,8 +35,7 @@ namespace corsika::process::EnergyLoss {
void Init() {}
corsika::process::EProcessReturn DoContinuous(corsika::setup::Stack::ParticleType&,
- corsika::setup::Trajectory&,
- corsika::setup::Stack&);
+ corsika::setup::Trajectory&);
corsika::units::si::LengthType MaxStepLength(corsika::setup::Stack::ParticleType&,
corsika::setup::Trajectory&);
diff --git a/Processes/HadronicElasticModel/HadronicElasticModel.cc b/Processes/HadronicElasticModel/HadronicElasticModel.cc
index 289968645c4c758d10f7230dbcf99babba40838c..33f162e41531bf4008f91b21bb98e8afc13a7abb 100644
--- a/Processes/HadronicElasticModel/HadronicElasticModel.cc
+++ b/Processes/HadronicElasticModel/HadronicElasticModel.cc
@@ -79,7 +79,7 @@ namespace corsika::process::HadronicElasticModel {
}
template <>
- process::EProcessReturn HadronicElasticInteraction::DoInteraction(Particle& p, Stack&) {
+ process::EProcessReturn HadronicElasticInteraction::DoInteraction(Particle& p) {
if (p.GetPID() != particles::Code::Proton) { return process::EProcessReturn::eOk; }
using namespace units::si;
diff --git a/Processes/HadronicElasticModel/HadronicElasticModel.h b/Processes/HadronicElasticModel/HadronicElasticModel.h
index fae8c97242bac32426eca9089ac8a8ff611b631e..4180c04c514d00089fd044ab39027b9128a853ff 100644
--- a/Processes/HadronicElasticModel/HadronicElasticModel.h
+++ b/Processes/HadronicElasticModel/HadronicElasticModel.h
@@ -59,8 +59,8 @@ namespace corsika::process::HadronicElasticModel {
template <typename Particle, typename Track>
corsika::units::si::GrammageType GetInteractionLength(Particle const& p, Track&);
- template <typename Particle, typename Stack>
- corsika::process::EProcessReturn DoInteraction(Particle& p, Stack&);
+ template <typename Particle>
+ corsika::process::EProcessReturn DoInteraction(Particle&);
};
} // namespace corsika::process::HadronicElasticModel
diff --git a/Processes/NullModel/NullModel.cc b/Processes/NullModel/NullModel.cc
index c99a464b4dd4eca06e288df524202ec8016bd3e3..802ee5b4a669f824ac94c717e0bd61b3a0cfd3a3 100644
--- a/Processes/NullModel/NullModel.cc
+++ b/Processes/NullModel/NullModel.cc
@@ -14,21 +14,23 @@
#include <corsika/setup/SetupTrajectory.h>
using namespace corsika;
-using namespace corsika::process::null_model;
+namespace corsika::process::null_model {
-void corsika::process::null_model::NullModel::Init() {}
+ void NullModel::Init() {}
-NullModel::NullModel(units::si::LengthType maxStepLength)
- : fMaxStepLength(maxStepLength) {}
+ NullModel::NullModel(units::si::LengthType maxStepLength)
+ : fMaxStepLength(maxStepLength) {}
-template <>
-process::EProcessReturn NullModel::DoContinuous(setup::Stack::ParticleType&,
- setup::Trajectory&, setup::Stack&) const {
- return process::EProcessReturn::eOk;
-}
+ template <>
+ process::EProcessReturn NullModel::DoContinuous(setup::Stack::ParticleType&,
+ setup::Trajectory&) const {
+ return process::EProcessReturn::eOk;
+ }
-template <>
-units::si::LengthType NullModel::MaxStepLength(setup::Stack::ParticleType&,
- setup::Trajectory&) const {
- return fMaxStepLength;
-}
+ template <>
+ units::si::LengthType NullModel::MaxStepLength(setup::Stack::ParticleType&,
+ setup::Trajectory&) const {
+ return fMaxStepLength;
+ }
+
+} // namespace corsika::process::null_model
diff --git a/Processes/NullModel/NullModel.h b/Processes/NullModel/NullModel.h
index 5579d69f8a1ec08348e242e4d5bb8d012a21feb5..d3a820da41ede79dfd40669e7da6d15faaa8cf79 100644
--- a/Processes/NullModel/NullModel.h
+++ b/Processes/NullModel/NullModel.h
@@ -12,11 +12,12 @@
#ifndef _Physics_NullModel_NullModel_h_
#define _Physics_NullModel_NullModel_h_
-#include <corsika/process/ContinuousProcess.h>
+#include <corsika/process/BaseProcess.h>
+#include <corsika/units/PhysicalUnits.h>
namespace corsika::process::null_model {
- class NullModel : public corsika::process::ContinuousProcess<NullModel> {
+ class NullModel : public corsika::process::BaseProcess<NullModel> {
corsika::units::si::LengthType const fMaxStepLength;
public:
@@ -25,8 +26,8 @@ namespace corsika::process::null_model {
void Init();
- template <typename Particle, typename Track, typename Stack>
- process::EProcessReturn DoContinuous(Particle&, Track&, Stack&) const;
+ template <typename Particle, typename Track>
+ process::EProcessReturn DoContinuous(Particle&, Track&) const;
template <typename Particle, typename Track>
corsika::units::si::LengthType MaxStepLength(Particle&, Track&) const;
diff --git a/Processes/NullModel/testNullModel.cc b/Processes/NullModel/testNullModel.cc
index ed639eca4cb54fa3f1f6a676a27c01eb37663b6c..930d5979fec7bbcc728212b74e944f7dd11599cf 100644
--- a/Processes/NullModel/testNullModel.cc
+++ b/Processes/NullModel/testNullModel.cc
@@ -39,23 +39,19 @@ TEST_CASE("NullModel", "[processes]") {
geometry::Trajectory<geometry::Line> track(line, 10_s);
setup::Stack stack;
- setup::Stack::ParticleType
- // auto
- particle = stack.AddParticle(
- std::tuple<corsika::particles::Code, corsika::units::si::HEPEnergyType,
- corsika::stack::MomentumVector, corsika::geometry::Point,
- corsika::units::si::TimeType>{
- particles::Code::Electron, 1.5_GeV,
- stack::MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
- geometry::Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s});
-
+ setup::Stack::ParticleType particle = stack.AddParticle(
+ std::tuple<particles::Code, units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
+ particles::Code::Electron, 100_GeV,
+ corsika::stack::MomentumVector(dummyCS, {0_GeV, 0_GeV, -1_GeV}),
+ geometry::Point(dummyCS, {0_m, 0_m, 10_km}), 0_ns});
SECTION("interface") {
NullModel model(10_m);
model.Init();
[[maybe_unused]] const process::EProcessReturn ret =
- model.DoContinuous(particle, track, stack);
+ model.DoContinuous(particle, track);
LengthType const length = model.MaxStepLength(particle, track);
CHECK((length / 10_m) == Approx(1));
diff --git a/Processes/Pythia/Decay.cc b/Processes/Pythia/Decay.cc
index 77241b106dff5d6f674f60ec9a2acfe763ca09be..222c68c8dd57d7fd57c453b62f53777e85ad062a 100644
--- a/Processes/Pythia/Decay.cc
+++ b/Processes/Pythia/Decay.cc
@@ -22,7 +22,8 @@ using std::vector;
using namespace corsika;
using namespace corsika::setup;
-using Particle = Stack::StackIterator; // ParticleType;
+using Projectile = corsika::setup::StackView::ParticleType;
+using Particle = corsika::setup::Stack::ParticleType;
using Track = Trajectory;
namespace corsika::process::pythia {
@@ -84,13 +85,13 @@ namespace corsika::process::pythia {
}
template <>
- void Decay::DoDecay(Particle& p, Stack&) {
+ void Decay::DoDecay(Projectile& vP) {
using geometry::Point;
using namespace units;
using namespace units::si;
- auto const decayPoint = p.GetPosition();
- auto const t0 = p.GetTime();
+ auto const decayPoint = vP.GetPosition();
+ auto const t0 = vP.GetTime();
// coordinate system, get global frame of reference
geometry::CoordinateSystem& rootCS =
@@ -103,17 +104,17 @@ namespace corsika::process::pythia {
event.reset();
// set particle unstable
- Decay::SetUnstable(p.GetPID());
+ Decay::SetUnstable(vP.GetPID());
// input particle PDG
- auto const pdgCode = static_cast<int>(particles::GetPDG(p.GetPID()));
+ auto const pdgCode = static_cast<int>(particles::GetPDG(vP.GetPID()));
- auto const pcomp = p.GetMomentum().GetComponents();
+ auto const pcomp = vP.GetMomentum().GetComponents();
double px = pcomp[0] / 1_GeV;
double py = pcomp[1] / 1_GeV;
double pz = pcomp[2] / 1_GeV;
- double en = p.GetEnergy() / 1_GeV;
- double m = particles::GetMass(p.GetPID()) / 1_GeV;
+ double en = vP.GetEnergy() / 1_GeV;
+ double m = particles::GetMass(vP.GetPID()) / 1_GeV;
// add particle to pythia stack
event.append(pdgCode, 1, 0, 0, px, py, pz, en, m);
@@ -138,17 +139,17 @@ namespace corsika::process::pythia {
cout << "particle: id=" << pyId << " momentum=" << pyP.GetComponents() / 1_GeV
<< " energy=" << pyEn << endl;
- p.AddSecondary(
+ vP.AddSecondary(
tuple<particles::Code, units::si::HEPEnergyType,
corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
pyId, pyEn, pyP, decayPoint, t0});
}
// set particle stable
- Decay::SetStable(p.GetPID());
+ Decay::SetStable(vP.GetPID());
// remove original particle from corsika stack
- p.Delete();
+ vP.Delete();
// if (fCount>10) throw std::runtime_error("stop here");
}
diff --git a/Processes/Pythia/Decay.h b/Processes/Pythia/Decay.h
index fc76137160177afabe0d6624b2bc1b708f545ab4..6506625b07a04ec076919bafcbf0e8c92395aaf4 100644
--- a/Processes/Pythia/Decay.h
+++ b/Processes/Pythia/Decay.h
@@ -35,11 +35,11 @@ namespace corsika::process {
void SetUnstable(const corsika::particles::Code);
void SetStable(const corsika::particles::Code);
- template <typename Particle>
- corsika::units::si::TimeType GetLifetime(Particle const& p);
+ template <typename TParticle>
+ corsika::units::si::TimeType GetLifetime(TParticle const&);
- template <typename Particle, typename Stack>
- void DoDecay(Particle& p, Stack&);
+ template <typename TProjectile>
+ void DoDecay(TProjectile&);
private:
Pythia8::Pythia fPythia;
diff --git a/Processes/Pythia/Interaction.cc b/Processes/Pythia/Interaction.cc
index 092563e52bb66cf8130604100ed6a3f7e5b1c18d..ae8104350086f5cfda94b174305d0f0366faa2ce 100644
--- a/Processes/Pythia/Interaction.cc
+++ b/Processes/Pythia/Interaction.cc
@@ -26,7 +26,8 @@ using std::tuple;
using namespace corsika;
using namespace corsika::setup;
-using Particle = Stack::StackIterator; // ParticleType;
+using Projectile = corsika::setup::StackView::ParticleType;
+using Particle = corsika::setup::Stack::ParticleType;
using Track = Trajectory;
namespace corsika::process::pythia {
@@ -224,8 +225,9 @@ namespace corsika::process::pythia {
i++;
cout << "Interaction: get interaction length for target: " << targetId << endl;
- [[maybe_unused]] auto const [productionCrossSection, elaCrossSection] =
+ auto const [productionCrossSection, elaCrossSection] =
GetCrossSection(corsikaBeamId, targetId, ECoM);
+ [[maybe_unused]] const auto& dummy_elaCrossSection = elaCrossSection;
cout << "Interaction: IntLength: pythia return (mb): "
<< productionCrossSection / 1_mb << endl
@@ -257,14 +259,14 @@ namespace corsika::process::pythia {
*/
template <>
- process::EProcessReturn Interaction::DoInteraction(Particle& p, Stack&) {
+ process::EProcessReturn Interaction::DoInteraction(Projectile& vP) {
using namespace units;
using namespace utl;
using namespace units::si;
using namespace geometry;
- const auto corsikaBeamId = p.GetPID();
+ const auto corsikaBeamId = vP.GetPID();
cout << "Pythia::Interaction: "
<< "DoInteraction: " << corsikaBeamId << " interaction? "
<< process::pythia::Interaction::CanInteract(corsikaBeamId) << endl;
@@ -280,8 +282,8 @@ namespace corsika::process::pythia {
RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
// position and time of interaction, not used in Sibyll
- Point pOrig = p.GetPosition();
- TimeType tOrig = p.GetTime();
+ Point pOrig = vP.GetPosition();
+ TimeType tOrig = vP.GetTime();
// define target
// FOR NOW: target is always at rest
@@ -290,8 +292,8 @@ namespace corsika::process::pythia {
const FourVector PtargLab(eTargetLab, pTargetLab);
// define projectile
- HEPEnergyType const eProjectileLab = p.GetEnergy();
- auto const pProjectileLab = p.GetMomentum();
+ HEPEnergyType const eProjectileLab = vP.GetEnergy();
+ auto const pProjectileLab = vP.GetMomentum();
cout << "Interaction: ebeam lab: " << eProjectileLab / 1_GeV << endl
<< "Interaction: pbeam lab: " << pProjectileLab.GetComponents() / 1_GeV
@@ -326,12 +328,12 @@ namespace corsika::process::pythia {
cout << "Interaction: time: " << tOrig << endl;
HEPEnergyType Etot = eProjectileLab + eTargetLab;
- MomentumVector Ptot = p.GetMomentum();
+ MomentumVector Ptot = vP.GetMomentum();
// invariant mass, i.e. cm. energy
HEPEnergyType Ecm = sqrt(Etot * Etot - Ptot.squaredNorm());
// sample target mass number
- const auto* currentNode = p.GetNode();
+ const auto* currentNode = vP.GetNode();
const auto& mediumComposition =
currentNode->GetModelProperties().GetNuclearComposition();
// get cross sections for target materials
@@ -345,8 +347,8 @@ namespace corsika::process::pythia {
for (size_t i = 0; i < compVec.size(); ++i) {
auto const targetId = compVec[i];
- [[maybe_unused]] const auto [sigProd, sigEla] =
- GetCrossSection(corsikaBeamId, targetId, Ecm);
+ const auto [sigProd, sigEla] = GetCrossSection(corsikaBeamId, targetId, Ecm);
+ [[maybe_unused]] const auto& dummy_sigEla = sigEla;
cross_section_of_components[i] = sigProd;
}
@@ -385,19 +387,19 @@ namespace corsika::process::pythia {
MomentumVector Plab_final(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
HEPEnergyType Elab_final = 0_GeV;
for (int i = 0; i < event.size(); ++i) {
- Pythia8::Particle& pp = event[i];
+ Pythia8::Particle& p8p = event[i];
// skip particles that have decayed in pythia
- if (!pp.isFinal()) continue;
+ if (!p8p.isFinal()) continue;
auto const pyId =
- particles::ConvertFromPDG(static_cast<particles::PDGCode>(pp.id()));
+ particles::ConvertFromPDG(static_cast<particles::PDGCode>(p8p.id()));
const MomentumVector pyPlab(
- rootCS, {pp.px() * 1_GeV, pp.py() * 1_GeV, pp.pz() * 1_GeV});
- HEPEnergyType const pyEn = pp.e() * 1_GeV;
+ rootCS, {p8p.px() * 1_GeV, p8p.py() * 1_GeV, p8p.pz() * 1_GeV});
+ HEPEnergyType const pyEn = p8p.e() * 1_GeV;
// add to corsika stack
- auto pnew = p.AddSecondary(
+ auto pnew = vP.AddSecondary(
tuple<particles::Code, units::si::HEPEnergyType, stack::MomentumVector,
geometry::Point, units::si::TimeType>{pyId, pyEn, pyPlab, pOrig,
tOrig});
@@ -410,7 +412,7 @@ namespace corsika::process::pythia {
<< ", Plab_final=" << (Plab_final / 1_GeV).GetComponents() << endl;
}
// delete current particle
- p.Delete();
+ vP.Delete();
}
return process::EProcessReturn::eOk;
}
diff --git a/Processes/Pythia/Interaction.h b/Processes/Pythia/Interaction.h
index 3e1d8813b22a55ffef913ea8298414a8d441ec7d..ac27825acf0561d0d9e7e8e8b66961bdfc674e2f 100644
--- a/Processes/Pythia/Interaction.h
+++ b/Processes/Pythia/Interaction.h
@@ -51,16 +51,16 @@ namespace corsika::process::pythia {
const corsika::particles::Code TargetId,
const corsika::units::si::HEPEnergyType CoMenergy);
- template <typename Particle, typename Track>
- corsika::units::si::GrammageType GetInteractionLength(Particle&, Track&);
+ template <typename TParticle, typename TTrack>
+ corsika::units::si::GrammageType GetInteractionLength(TParticle&, TTrack&);
/**
In this function PYTHIA is called to produce one event. The
event is copied (and boosted) into the shower lab frame.
*/
- template <typename Particle, typename Stack>
- corsika::process::EProcessReturn DoInteraction(Particle&, Stack&);
+ template <typename TProjctile>
+ corsika::process::EProcessReturn DoInteraction(TProjctile&);
private:
corsika::random::RNG& fRNG =
diff --git a/Processes/Pythia/testPythia.cc b/Processes/Pythia/testPythia.cc
index 1e6c57cca49e7fb53017776a8d4fb69790c9b278..3fb2df13ecd1c24e674997a3c3aa468164f21d06 100644
--- a/Processes/Pythia/testPythia.cc
+++ b/Processes/Pythia/testPythia.cc
@@ -145,11 +145,12 @@ TEST_CASE("pythia process") {
random::RNGManager::GetInstance().RegisterRandomStream("pythia");
- process::pythia::Decay model(particleList);
+ corsika::stack::SecondaryView view(particle);
+ auto projectile = view.GetProjectile();
+ process::pythia::Decay model(particleList);
model.Init();
- /*[[maybe_unused]] const process::EProcessReturn ret =*/model.DoDecay(particle,
- stack);
+ /*[[maybe_unused]] const process::EProcessReturn ret =*/model.DoDecay(projectile);
[[maybe_unused]] const TimeType time = model.GetLifetime(particle);
}
@@ -166,11 +167,12 @@ TEST_CASE("pythia process") {
corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
particles::Code::PiPlus, E0, plab, pos, 0_ns});
particle.SetNode(nodePtr);
- process::pythia::Interaction model;
+ corsika::stack::SecondaryView view(particle);
+ auto projectile = view.GetProjectile();
+ process::pythia::Interaction model;
model.Init();
- /*[[maybe_unused]] const process::EProcessReturn ret =*/model.DoInteraction(particle,
- stack);
+ [[maybe_unused]] const process::EProcessReturn ret = model.DoInteraction(projectile);
[[maybe_unused]] const GrammageType length =
model.GetInteractionLength(particle, track);
}
diff --git a/Processes/Sibyll/Decay.cc b/Processes/Sibyll/Decay.cc
index 3cd009c580d10845ea2219b59322dc1855a88295..7a865e92fa08d606c68c5bad782e0a70e826f101 100644
--- a/Processes/Sibyll/Decay.cc
+++ b/Processes/Sibyll/Decay.cc
@@ -24,7 +24,8 @@ using std::vector;
using namespace corsika;
using namespace corsika::setup;
-using Particle = Stack::StackIterator; // ParticleType;
+using Projectile = corsika::setup::StackView::ParticleType;
+using Particle = corsika::setup::Stack::ParticleType;
using Track = Trajectory;
namespace corsika::process::sibyll {
@@ -103,27 +104,27 @@ namespace corsika::process::sibyll {
}
template <>
- units::si::TimeType Decay::GetLifetime(Particle const& p) {
+ units::si::TimeType Decay::GetLifetime(Particle const& vP) {
using namespace units::si;
- HEPEnergyType E = p.GetEnergy();
- HEPMassType m = p.GetMass();
+ HEPEnergyType E = vP.GetEnergy();
+ HEPMassType m = vP.GetMass();
const double gamma = E / m;
- const TimeType t0 = particles::GetLifetime(p.GetPID());
+ const TimeType t0 = particles::GetLifetime(vP.GetPID());
auto const lifetime = gamma * t0;
const auto mkin =
- (E * E - p.GetMomentum().squaredNorm()); // delta_mass(p.GetMomentum(), E, m);
- cout << "Decay: code: " << p.GetPID() << endl;
+ (E * E - vP.GetMomentum().squaredNorm()); // delta_mass(vP.GetMomentum(), E, m);
+ cout << "Decay: code: " << vP.GetPID() << endl;
cout << "Decay: MinStep: t0: " << t0 << endl;
cout << "Decay: MinStep: energy: " << E / 1_GeV << " GeV" << endl;
- cout << "Decay: momentum: " << p.GetMomentum().GetComponents() / 1_GeV << " GeV"
+ cout << "Decay: momentum: " << vP.GetMomentum().GetComponents() / 1_GeV << " GeV"
<< endl;
cout << "Decay: momentum: shell mass-kin. inv. mass " << mkin / 1_GeV / 1_GeV << " "
<< m / 1_GeV * m / 1_GeV << endl;
- auto sib_id = process::sibyll::ConvertToSibyllRaw(p.GetPID());
+ auto sib_id = process::sibyll::ConvertToSibyllRaw(vP.GetPID());
cout << "Decay: sib mass: " << get_sibyll_mass2(sib_id) << endl;
cout << "Decay: MinStep: gamma: " << gamma << endl;
cout << "Decay: MinStep: tau: " << lifetime << endl;
@@ -132,23 +133,23 @@ namespace corsika::process::sibyll {
}
template <>
- void Decay::DoDecay(Particle& p, Stack&) {
+ void Decay::DoDecay(Projectile& vP) {
using geometry::Point;
using namespace units::si;
fCount++;
SibStack ss;
ss.Clear();
- const particles::Code pCode = p.GetPID();
+ const particles::Code pCode = vP.GetPID();
// copy particle to sibyll stack
- ss.AddParticle(process::sibyll::ConvertToSibyllRaw(pCode), p.GetEnergy(),
- p.GetMomentum(),
+ ss.AddParticle(process::sibyll::ConvertToSibyllRaw(pCode), vP.GetEnergy(),
+ vP.GetMomentum(),
// setting particle mass with Corsika values, may be inconsistent
// with sibyll internal values
particles::GetMass(pCode));
// remember position
- Point const decayPoint = p.GetPosition();
- TimeType const t0 = p.GetTime();
+ Point const decayPoint = vP.GetPosition();
+ TimeType const t0 = vP.GetTime();
// set all particles/hadrons unstable
// setHadronsUnstable();
SetUnstable(pCode);
@@ -166,7 +167,7 @@ namespace corsika::process::sibyll {
// FOR NOW: skip particles that have decayed in Sibyll, move to iterator?
if (psib.HasDecayed()) continue;
// add to corsika stack
- p.AddSecondary(
+ vP.AddSecondary(
tuple<particles::Code, units::si::HEPEnergyType, corsika::stack::MomentumVector,
geometry::Point, units::si::TimeType>{
process::sibyll::ConvertFromSibyll(psib.GetPID()), psib.GetEnergy(),
@@ -174,8 +175,6 @@ namespace corsika::process::sibyll {
}
// empty sibyll stack
ss.Clear();
- // remove original particle from corsika stack
- p.Delete();
}
} // namespace corsika::process::sibyll
diff --git a/Processes/Sibyll/Decay.h b/Processes/Sibyll/Decay.h
index 7e7a202e7c3172577c6a5662d09eb67bbb9433f2..5dc06bd0cfc6db1a59c74f04dd78971b3f01a748 100644
--- a/Processes/Sibyll/Decay.h
+++ b/Processes/Sibyll/Decay.h
@@ -37,10 +37,10 @@ namespace corsika::process {
void SetHadronsUnstable();
template <typename Particle>
- corsika::units::si::TimeType GetLifetime(Particle const& p);
+ corsika::units::si::TimeType GetLifetime(Particle const&);
- template <typename Particle, typename Stack>
- void DoDecay(Particle& p, Stack&);
+ template <typename Projectile>
+ void DoDecay(Projectile&);
};
} // namespace sibyll
} // namespace corsika::process
diff --git a/Processes/Sibyll/Interaction.cc b/Processes/Sibyll/Interaction.cc
index 333b6c9894a9f1e7e84a77fecaa283afe9afe047..ed550a012d63df6b89586198933ce10c4d01b8d2 100644
--- a/Processes/Sibyll/Interaction.cc
+++ b/Processes/Sibyll/Interaction.cc
@@ -27,8 +27,11 @@ using std::cout;
using std::endl;
using std::tuple;
-using SetupParticle = corsika::setup::Stack::StackIterator; // SetupParticleType;
-using Track = corsika::setup::Trajectory;
+using namespace corsika;
+using namespace corsika::setup;
+using Particle = Stack::StackIterator; // ParticleType;
+using Projectile = StackView::StackIterator; // ParticleType;
+using Track = Trajectory;
namespace corsika::process::sibyll {
@@ -115,8 +118,7 @@ namespace corsika::process::sibyll {
}
template <>
- units::si::GrammageType Interaction::GetInteractionLength(SetupParticle& p,
- Track&) const {
+ units::si::GrammageType Interaction::GetInteractionLength(Particle& vP, Track&) const {
using namespace units;
using namespace units::si;
@@ -126,7 +128,7 @@ namespace corsika::process::sibyll {
CoordinateSystem& rootCS =
RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
- const particles::Code corsikaBeamId = p.GetPID();
+ const particles::Code corsikaBeamId = vP.GetPID();
// beam particles for sibyll : 1, 2, 3 for p, pi, k
// read from cross section code table
@@ -136,9 +138,9 @@ namespace corsika::process::sibyll {
MomentumVector pTarget(rootCS, {0_GeV, 0_GeV, 0_GeV});
// total momentum and energy
- HEPEnergyType Elab = p.GetEnergy() + constants::nucleonMass;
+ HEPEnergyType Elab = vP.GetEnergy() + constants::nucleonMass;
MomentumVector pTotLab(rootCS, {0_GeV, 0_GeV, 0_GeV});
- pTotLab += p.GetMomentum();
+ pTotLab += vP.GetMomentum();
pTotLab += pTarget;
auto const pTotLabNorm = pTotLab.norm();
// calculate cm. energy
@@ -146,9 +148,9 @@ namespace corsika::process::sibyll {
(Elab + pTotLabNorm) * (Elab - pTotLabNorm)); // binomial for numerical accuracy
cout << "Interaction: LambdaInt: \n"
- << " input energy: " << p.GetEnergy() / 1_GeV << endl
+ << " input energy: " << vP.GetEnergy() / 1_GeV << endl
<< " beam can interact:" << kInteraction << endl
- << " beam pid:" << p.GetPID() << endl;
+ << " beam pid:" << vP.GetPID() << endl;
// TODO: move limits into variables
// FR: removed && Elab >= 8.5_GeV
@@ -160,7 +162,7 @@ namespace corsika::process::sibyll {
ideally as full particle object so that the four momenta
and the boosts can be defined..
*/
- auto const* currentNode = p.GetNode();
+ auto const* currentNode = vP.GetNode();
const auto mediumComposition =
currentNode->GetModelProperties().GetNuclearComposition();
// determine average interaction length
@@ -174,8 +176,9 @@ namespace corsika::process::sibyll {
i++;
cout << "Interaction: get interaction length for target: " << targetId << endl;
- [[maybe_unused]] auto const [productionCrossSection, elaCrossSection] =
+ auto const [productionCrossSection, elaCrossSection] =
GetCrossSection(corsikaBeamId, targetId, ECoM);
+ [[maybe_unused]] const auto& dummy_elaCX = elaCrossSection;
cout << "Interaction: "
<< " IntLength: sibyll return (mb): " << productionCrossSection / 1_mb
@@ -206,14 +209,14 @@ namespace corsika::process::sibyll {
*/
template <>
- process::EProcessReturn Interaction::DoInteraction(SetupParticle& p, setup::Stack&) {
+ process::EProcessReturn Interaction::DoInteraction(Projectile& vP) {
using namespace units;
using namespace utl;
using namespace units::si;
using namespace geometry;
- const auto corsikaBeamId = p.GetPID();
+ const auto corsikaBeamId = vP.GetPID();
cout << "ProcessSibyll: "
<< "DoInteraction: " << corsikaBeamId << " interaction? "
<< process::sibyll::CanInteract(corsikaBeamId) << endl;
@@ -228,8 +231,8 @@ namespace corsika::process::sibyll {
RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
// position and time of interaction, not used in Sibyll
- Point pOrig = p.GetPosition();
- TimeType tOrig = p.GetTime();
+ Point pOrig = vP.GetPosition();
+ TimeType tOrig = vP.GetTime();
// define target
// for Sibyll is always a single nucleon
@@ -239,8 +242,8 @@ namespace corsika::process::sibyll {
const FourVector PtargLab(eTargetLab, pTargetLab);
// define projectile
- HEPEnergyType const eProjectileLab = p.GetEnergy();
- auto const pProjectileLab = p.GetMomentum();
+ HEPEnergyType const eProjectileLab = vP.GetEnergy();
+ auto const pProjectileLab = vP.GetMomentum();
cout << "Interaction: ebeam lab: " << eProjectileLab / 1_GeV << endl
<< "Interaction: pbeam lab: " << pProjectileLab.GetComponents() / 1_GeV
@@ -275,12 +278,12 @@ namespace corsika::process::sibyll {
cout << "Interaction: time: " << tOrig << endl;
HEPEnergyType Etot = eProjectileLab + eTargetLab;
- MomentumVector Ptot = p.GetMomentum();
+ MomentumVector Ptot = vP.GetMomentum();
// invariant mass, i.e. cm. energy
HEPEnergyType Ecm = sqrt(Etot * Etot - Ptot.squaredNorm());
// sample target mass number
- auto const* currentNode = p.GetNode();
+ auto const* currentNode = vP.GetNode();
auto const& mediumComposition =
currentNode->GetModelProperties().GetNuclearComposition();
// get cross sections for target materials
@@ -294,8 +297,8 @@ namespace corsika::process::sibyll {
for (size_t i = 0; i < compVec.size(); ++i) {
auto const targetId = compVec[i];
- [[maybe_unsused]] const auto [sigProd, sigEla] =
- GetCrossSection(corsikaBeamId, targetId, Ecm);
+ const auto [sigProd, sigEla] = GetCrossSection(corsikaBeamId, targetId, Ecm);
+ [[maybe_unused]] const auto& dummy_sigEla = sigEla;
cross_section_of_components[i] = sigProd;
}
@@ -330,7 +333,6 @@ namespace corsika::process::sibyll {
<< " DoInteraction: should have dropped particle.. "
<< "THIS IS AN ERROR" << endl;
throw std::runtime_error("energy too low for SIBYLL");
- // p.Delete(); delete later... different process
} else {
fCount++;
// Sibyll does not know about units..
@@ -361,7 +363,7 @@ namespace corsika::process::sibyll {
auto const Plab = boost.fromCoM(FourVector(eCoM, pCoM));
// add to corsika stack
- auto pnew = p.AddSecondary(
+ auto pnew = vP.AddSecondary(
tuple<particles::Code, units::si::HEPEnergyType, stack::MomentumVector,
geometry::Point, units::si::TimeType>{
process::sibyll::ConvertFromSibyll(psib.GetPID()),
@@ -378,7 +380,7 @@ namespace corsika::process::sibyll {
}
}
// delete current particle
- p.Delete();
+ vP.Delete();
return process::EProcessReturn::eOk;
}
diff --git a/Processes/Sibyll/Interaction.h b/Processes/Sibyll/Interaction.h
index d2274d428e7819bb929b5a199e20ae2743da5a02..4e1899f685e9254c3207850c227328340a051e94 100644
--- a/Processes/Sibyll/Interaction.h
+++ b/Processes/Sibyll/Interaction.h
@@ -61,8 +61,8 @@ namespace corsika::process::sibyll {
event is copied (and boosted) into the shower lab frame.
*/
- template <typename Particle, typename Stack>
- corsika::process::EProcessReturn DoInteraction(Particle&, Stack&);
+ template <typename TProjectile>
+ corsika::process::EProcessReturn DoInteraction(TProjectile&);
private:
corsika::random::RNG& fRNG =
diff --git a/Processes/Sibyll/NuclearInteraction.cc b/Processes/Sibyll/NuclearInteraction.cc
index 52c585916e5b274f5700d2e8cde726a4c664117f..31abd9da350e1c80c24868eb2313e860b51bd821 100644
--- a/Processes/Sibyll/NuclearInteraction.cc
+++ b/Processes/Sibyll/NuclearInteraction.cc
@@ -29,10 +29,11 @@ using std::endl;
using std::tuple;
using std::vector;
-using corsika::setup::SetupEnvironment;
-using SetupStack = corsika::setup::Stack;
-using Particle = SetupStack::StackIterator; // ParticleType;
-using Track = corsika::setup::Trajectory;
+using namespace corsika;
+using namespace corsika::setup;
+using Particle = Stack::ParticleType; // StackIterator; // ParticleType;
+using Projectile = StackView::StackIterator; // StackView::ParticleType;
+using Track = Trajectory;
namespace corsika::process::sibyll {
@@ -177,15 +178,15 @@ namespace corsika::process::sibyll {
template <>
template <>
tuple<units::si::CrossSectionType, units::si::CrossSectionType>
- NuclearInteraction<SetupEnvironment>::GetCrossSection(Particle& p,
+ NuclearInteraction<SetupEnvironment>::GetCrossSection(Particle& vP,
const particles::Code TargetId) {
using namespace units::si;
- if (p.GetPID() != particles::Code::Nucleus)
+ if (vP.GetPID() != particles::Code::Nucleus)
throw std::runtime_error(
"NuclearInteraction: GetCrossSection: particle not a nucleus!");
- auto const iBeamA = p.GetNuclearA();
- HEPEnergyType LabEnergyPerNuc = p.GetEnergy() / iBeamA;
+ auto const iBeamA = vP.GetNuclearA();
+ HEPEnergyType LabEnergyPerNuc = vP.GetEnergy() / iBeamA;
cout << "NuclearInteraction: GetCrossSection: called with: beamNuclA= " << iBeamA
<< " TargetId= " << TargetId << " LabEnergyPerNuc= " << LabEnergyPerNuc / 1_GeV
<< endl;
@@ -216,7 +217,7 @@ namespace corsika::process::sibyll {
template <>
template <>
units::si::GrammageType NuclearInteraction<SetupEnvironment>::GetInteractionLength(
- Particle& p, Track&) {
+ Particle& vP, Track&) {
using namespace units;
using namespace units::si;
@@ -226,7 +227,7 @@ namespace corsika::process::sibyll {
CoordinateSystem& rootCS =
RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
- const particles::Code corsikaBeamId = p.GetPID();
+ const particles::Code corsikaBeamId = vP.GetPID();
if (!particles::IsNucleus(corsikaBeamId)) {
// no nuclear interaction
return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
@@ -243,12 +244,12 @@ namespace corsika::process::sibyll {
corsika::stack::MomentumVector pTarget(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
// total momentum and energy
- HEPEnergyType Elab = p.GetEnergy() + constants::nucleonMass;
- int const nuclA = p.GetNuclearA();
- auto const ElabNuc = p.GetEnergy() / nuclA;
+ HEPEnergyType Elab = vP.GetEnergy() + constants::nucleonMass;
+ int const nuclA = vP.GetNuclearA();
+ auto const ElabNuc = vP.GetEnergy() / nuclA;
corsika::stack::MomentumVector pTotLab(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
- pTotLab += p.GetMomentum();
+ pTotLab += vP.GetMomentum();
pTotLab += pTarget;
auto const pTotLabNorm = pTotLab.norm();
// calculate cm. energy
@@ -275,7 +276,7 @@ namespace corsika::process::sibyll {
ideally as full particle object so that the four momenta
and the boosts can be defined..
*/
- auto const* const currentNode = p.GetNode();
+ auto const* const currentNode = vP.GetNode();
auto const& mediumComposition =
currentNode->GetModelProperties().GetNuclearComposition();
// determine average interaction length
@@ -289,8 +290,9 @@ namespace corsika::process::sibyll {
i++;
cout << "NuclearInteraction: get interaction length for target: " << targetId
<< endl;
- [[maybe_unused]] auto const [productionCrossSection, elaCrossSection] =
- GetCrossSection(p, targetId);
+ auto const [productionCrossSection, elaCrossSection] =
+ GetCrossSection(vP, targetId);
+ [[maybe_unused]] auto& dummy_elaCrossSection = elaCrossSection;
cout << "NuclearInteraction: "
<< "IntLength: nuclib return (mb): " << productionCrossSection / 1_mb
@@ -317,7 +319,7 @@ namespace corsika::process::sibyll {
template <>
template <>
process::EProcessReturn NuclearInteraction<SetupEnvironment>::DoInteraction(
- Particle& p, SetupStack& s) {
+ Projectile& vP) {
// this routine superimposes different nucleon-nucleon interactions
// in a nucleus-nucleus interaction, based the SIBYLL routine SIBNUC
@@ -327,9 +329,9 @@ namespace corsika::process::sibyll {
using namespace units::si;
using namespace geometry;
- const auto ProjId = p.GetPID();
+ const auto ProjId = vP.GetPID();
// TODO: calculate projectile mass in nuclearStackExtension
- // const auto ProjMass = p.GetMass();
+ // const auto ProjMass = vP.GetMass();
cout << "NuclearInteraction: DoInteraction: called with:" << ProjId << endl;
if (!IsNucleus(ProjId)) {
@@ -346,8 +348,8 @@ namespace corsika::process::sibyll {
"should use NuclearStackExtension!");
auto const ProjMass =
- p.GetNuclearZ() * particles::Proton::GetMass() +
- (p.GetNuclearA() - p.GetNuclearZ()) * particles::Neutron::GetMass();
+ vP.GetNuclearZ() * particles::Proton::GetMass() +
+ (vP.GetNuclearA() - vP.GetNuclearZ()) * particles::Neutron::GetMass();
cout << "NuclearInteraction: projectile mass: " << ProjMass / 1_GeV << endl;
fCount++;
@@ -356,21 +358,21 @@ namespace corsika::process::sibyll {
RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
// position and time of interaction, not used in NUCLIB
- Point pOrig = p.GetPosition();
- TimeType tOrig = p.GetTime();
+ Point pOrig = vP.GetPosition();
+ TimeType tOrig = vP.GetTime();
cout << "Interaction: position of interaction: " << pOrig.GetCoordinates() << endl;
cout << "Interaction: time: " << tOrig << endl;
// projectile nucleon number
- const int kAProj = p.GetNuclearA(); // GetNucleusA(ProjId);
+ const int kAProj = vP.GetNuclearA(); // GetNucleusA(ProjId);
if (kAProj > GetMaxNucleusAProjectile())
throw std::runtime_error("Projectile nucleus too large for NUCLIB!");
// kinematics
// define projectile nucleus
- HEPEnergyType const eProjectileLab = p.GetEnergy();
- auto const pProjectileLab = p.GetMomentum();
+ HEPEnergyType const eProjectileLab = vP.GetEnergy();
+ auto const pProjectileLab = vP.GetMomentum();
const FourVector PprojLab(eProjectileLab, pProjectileLab);
cout << "NuclearInteraction: eProj lab: " << eProjectileLab / 1_GeV << endl
@@ -378,8 +380,8 @@ namespace corsika::process::sibyll {
<< endl;
// define projectile nucleon
- HEPEnergyType const eProjectileNucLab = p.GetEnergy() / kAProj;
- auto const pProjectileNucLab = p.GetMomentum() / kAProj;
+ HEPEnergyType const eProjectileNucLab = vP.GetEnergy() / kAProj;
+ auto const pProjectileNucLab = vP.GetMomentum() / kAProj;
const FourVector PprojNucLab(eProjectileNucLab, pProjectileNucLab);
cout << "NuclearInteraction: eProjNucleon lab: " << eProjectileNucLab / 1_GeV << endl
@@ -431,7 +433,7 @@ namespace corsika::process::sibyll {
//
// proton stand-in for nucleon
const auto beamId = particles::Proton::GetCode();
- auto const* const currentNode = p.GetNode();
+ auto const* const currentNode = vP.GetNode();
const auto& mediumComposition =
currentNode->GetModelProperties().GetNuclearComposition();
cout << "get nucleon-nucleus cross sections for target materials.." << endl;
@@ -556,15 +558,16 @@ namespace corsika::process::sibyll {
if (nuclA == 1)
// add nucleon
- p.AddSecondary(tuple<particles::Code, units::si::HEPEnergyType,
- stack::MomentumVector, geometry::Point, units::si::TimeType>{
- specCode, Plab.GetTimeLikeComponent(), Plab.GetSpaceLikeComponents(), pOrig,
- tOrig});
+ vP.AddSecondary(
+ tuple<particles::Code, units::si::HEPEnergyType, stack::MomentumVector,
+ geometry::Point, units::si::TimeType>{
+ specCode, Plab.GetTimeLikeComponent(), Plab.GetSpaceLikeComponents(),
+ pOrig, tOrig});
else
// add nucleus
- p.AddSecondary(tuple<particles::Code, units::si::HEPEnergyType,
- corsika::stack::MomentumVector, geometry::Point,
- units::si::TimeType, unsigned short, unsigned short>{
+ vP.AddSecondary(tuple<particles::Code, units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, geometry::Point,
+ units::si::TimeType, unsigned short, unsigned short>{
specCode, Plab.GetTimeLikeComponent(), Plab.GetSpaceLikeComponents(), pOrig,
tOrig, nuclA, nuclZ});
}
@@ -583,7 +586,7 @@ namespace corsika::process::sibyll {
const double mass_ratio = particles::GetMass(elaNucCode) / ProjMass;
auto const Plab = PprojLab * mass_ratio;
- p.AddSecondary(
+ vP.AddSecondary(
tuple<particles::Code, units::si::HEPEnergyType, corsika::stack::MomentumVector,
geometry::Point, units::si::TimeType>{
elaNucCode, Plab.GetTimeLikeComponent(), Plab.GetSpaceLikeComponents(),
@@ -597,19 +600,16 @@ namespace corsika::process::sibyll {
auto pCode = particles::Proton::GetCode();
// temporarily add to stack, will be removed after interaction in DoInteraction
cout << "inelastic interaction no. " << j << endl;
- auto inelasticNucleon = p.AddSecondary(
+ auto inelasticNucleon = vP.AddSecondary(
tuple<particles::Code, units::si::HEPEnergyType, corsika::stack::MomentumVector,
geometry::Point, units::si::TimeType>{
pCode, PprojNucLab.GetTimeLikeComponent(),
PprojNucLab.GetSpaceLikeComponents(), pOrig, tOrig});
// create inelastic interaction
cout << "calling HadronicInteraction..." << endl;
- fHadronicInteraction.DoInteraction(inelasticNucleon, s);
+ fHadronicInteraction.DoInteraction(inelasticNucleon);
}
- // delete parent particle
- p.Delete();
-
cout << "NuclearInteraction: DoInteraction: done" << endl;
return process::EProcessReturn::eOk;
diff --git a/Processes/Sibyll/NuclearInteraction.h b/Processes/Sibyll/NuclearInteraction.h
index b04bba538ec9cf2472f8d102fe4ec3bce070578d..6bee9ba7a7964e74112c2bf60e5cbd4b1c23c2fa 100644
--- a/Processes/Sibyll/NuclearInteraction.h
+++ b/Processes/Sibyll/NuclearInteraction.h
@@ -54,10 +54,10 @@ namespace corsika::process::sibyll {
GetCrossSection(Particle& p, const corsika::particles::Code TargetId);
template <typename Particle, typename Track>
- corsika::units::si::GrammageType GetInteractionLength(Particle& p, Track&);
+ corsika::units::si::GrammageType GetInteractionLength(Particle&, Track&);
- template <typename Particle, typename Stack>
- corsika::process::EProcessReturn DoInteraction(Particle& p, Stack& s);
+ template <typename Projectle>
+ corsika::process::EProcessReturn DoInteraction(Projectle&);
private:
TEnvironment const& fEnvironment;
diff --git a/Processes/Sibyll/testSibyll.cc b/Processes/Sibyll/testSibyll.cc
index 6acb2b4938282453a94997f69f9cdc01ad99e475..d5f5d9bc8e39ac775de3fb9d541a18b596941bad 100644
--- a/Processes/Sibyll/testSibyll.cc
+++ b/Processes/Sibyll/testSibyll.cc
@@ -119,12 +119,13 @@ TEST_CASE("SibyllInterface", "[processes]") {
corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
particles::Code::Proton, E0, plab, pos, 0_ns});
particle.SetNode(nodePtr);
+ corsika::stack::SecondaryView view(particle);
+ auto projectile = view.GetProjectile();
Interaction model;
model.Init();
- [[maybe_unused]] const process::EProcessReturn ret =
- model.DoInteraction(particle, stack);
+ [[maybe_unused]] const process::EProcessReturn ret = model.DoInteraction(projectile);
[[maybe_unused]] const GrammageType length =
model.GetInteractionLength(particle, track);
}
@@ -144,13 +145,14 @@ TEST_CASE("SibyllInterface", "[processes]") {
units::si::TimeType, unsigned short, unsigned short>{
particles::Code::Nucleus, E0, plab, pos, 0_ns, 4, 2});
particle.SetNode(nodePtr);
+ corsika::stack::SecondaryView view(particle);
+ auto projectile = view.GetProjectile();
Interaction hmodel;
NuclearInteraction model(hmodel, env);
model.Init();
- [[maybe_unused]] const process::EProcessReturn ret =
- model.DoInteraction(particle, stack);
+ [[maybe_unused]] const process::EProcessReturn ret = model.DoInteraction(projectile);
[[maybe_unused]] const GrammageType length =
model.GetInteractionLength(particle, track);
}
@@ -167,6 +169,8 @@ TEST_CASE("SibyllInterface", "[processes]") {
std::tuple<particles::Code, units::si::HEPEnergyType,
corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
particles::Code::Proton, E0, plab, pos, 0_ns});
+ corsika::stack::SecondaryView view(particle);
+ auto projectile = view.GetProjectile();
const std::vector<particles::Code> particleList = {
particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus,
@@ -175,8 +179,7 @@ TEST_CASE("SibyllInterface", "[processes]") {
Decay model(particleList);
model.Init();
- /*[[maybe_unused]] const process::EProcessReturn ret =*/model.DoDecay(particle,
- stack);
+ /*[[maybe_unused]] const process::EProcessReturn ret =*/model.DoDecay(projectile);
[[maybe_unused]] const TimeType time = model.GetLifetime(particle);
}
}
diff --git a/Processes/StackInspector/CMakeLists.txt b/Processes/StackInspector/CMakeLists.txt
index f00f44c2ab6b5513630f5344545f988b8e78dc3f..f8297639672f215e7a34e2330bd8e8387454d185 100644
--- a/Processes/StackInspector/CMakeLists.txt
+++ b/Processes/StackInspector/CMakeLists.txt
@@ -57,7 +57,6 @@ add_executable (testStackInspector testStackInspector.cc)
target_link_libraries (
testStackInspector
ProcessStackInspector
- CORSIKAsetup
CORSIKAgeometry
CORSIKAunits
CORSIKAthirdparty # for catch2
diff --git a/Processes/StackInspector/StackInspector.cc b/Processes/StackInspector/StackInspector.cc
index a22742165bec0a6f4a9bfe3b06e33e25ca2e9c3c..814763d5cbfcbc0bd4201db70f220dc4e8ea123f 100644
--- a/Processes/StackInspector/StackInspector.cc
+++ b/Processes/StackInspector/StackInspector.cc
@@ -16,7 +16,6 @@
#include <corsika/logging/Logger.h>
-#include <corsika/cascade/testCascade.h>
#include <corsika/setup/SetupTrajectory.h>
#include <iostream>
@@ -28,23 +27,22 @@ using namespace corsika::particles;
using namespace corsika::units::si;
using namespace corsika::process::stack_inspector;
-template <typename Stack>
-StackInspector<Stack>::StackInspector(const bool aReport)
- : fReport(aReport)
+template <typename TStack>
+StackInspector<TStack>::StackInspector(const int nStep, const bool aReport)
+ : StackProcess<StackInspector<TStack>>(nStep)
+ , fReport(aReport)
, fCountStep(0) {}
-template <typename Stack>
-StackInspector<Stack>::~StackInspector() {}
-
-template <typename Stack>
-process::EProcessReturn StackInspector<Stack>::DoContinuous(Particle&, setup::Trajectory&,
- Stack& s) {
+template <typename TStack>
+StackInspector<TStack>::~StackInspector() {}
+template <typename TStack>
+process::EProcessReturn StackInspector<TStack>::DoStack(TStack& vS) {
if (!fReport) return process::EProcessReturn::eOk;
[[maybe_unused]] int i = 0;
HEPEnergyType Etot = 0_GeV;
- for (auto& iterP : s) {
+ for (auto& iterP : vS) {
HEPEnergyType E = iterP.GetEnergy();
Etot += E;
geometry::CoordinateSystem& rootCS = geometry::RootCoordinateSystem::GetInstance()
@@ -53,27 +51,21 @@ process::EProcessReturn StackInspector<Stack>::DoContinuous(Particle&, setup::Tr
cout << "StackInspector: i=" << setw(5) << fixed << (i++) << ", id=" << setw(30)
<< iterP.GetPID() << " E=" << setw(15) << scientific << (E / 1_GeV) << " GeV, "
<< " pos=" << pos << " node = " << iterP.GetNode();
- // if (iterP.GetPID()==Code::Nucleus)
- // cout << " nuc_ref=" << iterP.GetNucleusRef();
+ if (iterP.GetPID() == Code::Nucleus) cout << " nuc_ref=" << iterP.GetNucleusRef();
cout << endl;
}
fCountStep++;
- cout << "StackInspector: nStep=" << fCountStep << " stackSize=" << s.GetSize()
+ cout << "StackInspector: nStep=" << fCountStep << " stackSize=" << vS.GetSize()
<< " Estack=" << Etot / 1_GeV << " GeV" << endl;
return process::EProcessReturn::eOk;
}
-template <typename Stack>
-corsika::units::si::LengthType StackInspector<Stack>::MaxStepLength(Particle&,
- setup::Trajectory&) {
- return std::numeric_limits<double>::infinity() * meter;
-}
-
-template <typename Stack>
-void StackInspector<Stack>::Init() {
+template <typename TStack>
+void StackInspector<TStack>::Init() {
fCountStep = 0;
}
+#include <corsika/cascade/testCascade.h>
#include <corsika/setup/SetupStack.h>
template class process::stack_inspector::StackInspector<setup::Stack>;
diff --git a/Processes/StackInspector/StackInspector.h b/Processes/StackInspector/StackInspector.h
index 9488767d780c9b2747ace740a14b21aa7da341a4..57a929dd2362b302d43d05175878a257c5c8a832 100644
--- a/Processes/StackInspector/StackInspector.h
+++ b/Processes/StackInspector/StackInspector.h
@@ -12,7 +12,7 @@
#ifndef _Physics_StackInspector_StackInspector_h_
#define _Physics_StackInspector_StackInspector_h_
-#include <corsika/process/ContinuousProcess.h>
+#include <corsika/process/StackProcess.h>
#include <corsika/setup/SetupTrajectory.h>
#include <corsika/units/PhysicalUnits.h>
@@ -20,20 +20,17 @@ namespace corsika::process {
namespace stack_inspector {
- template <typename Stack>
- class StackInspector
- : public corsika::process::ContinuousProcess<StackInspector<Stack>> {
+ template <typename TStack>
+ class StackInspector : public corsika::process::StackProcess<StackInspector<TStack>> {
- typedef typename Stack::ParticleType Particle;
+ typedef typename TStack::ParticleType Particle;
public:
- StackInspector(const bool aReport);
+ StackInspector(const int nStep, const bool aReport);
~StackInspector();
void Init();
- EProcessReturn DoContinuous(Particle&, corsika::setup::Trajectory&, Stack& s);
- corsika::units::si::LengthType MaxStepLength(Particle&,
- corsika::setup::Trajectory&);
+ EProcessReturn DoStack(TStack&);
private:
bool fReport;
diff --git a/Processes/StackInspector/testStackInspector.cc b/Processes/StackInspector/testStackInspector.cc
index 11e9cff3e4eae58513694dd3b64ba6a6ea64bead..832bea00a8074214ef9f3a235a6a574f7286cf02 100644
--- a/Processes/StackInspector/testStackInspector.cc
+++ b/Processes/StackInspector/testStackInspector.cc
@@ -21,12 +21,12 @@
#include <corsika/units/PhysicalUnits.h>
-#include <corsika/setup/SetupStack.h>
-#include <corsika/setup/SetupTrajectory.h>
+#include <corsika/cascade/testCascade.h>
using namespace corsika::units::si;
using namespace corsika::process::stack_inspector;
using namespace corsika;
+using namespace corsika::geometry;
TEST_CASE("StackInspector", "[processes]") {
@@ -38,21 +38,21 @@ TEST_CASE("StackInspector", "[processes]") {
geometry::Line line(origin, v);
geometry::Trajectory<geometry::Line> track(line, 10_s);
- setup::Stack stack;
- auto particle = stack.AddParticle(
+ TestCascadeStack stack;
+ stack.Clear();
+ HEPEnergyType E0 = 100_GeV;
+ stack.AddParticle(
std::tuple<particles::Code, units::si::HEPEnergyType,
corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
- particles::Code::Electron, 10_GeV,
+ particles::Code::Electron, E0,
corsika::stack::MomentumVector(rootCS, {0_GeV, 0_GeV, -1_GeV}),
- geometry::Point(rootCS, {0_m, 0_m, 10_km}), 0_ns});
+ Point(rootCS, {0_m, 0_m, 10_km}), 0_ns});
SECTION("interface") {
- StackInspector<setup::Stack> model(true);
+ StackInspector<TestCascadeStack> model(1, true);
model.Init();
- [[maybe_unused]] const process::EProcessReturn ret =
- model.DoContinuous(particle, track, stack);
- [[maybe_unused]] const LengthType length = model.MaxStepLength(particle, track);
+ [[maybe_unused]] const process::EProcessReturn ret = model.DoStack(stack);
}
}
diff --git a/Processes/TrackWriter/TrackWriter.cc b/Processes/TrackWriter/TrackWriter.cc
index b58a40fbbeca28b8e295ea246901253d674bef8a..5ea8588d9f19aca5f430aa7e749e7d3a51005e59 100644
--- a/Processes/TrackWriter/TrackWriter.cc
+++ b/Processes/TrackWriter/TrackWriter.cc
@@ -33,13 +33,13 @@ namespace corsika::process::TrackWriter {
}
template <>
- process::EProcessReturn TrackWriter::DoContinuous(Particle& p, Track& t, Stack&) {
+ process::EProcessReturn TrackWriter::DoContinuous(Particle& vP, Track& vT) {
using namespace units::si;
- auto const start = t.GetPosition(0).GetCoordinates();
- auto const delta = t.GetPosition(1).GetCoordinates() - start;
- auto const& name = particles::GetName(p.GetPID());
+ auto const start = vT.GetPosition(0).GetCoordinates();
+ auto const delta = vT.GetPosition(1).GetCoordinates() - start;
+ auto const& name = particles::GetName(vP.GetPID());
- fFile << name << " " << p.GetEnergy() / 1_eV << ' ' << start[0] / 1_m << ' '
+ fFile << name << " " << vP.GetEnergy() / 1_eV << ' ' << start[0] / 1_m << ' '
<< start[1] / 1_m << ' ' << start[2] / 1_m << " " << delta[0] / 1_m << ' '
<< delta[1] / 1_m << ' ' << delta[2] / 1_m << '\n';
diff --git a/Processes/TrackWriter/TrackWriter.h b/Processes/TrackWriter/TrackWriter.h
index fa2430569240b159f0c8cb94dd06bc314472be05..51fd5a157da9955aed9a06c69a1b6270626cfb5d 100644
--- a/Processes/TrackWriter/TrackWriter.h
+++ b/Processes/TrackWriter/TrackWriter.h
@@ -28,8 +28,8 @@ namespace corsika::process::TrackWriter {
void Init();
- template <typename Particle, typename Track, typename Stack>
- corsika::process::EProcessReturn DoContinuous(Particle& p, Track& t, Stack&);
+ template <typename Particle, typename Track>
+ corsika::process::EProcessReturn DoContinuous(Particle&, Track&);
template <typename Particle, typename Track>
corsika::units::si::LengthType MaxStepLength(Particle&, Track&);
diff --git a/Processes/TrackingLine/TrackingLine.h b/Processes/TrackingLine/TrackingLine.h
index e847a3fd7ebef23583c632cc5b71af474c883428..a801a2ee2baa206424211151c6a7fb9a616d0b70 100644
--- a/Processes/TrackingLine/TrackingLine.h
+++ b/Processes/TrackingLine/TrackingLine.h
@@ -105,7 +105,8 @@ namespace corsika::process {
currentLogicalVolumeNode->GetVolume());
// for the moment we are a bit bold here and assume
// everything is a sphere, crashes with exception if not
- auto const [t1, t2] = *TimeOfIntersection(line, sphere);
+ [[maybe_unused]] auto const [t1, t2] = *TimeOfIntersection(line, sphere);
+ [[maybe_unused]] auto dummy_t1 = t1;
intersections.emplace_back(t2, currentLogicalVolumeNode->GetParent());
}
diff --git a/Processes/TrackingLine/testTrackingLine.cc b/Processes/TrackingLine/testTrackingLine.cc
index fb213f9d8e3729015a8ad2f19c092e3d476bf9dc..1acff46eec0d810138203856b0bb6900e153ce87 100644
--- a/Processes/TrackingLine/testTrackingLine.cc
+++ b/Processes/TrackingLine/testTrackingLine.cc
@@ -91,6 +91,8 @@ TEST_CASE("TrackingLine") {
Line line(origin, v);
auto const [traj, geomMaxLength, nextVol] = tracking.GetTrack(p);
+ [[maybe_unused]] auto dummy_geomMaxLength = geomMaxLength;
+ [[maybe_unused]] auto dummy_nextVol = nextVol;
REQUIRE((traj.GetPosition(1.) - Point(cs, 0_m, 0_m, radius))
.GetComponents(cs)
diff --git a/Setup/SetupStack.h b/Setup/SetupStack.h
index dfcfa757dc2d0ad183a1ae9919bd603fbe8bb250..a59ca2aa3815b5a956e17a07faabf61e2834d972 100644
--- a/Setup/SetupStack.h
+++ b/Setup/SetupStack.h
@@ -29,6 +29,12 @@
// definition of stack-data object to store geometry information
template <typename TEnvType>
+
+/**
+ * @class GeometryData
+ *
+ * definition of stack-data object to store geometry information
+ */
class GeometryData {
public:
@@ -57,8 +63,14 @@ private:
std::vector<const BaseNodeType*> fNode;
};
+/**
+ * @class GeometryDataInterface
+ *
+ * corresponding defintion of a stack-readout object, the iteractor
+ * dereference operator will deliver access to these function
// defintion of a stack-readout object, the iteractor dereference
// operator will deliver access to these function
+ */
template <typename T, typename TEnvType>
class GeometryDataInterface : public T {
@@ -120,6 +132,25 @@ namespace corsika::setup {
// this is the REAL stack we use:
using Stack = detail::StackWithGeometry;
+ /*
+ See Issue 161
+
+ unfortunately clang does not support this in the same way (yet) as
+ gcc, so we have to distinguish here. If clang cataches up, we
+ could remove the clang branch here and also in
+ corsika::Cascade. The gcc code is much more generic and
+ universal. If we could do the gcc version, we won't had to define
+ StackView globally, we could do it with MakeView whereever it is
+ actually needed. Keep an eye on this!
+ */
+#if defined(__clang__)
+ using StackView =
+ corsika::stack::SecondaryView<typename corsika::setup::Stack::StackImpl,
+ corsika::setup::detail::StackWithGeometryInterface>;
+#elif defined(__GNUC__) || defined(__GNUG__)
+ using StackView = corsika::stack::MakeView<corsika::setup::Stack>::type;
+#endif
+
} // namespace corsika::setup
#endif
diff --git a/Stack/NuclearStackExtension/CMakeLists.txt b/Stack/NuclearStackExtension/CMakeLists.txt
index f2f34cfcfa4cf6108651282b464c476de9c125a5..6032c0b90218790d6622a0474af371c62687b711 100644
--- a/Stack/NuclearStackExtension/CMakeLists.txt
+++ b/Stack/NuclearStackExtension/CMakeLists.txt
@@ -40,8 +40,7 @@ target_link_libraries (
testNuclearStackExtension
SuperStupidStack
NuclearStackExtension
- # CORSIKAutls
- ProcessStackInspector
+ CORSIKAparticles
CORSIKAgeometry
CORSIKAunits
CORSIKAthirdparty # for catch2
diff --git a/Stack/NuclearStackExtension/NuclearStackExtension.h b/Stack/NuclearStackExtension/NuclearStackExtension.h
index 6c1fd2b801120bbe0c22c6c76d22ccadc652f4af..40eb5e0e672883e58cc207324e99d11afdd6be18 100644
--- a/Stack/NuclearStackExtension/NuclearStackExtension.h
+++ b/Stack/NuclearStackExtension/NuclearStackExtension.h
@@ -178,6 +178,8 @@ namespace corsika::stack {
return InnerParticleInterface<StackIteratorInterface>::GetChargeNumber();
}
+ int GetNucleusRef() const { return GetStackData().GetNucleusRef(GetIndex()); }
+
protected:
void SetNucleusRef(const int vR) { GetStackData().SetNucleusRef(GetIndex(), vR); }
};
diff --git a/Stack/NuclearStackExtension/testNuclearStackExtension.cc b/Stack/NuclearStackExtension/testNuclearStackExtension.cc
index 5185c3538f7b915dd825b8fe1ec98c751903c2a2..c92e5471fdf53e74aa6fdf3e0312e8f0492ae07e 100644
--- a/Stack/NuclearStackExtension/testNuclearStackExtension.cc
+++ b/Stack/NuclearStackExtension/testNuclearStackExtension.cc
@@ -29,10 +29,11 @@ using namespace corsika::units::si;
// this is an auxiliary help typedef, which I don't know how to put
// into NuclearStackExtension.h where it belongs...
template <typename StackIter>
-using ExtendedParticleInterfaceType = stack::nuclear_extension::NuclearParticleInterface<
- stack::super_stupid::SuperStupidStack::PIType, StackIter>;
+using ExtendedParticleInterfaceType =
+ corsika::stack::nuclear_extension::NuclearParticleInterface<
+ corsika::stack::super_stupid::SuperStupidStack::template PIType, StackIter>;
-using ExtStack = NuclearStackExtension<stack::super_stupid::SuperStupidStack,
+using ExtStack = NuclearStackExtension<corsika::stack::super_stupid::SuperStupidStack,
ExtendedParticleInterfaceType>;
#include <iostream>
diff --git a/Stack/SuperStupidStack/CMakeLists.txt b/Stack/SuperStupidStack/CMakeLists.txt
index de96527d4a0fb757ac3e79b61ef4ff7481981fa1..d787f0d033ce0977a96c875412685c4dc1f1d38c 100644
--- a/Stack/SuperStupidStack/CMakeLists.txt
+++ b/Stack/SuperStupidStack/CMakeLists.txt
@@ -38,9 +38,8 @@ add_executable (
target_link_libraries (
testSuperStupidStack
-# CORSIKAutls
- ProcessStackInspector
CORSIKAgeometry
+ CORSIKAparticles
CORSIKAunits
CORSIKAthirdparty # for catch2
)
diff --git a/Stack/SuperStupidStack/SuperStupidStack.h b/Stack/SuperStupidStack/SuperStupidStack.h
index 1610aede7ec0cad53a8e5e92d14ffba5ce1ab973..0b67b95f184bba24358e10980ededdf1dc86c551 100644
--- a/Stack/SuperStupidStack/SuperStupidStack.h
+++ b/Stack/SuperStupidStack/SuperStupidStack.h
@@ -39,6 +39,8 @@ namespace corsika::stack {
protected:
using corsika::stack::ParticleBase<StackIteratorInterface>::GetStack;
using corsika::stack::ParticleBase<StackIteratorInterface>::GetStackData;
+
+ public:
using corsika::stack::ParticleBase<StackIteratorInterface>::GetIndex;
public: