diff --git a/Documentation/Examples/CMakeLists.txt b/Documentation/Examples/CMakeLists.txt
index 131476a3a7c133e7d1b093a05a67ca830494be39..ffc14426f790c8c760d77fdc056d782836bc1a94 100644
--- a/Documentation/Examples/CMakeLists.txt
+++ b/Documentation/Examples/CMakeLists.txt
@@ -31,7 +31,7 @@ target_link_libraries (cascade_example SuperStupidStack CORSIKAunits
ProcessPythia
CORSIKAcascade
ProcessEnergyLoss
-# ProcessStackInspector
+ ProcessStackInspector
ProcessTrackWriter
ProcessTrackingLine
CORSIKAprocesses
@@ -71,7 +71,6 @@ target_link_libraries (cascade_proton_example SuperStupidStack CORSIKAunits
ProcessPythia
CORSIKAcascade
ProcessEnergyLoss
-# ProcessStackInspector
ProcessTrackWriter
ProcessTrackingLine
ProcessHadronicElasticModel
@@ -82,6 +81,7 @@ target_link_libraries (cascade_proton_example SuperStupidStack CORSIKAunits
CORSIKAenvironment
CORSIKAprocesssequence
)
+
install (TARGETS cascade_proton_example DESTINATION share/examples)
CORSIKA_ADD_TEST (cascade_proton_example)
@@ -94,7 +94,6 @@ target_link_libraries (vertical_EAS SuperStupidStack CORSIKAunits
ProcessPythia
CORSIKAcascade
ProcessEnergyLoss
-# ProcessStackInspector
ProcessTrackWriter
ProcessTrackingLine
ProcessHadronicElasticModel
@@ -106,7 +105,7 @@ target_link_libraries (vertical_EAS SuperStupidStack CORSIKAunits
CORSIKAprocesssequence
)
install (TARGETS vertical_EAS DESTINATION share/examples)
-CORSIKA_ADD_TEST (vertical_EAS)
+# 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
diff --git a/Documentation/Examples/boundary_example.cc b/Documentation/Examples/boundary_example.cc
index 0c767160e08123292379ad330941c2f37479b627..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());
@@ -282,15 +273,16 @@ int main() {
random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
process::sibyll::Interaction sibyll;
- process::sibyll::NuclearInteraction sibyllNuc(sibyll);
- process::sibyll::Decay decay;
+ process::sibyll::Decay decay{{particles::Code::PiPlus, particles::Code::PiMinus,
+ particles::Code::KPlus, particles::Code::KMinus,
+ particles::Code::K0Long, particles::Code::K0Short}};
ProcessCut cut(20_GeV);
process::TrackWriter::TrackWriter trackWriter("tracks.dat");
MyBoundaryCrossingProcess<true> boundaryCrossing("crossings.dat");
// assemble all processes into an ordered process list
- auto sequence = sibyll << sibyllNuc << decay << cut << boundaryCrossing << trackWriter;
+ auto sequence = sibyll << decay << cut << boundaryCrossing << trackWriter;
// setup particle stack, and add primary particles
setup::Stack stack;
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index aa70f3d62dbb880951985bc54754c2017a8ec0e9..66b1aa0b7f51dd51e0e6eff8790b86f448d59913 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -12,6 +12,7 @@
#include <corsika/cascade/Cascade.h>
#include <corsika/process/ProcessSequence.h>
#include <corsika/process/energy_loss/EnergyLoss.h>
+#include <corsika/process/stack_inspector/StackInspector.h>
#include <corsika/process/tracking_line/TrackingLine.h>
#include <corsika/setup/SetupEnvironment.h>
@@ -40,7 +41,6 @@
#include <iostream>
#include <limits>
-#include <typeinfo>
using namespace corsika;
using namespace corsika::process;
@@ -238,7 +238,7 @@ public:
int main() {
const LengthType height_atmosphere = 112.8_km;
-
+
feenableexcept(FE_INVALID);
// initialize random number sequence(s)
random::RNGManager::GetInstance().RegisterRandomStream("cascade");
@@ -273,8 +273,8 @@ int main() {
universe.AddChild(std::move(outerMedium));
// setup processes, decays and interactions
- tracking_line::TrackingLine<setup::Stack, setup::Trajectory> tracking(env);
- // stack_inspector::StackInspector<setup::Stack> stackInspect(true);
+ 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,
@@ -282,17 +282,18 @@ int main() {
random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
random::RNGManager::GetInstance().RegisterRandomStream("pythia");
- process::sibyll::Interaction sibyll(env);
- process::sibyll::NuclearInteraction sibyllNuc(env, sibyll);
+ process::sibyll::Interaction sibyll;
+ process::sibyll::NuclearInteraction sibyllNuc(sibyll, env);
process::sibyll::Decay decay(trackedHadrons);
ProcessCut cut(20_GeV);
- ObservationLevel obsLevel(height_atmosphere - 2000_m);//1400_m);
+ ObservationLevel obsLevel(height_atmosphere - 2000_m); // 1400_m);
process::TrackWriter::TrackWriter trackWriter("tracks.dat");
process::EnergyLoss::EnergyLoss eLoss;
// assemble all processes into an ordered process list
- auto sequence = sibyll << sibyllNuc << decay << eLoss << cut << trackWriter; // << obsLevel
+ auto sequence = stackInspect << sibyll << sibyllNuc << decay << eLoss << cut
+ << trackWriter;
// setup particle stack, and add primary particle
setup::Stack stack;
@@ -301,7 +302,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.;
diff --git a/Documentation/Examples/cascade_proton_example.cc b/Documentation/Examples/cascade_proton_example.cc
index 7cbfda4c0c4e1379418b8f39157e143ece4a335c..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,27 +219,26 @@ int main() {
random::RNGManager::GetInstance().RegisterRandomStream("cascade");
// setup environment, geometry
- environment::Environment env;
+ using EnvType = Environment<setup::IEnvironmentModel>;
+ EnvType env;
auto& universe = *(env.GetUniverse());
- auto theMedium = environment::Environment::CreateNode<Sphere>(
- Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m},
- 1_km * std::numeric_limits<double>::infinity());
+ auto theMedium =
+ 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));
const CoordinateSystem& rootCS = env.GetCoordinateSystem();
// setup processes, decays and interactions
- tracking_line::TrackingLine<setup::Stack, setup::Trajectory> tracking(env);
- stack_inspector::StackInspector<setup::Stack> p0(true);
+ tracking_line::TrackingLine tracking;
const std::vector<particles::Code> trackedHadrons = {
particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus,
@@ -258,7 +247,7 @@ int main() {
random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
random::RNGManager::GetInstance().RegisterRandomStream("pythia");
// process::sibyll::Interaction sibyll(env);
- process::pythia::Interaction pythia(env);
+ process::pythia::Interaction pythia;
// process::sibyll::NuclearInteraction sibyllNuc(env, sibyll);
// process::sibyll::Decay decay(trackedHadrons);
process::pythia::Decay decay(trackedHadrons);
@@ -271,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 = p0 << 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/vertical_EAS.cc b/Documentation/Examples/vertical_EAS.cc
index 6c4036debdcfe1d46f5350321f06016ebdf7cd16..efeb157b75b066746aa0e7e0c1ae2506739e866a 100644
--- a/Documentation/Examples/vertical_EAS.cc
+++ b/Documentation/Examples/vertical_EAS.cc
@@ -13,7 +13,6 @@
#include <corsika/process/ProcessSequence.h>
#include <corsika/process/energy_loss/EnergyLoss.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>
@@ -245,16 +244,17 @@ int main() {
random::RNGManager::GetInstance().RegisterRandomStream("cascade");
// setup environment, geometry
- environment::Environment env;
+ using EnvType = Environment<setup::IEnvironmentModel>;
+ EnvType env;
auto& universe = *(env.GetUniverse());
- auto theMedium = environment::Environment::CreateNode<Sphere>(
- Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m},
- 1_km * std::numeric_limits<double>::infinity());
+ 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 = environment::HomogeneousMedium<environment::IMediumModel>;
+ using MyHomogeneousModel = HomogeneousMedium<environment::IMediumModel>;
theMedium->SetModelProperties<MyHomogeneousModel>(
1_kg / (1_m * 1_m * 1_m),
environment::NuclearComposition(
@@ -267,16 +267,15 @@ int main() {
const CoordinateSystem& rootCS = env.GetCoordinateSystem();
// setup processes, decays and interactions
- tracking_line::TrackingLine<setup::Stack, setup::Trajectory> tracking(env);
- // stack_inspector::StackInspector<setup::Stack> stackInspect(true);
+ 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(env);
- process::sibyll::NuclearInteraction sibyllNuc(env, sibyll);
+ 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);
@@ -291,11 +290,8 @@ int main() {
process::EnergyLoss::EnergyLoss eLoss;
// assemble all processes into an ordered process list
- // auto sequence = stackInspect << sibyll << decay << hadronicElastic << cut <<
- // trackWriter; auto sequence = stackInspect << sibyll << sibyllNuc << decay << eLoss <<
// cut << trackWriter;
auto sequence = sibyll << sibyllNuc << decay << eLoss << cut << obsLevel;
- // auto sequence = stackInspect << sibyll << sibyllNuc << decay << eLoss << cut;
// cout << "decltype(sequence)=" << type_id_with_cvr<decltype(sequence)>().pretty_name()
// << "\n";
diff --git a/Environment/FlatAtmosphere/FlatAtmosphere.h b/Environment/FlatAtmosphere/FlatAtmosphere.h
index 6839def48c2f6fddfafd0beb5fd4ca46058ebd09..4227e0de4f5490ee2ac5a90857d68c5653e865ec 100644
--- a/Environment/FlatAtmosphere/FlatAtmosphere.h
+++ b/Environment/FlatAtmosphere/FlatAtmosphere.h
@@ -8,3 +8,5 @@
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
+
+
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 7c32076369e63e2982e8fe744a2a23bf2395a268..7a8b2e785f134ba59546e5622366a03c8f7554a4 100644
--- a/Environment/VolumeTreeNode.h
+++ b/Environment/VolumeTreeNode.h
@@ -108,7 +108,7 @@ namespace corsika::environment {
auto const& GetModelProperties() const { return *fModelProperties; }
- auto const& HasModelProperties() const { return fModelProperties.get() != 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 2903d814ce8ad3eec8ba995519ac4745faeb644d..cce178637582bb48a8b17d1017bc489802c192aa 100644
--- a/Framework/Cascade/CMakeLists.txt
+++ b/Framework/Cascade/CMakeLists.txt
@@ -43,7 +43,7 @@ target_link_libraries (
CORSIKArandom
ProcessSibyll
CORSIKAcascade
-# ProcessStackInspector
+ ProcessStackInspector
ProcessTrackingLine
ProcessNullModel
CORSIKAstackinterface
diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h
index d3ee524313d5af50c4bbb50ce34e00cd7a6218b7..75f29e4ed51bc70d7d31d4917ef4c186baee4bef 100644
--- a/Framework/Cascade/Cascade.h
+++ b/Framework/Cascade/Cascade.h
@@ -36,6 +36,9 @@
#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.
*/
@@ -54,6 +57,7 @@ namespace corsika::cascade {
* <code>auto GetTrack(Particle const& p)</auto>,
* with the return type <code>geometry::Trajectory<corsika::geometry::Line>
* </code>
+ *
* <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>
@@ -61,7 +65,12 @@ namespace corsika::cascade {
*
*/
- template <typename TTracking, typename TProcessList, typename TStack>
+ 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 TStack::ParticleType;
using VolumeTreeNode =
@@ -115,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
@@ -139,6 +149,7 @@ namespace corsika::cascade {
// determine geometric tracking
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 =
@@ -151,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
@@ -217,25 +228,17 @@ namespace corsika::cascade {
important to use projectle (and not vParticle) for Interaction,
and Decay!
*/
- setup::StackView secondaries(vParticle);
- [[maybe_unused]] auto projectile = secondaries.GetProjectile();
- /*
- 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!
- */
- setup::StackView secondaries(vParticle);
+ TStackView secondaries(vParticle);
[[maybe_unused]] auto projectile = secondaries.GetProjectile();
- if (min_distance < distance_max) { // interaction to happen within geometric limit
- // check whether decay or interaction limits this step
-
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;
@@ -255,26 +258,31 @@ namespace corsika::cascade {
random::UniformRealDistribution<InverseTimeType> uniDist(actual_decay_time);
const auto sample_process = uniDist(fRNG);
InverseTimeType inv_decay_count = 0 / second;
- fProcessSequence.SelectDecay(vParticle, projectile, 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(); // last thing in Step function
+
+ 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);
}
}
diff --git a/Framework/Cascade/testCascade.cc b/Framework/Cascade/testCascade.cc
index 1a46850c65c81570acd42bc50ec88d4f58344351..3f95c7a3a63fa6766dd6ed2e26bc5a38e2640a44 100644
--- a/Framework/Cascade/testCascade.cc
+++ b/Framework/Cascade/testCascade.cc
@@ -9,8 +9,6 @@
* the license.
*/
-#include <limits>
-
#include <corsika/cascade/testCascade.h>
#include <corsika/cascade/Cascade.h>
@@ -29,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>
@@ -44,6 +38,7 @@ using namespace corsika::units::si;
using namespace corsika::geometry;
#include <iostream>
+#include <limits>
using namespace std;
auto MakeDummyEnv() {
@@ -52,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.}));
@@ -65,47 +60,75 @@ auto MakeDummyEnv() {
return env;
}
-class ProcessSplit : public process::ContinuousProcess<ProcessSplit> {
+class ProcessSplit : public process::InteractionProcess<ProcessSplit> {
- int fCount = 0;
int fCalls = 0;
- HEPEnergyType fEcrit;
+ GrammageType fX0;
public:
- ProcessSplit(HEPEnergyType e)
- : fEcrit(e) {}
+ ProcessSplit(GrammageType const X0)
+ : fX0(X0) {}
template <typename Particle, typename Track>
- LengthType MaxStepLength(Particle&, Track&) const {
- return 1_m;
+ corsika::units::si::GrammageType GetInteractionLength(Particle&, Track&) const {
+ return fX0;
}
- template <typename Particle, typename Track>
- EProcessReturn DoContinuous(Particle& p, Track&) {
+ 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:
+ ProcessCut(HEPEnergyType e)
+ : fEcrit(e) {}
+
+ 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]") {
@@ -118,12 +141,16 @@ TEST_CASE("Cascade", "[Cascade]") {
stack_inspector::StackInspector<TestCascadeStack> stackInspect(true);
null_model::NullModel nullModel;
+ const GrammageType X0 = 20_g / square(1_cm);
const HEPEnergyType Ecrit = 85_MeV;
- ProcessSplit p1(Ecrit);
- auto sequence = nullModel /* << stackInspect*/ << 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();
@@ -138,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/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 ac49b2b38045858ab2eb50377145ceda457ebff6..33aad8236f6b2943d083fd95addcea1698700d43 100644
--- a/Framework/ProcessSequence/CMakeLists.txt
+++ b/Framework/ProcessSequence/CMakeLists.txt
@@ -15,6 +15,7 @@ set (
ContinuousProcess.h
SecondariesProcess.h
InteractionProcess.h
+ StackProcess.h
DecayProcess.h
ProcessSequence.h
ProcessReturn.h
diff --git a/Framework/ProcessSequence/ProcessSequence.h b/Framework/ProcessSequence/ProcessSequence.h
index 51380c15bd89896cf2aa7c240e92ec44111fe101..ee154d8ae2810c35aa7b12b4c8ba25162504e3de 100644
--- a/Framework/ProcessSequence/ProcessSequence.h
+++ b/Framework/ProcessSequence/ProcessSequence.h
@@ -19,6 +19,7 @@
#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>
@@ -118,16 +119,16 @@ namespace corsika::process {
return ret;
}
- template <typename TSecondaries>
- EProcessReturn DoSecondaries(TSecondaries& vS) {
+ template <typename TStack>
+ EProcessReturn DoStack(TStack& vS) {
EProcessReturn ret = EProcessReturn::eOk;
- if constexpr (std::is_base_of<SecondariesProcess<T1type>, T1type>::value ||
+ if constexpr (std::is_base_of<StackProcess<T1type>, T1type>::value ||
is_process_sequence<T1>::value) {
- ret |= A.DoSecondaries(vS);
+ ret |= A.DoStack(vS);
}
- if constexpr (std::is_base_of<SecondariesProcess<T2type>, T2type>::value ||
+ if constexpr (std::is_base_of<StackProcess<T2type>, T2type>::value ||
is_process_sequence<T2>::value) {
- ret |= B.DoSecondaries(vS);
+ ret |= B.DoStack(vS);
}
return ret;
}
diff --git a/Framework/ProcessSequence/StackProcess.h b/Framework/ProcessSequence/StackProcess.h
new file mode 100644
index 0000000000000000000000000000000000000000..0b8775d716154d97bf6e7bea646fcbbae52fdcdd
--- /dev/null
+++ b/Framework/ProcessSequence/StackProcess.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_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>
+ struct StackProcess {
+
+ 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&);
+ };
+
+ // 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/StackInterface/CombinedStack.h b/Framework/StackInterface/CombinedStack.h
index 7633d383bef2eb1a6131d10e3caf282d606e1adf..8b987d5100732a7a339a23d6fadd3530870bb2ef 100644
--- a/Framework/StackInterface/CombinedStack.h
+++ b/Framework/StackInterface/CombinedStack.h
@@ -41,6 +41,10 @@ namespace corsika::stack {
class CombinedParticleInterface
: 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 PI_C =
CombinedParticleInterface<ParticleInterfaceA, ParticleInterfaceB, StackIterator>;
using PI_A = ParticleInterfaceA<StackIterator>;
diff --git a/Framework/StackInterface/ParticleBase.h b/Framework/StackInterface/ParticleBase.h
index e3504fbeff43ca80f584b179db7eb8bd907609c8..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
@@ -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 3e781f8e6c70edbca333c5ffd53d6f630961678d..563de91539ee39fd32730eb1baae6144f0e8d634 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 {
@@ -53,7 +54,7 @@ namespace corsika::stack {
*/
template <typename StackDataType, template <typename> typename ParticleInterface>
- class SecondaryView : public Stack<StackDataType&, ParticleInterface> {
+ class SecondaryView : public corsika::stack::Stack<StackDataType&, ParticleInterface> {
using ViewType = SecondaryView<StackDataType, ParticleInterface>;
@@ -100,9 +101,13 @@ namespace corsika::stack {
* new stack.
*/
template <typename... Args>
- SecondaryView(Args... args);
+ SecondaryView(Args... args) = delete;
public:
+ /**
+ 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()) {}
@@ -207,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..a666e0ff02817164564576245b990348f1dfb48c 100644
--- a/Framework/StackInterface/Stack.h
+++ b/Framework/StackInterface/Stack.h
@@ -13,12 +13,11 @@
#define _include_Stack_h__
#include <corsika/stack/StackIteratorInterface.h>
+#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 +84,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 4e5cc6d059967d8d64ffdb51a2474052e7bd0d2b..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>
@@ -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) {}
diff --git a/Framework/StackInterface/testCombinedStack.cc b/Framework/StackInterface/testCombinedStack.cc
index 9ae2f8ad60c8cee4028466e977edb37b6e7f508f..fbe4a031c928da4ea5a6c326ad1a98c903933972 100644
--- a/Framework/StackInterface/testCombinedStack.cc
+++ b/Framework/StackInterface/testCombinedStack.cc
@@ -335,20 +335,18 @@ using StackTest2 = CombinedStack<typename StackTest::StackImpl, TestStackData3,
#if defined(__clang__)
using StackTestView = SecondaryView<TestStackData, TestParticleInterface>;
#elif 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>;
-};
-using StackTestView = MakeView<StackTest2>::type;
+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();
diff --git a/Framework/StackInterface/testSecondaryView.cc b/Framework/StackInterface/testSecondaryView.cc
index 6f13ceecd1d5f88171a1b6e0574b5e1aaaab2c93..2d2064b64e2e09c96ec5f77ebbf25d50f5d60e9e 100644
--- a/Framework/StackInterface/testSecondaryView.cc
+++ b/Framework/StackInterface/testSecondaryView.cc
@@ -44,13 +44,11 @@ typedef Stack<TestStackData, TestParticleInterface> StackTest;
#if defined(__clang__)
using StackTestView = SecondaryView<TestStackData, TestParticleInterface>;
#elif 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>;
-};
using StackTestView = MakeView<StackTest>::type;
#endif
+using Particle = typename StackTest::ParticleType;
+
TEST_CASE("SecondaryStack", "[stack]") {
// helper function for sum over stack data
diff --git a/Main/Environment.h b/Main/Environment.h
index 6839def48c2f6fddfafd0beb5fd4ca46058ebd09..4227e0de4f5490ee2ac5a90857d68c5653e865ec 100644
--- a/Main/Environment.h
+++ b/Main/Environment.h
@@ -8,3 +8,5 @@
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
+
+
diff --git a/Main/Stack.h b/Main/Stack.h
index 6839def48c2f6fddfafd0beb5fd4ca46058ebd09..4227e0de4f5490ee2ac5a90857d68c5653e865ec 100644
--- a/Main/Stack.h
+++ b/Main/Stack.h
@@ -8,3 +8,5 @@
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
+
+
diff --git a/Processes/CMakeLists.txt b/Processes/CMakeLists.txt
index cddec6bd425a19f361346a7e50b6da843a0167c9..d6162179ca3440e788a40a9a6ab0c18edb7d44f6 100644
--- a/Processes/CMakeLists.txt
+++ b/Processes/CMakeLists.txt
@@ -4,7 +4,7 @@ add_subdirectory (Sibyll)
if (PYTHIA8_FOUND)
add_subdirectory (Pythia)
endif (PYTHIA8_FOUND)
-# add_subdirectory (StackInspector)
+add_subdirectory (StackInspector)
add_subdirectory (TrackWriter)
add_subdirectory (TrackingLine)
add_subdirectory (HadronicElasticModel)
diff --git a/Processes/EnergyLoss/EnergyLoss.cc b/Processes/EnergyLoss/EnergyLoss.cc
index 094eec177907048fabd4b70c10bc9b54f49941a8..dde5cf442c175d167cd14f605a04f4e5fd288736 100644
--- a/Processes/EnergyLoss/EnergyLoss.cc
+++ b/Processes/EnergyLoss/EnergyLoss.cc
@@ -125,8 +125,8 @@ namespace corsika::process::EnergyLoss {
// Barkas correction O(Z3) higher-order Born approximation
// see Appl. Phys. 85 (1999) 1249
- //double A = 1;
- //if (p.GetPID() == particles::Code::Nucleus) A = p.GetNuclearA();
+ // double A = 1;
+ // if (p.GetPID() == particles::Code::Nucleus) A = p.GetNuclearA();
// double const Erel = (p.GetEnergy()-p.GetMass()) / A / 1_keV;
// double const Llow = 0.01 * Erel;
// double const Lhigh = 1.5/pow(Erel, 0.4) + 45000./Zmat * pow(Erel, 1.6);
diff --git a/Processes/NullModel/NullModel.h b/Processes/NullModel/NullModel.h
index db157b55fefdc1ed584eb7f7437bc45d9d2986b0..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:
diff --git a/Processes/NullModel/testNullModel.cc b/Processes/NullModel/testNullModel.cc
index 4ab4c0422d5638703e5dff5d6e81f9eabaa14819..930d5979fec7bbcc728212b74e944f7dd11599cf 100644
--- a/Processes/NullModel/testNullModel.cc
+++ b/Processes/NullModel/testNullModel.cc
@@ -39,16 +39,12 @@ 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);
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 66c445dd037b5987f10cf96b3e406d5dc8a8a04e..ae8104350086f5cfda94b174305d0f0366faa2ce 100644
--- a/Processes/Pythia/Interaction.cc
+++ b/Processes/Pythia/Interaction.cc
@@ -26,16 +26,14 @@ 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 {
typedef corsika::geometry::Vector<corsika::units::si::hepmomentum_d> MomentumVector;
- Interaction::Interaction(environment::Environment const& env)
- : fEnvironment(env) {}
-
Interaction::~Interaction() { cout << "Pythia::Interaction n=" << fCount << endl; }
void Interaction::Init() {
@@ -78,7 +76,7 @@ namespace corsika::process::pythia {
}
void Interaction::SetParticleListStable(
- const std::vector<particles::Code> particleList) {
+ std::vector<particles::Code> const& particleList) {
for (auto p : particleList) Interaction::SetStable(p);
}
@@ -213,8 +211,7 @@ namespace corsika::process::pythia {
ideally as full particle object so that the four momenta
and the boosts can be defined..
*/
- const auto currentNode =
- fEnvironment.GetUniverse()->GetContainingNode(p.GetPosition());
+ const auto* currentNode = p.GetNode();
const auto mediumComposition =
currentNode->GetModelProperties().GetNuclearComposition();
// determine average interaction length
@@ -222,7 +219,7 @@ namespace corsika::process::pythia {
int i = -1;
si::CrossSectionType weightedProdCrossSection = 0_mb;
// get weights of components from environment/medium
- const auto w = mediumComposition.GetFractions();
+ const auto& w = mediumComposition.GetFractions();
// loop over components in medium
for (auto const targetId : mediumComposition.GetComponents()) {
i++;
@@ -230,8 +227,7 @@ namespace corsika::process::pythia {
auto const [productionCrossSection, elaCrossSection] =
GetCrossSection(corsikaBeamId, targetId, ECoM);
- [[maybe_unused]] auto elaCrossSectionCopy =
- elaCrossSection; // ONLY TO AVOID COMPILER WARNING
+ [[maybe_unused]] const auto& dummy_elaCrossSection = elaCrossSection;
cout << "Interaction: IntLength: pythia return (mb): "
<< productionCrossSection / 1_mb << endl
@@ -263,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;
@@ -286,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
@@ -296,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
@@ -332,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 = fEnvironment.GetUniverse()->GetContainingNode(pOrig);
+ const auto* currentNode = vP.GetNode();
const auto& mediumComposition =
currentNode->GetModelProperties().GetNuclearComposition();
// get cross sections for target materials
@@ -352,9 +348,8 @@ namespace corsika::process::pythia {
for (size_t i = 0; i < compVec.size(); ++i) {
auto const targetId = compVec[i];
const auto [sigProd, sigEla] = GetCrossSection(corsikaBeamId, targetId, Ecm);
+ [[maybe_unused]] const auto& dummy_sigEla = sigEla;
cross_section_of_components[i] = sigProd;
- [[maybe_unused]] auto sigElaCopy =
- sigEla; // to avoid not used warning in array binding
}
const auto corsikaTargetId =
@@ -392,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});
@@ -417,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 7476d09002c2015728cf164e0432f8a73e316026..ac27825acf0561d0d9e7e8e8b66961bdfc674e2f 100644
--- a/Processes/Pythia/Interaction.h
+++ b/Processes/Pythia/Interaction.h
@@ -19,10 +19,6 @@
#include <corsika/units/PhysicalUnits.h>
#include <tuple>
-namespace corsika::environment {
- class Environment;
-}
-
namespace corsika::process::pythia {
class Interaction : public corsika::process::InteractionProcess<Interaction> {
@@ -31,12 +27,12 @@ namespace corsika::process::pythia {
bool fInitialized = false;
public:
- Interaction(corsika::environment::Environment const& env);
+ Interaction() {}
~Interaction();
void Init();
- void SetParticleListStable(const std::vector<particles::Code>);
+ void SetParticleListStable(std::vector<particles::Code> const&);
void SetUnstable(const corsika::particles::Code);
void SetStable(const corsika::particles::Code);
@@ -55,19 +51,18 @@ 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::environment::Environment const& fEnvironment;
corsika::random::RNG& fRNG =
corsika::random::RNGManager::GetInstance().GetRandomStream("pythia");
Pythia8::Pythia fPythia;
diff --git a/Processes/Pythia/testPythia.cc b/Processes/Pythia/testPythia.cc
index a86e4de716cabe69e611ffd142c7065ad7bff40f..3fb2df13ecd1c24e674997a3c3aa468164f21d06 100644
--- a/Processes/Pythia/testPythia.cc
+++ b/Processes/Pythia/testPythia.cc
@@ -98,12 +98,15 @@ using namespace corsika::units::si;
TEST_CASE("pythia process") {
// setup environment, geometry
- environment::Environment env;
+ environment::Environment<environment::IMediumModel> env;
auto& universe = *(env.GetUniverse());
- auto theMedium = environment::Environment::CreateNode<geometry::Sphere>(
- geometry::Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m},
- 1_km * std::numeric_limits<double>::infinity());
+ geometry::CoordinateSystem const& cs = env.GetCoordinateSystem();
+
+ auto theMedium =
+ environment::Environment<environment::IMediumModel>::CreateNode<geometry::Sphere>(
+ geometry::Point{cs, 0_m, 0_m, 0_m},
+ 1_km * std::numeric_limits<double>::infinity());
using MyHomogeneousModel = environment::HomogeneousMedium<environment::IMediumModel>;
theMedium->SetModelProperties<MyHomogeneousModel>(
@@ -112,10 +115,9 @@ TEST_CASE("pythia process") {
std::vector<particles::Code>{particles::Code::Hydrogen},
std::vector<float>{1.}));
+ auto const* nodePtr = theMedium.get(); // save the medium for later use before moving it
universe.AddChild(std::move(theMedium));
- const geometry::CoordinateSystem& cs = env.GetCoordinateSystem();
-
geometry::Point const origin(cs, {0_m, 0_m, 0_m});
geometry::Vector<units::si::SpeedType::dimension_type> v(cs, 0_m / second, 0_m / second,
1_m / second);
@@ -143,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);
}
@@ -163,12 +166,13 @@ TEST_CASE("pythia process") {
std::tuple<particles::Code, units::si::HEPEnergyType,
corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
particles::Code::PiPlus, E0, plab, pos, 0_ns});
+ particle.SetNode(nodePtr);
+ corsika::stack::SecondaryView view(particle);
+ auto projectile = view.GetProjectile();
- process::pythia::Interaction model(env);
-
+ 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/Interaction.cc b/Processes/Sibyll/Interaction.cc
index 83348508b37938b07eff8e664322e780930ba23c..ed550a012d63df6b89586198933ce10c4d01b8d2 100644
--- a/Processes/Sibyll/Interaction.cc
+++ b/Processes/Sibyll/Interaction.cc
@@ -53,7 +53,7 @@ namespace corsika::process::sibyll {
if (fInternalDecays) {
// define which particles are passed to corsika, i.e. which particles make it into
// history even very shortlived particles like charm or pi0 are of interest here
- const std::vector<particles::Code> HadronsWeWantTrackedByCorsika = {
+ const std::vector<particles::Code> hadronsWeWantTrackedByCorsika = {
particles::Code::PiPlus, particles::Code::PiMinus,
particles::Code::Pi0, particles::Code::KMinus,
particles::Code::KPlus, particles::Code::K0Long,
@@ -64,7 +64,7 @@ namespace corsika::process::sibyll {
particles::Code::DMinus, particles::Code::D0,
particles::Code::D0Bar};
- Interaction::SetParticleListStable(HadronsWeWantTrackedByCorsika);
+ Interaction::SetParticleListStable(hadronsWeWantTrackedByCorsika);
}
fInitialized = true;
@@ -72,7 +72,7 @@ namespace corsika::process::sibyll {
}
void Interaction::SetParticleListStable(
- const std::vector<particles::Code> particleList) {
+ std::vector<particles::Code> const& particleList) {
for (auto p : particleList) Interaction::SetStable(p);
}
@@ -91,7 +91,7 @@ namespace corsika::process::sibyll {
tuple<units::si::CrossSectionType, units::si::CrossSectionType>
Interaction::GetCrossSection(const particles::Code BeamId,
const particles::Code TargetId,
- const units::si::HEPEnergyType CoMenergy) {
+ const units::si::HEPEnergyType CoMenergy) const {
using namespace units::si;
double sigProd, sigEla, dummy, dum1, dum3, dum4;
double dumdif[3];
@@ -118,7 +118,7 @@ namespace corsika::process::sibyll {
}
template <>
- units::si::GrammageType Interaction::GetInteractionLength(Particle& p, Track&) {
+ units::si::GrammageType Interaction::GetInteractionLength(Particle& vP, Track&) const {
using namespace units;
using namespace units::si;
@@ -128,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
@@ -138,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
@@ -148,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
@@ -162,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
@@ -176,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
@@ -208,14 +209,14 @@ namespace corsika::process::sibyll {
*/
template <>
- process::EProcessReturn Interaction::DoInteraction(Projectile& p) {
+ 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;
@@ -230,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
@@ -241,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
@@ -277,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
@@ -296,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;
}
@@ -362,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()),
@@ -379,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 4d70a0afb8bc53cb07c1ad09909bce103419dc47..4e1899f685e9254c3207850c227328340a051e94 100644
--- a/Processes/Sibyll/Interaction.h
+++ b/Processes/Sibyll/Interaction.h
@@ -32,18 +32,18 @@ namespace corsika::process::sibyll {
void Init();
- void SetParticleListStable(const std::vector<particles::Code>);
+ void SetParticleListStable(std::vector<particles::Code> const&);
void SetUnstable(const corsika::particles::Code);
void SetStable(const corsika::particles::Code);
bool WasInitialized() { return fInitialized; }
- bool IsValidCoMEnergy(corsika::units::si::HEPEnergyType ecm) {
+ bool IsValidCoMEnergy(corsika::units::si::HEPEnergyType ecm) const {
return (fMinEnergyCoM <= ecm) && (ecm <= fMaxEnergyCoM);
}
- int GetMaxTargetMassNumber() { return fMaxTargetMassNumber; }
- corsika::units::si::HEPEnergyType GetMinEnergyCoM() { return fMinEnergyCoM; }
- corsika::units::si::HEPEnergyType GetMaxEnergyCoM() { return fMaxEnergyCoM; }
- bool IsValidTarget(corsika::particles::Code TargetId) {
+ int GetMaxTargetMassNumber() const { return fMaxTargetMassNumber; }
+ corsika::units::si::HEPEnergyType GetMinEnergyCoM() const { return fMinEnergyCoM; }
+ corsika::units::si::HEPEnergyType GetMaxEnergyCoM() const { return fMaxEnergyCoM; }
+ bool IsValidTarget(corsika::particles::Code TargetId) const {
return (corsika::particles::GetNucleusA(TargetId) < fMaxTargetMassNumber) &&
corsika::particles::IsNucleus(TargetId);
}
@@ -51,18 +51,18 @@ namespace corsika::process::sibyll {
std::tuple<corsika::units::si::CrossSectionType, corsika::units::si::CrossSectionType>
GetCrossSection(const corsika::particles::Code BeamId,
const corsika::particles::Code TargetId,
- const corsika::units::si::HEPEnergyType CoMenergy);
+ const corsika::units::si::HEPEnergyType CoMenergy) const;
template <typename Particle, typename Track>
- corsika::units::si::GrammageType GetInteractionLength(Particle&, Track&);
+ corsika::units::si::GrammageType GetInteractionLength(Particle&, Track&) const;
/**
In this function SIBYLL is called to produce one event. The
event is copied (and boosted) into the shower lab frame.
*/
- template <typename Particle>
- corsika::process::EProcessReturn DoInteraction(Particle&);
+ 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 da85b0f3da83729de49e21d921fce068b4990cea..31abd9da350e1c80c24868eb2313e860b51bd821 100644
--- a/Processes/Sibyll/NuclearInteraction.cc
+++ b/Processes/Sibyll/NuclearInteraction.cc
@@ -37,37 +37,42 @@ using Track = Trajectory;
namespace corsika::process::sibyll {
- NuclearInteraction::NuclearInteraction(process::sibyll::Interaction& hadint)
- : fHadronicInteraction(hadint) {}
+ template <>
+ NuclearInteraction<SetupEnvironment>::NuclearInteraction(
+ process::sibyll::Interaction& hadint, SetupEnvironment const& env)
+ : fEnvironment(env)
+ , fHadronicInteraction(hadint) {}
- NuclearInteraction::~NuclearInteraction() {
+ template <>
+ NuclearInteraction<SetupEnvironment>::~NuclearInteraction() {
cout << "Nuclib::NuclearInteraction n=" << fCount << " Nnuc=" << fNucCount << endl;
- fTargetComponentsIndex.clear();
}
- void NuclearInteraction::Init() {
-
- using random::RNGManager;
-
- // initialize hadronic interaction module
- // TODO: safe to run multiple initializations?
- if (!fHadronicInteraction.WasInitialized()) fHadronicInteraction.Init();
-
- // check compatibility of energy ranges, someone could try to use low-energy model..
- if (!fHadronicInteraction.IsValidCoMEnergy(GetMinEnergyPerNucleonCoM()) ||
- !fHadronicInteraction.IsValidCoMEnergy(GetMaxEnergyPerNucleonCoM()))
- throw std::runtime_error(
- "NuclearInteraction: hadronic interaction model incompatible!");
-
- // initialize nuclib
- // TODO: make sure this does not overlap with sibyll
- nuc_nuc_ini_();
+ template <>
+ void NuclearInteraction<SetupEnvironment>::PrintCrossSectionTable(
+ corsika::particles::Code pCode) {
+ using namespace corsika::particles;
+ const int k = fTargetComponentsIndex.at(pCode);
+ Code pNuclei[] = {Code::Helium, Code::Lithium7, Code::Oxygen,
+ Code::Neon, Code::Argon, Code::Iron};
+ cout << "en/A ";
+ for (auto& j : pNuclei) cout << std::setw(9) << j;
+ cout << endl;
- // initialize cross sections
- InitializeNuclearCrossSections();
+ // loop over energy bins
+ for (int i = 0; i < GetNEnergyBins(); ++i) {
+ cout << " " << i << " ";
+ for (auto& n : pNuclei) {
+ auto const j = GetNucleusA(n);
+ cout << " " << std::setprecision(5) << std::setw(8)
+ << cnucsignuc_.sigma[j - 1][k][i];
+ }
+ cout << endl;
+ }
}
- void NuclearInteraction::InitializeNuclearCrossSections() {
+ template <>
+ void NuclearInteraction<SetupEnvironment>::InitializeNuclearCrossSections() {
using namespace corsika::particles;
using namespace units::si;
@@ -76,9 +81,9 @@ namespace corsika::process::sibyll {
auto const allElementsInUniverse = std::invoke([&]() {
std::set<particles::Code> allElementsInUniverse;
auto collectElements = [&](auto& vtn) {
- if (auto const mp = vtn.GetModelPropertiesPtr();
- mp != nullptr) { // do not query Universe it self, it has no ModelProperties
- auto const& comp = mp->GetNuclearComposition().GetComponents();
+ if (vtn.HasModelProperties()) {
+ auto const& comp =
+ vtn.GetModelProperties().GetNuclearComposition().GetComponents();
for (auto const c : comp) allElementsInUniverse.insert(c);
}
};
@@ -113,10 +118,10 @@ namespace corsika::process::sibyll {
const double dsig = siginel / 1_mb;
const double dsigela = sigela / 1_mb;
// loop over projectiles, mass numbers from 2 to fMaxNucleusAProjectile
- for (int j = 1; j < fMaxNucleusAProjectile; ++j) {
+ for (int j = 1; j < gMaxNucleusAProjectile; ++j) {
const int jj = j + 1;
double sig_out, dsig_out, sigqe_out, dsigqe_out;
- sigma_mc_(jj, ib, dsig, dsigela, fNSample, sig_out, dsig_out, sigqe_out,
+ sigma_mc_(jj, ib, dsig, dsigela, gNSample, sig_out, dsig_out, sigqe_out,
dsigqe_out);
// write to table
cnucsignuc_.sigma[j][k][i] = sig_out;
@@ -132,28 +137,28 @@ namespace corsika::process::sibyll {
}
}
- void NuclearInteraction::PrintCrossSectionTable(corsika::particles::Code pCode) {
- using namespace corsika::particles;
- const int k = fTargetComponentsIndex.at(pCode);
- Code pNuclei[] = {Code::Helium, Code::Lithium7, Code::Oxygen,
- Code::Neon, Code::Argon, Code::Iron};
- cout << "en/A ";
- for (auto& j : pNuclei) cout << std::setw(9) << j;
- cout << endl;
+ template <>
+ void NuclearInteraction<SetupEnvironment>::Init() {
+ // initialize hadronic interaction module
+ // TODO: safe to run multiple initializations?
+ if (!fHadronicInteraction.WasInitialized()) fHadronicInteraction.Init();
- // loop over energy bins
- for (int i = 0; i < GetNEnergyBins(); ++i) {
- cout << " " << i << " ";
- for (auto& n : pNuclei) {
- auto const j = GetNucleusA(n);
- cout << " " << std::setprecision(5) << std::setw(8)
- << cnucsignuc_.sigma[j - 1][k][i];
- }
- cout << endl;
- }
+ // check compatibility of energy ranges, someone could try to use low-energy model..
+ if (!fHadronicInteraction.IsValidCoMEnergy(GetMinEnergyPerNucleonCoM()) ||
+ !fHadronicInteraction.IsValidCoMEnergy(GetMaxEnergyPerNucleonCoM()))
+ throw std::runtime_error(
+ "NuclearInteraction: hadronic interaction model incompatible!");
+
+ // initialize nuclib
+ // TODO: make sure this does not overlap with sibyll
+ nuc_nuc_ini_();
+
+ // initialize cross sections
+ InitializeNuclearCrossSections();
}
- units::si::CrossSectionType NuclearInteraction::ReadCrossSectionTable(
+ template <>
+ units::si::CrossSectionType NuclearInteraction<SetupEnvironment>::ReadCrossSectionTable(
const int ia, particles::Code pTarget, units::si::HEPEnergyType elabnuc) {
using namespace corsika::particles;
using namespace units::si;
@@ -171,15 +176,17 @@ namespace corsika::process::sibyll {
// TODO: remove elastic cross section?
template <>
+ template <>
tuple<units::si::CrossSectionType, units::si::CrossSectionType>
- NuclearInteraction::GetCrossSection(Particle& p, const particles::Code TargetId) {
+ 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;
@@ -208,7 +215,9 @@ namespace corsika::process::sibyll {
}
template <>
- units::si::GrammageType NuclearInteraction::GetInteractionLength(Particle& p, Track&) {
+ template <>
+ units::si::GrammageType NuclearInteraction<SetupEnvironment>::GetInteractionLength(
+ Particle& vP, Track&) {
using namespace units;
using namespace units::si;
@@ -218,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);
@@ -235,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
@@ -258,8 +267,8 @@ namespace corsika::process::sibyll {
// energy limits
// TODO: values depend on hadronic interaction model !! this is sibyll specific
- if (ElabNuc >= 8.5_GeV && ECoMNN >= fMinEnergyPerNucleonCoM &&
- ECoMNN < fMaxEnergyPerNucleonCoM) {
+ if (ElabNuc >= 8.5_GeV && ECoMNN >= gMinEnergyPerNucleonCoM &&
+ ECoMNN < gMaxEnergyPerNucleonCoM) {
// get target from environment
/*
@@ -267,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
@@ -275,15 +284,15 @@ namespace corsika::process::sibyll {
int i = -1;
si::CrossSectionType weightedProdCrossSection = 0_mb;
// get weights of components from environment/medium
- const auto w = mediumComposition.GetFractions();
+ const auto& w = mediumComposition.GetFractions();
// loop over components in medium
for (auto const targetId : mediumComposition.GetComponents()) {
i++;
cout << "NuclearInteraction: get interaction length for target: " << targetId
<< endl;
auto const [productionCrossSection, elaCrossSection] =
- GetCrossSection(p, targetId);
- [[maybe_unused]] auto elaCrossSectionCopy = elaCrossSection;
+ GetCrossSection(vP, targetId);
+ [[maybe_unused]] auto& dummy_elaCrossSection = elaCrossSection;
cout << "NuclearInteraction: "
<< "IntLength: nuclib return (mb): " << productionCrossSection / 1_mb
@@ -298,7 +307,7 @@ namespace corsika::process::sibyll {
GrammageType const int_length = mediumComposition.GetAverageMassNumber() *
units::constants::u / weightedProdCrossSection;
cout << "NuclearInteraction: "
- << "interaction length (g/cm2): " << int_length / (0.001_kg) * 1_cm * 1_cm
+ << "interaction length (g/cm2): " << int_length * (1_cm * 1_cm / (0.001_kg))
<< endl;
return int_length;
@@ -308,7 +317,9 @@ namespace corsika::process::sibyll {
}
template <>
- process::EProcessReturn NuclearInteraction::DoInteraction(Projectile& p) {
+ template <>
+ process::EProcessReturn NuclearInteraction<SetupEnvironment>::DoInteraction(
+ Projectile& vP) {
// this routine superimposes different nucleon-nucleon interactions
// in a nucleus-nucleus interaction, based the SIBYLL routine SIBNUC
@@ -318,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)) {
@@ -337,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++;
@@ -347,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
@@ -369,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
@@ -422,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;
@@ -547,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});
}
@@ -574,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(),
@@ -588,7 +600,7 @@ 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(),
diff --git a/Processes/Sibyll/NuclearInteraction.h b/Processes/Sibyll/NuclearInteraction.h
index 520265fbb692a34571f5e8baa23de3d5e9c4372e..6bee9ba7a7964e74112c2bf60e5cbd4b1c23c2fa 100644
--- a/Processes/Sibyll/NuclearInteraction.h
+++ b/Processes/Sibyll/NuclearInteraction.h
@@ -24,15 +24,15 @@ namespace corsika::process::sibyll {
*
*
**/
-
+ template <class TEnvironment>
class NuclearInteraction
- : public corsika::process::InteractionProcess<NuclearInteraction> {
+ : public corsika::process::InteractionProcess<NuclearInteraction<TEnvironment>> {
int fCount = 0;
int fNucCount = 0;
public:
- NuclearInteraction(corsika::process::sibyll::Interaction& hadint);
+ NuclearInteraction(corsika::process::sibyll::Interaction&, TEnvironment const&);
~NuclearInteraction();
void Init();
void InitializeNuclearCrossSections();
@@ -40,38 +40,41 @@ namespace corsika::process::sibyll {
corsika::units::si::CrossSectionType ReadCrossSectionTable(
const int, corsika::particles::Code, corsika::units::si::HEPEnergyType);
corsika::units::si::HEPEnergyType GetMinEnergyPerNucleonCoM() {
- return fMinEnergyPerNucleonCoM;
+ return gMinEnergyPerNucleonCoM;
}
corsika::units::si::HEPEnergyType GetMaxEnergyPerNucleonCoM() {
- return fMaxEnergyPerNucleonCoM;
+ return gMaxEnergyPerNucleonCoM;
}
- int GetMaxNucleusAProjectile() { return fMaxNucleusAProjectile; }
- int GetMaxNFragments() { return fMaxNFragments; }
- int GetNEnergyBins() { return fNEnBins; }
+ int constexpr GetMaxNucleusAProjectile() { return gMaxNucleusAProjectile; }
+ int constexpr GetMaxNFragments() { return gMaxNFragments; }
+ int constexpr GetNEnergyBins() { return gNEnBins; }
+
template <typename Particle>
std::tuple<corsika::units::si::CrossSectionType, corsika::units::si::CrossSectionType>
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 Projectle>
- corsika::process::EProcessReturn DoInteraction(Projectle& p);
+ corsika::process::EProcessReturn DoInteraction(Projectle&);
private:
+ TEnvironment const& fEnvironment;
corsika::process::sibyll::Interaction& fHadronicInteraction;
std::map<corsika::particles::Code, int> fTargetComponentsIndex;
corsika::random::RNG& fRNG =
corsika::random::RNGManager::GetInstance().GetRandomStream("s_rndm");
- const int fNSample = 500; // number of samples in MC estimation of cross section
- const int fMaxNucleusAProjectile = 56;
- const int fNEnBins = 6;
- const int fMaxNFragments = 60;
+ static constexpr int gNSample =
+ 500; // number of samples in MC estimation of cross section
+ static constexpr int gMaxNucleusAProjectile = 56;
+ static constexpr int gNEnBins = 6;
+ static constexpr int gMaxNFragments = 60;
// energy limits defined by table used for cross section in signuc.f
// 10**1 GeV to 10**6 GeV
- const corsika::units::si::HEPEnergyType fMinEnergyPerNucleonCoM =
+ static constexpr corsika::units::si::HEPEnergyType gMinEnergyPerNucleonCoM =
10. * 1e9 * corsika::units::si::electronvolt;
- const corsika::units::si::HEPEnergyType fMaxEnergyPerNucleonCoM =
+ static constexpr corsika::units::si::HEPEnergyType gMaxEnergyPerNucleonCoM =
1.e6 * 1e9 * corsika::units::si::electronvolt;
};
diff --git a/Processes/Sibyll/testSibyll.cc b/Processes/Sibyll/testSibyll.cc
index ce096994b7829547f28e5f88bf53ff623108daab..d5f5d9bc8e39ac775de3fb9d541a18b596941bad 100644
--- a/Processes/Sibyll/testSibyll.cc
+++ b/Processes/Sibyll/testSibyll.cc
@@ -119,7 +119,7 @@ TEST_CASE("SibyllInterface", "[processes]") {
corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
particles::Code::Proton, E0, plab, pos, 0_ns});
particle.SetNode(nodePtr);
- stack::SecondaryView view(particle);
+ corsika::stack::SecondaryView view(particle);
auto projectile = view.GetProjectile();
Interaction model;
@@ -145,11 +145,11 @@ TEST_CASE("SibyllInterface", "[processes]") {
units::si::TimeType, unsigned short, unsigned short>{
particles::Code::Nucleus, E0, plab, pos, 0_ns, 4, 2});
particle.SetNode(nodePtr);
- stack::SecondaryView view(particle);
+ corsika::stack::SecondaryView view(particle);
auto projectile = view.GetProjectile();
Interaction hmodel;
- NuclearInteraction model(hmodel);
+ NuclearInteraction model(hmodel, env);
model.Init();
[[maybe_unused]] const process::EProcessReturn ret = model.DoInteraction(projectile);
@@ -169,7 +169,7 @@ 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});
- stack::SecondaryView view(particle);
+ corsika::stack::SecondaryView view(particle);
auto projectile = view.GetProjectile();
const std::vector<particles::Code> particleList = {
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 eb0c9e11f871dfc22f5b2a27f46ef8f0e6051d3d..f192926ca73fa2032d05204a80a9250317504f9b 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>
@@ -37,14 +36,12 @@ template <typename Stack>
StackInspector<Stack>::~StackInspector() {}
template <typename Stack>
-process::EProcessReturn StackInspector<Stack>::DoContinuous(Particle&,
- setup::Trajectory&) {
-
+process::EProcessReturn StackInspector<Stack>::DoStack(Stack& 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 +50,21 @@ process::EProcessReturn StackInspector<Stack>::DoContinuous(Particle&,
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() {
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 c4262c1754ddb66cdd55ddc3689e507541da7a7c..4674a591bd1a0ecec1adf5cf80030dd8d79b91a8 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>
@@ -21,8 +21,7 @@ namespace corsika::process {
namespace stack_inspector {
template <typename Stack>
- class StackInspector
- : public corsika::process::ContinuousProcess<StackInspector<Stack>> {
+ class StackInspector : public corsika::process::StackProcess<StackInspector<Stack>> {
typedef typename Stack::ParticleType Particle;
@@ -31,9 +30,7 @@ namespace corsika::process {
~StackInspector();
void Init();
- EProcessReturn DoContinuous(Particle&, corsika::setup::Trajectory&);
- corsika::units::si::LengthType MaxStepLength(Particle&,
- corsika::setup::Trajectory&);
+ EProcessReturn DoStack(Stack&);
private:
bool fReport;
diff --git a/Processes/StackInspector/testStackInspector.cc b/Processes/StackInspector/testStackInspector.cc
index 11e9cff3e4eae58513694dd3b64ba6a6ea64bead..d1cc10739c09102d7138a69ba06528fa503ea81b 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(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/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 2dde8fb4fddddf27689c517449bb8ff23f4e94c1..a59ca2aa3815b5a956e17a07faabf61e2834d972 100644
--- a/Setup/SetupStack.h
+++ b/Setup/SetupStack.h
@@ -68,6 +68,8 @@ private:
*
* 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 {
@@ -146,11 +148,7 @@ namespace corsika::setup {
corsika::stack::SecondaryView<typename corsika::setup::Stack::StackImpl,
corsika::setup::detail::StackWithGeometryInterface>;
#elif 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>;
- };
- using StackView = MakeView<corsika::setup::Stack>::type;
+ using StackView = corsika::stack::MakeView<corsika::setup::Stack>::type;
#endif
} // namespace corsika::setup
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/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: