From 93d9056eef3cf04bb2602891614fae0471948ce9 Mon Sep 17 00:00:00 2001
From: ralfulrich <ralf.ulrich@kit.edu>
Date: Sat, 9 Feb 2019 13:37:02 +0100
Subject: [PATCH] added VolumeTreeNode to Stack
---
Documentation/Examples/cascade_example.cc | 46 ++++----
Documentation/Examples/stack_example.cc | 11 +-
Framework/Cascade/Cascade.h | 7 ++
Framework/Cascade/testCascade.cc | 52 +++------
Framework/Particles/ParticleProperties.cc | 20 ++--
Framework/Particles/ParticleProperties.h | 5 +-
Framework/Particles/testParticles.cc | 10 +-
Framework/StackInterface/CombinedStack.h | 4 +-
Framework/StackInterface/testCombinedStack.cc | 12 +-
Framework/StackInterface/testSecondaryView.cc | 1 +
.../StackInterface/testStackInterface.cc | 1 +
Framework/Units/PhysicalConstants.h | 4 +-
Processes/NullModel/testNullModel.cc | 15 +--
Processes/Sibyll/Decay.h | 26 +++--
Processes/Sibyll/Interaction.h | 10 +-
Processes/Sibyll/NuclearInteraction.h | 39 +++++--
Processes/Sibyll/testSibyll.cc | 68 ++++++-----
.../StackInspector/testStackInspector.cc | 12 +-
SCIENTIFIC_AUTHORS | 3 +-
Setup/SetupStack.h | 107 +++++++++++++++--
.../NuclearStackExtension.h | 109 ++++++++++++------
.../testNuclearStackExtension.cc | 105 ++++++++++-------
Stack/SuperStupidStack/SuperStupidStack.h | 49 +++++---
.../SuperStupidStack/testSuperStupidStack.cc | 21 ++--
24 files changed, 478 insertions(+), 259 deletions(-)
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index 67b6f4d91..82879ccb3 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -55,7 +55,7 @@ using namespace corsika::environment;
using namespace std;
using namespace corsika::units::si;
-class ProcessCut : public corsika::process::ContinuousProcess<ProcessCut> {
+class ProcessCut : public process::ContinuousProcess<ProcessCut> {
HEPEnergyType fECut;
@@ -72,7 +72,7 @@ public:
template <typename Particle>
bool isBelowEnergyCut(Particle& p) const {
// nuclei
- if (p.GetPID() == corsika::particles::Code::Nucleus) {
+ if (p.GetPID() == particles::Code::Nucleus) {
auto const ElabNuc = p.GetEnergy() / p.GetNuclearA();
auto const EcmNN = sqrt(2. * ElabNuc * 0.93827_GeV);
if (ElabNuc < fECut || EcmNN < 10_GeV)
@@ -223,25 +223,24 @@ public:
int main() {
feenableexcept(FE_INVALID);
// initialize random number sequence(s)
- corsika::random::RNGManager::GetInstance().RegisterRandomStream("cascade");
+ random::RNGManager::GetInstance().RegisterRandomStream("cascade");
// setup environment, geometry
- corsika::environment::Environment env;
+ environment::Environment env;
auto& universe = *(env.GetUniverse());
- auto theMedium = corsika::environment::Environment::CreateNode<Sphere>(
+ auto theMedium = environment::Environment::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 =
- corsika::environment::HomogeneousMedium<corsika::environment::IMediumModel>;
+ using MyHomogeneousModel = environment::HomogeneousMedium<environment::IMediumModel>;
theMedium->SetModelProperties<MyHomogeneousModel>(
1_kg / (1_m * 1_m * 1_m),
- corsika::environment::NuclearComposition(
- std::vector<corsika::particles::Code>{corsika::particles::Code::Nitrogen,
- corsika::particles::Code::Oxygen},
+ environment::NuclearComposition(
+ std::vector<particles::Code>{particles::Code::Nitrogen,
+ particles::Code::Oxygen},
std::vector<float>{(float)1. - fox, fox}));
universe.AddChild(std::move(theMedium));
@@ -252,17 +251,17 @@ int main() {
tracking_line::TrackingLine<setup::Stack> tracking(env);
stack_inspector::StackInspector<setup::Stack> p0(true);
- corsika::random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
- corsika::process::sibyll::Interaction sibyll(env);
- corsika::process::sibyll::NuclearInteraction sibyllNuc(env, sibyll);
- corsika::process::sibyll::Decay decay;
+ random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
+ process::sibyll::Interaction sibyll(env);
+ process::sibyll::NuclearInteraction sibyllNuc(env, sibyll);
+ process::sibyll::Decay decay;
ProcessCut cut(20_GeV);
- // corsika::random::RNGManager::GetInstance().RegisterRandomStream("HadronicElasticModel");
- // corsika::process::HadronicElasticModel::HadronicElasticInteraction
+ // random::RNGManager::GetInstance().RegisterRandomStream("HadronicElasticModel");
+ // process::HadronicElasticModel::HadronicElasticInteraction
// hadronicElastic(env);
- corsika::process::TrackWriter::TrackWriter trackWriter("tracks.dat");
+ process::TrackWriter::TrackWriter trackWriter("tracks.dat");
// assemble all processes into an ordered process list
// auto sequence = p0 << sibyll << decay << hadronicElastic << cut << trackWriter;
@@ -277,8 +276,8 @@ int main() {
const Code beamCode = Code::Nucleus;
const int nuclA = 56;
const int nuclZ = int(nuclA / 2.15 + 0.7);
- const HEPMassType mass = corsika::particles::Proton::GetMass() * nuclZ +
- (nuclA - nuclZ) * corsika::particles::Neutron::GetMass();
+ const HEPMassType mass = particles::Proton::GetMass() * nuclZ +
+ (nuclA - nuclZ) * particles::Neutron::GetMass();
const HEPEnergyType E0 =
nuclA *
100_GeV; // 1_PeV crashes with bad COMboost in second interaction (crash later)
@@ -296,16 +295,19 @@ int main() {
};
auto const [px, py, pz] =
momentumComponents(theta / 180. * M_PI, phi / 180. * M_PI, P0);
- auto plab = stack::super_stupid::MomentumVector(rootCS, {px, py, pz});
+ auto plab = corsika::stack::MomentumVector(rootCS, {px, py, pz});
cout << "input particle: " << beamCode << endl;
cout << "input angles: theta=" << theta << " phi=" << phi << endl;
cout << "input momentum: " << plab.GetComponents() / 1_GeV << endl;
Point pos(rootCS, 0_m, 0_m, 0_m);
- stack.AddParticle(beamCode, E0, plab, pos, 0_ns, nuclA, nuclZ);
+ stack.AddParticle(std::tuple<particles::Code, units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, geometry::Point,
+ units::si::TimeType, unsigned short, unsigned short>{
+ beamCode, E0, plab, pos, 0_ns, nuclA, nuclZ});
}
// define air shower object, run simulation
- corsika::cascade::Cascade EAS(env, tracking, sequence, stack);
+ cascade::Cascade EAS(env, tracking, sequence, stack);
EAS.Init();
EAS.Run();
diff --git a/Documentation/Examples/stack_example.cc b/Documentation/Examples/stack_example.cc
index 3431db459..85c68edb6 100644
--- a/Documentation/Examples/stack_example.cc
+++ b/Documentation/Examples/stack_example.cc
@@ -29,9 +29,12 @@ void fill(corsika::stack::super_stupid::SuperStupidStack& s) {
const geometry::CoordinateSystem& rootCS =
geometry::RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
for (int i = 0; i < 11; ++i) {
- s.AddParticle(corsika::particles::Code::Electron, 1.5_GeV * i,
- stack::super_stupid::MomentumVector(rootCS, {0_GeV, 0_GeV, 1_GeV}),
- geometry::Point(rootCS, 0_m, 0_m, 0_m), 0_ns);
+ s.AddParticle(
+ std::tuple<particles::Code, units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
+ particles::Code::Electron, 1.5_GeV * i,
+ corsika::stack::MomentumVector(rootCS, {0_GeV, 0_GeV, 1_GeV}),
+ geometry::Point(rootCS, 0_m, 0_m, 0_m), 0_ns});
}
}
@@ -43,7 +46,7 @@ void read(corsika::stack::super_stupid::SuperStupidStack& s) {
for (auto& p : s) {
total_energy += p.GetEnergy();
// particles are electrons with 1.5 GeV energy times i
- assert(p.GetPID() == corsika::particles::Code::Electron);
+ assert(p.GetPID() == particles::Code::Electron);
assert(p.GetEnergy() == 1.5_GeV * (i++));
}
}
diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h
index d989fba6e..03204dde7 100644
--- a/Framework/Cascade/Cascade.h
+++ b/Framework/Cascade/Cascade.h
@@ -167,6 +167,13 @@ namespace corsika::cascade {
step.LimitEndTo(min_distance);
+ // particle.GetNode(); // previous VolumeNode
+ particle.SetNode(
+ currentNode); // NOTE @Max : here we need to distinguish: IF particle step is
+ // limited by tracking (via fTracking.GetTrack()), THEN we need
+ // to check/update VolumeNodes. In all other cases it is
+ // guaranteed that we are still in the same volume
+
// apply all continuous processes on particle + track
corsika::process::EProcessReturn status =
fProcessSequence.DoContinuous(particle, step, fStack);
diff --git a/Framework/Cascade/testCascade.cc b/Framework/Cascade/testCascade.cc
index 7abf985ad..91505ad8e 100644
--- a/Framework/Cascade/testCascade.cc
+++ b/Framework/Cascade/testCascade.cc
@@ -19,8 +19,6 @@
#include <corsika/process/stack_inspector/StackInspector.h>
#include <corsika/process/tracking_line/TrackingLine.h>
-#include <corsika/stack/super_stupid/SuperStupidStack.h>
-
#include <corsika/particles/ParticleProperties.h>
#include <corsika/geometry/Point.h>
@@ -42,34 +40,32 @@ using corsika::setup::Trajectory;
using namespace corsika;
using namespace corsika::process;
using namespace corsika::units;
+using namespace corsika::units::si;
using namespace corsika::geometry;
#include <iostream>
using namespace std;
-using namespace corsika::units::si;
-corsika::environment::Environment MakeDummyEnv() {
- corsika::environment::Environment env; // dummy environment
+environment::Environment MakeDummyEnv() {
+ environment::Environment env; // dummy environment
auto& universe = *(env.GetUniverse());
- auto theMedium = corsika::environment::Environment::CreateNode<Sphere>(
+ auto theMedium = environment::Environment::CreateNode<Sphere>(
Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m},
1_km * std::numeric_limits<double>::infinity());
- using MyHomogeneousModel =
- corsika::environment::HomogeneousMedium<corsika::environment::IMediumModel>;
+ using MyHomogeneousModel = environment::HomogeneousMedium<environment::IMediumModel>;
theMedium->SetModelProperties<MyHomogeneousModel>(
1_g / (1_m * 1_m * 1_m),
- corsika::environment::NuclearComposition(
- std::vector<corsika::particles::Code>{corsika::particles::Code::Proton},
- std::vector<float>{1.}));
+ environment::NuclearComposition(
+ std::vector<particles::Code>{particles::Code::Proton}, std::vector<float>{1.}));
universe.AddChild(std::move(theMedium));
return env;
}
-class ProcessSplit : public corsika::process::ContinuousProcess<ProcessSplit> {
+class ProcessSplit : public process::ContinuousProcess<ProcessSplit> {
int fCount = 0;
int fCalls = 0;
@@ -93,7 +89,10 @@ public:
fCount++;
} else {
p.SetEnergy(E / 2);
- p.AddSecondary(p.GetPID(), E / 2, p.GetMomentum(), p.GetPosition(), p.GetTime());
+ 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()});
}
return EProcessReturn::eOk;
}
@@ -110,7 +109,7 @@ private:
};
TEST_CASE("Cascade", "[Cascade]") {
- corsika::random::RNGManager& rmng = corsika::random::RNGManager::GetInstance();
+ random::RNGManager& rmng = random::RNGManager::GetInstance();
rmng.RegisterRandomStream("cascade");
auto env = MakeDummyEnv();
@@ -123,34 +122,21 @@ TEST_CASE("Cascade", "[Cascade]") {
auto sequence = p0 << p1;
setup::Stack stack;
- corsika::cascade::Cascade EAS(env, tracking, sequence, stack);
+ cascade::Cascade EAS(env, tracking, sequence, stack);
CoordinateSystem const& rootCS =
RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
stack.Clear();
HEPEnergyType E0 = 100_GeV;
stack.AddParticle(
- particles::Code::Electron, E0,
- corsika::stack::super_stupid::MomentumVector(rootCS, {0_GeV, 0_GeV, -1_GeV}),
- Point(rootCS, {0_m, 0_m, 10_km}), 0_ns);
+ std::tuple<particles::Code, units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
+ particles::Code::Electron, E0,
+ corsika::stack::MomentumVector(rootCS, {0_GeV, 0_GeV, -1_GeV}),
+ Point(rootCS, {0_m, 0_m, 10_km}), 0_ns});
EAS.Init();
EAS.Run();
CHECK(p1.GetCount() == 2048);
CHECK(p1.GetCalls() == 4095);
-
- /*
- SECTION("sectionTwo") {
- for (int i = 0; i < 0; ++i) {
- stack.Clear();
- auto particle = stack.NewParticle();
- HEPEnergyType E0 = 100_GeV * pow(10, i);
- particle.SetEnergy(E0);
- EAS.Init();
- EAS.Run();
-
- // cout << "Result: E0=" << E0 / 1_GeV << "GeV, count=" << p1.GetCount() << endl;
- }
- }
- */
}
diff --git a/Framework/Particles/ParticleProperties.cc b/Framework/Particles/ParticleProperties.cc
index d59846fcf..a1a3e8cf9 100644
--- a/Framework/Particles/ParticleProperties.cc
+++ b/Framework/Particles/ParticleProperties.cc
@@ -27,17 +27,17 @@ namespace corsika::particles {
std::ostream& operator<<(std::ostream& stream, corsika::particles::Code const p) {
return stream << corsika::particles::GetName(p);
}
-
+
Code ConvertFromPDG(PDGCode p) {
- static_assert(detail::conversionArray.size() % 2 == 1);
- // this will fail, for the strange case where the maxPDG is negative...
- unsigned int constexpr maxPDG{(detail::conversionArray.size() - 1) >> 1};
- auto k = static_cast<PDGCodeType>(p);
- if ((unsigned int)abs(k) <= maxPDG) {
- return detail::conversionArray[k + maxPDG];
- } else {
- return detail::conversionMap.at(p);
- }
+ static_assert(detail::conversionArray.size() % 2 == 1);
+ // this will fail, for the strange case where the maxPDG is negative...
+ unsigned int constexpr maxPDG{(detail::conversionArray.size() - 1) >> 1};
+ auto k = static_cast<PDGCodeType>(p);
+ if ((unsigned int)abs(k) <= maxPDG) {
+ return detail::conversionArray[k + maxPDG];
+ } else {
+ return detail::conversionMap.at(p);
+ }
}
} // namespace corsika::particles
diff --git a/Framework/Particles/ParticleProperties.h b/Framework/Particles/ParticleProperties.h
index 7b2b6f7e5..dcbbe6398 100644
--- a/Framework/Particles/ParticleProperties.h
+++ b/Framework/Particles/ParticleProperties.h
@@ -90,8 +90,7 @@ namespace corsika::particles {
}
corsika::units::si::TimeType constexpr GetLifetime(Code const p) {
- return detail::lifetime[static_cast<CodeIntType>(p)] *
- corsika::units::si::second;
+ return detail::lifetime[static_cast<CodeIntType>(p)] * corsika::units::si::second;
}
bool constexpr IsNucleus(Code const p) {
@@ -111,7 +110,7 @@ namespace corsika::particles {
**/
std::ostream& operator<<(std::ostream& stream, corsika::particles::Code const p);
-
+
Code ConvertFromPDG(PDGCode);
} // namespace corsika::particles
diff --git a/Framework/Particles/testParticles.cc b/Framework/Particles/testParticles.cc
index cba91b476..8bb73e7ac 100644
--- a/Framework/Particles/testParticles.cc
+++ b/Framework/Particles/testParticles.cc
@@ -34,8 +34,10 @@ TEST_CASE("ParticleProperties", "[Particles]") {
SECTION("Masses") {
REQUIRE(Electron::GetMass() / (511_keV) == Approx(1));
REQUIRE(Electron::GetMass() / GetMass(Code::Electron) == Approx(1));
-
- REQUIRE((Proton::GetMass() + Neutron::GetMass()) / corsika::units::constants::nucleonMass == Approx(2));
+
+ REQUIRE((Proton::GetMass() + Neutron::GetMass()) /
+ corsika::units::constants::nucleonMass ==
+ Approx(2));
}
SECTION("Charges") {
@@ -57,14 +59,14 @@ TEST_CASE("ParticleProperties", "[Particles]") {
REQUIRE(GetPDG(Code::NuMu) == PDGCode::NuMu);
REQUIRE(GetPDG(Code::NuE) == PDGCode::NuE);
REQUIRE(GetPDG(Code::MuMinus) == PDGCode::MuMinus);
-
+
REQUIRE(static_cast<int>(GetPDG(Code::PiPlus)) == 211);
REQUIRE(static_cast<int>(GetPDG(Code::DPlus)) == 411);
REQUIRE(static_cast<int>(GetPDG(Code::NuMu)) == 14);
REQUIRE(static_cast<int>(GetPDG(Code::NuEBar)) == -12);
REQUIRE(static_cast<int>(GetPDG(Code::MuMinus)) == 13);
}
-
+
SECTION("Conversion PDG -> internal") {
REQUIRE(ConvertFromPDG(PDGCode::KStarMinus) == Code::KStarMinus);
REQUIRE(ConvertFromPDG(PDGCode::MuPlus) == Code::MuPlus);
diff --git a/Framework/StackInterface/CombinedStack.h b/Framework/StackInterface/CombinedStack.h
index 95085c66a..44a29a628 100644
--- a/Framework/StackInterface/CombinedStack.h
+++ b/Framework/StackInterface/CombinedStack.h
@@ -51,8 +51,8 @@ namespace corsika::stack {
template <typename... Args1>
void SetParticleData(C& p, const std::tuple<Args1...> vA) {
// static_assert(MT<I>::has_not, "error");
- I::SetParticleData(static_cast<I&>(p), vA);
- T::SetParticleData(static_cast<T&>(p));
+ I::SetParticleData(static_cast<I&>(p), vA); // original stack
+ T::SetParticleData(static_cast<T&>(p)); // addon stack
}
template <typename... Args1, typename... Args2>
void SetParticleData(C& p, const std::tuple<Args1...> vA,
diff --git a/Framework/StackInterface/testCombinedStack.cc b/Framework/StackInterface/testCombinedStack.cc
index 5cd65446e..63e984b1b 100644
--- a/Framework/StackInterface/testCombinedStack.cc
+++ b/Framework/StackInterface/testCombinedStack.cc
@@ -27,6 +27,7 @@ using boost::typeindex::type_id_with_cvr;
// cpp file
#include <catch2/catch.hpp>
+using namespace corsika;
using namespace corsika::stack;
using namespace std;
@@ -64,8 +65,7 @@ private:
// defintion of a stack-readout object, the iteractor dereference
// operator will deliver access to these function
template <typename T>
-class TestParticleInterface2
- : public T { // corsika::stack::ParticleBaseAdd<T> {//public T {
+class TestParticleInterface2 : public T {
public:
using T::GetIndex;
@@ -254,8 +254,8 @@ public:
// combined stack
template <typename StackIter>
using CombinedTestInterfaceType2 =
- corsika::stack::CombinedParticleInterface<CombinedTestInterfaceType,
- TestParticleInterface3, StackIter>;
+ corsika::stack::CombinedParticleInterface<StackTest::PIType, TestParticleInterface3,
+ StackIter>;
using StackTest2 = CombinedStack<typename StackTest::StackImpl, TestStackData3,
CombinedTestInterfaceType2>;
@@ -270,8 +270,8 @@ TEST_CASE("Combined Stack - multi", "[stack]") {
auto p1 = s.AddParticle(std::tuple{9.9});
auto p2 = s.AddParticle(std::tuple{8.8}, std::tuple{0.1});
p2.SetData2(0.1); // not clear why this is needed, need to check
- // SetParticleData workflow for more complicated
- // settings
+ // SetParticleData workflow for more complicated
+ // settings
// auto p3 = s.AddParticle( std::tuple {8.8}, std::tuple{1.}, std::tuple{0.1} );
p1.SetData3(20.2);
p2.SetData3(10.3);
diff --git a/Framework/StackInterface/testSecondaryView.cc b/Framework/StackInterface/testSecondaryView.cc
index ebd2182a7..df8521386 100644
--- a/Framework/StackInterface/testSecondaryView.cc
+++ b/Framework/StackInterface/testSecondaryView.cc
@@ -27,6 +27,7 @@ using boost::typeindex::type_id_with_cvr;
// cpp file
#include <catch2/catch.hpp>
+using namespace corsika;
using namespace corsika::stack;
using namespace std;
diff --git a/Framework/StackInterface/testStackInterface.cc b/Framework/StackInterface/testStackInterface.cc
index 7b80b3019..448aab927 100644
--- a/Framework/StackInterface/testStackInterface.cc
+++ b/Framework/StackInterface/testStackInterface.cc
@@ -28,6 +28,7 @@ using boost::typeindex::type_id_with_cvr;
// cpp file
#include <catch2/catch.hpp>
+using namespace corsika;
using namespace corsika::stack;
using namespace std;
diff --git a/Framework/Units/PhysicalConstants.h b/Framework/Units/PhysicalConstants.h
index 1bf9872a0..79dd0e6de 100644
--- a/Framework/Units/PhysicalConstants.h
+++ b/Framework/Units/PhysicalConstants.h
@@ -61,9 +61,9 @@ namespace corsika::units::constants {
// unified atomic mass unit
constexpr quantity<mass_d> u{Rep(1.6605402e-27L) * kilogram};
-
+
auto constexpr nucleonMass = 0.5 * (0.93827 + 0.93957) * 1e9 * electronvolt;
-
+
// etc.
} // namespace corsika::units::constants
diff --git a/Processes/NullModel/testNullModel.cc b/Processes/NullModel/testNullModel.cc
index 6bf7d99da..9be839d9c 100644
--- a/Processes/NullModel/testNullModel.cc
+++ b/Processes/NullModel/testNullModel.cc
@@ -21,8 +21,6 @@
#include <corsika/units/PhysicalUnits.h>
-#include <corsika/stack/super_stupid/SuperStupidStack.h>
-
#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h>
@@ -35,16 +33,19 @@ TEST_CASE("NullModel", "[processes]") {
auto const& dummyCS =
geometry::RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
geometry::Point const origin(dummyCS, {0_m, 0_m, 0_m});
- geometry::Vector<corsika::units::si::SpeedType::dimension_type> v(
- dummyCS, 0_m / second, 0_m / second, 1_m / second);
+ geometry::Vector<units::si::SpeedType::dimension_type> v(dummyCS, 0_m / second,
+ 0_m / second, 1_m / second);
geometry::Line line(origin, v);
geometry::Trajectory<geometry::Line> track(line, 10_s);
setup::Stack stack;
auto particle = stack.AddParticle(
- particles::Code::Electron, 1.5_GeV,
- stack::super_stupid::MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
- geometry::Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s);
+ 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});
SECTION("interface") {
diff --git a/Processes/Sibyll/Decay.h b/Processes/Sibyll/Decay.h
index eddae6138..24ccffcb8 100644
--- a/Processes/Sibyll/Decay.h
+++ b/Processes/Sibyll/Decay.h
@@ -171,13 +171,17 @@ namespace corsika::process {
ss.Clear();
const corsika::particles::Code pCode = p.GetPID();
// copy particle to sibyll stack
- ss.AddParticle(process::sibyll::ConvertToSibyllRaw(pCode), p.GetEnergy(),
- p.GetMomentum(),
- // setting particle mass with Corsika values, may be inconsistent
- // with sibyll internal values
- // TODO: #warning setting particle mass with Corsika values, may be
- // inconsistent with sibyll internal values
- corsika::particles::GetMass(pCode));
+ ss.AddParticle(
+ // std::tuple<corsika::particles::Code, corsika::units::si::HEPEnergyType,
+ // corsika::stack::MomentumVector, corsika::geometry::Point,
+ // corsika::units::si::TimeType>{
+ corsika::process::sibyll::ConvertToSibyllRaw(pCode), p.GetEnergy(),
+ p.GetMomentum(),
+ // setting particle mass with Corsika values, may be inconsistent
+ // with sibyll internal values
+ // TODO: #warning setting particle mass with Corsika values, may be
+ // inconsistent with sibyll internal values
+ corsika::particles::GetMass(pCode));
// remember position
Point const decayPoint = p.GetPosition();
TimeType const t0 = p.GetTime();
@@ -200,8 +204,12 @@ namespace corsika::process {
// FOR NOW: skip particles that have decayed in Sibyll, move to iterator?
if (psib.HasDecayed()) continue;
// add to corsika stack
- p.AddSecondary(process::sibyll::ConvertFromSibyll(psib.GetPID()),
- psib.GetEnergy(), psib.GetMomentum(), decayPoint, t0);
+ p.AddSecondary(
+ std::tuple<corsika::particles::Code, corsika::units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, corsika::geometry::Point,
+ corsika::units::si::TimeType>{
+ corsika::process::sibyll::ConvertFromSibyll(psib.GetPID()),
+ psib.GetEnergy(), psib.GetMomentum(), decayPoint, t0});
}
// empty sibyll stack
ss.Clear();
diff --git a/Processes/Sibyll/Interaction.h b/Processes/Sibyll/Interaction.h
index 69755f00e..6a85592fe 100644
--- a/Processes/Sibyll/Interaction.h
+++ b/Processes/Sibyll/Interaction.h
@@ -355,9 +355,13 @@ namespace corsika::process::sibyll {
auto const Plab = boost.fromCoM(FourVector(eCoM, pCoM));
// add to corsika stack
- auto pnew = p.AddSecondary(process::sibyll::ConvertFromSibyll(psib.GetPID()),
- Plab.GetTimeLikeComponent(),
- Plab.GetSpaceLikeComponents(), pOrig, tOrig);
+ auto pnew =
+ p.AddSecondary(std::tuple<particles::Code, units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, geometry::Point,
+ units::si::TimeType>{
+ process::sibyll::ConvertFromSibyll(psib.GetPID()),
+ Plab.GetTimeLikeComponent(), Plab.GetSpaceLikeComponents(), pOrig,
+ tOrig});
Plab_final += pnew.GetMomentum();
Elab_final += pnew.GetEnergy();
diff --git a/Processes/Sibyll/NuclearInteraction.h b/Processes/Sibyll/NuclearInteraction.h
index 7b7492d4c..e050b8c7e 100644
--- a/Processes/Sibyll/NuclearInteraction.h
+++ b/Processes/Sibyll/NuclearInteraction.h
@@ -251,10 +251,10 @@ namespace corsika::process::sibyll {
"should use NuclearStackExtension!");
auto const ProjMass =
- p.GetNuclearZ() * corsika::particles::Proton::GetMass() +
- (p.GetNuclearA() - p.GetNuclearZ()) * corsika::particles::Neutron::GetMass();
+ p.GetNuclearZ() * corsika::particles::Proton::GetMass() +
+ (p.GetNuclearA() - p.GetNuclearZ()) * corsika::particles::Neutron::GetMass();
cout << "NuclearInteraction: projectile mass: " << ProjMass / 1_GeV << endl;
-
+
fCount++;
const CoordinateSystem& rootCS =
@@ -469,12 +469,20 @@ namespace corsika::process::sibyll {
if (nuclA == 1)
// add nucleon
- p.AddSecondary(specCode, Plab.GetTimeLikeComponent(),
- Plab.GetSpaceLikeComponents(), pOrig, tOrig);
+ p.AddSecondary(
+ std::tuple<corsika::particles::Code, corsika::units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, corsika::geometry::Point,
+ corsika::units::si::TimeType>{
+ specCode, Plab.GetTimeLikeComponent(), Plab.GetSpaceLikeComponents(),
+ pOrig, tOrig});
else
// add nucleus
- p.AddSecondary(specCode, Plab.GetTimeLikeComponent(),
- Plab.GetSpaceLikeComponents(), pOrig, tOrig, nuclA, nuclZ);
+ p.AddSecondary(
+ std::tuple<corsika::particles::Code, corsika::units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, corsika::geometry::Point,
+ corsika::units::si::TimeType, unsigned short, unsigned short>{
+ specCode, Plab.GetTimeLikeComponent(), Plab.GetSpaceLikeComponents(),
+ pOrig, tOrig, nuclA, nuclZ});
}
// add elastic nucleons to corsika stack
@@ -491,8 +499,12 @@ namespace corsika::process::sibyll {
const double mass_ratio = corsika::particles::GetMass(elaNucCode) / ProjMass;
auto const Plab = PprojLab * mass_ratio;
- p.AddSecondary(elaNucCode, Plab.GetTimeLikeComponent(),
- Plab.GetSpaceLikeComponents(), pOrig, tOrig);
+ p.AddSecondary(
+ std::tuple<corsika::particles::Code, corsika::units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, corsika::geometry::Point,
+ corsika::units::si::TimeType>{
+ elaNucCode, Plab.GetTimeLikeComponent(), Plab.GetSpaceLikeComponents(),
+ pOrig, tOrig});
}
// add inelastic interactions
@@ -502,9 +514,12 @@ namespace corsika::process::sibyll {
auto pCode = corsika::particles::Proton::GetCode();
// temporarily add to stack, will be removed after interaction in DoInteraction
cout << "inelastic interaction no. " << j << endl;
- auto inelasticNucleon =
- p.AddSecondary(pCode, PprojNucLab.GetTimeLikeComponent(),
- PprojNucLab.GetSpaceLikeComponents(), pOrig, tOrig);
+ auto inelasticNucleon = p.AddSecondary(
+ std::tuple<corsika::particles::Code, corsika::units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, corsika::geometry::Point,
+ corsika::units::si::TimeType>{
+ pCode, PprojNucLab.GetTimeLikeComponent(),
+ PprojNucLab.GetSpaceLikeComponents(), pOrig, tOrig});
// create inelastic interaction
cout << "calling HadronicInteraction..." << endl;
fHadronicInteraction.DoInteraction(inelasticNucleon, s);
diff --git a/Processes/Sibyll/testSibyll.cc b/Processes/Sibyll/testSibyll.cc
index 8f53c6e12..93a2840e0 100644
--- a/Processes/Sibyll/testSibyll.cc
+++ b/Processes/Sibyll/testSibyll.cc
@@ -31,39 +31,37 @@ using namespace corsika::process::sibyll;
TEST_CASE("Sibyll", "[processes]") {
SECTION("Sibyll -> Corsika") {
- REQUIRE(corsika::particles::Electron::GetCode() ==
+ REQUIRE(particles::Electron::GetCode() ==
process::sibyll::ConvertFromSibyll(process::sibyll::SibyllCode::Electron));
}
SECTION("Corsika -> Sibyll") {
- REQUIRE(process::sibyll::ConvertToSibyll(corsika::particles::Electron::GetCode()) ==
+ REQUIRE(process::sibyll::ConvertToSibyll(particles::Electron::GetCode()) ==
process::sibyll::SibyllCode::Electron);
- REQUIRE(process::sibyll::ConvertToSibyllRaw(corsika::particles::Proton::GetCode()) ==
- 13);
+ REQUIRE(process::sibyll::ConvertToSibyllRaw(particles::Proton::GetCode()) == 13);
}
SECTION("KnownBySibyll") {
- REQUIRE(process::sibyll::KnownBySibyll(corsika::particles::Electron::GetCode()));
+ REQUIRE(process::sibyll::KnownBySibyll(particles::Electron::GetCode()));
- REQUIRE_FALSE(
- process::sibyll::KnownBySibyll(corsika::particles::XiPrimeC0::GetCode()));
+ REQUIRE_FALSE(process::sibyll::KnownBySibyll(particles::XiPrimeC0::GetCode()));
}
SECTION("canInteractInSibyll") {
- REQUIRE(process::sibyll::CanInteract(corsika::particles::Proton::GetCode()));
- REQUIRE(process::sibyll::CanInteract(corsika::particles::Code::XiCPlus));
+ REQUIRE(process::sibyll::CanInteract(particles::Proton::GetCode()));
+ REQUIRE(process::sibyll::CanInteract(particles::Code::XiCPlus));
- REQUIRE_FALSE(process::sibyll::CanInteract(corsika::particles::Electron::GetCode()));
- REQUIRE_FALSE(process::sibyll::CanInteract(corsika::particles::SigmaC0::GetCode()));
+ REQUIRE_FALSE(process::sibyll::CanInteract(particles::Electron::GetCode()));
+ REQUIRE_FALSE(process::sibyll::CanInteract(particles::SigmaC0::GetCode()));
}
SECTION("cross-section type") {
- REQUIRE(process::sibyll::GetSibyllXSCode(corsika::particles::Code::Electron) == 0);
- REQUIRE(process::sibyll::GetSibyllXSCode(corsika::particles::Code::K0Long) == 3);
- REQUIRE(process::sibyll::GetSibyllXSCode(corsika::particles::Code::SigmaPlus) == 1);
- REQUIRE(process::sibyll::GetSibyllXSCode(corsika::particles::Code::PiMinus) == 2);
+ REQUIRE(process::sibyll::GetSibyllXSCode(particles::Code::Electron) == 0);
+ REQUIRE(process::sibyll::GetSibyllXSCode(particles::Code::K0Long) == 3);
+ REQUIRE(process::sibyll::GetSibyllXSCode(particles::Code::SigmaPlus) == 1);
+ REQUIRE(process::sibyll::GetSibyllXSCode(particles::Code::PiMinus) == 2);
}
}
@@ -87,32 +85,30 @@ using namespace corsika::units;
TEST_CASE("SibyllInterface", "[processes]") {
// setup environment, geometry
- corsika::environment::Environment env;
+ environment::Environment env;
auto& universe = *(env.GetUniverse());
- auto theMedium = corsika::environment::Environment::CreateNode<geometry::Sphere>(
+ auto theMedium = environment::Environment::CreateNode<geometry::Sphere>(
geometry::Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m},
1_km * std::numeric_limits<double>::infinity());
- using MyHomogeneousModel =
- corsika::environment::HomogeneousMedium<corsika::environment::IMediumModel>;
+ using MyHomogeneousModel = environment::HomogeneousMedium<environment::IMediumModel>;
theMedium->SetModelProperties<MyHomogeneousModel>(
1_kg / (1_m * 1_m * 1_m),
- corsika::environment::NuclearComposition(
- std::vector<corsika::particles::Code>{corsika::particles::Code::Oxygen},
- std::vector<float>{1.}));
+ environment::NuclearComposition(
+ std::vector<particles::Code>{particles::Code::Oxygen}, std::vector<float>{1.}));
universe.AddChild(std::move(theMedium));
const geometry::CoordinateSystem& cs = env.GetCoordinateSystem();
geometry::Point const origin(cs, {0_m, 0_m, 0_m});
- geometry::Vector<corsika::units::si::SpeedType::dimension_type> v(
- cs, 0_m / second, 0_m / second, 1_m / second);
+ geometry::Vector<units::si::SpeedType::dimension_type> v(cs, 0_m / second, 0_m / second,
+ 1_m / second);
geometry::Line line(origin, v);
geometry::Trajectory<geometry::Line> track(line, 10_s);
- corsika::random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
+ random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
SECTION("InteractionInterface") {
@@ -120,9 +116,12 @@ TEST_CASE("SibyllInterface", "[processes]") {
const HEPEnergyType E0 = 100_GeV;
HEPMomentumType P0 =
sqrt(E0 * E0 - particles::Proton::GetMass() * particles::Proton::GetMass());
- auto plab = stack::super_stupid::MomentumVector(cs, {0_GeV, 0_GeV, -P0});
+ auto plab = corsika::stack::MomentumVector(cs, {0_GeV, 0_GeV, -P0});
geometry::Point pos(cs, 0_m, 0_m, 0_m);
- auto particle = stack.AddParticle(particles::Code::Proton, E0, plab, pos, 0_ns);
+ auto particle = stack.AddParticle(
+ std::tuple<particles::Code, units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
+ particles::Code::Proton, E0, plab, pos, 0_ns});
Interaction model(env);
@@ -139,10 +138,14 @@ TEST_CASE("SibyllInterface", "[processes]") {
const HEPEnergyType E0 = 400_GeV;
HEPMomentumType P0 =
sqrt(E0 * E0 - particles::Proton::GetMass() * particles::Proton::GetMass());
- auto plab = stack::super_stupid::MomentumVector(cs, {0_GeV, 0_GeV, -P0});
+ auto plab = corsika::stack::MomentumVector(cs, {0_GeV, 0_GeV, -P0});
geometry::Point pos(cs, 0_m, 0_m, 0_m);
- auto particle = stack.AddParticle(particles::Code::Nucleus, E0, plab, pos, 0_ns, 4, 2);
+ auto particle =
+ stack.AddParticle(std::tuple<particles::Code, units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, geometry::Point,
+ units::si::TimeType, unsigned short, unsigned short>{
+ particles::Code::Nucleus, E0, plab, pos, 0_ns, 4, 2});
Interaction hmodel(env);
NuclearInteraction model(env, hmodel);
@@ -160,9 +163,12 @@ TEST_CASE("SibyllInterface", "[processes]") {
const HEPEnergyType E0 = 10_GeV;
HEPMomentumType P0 =
sqrt(E0 * E0 - particles::Proton::GetMass() * particles::Proton::GetMass());
- auto plab = stack::super_stupid::MomentumVector(cs, {0_GeV, 0_GeV, -P0});
+ auto plab = corsika::stack::MomentumVector(cs, {0_GeV, 0_GeV, -P0});
geometry::Point pos(cs, 0_m, 0_m, 0_m);
- auto particle = stack.AddParticle(particles::Code::Proton, E0, plab, pos, 0_ns);
+ auto particle = stack.AddParticle(
+ std::tuple<particles::Code, units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
+ particles::Code::Proton, E0, plab, pos, 0_ns});
Decay model;
diff --git a/Processes/StackInspector/testStackInspector.cc b/Processes/StackInspector/testStackInspector.cc
index 784be88f9..11e9cff3e 100644
--- a/Processes/StackInspector/testStackInspector.cc
+++ b/Processes/StackInspector/testStackInspector.cc
@@ -33,16 +33,18 @@ TEST_CASE("StackInspector", "[processes]") {
auto const& rootCS =
geometry::RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
geometry::Point const origin(rootCS, {0_m, 0_m, 0_m});
- geometry::Vector<corsika::units::si::SpeedType::dimension_type> v(
- rootCS, 0_m / second, 0_m / second, 1_m / second);
+ geometry::Vector<units::si::SpeedType::dimension_type> v(rootCS, 0_m / second,
+ 0_m / second, 1_m / second);
geometry::Line line(origin, v);
geometry::Trajectory<geometry::Line> track(line, 10_s);
setup::Stack stack;
auto particle = stack.AddParticle(
- particles::Code::Electron, 10_GeV,
- corsika::stack::super_stupid::MomentumVector(rootCS, {0_GeV, 0_GeV, -1_GeV}),
- geometry::Point(rootCS, {0_m, 0_m, 10_km}), 0_ns);
+ std::tuple<particles::Code, units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
+ particles::Code::Electron, 10_GeV,
+ corsika::stack::MomentumVector(rootCS, {0_GeV, 0_GeV, -1_GeV}),
+ geometry::Point(rootCS, {0_m, 0_m, 10_km}), 0_ns});
SECTION("interface") {
diff --git a/SCIENTIFIC_AUTHORS b/SCIENTIFIC_AUTHORS
index 2e7676f37..823eb4b4e 100644
--- a/SCIENTIFIC_AUTHORS
+++ b/SCIENTIFIC_AUTHORS
@@ -1,2 +1,3 @@
-(empty...tbd)
+See https://gitlab.ikp.kit.edu/AirShowerPhysics/corsika/wikis/ICRC2019
+
diff --git a/Setup/SetupStack.h b/Setup/SetupStack.h
index fdfde3d8b..c50b6f400 100644
--- a/Setup/SetupStack.h
+++ b/Setup/SetupStack.h
@@ -12,22 +12,109 @@
#ifndef _corsika_setup_setupstack_h_
#define _corsika_setup_setupstack_h_
-#include <corsika/stack/nuclear_extension/NuclearStackExtension.h>
+// the basic particle data stack:
#include <corsika/stack/super_stupid/SuperStupidStack.h>
-// 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 =
- corsika::stack::nuclear_extension::NuclearParticleInterface<
- corsika::stack::super_stupid::SuperStupidStack::PIType, StackIter>;
+// extension with nuclear data for Code::Nucleus
+#include <corsika/stack/nuclear_extension/NuclearStackExtension.h>
+
+// extension with geometry information for tracking
+#include <corsika/environment/Environment.h>
+#include <corsika/stack/CombinedStack.h>
+
+#include <tuple>
+#include <vector>
+
+// definition of stack-data object to store geometry information
+class GeometryData {
+
+public:
+ // these functions are needed for the Stack interface
+ void Init() {}
+ void Clear() { fNode.clear(); }
+ unsigned int GetSize() const { return fNode.size(); }
+ unsigned int GetCapacity() const { return fNode.size(); }
+ void Copy(const int i1, const int i2) { fNode[i2] = fNode[i1]; }
+ void Swap(const int i1, const int i2) { std::swap(fNode[i1], fNode[i2]); }
+
+ // custom data access function
+ void SetNode(const int i, const corsika::environment::BaseNodeType* v) { fNode[i] = v; }
+ const corsika::environment::BaseNodeType* GetNode(const int i) const {
+ return fNode[i];
+ }
+
+ // these functions are also needed by the Stack interface
+ void IncrementSize() { fNode.push_back(nullptr); }
+ void DecrementSize() {
+ if (fNode.size() > 0) { fNode.pop_back(); }
+ }
+
+ // custom private data section
+private:
+ std::vector<const corsika::environment::BaseNodeType*> fNode;
+};
+
+// defintion of a stack-readout object, the iteractor dereference
+// operator will deliver access to these function
+template <typename T>
+class GeometryDataInterface : public T {
+
+public:
+ using T::GetIndex;
+ using T::GetStackData;
+ using T::SetParticleData;
+
+ // default version for particle-creation from input data
+ void SetParticleData(const std::tuple<const corsika::environment::BaseNodeType*> v) {
+ SetNode(std::get<0>(v));
+ }
+ void SetParticleData(
+ GeometryDataInterface<T>& parent,
+ const std::tuple<const corsika::environment::BaseNodeType*>) { // v = {nullptr}) {
+ SetNode(parent.GetNode());
+ }
+ void SetParticleData() { SetNode(nullptr); }
+ void SetParticleData(GeometryDataInterface<T>& parent) { SetNode(parent.GetNode()); }
+ void SetNode(const corsika::environment::BaseNodeType* v) {
+ GetStackData().SetNode(GetIndex(), v);
+ }
+ const corsika::environment::BaseNodeType* GetNode() const {
+ return GetStackData().GetNode(GetIndex());
+ }
+};
namespace corsika::setup {
- using Stack = corsika::stack::nuclear_extension::NuclearStackExtension<
- corsika::stack::super_stupid::SuperStupidStack, ExtendedParticleInterfaceType>;
+ namespace detail {
+
+ //
+ // 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 =
+ corsika::stack::nuclear_extension::NuclearParticleInterface<
+ corsika::stack::super_stupid::SuperStupidStack::PIType, StackIter>;
+ //
+
+ // the particle data stack with extra nuclear information:
+ using ParticleDataStack = corsika::stack::nuclear_extension::NuclearStackExtension<
+ corsika::stack::super_stupid::SuperStupidStack, ExtendedParticleInterfaceType>;
+
+ // combine particle data stack with geometry information for tracking
+ template <typename StackIter>
+ using StackWithGeometryInterface =
+ corsika::stack::CombinedParticleInterface<ParticleDataStack::PIType,
+ GeometryDataInterface, StackIter>;
+
+ using StackWithGeometry =
+ corsika::stack::CombinedStack<typename ParticleDataStack::StackImpl, GeometryData,
+ StackWithGeometryInterface>;
+
+ } // namespace detail
+
+ // this is the REAL stack we use:
+ using Stack = detail::StackWithGeometry;
- // typedef corsika::stack::super_stupid::SuperStupidStack Stack;
} // namespace corsika::setup
#endif
diff --git a/Stack/NuclearStackExtension/NuclearStackExtension.h b/Stack/NuclearStackExtension/NuclearStackExtension.h
index 355e9821d..11923aba8 100644
--- a/Stack/NuclearStackExtension/NuclearStackExtension.h
+++ b/Stack/NuclearStackExtension/NuclearStackExtension.h
@@ -20,6 +20,7 @@
#include <corsika/geometry/Vector.h>
#include <algorithm>
+#include <tuple>
#include <vector>
namespace corsika::stack {
@@ -39,9 +40,9 @@ namespace corsika::stack {
*
*/
- namespace nuclear_extension {
+ typedef corsika::geometry::Vector<corsika::units::si::hepmomentum_d> MomentumVector;
- typedef corsika::geometry::Vector<corsika::units::si::hepmomentum_d> MomentumVector;
+ namespace nuclear_extension {
/**
* @class NuclearParticleInterface
@@ -54,46 +55,88 @@ namespace corsika::stack {
class NuclearParticleInterface
: public InnerParticleInterface<StackIteratorInterface> {
- // public:
- // using ExtendedParticleInterface =
- // NuclearParticleInterface<InnerParticleInterface, StackIteratorInterface>;
-
protected:
using InnerParticleInterface<StackIteratorInterface>::GetStackData;
using InnerParticleInterface<StackIteratorInterface>::GetIndex;
public:
- void SetParticleData(const corsika::particles::Code vDataPID,
- const corsika::units::si::HEPEnergyType vDataE,
- const MomentumVector& vMomentum,
- const corsika::geometry::Point& vPosition,
- const corsika::units::si::TimeType vTime,
- const unsigned short vA = 0, const unsigned short vZ = 0) {
- if (vDataPID == corsika::particles::Code::Nucleus) {
- if (vA == 0 || vZ == 0) {
- std::ostringstream err;
- err << "NuclearStackExtension: no A and Z specified for new Nucleus!";
- throw std::runtime_error(err.str());
- }
- SetNucleusRef(
- GetStackData().GetNucleusNextRef()); // store this nucleus data ref
- SetNuclearA(vA);
- SetNuclearZ(vZ);
- } else {
- SetNucleusRef(-1); // this is not a nucleus
+ void SetParticleData(
+ const std::tuple<corsika::particles::Code, corsika::units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, corsika::geometry::Point,
+ corsika::units::si::TimeType>& v) {
+ if (std::get<0>(v) == corsika::particles::Code::Nucleus) {
+ std::ostringstream err;
+ err << "NuclearStackExtension: no A and Z specified for new Nucleus!";
+ throw std::runtime_error(err.str());
+ }
+ InnerParticleInterface<StackIteratorInterface>::SetParticleData(v);
+ SetNucleusRef(-1); // this is not a nucleus
+ }
+
+ void SetParticleData(
+ const std::tuple<corsika::particles::Code, corsika::units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, corsika::geometry::Point,
+ corsika::units::si::TimeType, unsigned short, unsigned short>&
+ v) {
+ const unsigned short A = std::get<5>(v);
+ const unsigned short Z = std::get<6>(v);
+ if (std::get<0>(v) != corsika::particles::Code::Nucleus || A == 0 || Z == 0) {
+ std::ostringstream err;
+ err << "NuclearStackExtension: no A and Z specified for new Nucleus!";
+ throw std::runtime_error(err.str());
}
+ SetNucleusRef(GetStackData().GetNucleusNextRef()); // store this nucleus data ref
+ SetNuclearA(A);
+ SetNuclearZ(Z);
InnerParticleInterface<StackIteratorInterface>::SetParticleData(
- vDataPID, vDataE, vMomentum, vPosition, vTime);
+ std::tuple<corsika::particles::Code, corsika::units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, corsika::geometry::Point,
+ corsika::units::si::TimeType>{std::get<0>(v), std::get<1>(v),
+ std::get<2>(v), std::get<3>(v),
+ std::get<4>(v)});
}
- void SetParticleData(InnerParticleInterface<StackIteratorInterface>&,
- const corsika::particles::Code vDataPID,
- const corsika::units::si::HEPEnergyType vDataE,
- const MomentumVector& vMomentum,
- const corsika::geometry::Point& vPosition,
- const corsika::units::si::TimeType vTime,
- const unsigned short vA = 0, const unsigned short vZ = 0) {
- SetParticleData(vDataPID, vDataE, vMomentum, vPosition, vTime, vA, vZ);
+ void SetParticleData(
+ InnerParticleInterface<StackIteratorInterface>& p,
+ const std::tuple<corsika::particles::Code, corsika::units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, corsika::geometry::Point,
+ corsika::units::si::TimeType>& v) {
+ if (std::get<0>(v) == corsika::particles::Code::Nucleus) {
+ std::ostringstream err;
+ err << "NuclearStackExtension: no A and Z specified for new Nucleus!";
+ throw std::runtime_error(err.str());
+ }
+ InnerParticleInterface<StackIteratorInterface>::SetParticleData(
+ p, std::tuple<corsika::particles::Code, corsika::units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, corsika::geometry::Point,
+ corsika::units::si::TimeType>{std::get<0>(v), std::get<1>(v),
+ std::get<2>(v), std::get<3>(v),
+ std::get<4>(v)});
+ SetNucleusRef(-1); // this is not a nucleus
+ }
+
+ void SetParticleData(
+ InnerParticleInterface<StackIteratorInterface>& p,
+ const std::tuple<corsika::particles::Code, corsika::units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, corsika::geometry::Point,
+ corsika::units::si::TimeType, unsigned short, unsigned short>&
+ v) {
+ const unsigned short A = std::get<5>(v);
+ const unsigned short Z = std::get<6>(v);
+ if (std::get<0>(v) != corsika::particles::Code::Nucleus || A == 0 || Z == 0) {
+ std::ostringstream err;
+ err << "NuclearStackExtension: no A and Z specified for new Nucleus!";
+ throw std::runtime_error(err.str());
+ }
+ SetNucleusRef(GetStackData().GetNucleusNextRef()); // store this nucleus data ref
+ SetNuclearA(A);
+ SetNuclearZ(Z);
+ InnerParticleInterface<StackIteratorInterface>::SetParticleData(
+ p, std::tuple<corsika::particles::Code, corsika::units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, corsika::geometry::Point,
+ corsika::units::si::TimeType>{std::get<0>(v), std::get<1>(v),
+ std::get<2>(v), std::get<3>(v),
+ std::get<4>(v)});
}
/**
diff --git a/Stack/NuclearStackExtension/testNuclearStackExtension.cc b/Stack/NuclearStackExtension/testNuclearStackExtension.cc
index 059837949..5185c3538 100644
--- a/Stack/NuclearStackExtension/testNuclearStackExtension.cc
+++ b/Stack/NuclearStackExtension/testNuclearStackExtension.cc
@@ -17,6 +17,8 @@
#include <boost/type_index.hpp>
using boost::typeindex::type_id_with_cvr;
+using namespace corsika;
+using namespace corsika::stack::nuclear_extension;
using namespace corsika::geometry;
using namespace corsika::units::si;
@@ -24,17 +26,13 @@ using namespace corsika::units::si;
// cpp file
#include <catch2/catch.hpp>
-using namespace corsika;
-using namespace corsika::stack::nuclear_extension;
-
// 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 =
- corsika::stack::nuclear_extension::NuclearParticleInterface<
- corsika::stack::super_stupid::SuperStupidStack::PIType, StackIter>;
+using ExtendedParticleInterfaceType = stack::nuclear_extension::NuclearParticleInterface<
+ stack::super_stupid::SuperStupidStack::PIType, StackIter>;
-using ExtStack = NuclearStackExtension<corsika::stack::super_stupid::SuperStupidStack,
+using ExtStack = NuclearStackExtension<stack::super_stupid::SuperStupidStack,
ExtendedParticleInterfaceType>;
#include <iostream>
@@ -45,15 +43,16 @@ TEST_CASE("NuclearStackExtension", "[stack]") {
geometry::CoordinateSystem& dummyCS =
geometry::RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
- // cout << "ParticleType=" << type_id_with_cvr<ParticleType>().pretty_name() << endl;
-
SECTION("write non nucleus") {
NuclearStackExtension<corsika::stack::super_stupid::SuperStupidStack,
ExtendedParticleInterfaceType>
s;
- s.AddParticle(particles::Code::Electron, 1.5_GeV,
- MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
- Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s);
+ s.AddParticle(
+ std::tuple<particles::Code, units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
+ particles::Code::Electron, 1.5_GeV,
+ corsika::stack::MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
+ Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s});
REQUIRE(s.GetSize() == 1);
}
@@ -61,24 +60,34 @@ TEST_CASE("NuclearStackExtension", "[stack]") {
NuclearStackExtension<corsika::stack::super_stupid::SuperStupidStack,
ExtendedParticleInterfaceType>
s;
- s.AddParticle(particles::Code::Nucleus, 1.5_GeV,
- MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
- Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s, 10, 10);
+ s.AddParticle(std::tuple<particles::Code, units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, geometry::Point,
+ units::si::TimeType, unsigned short, unsigned short>{
+ particles::Code::Nucleus, 1.5_GeV,
+ corsika::stack::MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
+ Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s, 10, 10});
REQUIRE(s.GetSize() == 1);
}
SECTION("write invalid nucleus") {
ExtStack s;
- REQUIRE_THROWS(s.AddParticle(
- particles::Code::Nucleus, 1.5_GeV, MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
- Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s, 0, 0));
+ REQUIRE_THROWS(
+ s.AddParticle(std::tuple<particles::Code, units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, geometry::Point,
+ units::si::TimeType, unsigned short, unsigned short>{
+ particles::Code::Nucleus, 1.5_GeV,
+ corsika::stack::MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
+ Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s, 0, 0}));
}
SECTION("read non nucleus") {
ExtStack s;
- s.AddParticle(particles::Code::Electron, 1.5_GeV,
- MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
- Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s);
+ s.AddParticle(
+ std::tuple<particles::Code, units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
+ particles::Code::Electron, 1.5_GeV,
+ corsika::stack::MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
+ Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s});
const auto pout = s.GetNextParticle();
REQUIRE(pout.GetPID() == particles::Code::Electron);
REQUIRE(pout.GetEnergy() == 1.5_GeV);
@@ -87,9 +96,12 @@ TEST_CASE("NuclearStackExtension", "[stack]") {
SECTION("read nucleus") {
ExtStack s;
- s.AddParticle(particles::Code::Nucleus, 1.5_GeV,
- MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
- Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s, 10, 9);
+ s.AddParticle(std::tuple<particles::Code, units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, geometry::Point,
+ units::si::TimeType, unsigned short, unsigned short>{
+ particles::Code::Nucleus, 1.5_GeV,
+ corsika::stack::MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
+ Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s, 10, 9});
const auto pout = s.GetNextParticle();
REQUIRE(pout.GetPID() == particles::Code::Nucleus);
REQUIRE(pout.GetEnergy() == 1.5_GeV);
@@ -100,9 +112,12 @@ TEST_CASE("NuclearStackExtension", "[stack]") {
SECTION("read invalid nucleus") {
ExtStack s;
- s.AddParticle(particles::Code::Electron, 1.5_GeV,
- MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
- Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s);
+ s.AddParticle(
+ std::tuple<particles::Code, units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
+ particles::Code::Electron, 1.5_GeV,
+ corsika::stack::MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
+ Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s});
const auto pout = s.GetNextParticle();
REQUIRE_THROWS(pout.GetNuclearA());
REQUIRE_THROWS(pout.GetNuclearZ());
@@ -114,13 +129,19 @@ TEST_CASE("NuclearStackExtension", "[stack]") {
// add 99 particles, each 10th particle is a nucleus with A=i and Z=A/2!
for (int i = 0; i < 99; ++i) {
if ((i + 1) % 10 == 0) {
- s.AddParticle(particles::Code::Nucleus, 1.5_GeV,
- MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
- Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s, i, i / 2);
+ s.AddParticle(std::tuple<particles::Code, units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, geometry::Point,
+ units::si::TimeType, unsigned short, unsigned short>{
+ particles::Code::Nucleus, 1.5_GeV,
+ corsika::stack::MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
+ Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s, i, i / 2});
} else {
- s.AddParticle(particles::Code::Electron, 1.5_GeV,
- MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
- Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s);
+ s.AddParticle(std::tuple<particles::Code, units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, geometry::Point,
+ units::si::TimeType>{
+ particles::Code::Electron, 1.5_GeV,
+ corsika::stack::MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
+ Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s});
}
}
@@ -135,13 +156,19 @@ TEST_CASE("NuclearStackExtension", "[stack]") {
// add 99 particles, each 10th particle is a nucleus with A=i and Z=A/2!
for (int i = 0; i < 99; ++i) {
if ((i + 1) % 10 == 0) {
- s.AddParticle(particles::Code::Nucleus, i * 15_GeV,
- MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
- Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s, i, i / 2);
+ s.AddParticle(std::tuple<particles::Code, units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, geometry::Point,
+ units::si::TimeType, unsigned short, unsigned short>{
+ particles::Code::Nucleus, i * 15_GeV,
+ corsika::stack::MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
+ Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s, i, i / 2});
} else {
- s.AddParticle(particles::Code::Electron, i * 1.5_GeV,
- MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
- Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s);
+ s.AddParticle(std::tuple<particles::Code, units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, geometry::Point,
+ units::si::TimeType>{
+ particles::Code::Electron, i * 1.5_GeV,
+ corsika::stack::MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
+ Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s});
}
}
diff --git a/Stack/SuperStupidStack/SuperStupidStack.h b/Stack/SuperStupidStack/SuperStupidStack.h
index cb63b0336..365a45aa9 100644
--- a/Stack/SuperStupidStack/SuperStupidStack.h
+++ b/Stack/SuperStupidStack/SuperStupidStack.h
@@ -25,9 +25,9 @@
namespace corsika::stack {
- namespace super_stupid {
+ typedef corsika::geometry::Vector<corsika::units::si::hepmomentum_d> MomentumVector;
- typedef corsika::geometry::Vector<corsika::units::si::hepmomentum_d> MomentumVector;
+ namespace super_stupid {
/**
* Example of a particle object on the stack.
@@ -42,19 +42,36 @@ namespace corsika::stack {
using corsika::stack::ParticleBase<StackIteratorInterface>::GetIndex;
public:
- void SetParticleData(const corsika::particles::Code vDataPID,
- const corsika::units::si::HEPEnergyType vDataE,
- const MomentumVector& vMomentum,
- const corsika::geometry::Point& vPosition,
- const corsika::units::si::TimeType vTime) {
- SetPID(vDataPID);
- SetEnergy(vDataE);
- SetMomentum(vMomentum);
- SetPosition(vPosition);
- SetTime(vTime);
- }
-
- void SetParticleData(ParticleInterface<StackIteratorInterface>& /*parent*/,
+ void SetParticleData(
+ const std::tuple<corsika::particles::Code, corsika::units::si::HEPEnergyType,
+ MomentumVector, corsika::geometry::Point,
+ corsika::units::si::TimeType>& v) {
+ SetPID(std::get<0>(v));
+ SetEnergy(std::get<1>(v));
+ SetMomentum(std::get<2>(v));
+ SetPosition(std::get<3>(v));
+ SetTime(std::get<4>(v));
+ }
+ /*
+ void SetParticleData(const corsika::particles::Code vDataPID,
+ const corsika::units::si::HEPEnergyType vDataE,
+ const MomentumVector& vMomentum,
+ const corsika::geometry::Point& vPosition,
+ const corsika::units::si::TimeType vTime) {
+ }*/
+
+ void SetParticleData(
+ ParticleInterface<StackIteratorInterface>&,
+ const std::tuple<corsika::particles::Code, corsika::units::si::HEPEnergyType,
+ MomentumVector, corsika::geometry::Point,
+ corsika::units::si::TimeType>& v) {
+ SetPID(std::get<0>(v));
+ SetEnergy(std::get<1>(v));
+ SetMomentum(std::get<2>(v));
+ SetPosition(std::get<3>(v));
+ SetTime(std::get<4>(v));
+ }
+ /* void SetParticleData(ParticleInterface<StackIteratorInterface>&,
const corsika::particles::Code vDataPID,
const corsika::units::si::HEPEnergyType vDataE,
const MomentumVector& vMomentum,
@@ -65,7 +82,7 @@ namespace corsika::stack {
SetMomentum(vMomentum);
SetPosition(vPosition);
SetTime(vTime);
- }
+ }*/
/// individual setters
void SetPID(const corsika::particles::Code id) {
diff --git a/Stack/SuperStupidStack/testSuperStupidStack.cc b/Stack/SuperStupidStack/testSuperStupidStack.cc
index 0e3628b46..3e9aa8466 100644
--- a/Stack/SuperStupidStack/testSuperStupidStack.cc
+++ b/Stack/SuperStupidStack/testSuperStupidStack.cc
@@ -34,16 +34,19 @@ TEST_CASE("SuperStupidStack", "[stack]") {
SECTION("read+write") {
SuperStupidStack s;
- s.AddParticle(particles::Code::Electron, 1.5_GeV,
- MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
- Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s);
+ s.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,
+ corsika::stack::MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
+ Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s});
// read
REQUIRE(s.GetSize() == 1);
auto pout = s.GetNextParticle();
REQUIRE(pout.GetPID() == particles::Code::Electron);
REQUIRE(pout.GetEnergy() == 1.5_GeV);
- // REQUIRE(pout.GetMomentum() == stack::super_stupid::MomentumVector(dummyCS, {1_GeV,
+ // REQUIRE(pout.GetMomentum() == stack::MomentumVector(dummyCS, {1_GeV,
// 1_GeV, 1_GeV})); REQUIRE(pout.GetPosition() == Point(dummyCS, {1 * meter, 1 *
// meter, 1 * meter}));
REQUIRE(pout.GetTime() == 100_s);
@@ -53,9 +56,13 @@ TEST_CASE("SuperStupidStack", "[stack]") {
SuperStupidStack s;
for (int i = 0; i < 99; ++i)
- s.AddParticle(particles::Code::Electron, 1.5_GeV,
- MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
- Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s);
+ s.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,
+ corsika::stack::MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
+ Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s});
REQUIRE(s.GetSize() == 99);
--
GitLab