diff --git a/corsika/detail/framework/geometry/Plane.inl b/corsika/detail/framework/geometry/Plane.inl
index f2170b48ea1cebd0076f13bee6fcc513c06ebf3f..d1a41c3eafebcdf6a8f62c15a6734030a39b1163 100644
--- a/corsika/detail/framework/geometry/Plane.inl
+++ b/corsika/detail/framework/geometry/Plane.inl
@@ -22,7 +22,7 @@ namespace corsika {
return fNormal.dot(vP - fCenter) > corsika::units::si::LengthType::zero();
}
- units::si::LengthType Plane::DistanceTo(geometry::Point const& vP) const
+ units::si::LengthType Plane::DistanceTo(corsika::Point const& vP) const
{
return (fNormal * (vP - fCenter).dot(fNormal)).norm();
}
diff --git a/corsika/detail/media/Universe.inl b/corsika/detail/media/Universe.inl
index d1f16c4e5f0773bf75192bc2cc520c99568ea8a5..c48d97b887c2b8a35c3c081612c0a329d928aa28 100644
--- a/corsika/detail/media/Universe.inl
+++ b/corsika/detail/media/Universe.inl
@@ -14,18 +14,12 @@
namespace corsika {
- Universe::Universe(corsika::CoordinateSystem const& pCS):
- corsika::Sphere(
- corsika::Point{pCS,
- 0 * corsika::units::si::meter,
- 0 * corsika::units::si::meter,
- 0 * corsika::units::si::meter},
- corsika::units::si::meter * std::numeric_limits<double>::infinity())
- {}
+ Universe::Universe(corsika::CoordinateSystem const& pCS)
+ : corsika::Sphere(
+ corsika::Point{pCS, 0 * corsika::units::si::meter,
+ 0 * corsika::units::si::meter, 0 * corsika::units::si::meter},
+ corsika::units::si::meter * std::numeric_limits<double>::infinity()) {}
- bool Universe::Contains(corsika::Point const&) const
- {
- return true;
- }
+ bool Universe::Contains(corsika::Point const&) const { return true; }
-}
\ No newline at end of file
+} // namespace corsika
\ No newline at end of file
diff --git a/corsika/detail/modules/null_model/NullModel.inl b/corsika/detail/modules/null_model/NullModel.inl
new file mode 100644
index 0000000000000000000000000000000000000000..87115f4e5748ca959e915d46d67e1423abb0a96a
--- /dev/null
+++ b/corsika/detail/modules/null_model/NullModel.inl
@@ -0,0 +1,36 @@
+/*
+ * (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.
+ */
+
+#pragma once
+
+#include <corsika/modules/null_model/NullModel.hpp>
+#include <corsika/setup/SetupStack.hpp>
+#include <corsika/setup/SetupTrajectory.hpp>
+
+namespace corsika::null_model {
+
+ void NullModel::Init() {}
+
+ NullModel::NullModel(corsika::units::si::LengthType maxStepLength)
+ : fMaxStepLength(maxStepLength) {}
+
+ template <>
+ corsika::EProcessReturn NullModel::DoContinuous(corsika::setup::Stack::ParticleType&,
+ corsika::setup::Trajectory&) const {
+ return corsika::EProcessReturn::eOk;
+ }
+
+ template <>
+ units::si::LengthType NullModel::MaxStepLength(corsika::setup::Stack::ParticleType&,
+ corsika::setup::Trajectory&) const {
+ return fMaxStepLength;
+ }
+
+} // namespace corsika::null_model
diff --git a/corsika/detail/modules/observation_plane/ObservationPlane.inl b/corsika/detail/modules/observation_plane/ObservationPlane.inl
new file mode 100644
index 0000000000000000000000000000000000000000..87a7a0410fa70edd6bf68761d49550a20c2eae1b
--- /dev/null
+++ b/corsika/detail/modules/observation_plane/ObservationPlane.inl
@@ -0,0 +1,68 @@
+/*
+ * (c) Copyright 2019 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#include <corsika/modules/observation_plane/ObservationPlane.hpp>
+
+#include <fstream>
+
+namespace corsika::observation_plane {
+
+ObservationPlane::ObservationPlane(corsika::Plane const& obsPlane,
+ std::string const& filename, bool deleteOnHit)
+ : plane_(obsPlane)
+ , outputStream_(filename)
+ , deleteOnHit_(deleteOnHit) {
+ outputStream_ << "#PDG code, energy / eV, distance to center / m" << std::endl;
+}
+
+corsika::EProcessReturn ObservationPlane::DoContinuous(
+ corsika::setup::Stack::ParticleType const& particle, corsika::setup::Trajectory const& trajectory) {
+ using namespace corsika::units::si;
+
+ TimeType const timeOfIntersection =
+ (plane_.GetCenter() - trajectory.GetR0()).dot(plane_.GetNormal()) /
+ trajectory.GetV0().dot(plane_.GetNormal());
+
+ if (timeOfIntersection < TimeType::zero()) { return corsika::EProcessReturn::eOk; }
+
+ if (plane_.IsAbove(trajectory.GetR0()) == plane_.IsAbove(trajectory.GetPosition(1))) {
+ return corsika::EProcessReturn::eOk;
+ }
+
+ outputStream_ << static_cast<int>(corsika::GetPDG(particle.GetPID())) << ' '
+ << particle.GetEnergy() * (1 / 1_eV) << ' '
+ << (trajectory.GetPosition(1) - plane_.GetCenter()).norm() / 1_m
+ << std::endl;
+
+ if (deleteOnHit_) {
+ return corsika::EProcessReturn::eParticleAbsorbed;
+ } else {
+ return corsika::EProcessReturn::eOk;
+ }
+}
+
+corsika::LengthType ObservationPlane::MaxStepLength(corsika::setup::Stack::ParticleType const&,
+ corsika::setup::Trajectory const& trajectory) {
+
+ using namespace corsika::units::si;
+
+TimeType const timeOfIntersection =
+ (plane_.GetCenter() - trajectory.GetR0()).dot(plane_.GetNormal()) /
+ trajectory.GetV0().dot(plane_.GetNormal());
+
+ if (timeOfIntersection < TimeType::zero()) {
+ return std::numeric_limits<double>::infinity() * 1_m;
+ }
+
+ auto const pointOfIntersection = trajectory.GetPosition(timeOfIntersection);
+ return (trajectory.GetR0() - pointOfIntersection).norm() * 1.0001;
+}
+
+}
\ No newline at end of file
diff --git a/corsika/framework/geometry/Plane.hpp b/corsika/framework/geometry/Plane.hpp
index 6a0d04f4c0c8e776b1ef6a8e324978897be5c0ae..fed4dab83983d44d94f940124893a2f8ef8b5df9 100644
--- a/corsika/framework/geometry/Plane.hpp
+++ b/corsika/framework/geometry/Plane.hpp
@@ -29,7 +29,7 @@ namespace corsika {
bool IsAbove(Point const& vP) const ;
- units::si::LengthType DistanceTo(geometry::Point const& vP) const ;
+ units::si::LengthType DistanceTo(corsika::Point const& vP) const ;
Point const& GetCenter() const ;
diff --git a/corsika/media/Universe.hpp b/corsika/media/Universe.hpp
index 7241b794b05ed91f06595f39eb9cb27e70cdd2dd..fe61ba0feeb13530e5ccaca1b9666c36e8cc5103 100644
--- a/corsika/media/Universe.hpp
+++ b/corsika/media/Universe.hpp
@@ -17,7 +17,7 @@ namespace corsika {
struct Universe : public corsika::Sphere
{
- Universe(corsika::CoordinateSystem const& pCS);
+ inline Universe(corsika::CoordinateSystem const& pCS);
inline bool Contains(corsika::Point const&) const override;
};
diff --git a/corsika/modules/null_model/NullModel.hpp b/corsika/modules/null_model/NullModel.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..dd754b734a099008ae227d3f77d432c28ffbe879
--- /dev/null
+++ b/corsika/modules/null_model/NullModel.hpp
@@ -0,0 +1,36 @@
+/*
+ * (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.
+ */
+
+#pragma once
+
+#include <corsika/framework/sequence/BaseProcess.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+
+namespace corsika::null_model {
+
+ class NullModel : public corsika::BaseProcess<NullModel> {
+ corsika::units::si::LengthType const fMaxStepLength;
+
+ public:
+ NullModel(corsika::units::si::LengthType maxStepLength =
+ corsika::units::si::meter * std::numeric_limits<double>::infinity());
+
+ void Init();
+
+ template <typename Particle, typename Track>
+ corsika::EProcessReturn DoContinuous(Particle&, Track&) const;
+
+ template <typename Particle, typename Track>
+ corsika::units::si::LengthType MaxStepLength(Particle&, Track&) const;
+ };
+
+} // namespace corsika::process::null_model
+
+#include <corsika/detail/modules/null_model/NullModel.inl>
diff --git a/corsika/modules/observation_plane/ObservationPlane.hpp b/corsika/modules/observation_plane/ObservationPlane.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..dc3410647ec1ca8f8605542a5d562c9760becf03
--- /dev/null
+++ b/corsika/modules/observation_plane/ObservationPlane.hpp
@@ -0,0 +1,49 @@
+/*
+ * (c) Copyright 2019 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.
+ */
+
+#pragma once
+
+#include <corsika/framework/geometry/Plane.hpp>
+#include <corsika/framework/sequence/ContinuousProcess.hpp>
+#include <corsika/setup/SetupStack.hpp>
+#include <corsika/setup/SetupTrajectory.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+
+#include <fstream>
+
+namespace corsika::observation_plane {
+
+ /**
+ * The ObservationPlane writes PDG codes, energies, and distances of particles to the
+ * central point of the plane into its output file. The particles are considered
+ * "absorbed" afterwards.
+ */
+ class ObservationPlane : public corsika::ContinuousProcess<ObservationPlane> {
+
+ public:
+ ObservationPlane(corsika::Plane const&, std::string const&, bool = true);
+ void Init() {}
+
+ corsika::EProcessReturn DoContinuous(
+ corsika::setup::Stack::ParticleType const& vParticle,
+ corsika::setup::Trajectory const& vTrajectory);
+
+ corsika::units::si::LengthType MaxStepLength(
+ corsika::setup::Stack::ParticleType const&,
+ corsika::setup::Trajectory const& vTrajectory);
+
+ private:
+ corsika::Plane const plane_;
+ std::ofstream outputStream_;
+ bool const deleteOnHit_;
+ };
+} // namespace corsika::process::observation_plane
+
+#include <corsika/detail/modules/observation_plane/ObservationPlane.inl>
diff --git a/corsika/modules/sibyll/Random.hpp b/corsika/modules/sibyll/Random.hpp
index 3d82be78bc0100e40e41c867a04f3a48cc6c6eb0..8c0ddcf91523a0e392bc4a1c4e636ec6f271e05f 100644
--- a/corsika/modules/sibyll/Random.hpp
+++ b/corsika/modules/sibyll/Random.hpp
@@ -13,10 +13,9 @@ namespace sibyll {
double rndm_interface() {
static corsika::RNG& rng =
- corsika::RNGManager::GetInstance().GetRandomStream("s_rndm");
+ corsika::RNGManager::GetInstance().GetRandomStream("s_rndm");
std::uniform_real_distribution<double> dist;
return dist(rng);
}
}
-
diff --git a/tests/modules/CMakeLists.txt b/tests/modules/CMakeLists.txt
index f4c2fea50aeeea3e3e8d54bcb295f89a2287a19e..e4449b878d8affe9577cbc2723eba46304c74289 100644
--- a/tests/modules/CMakeLists.txt
+++ b/tests/modules/CMakeLists.txt
@@ -2,9 +2,18 @@
set (test_modules_sources
TestMain.cpp
testSibyll.cpp
+ testNullModel.cpp
+ testObservationPlane.cpp
+ #testParticleCut.cpp
+ #testPythia.cpp
+ #testQGSJetII.cpp
+ #testSibyll.cpp
+ #testStackInspector.cpp
+ #testSwitchProcess.cpp
+ #testTrackingLine.cpp
+ #testUrQMD.cpp
)
-
add_executable (testModules ${test_modules_sources})
target_link_libraries (testModules CORSIKA8 Sibyll_static Catch2)
diff --git a/tests/modules/testNullModel.cpp b/tests/modules/testNullModel.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..349f193470f7dd96286413afd14fe4354a91872e
--- /dev/null
+++ b/tests/modules/testNullModel.cpp
@@ -0,0 +1,57 @@
+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#include <catch2/catch.hpp>
+
+#include <corsika/modules/null_model/NullModel.hpp>
+
+#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/framework/geometry/RootCoordinateSystem.hpp>
+#include <corsika/framework/geometry/Vector.hpp>
+
+#include <corsika/framework/core/PhysicalUnits.hpp>
+
+#include <corsika/setup/SetupStack.hpp>
+#include <corsika/setup/SetupTrajectory.hpp>
+
+using namespace corsika::units::si;
+using namespace corsika::null_model;
+using namespace corsika;
+
+TEST_CASE("NullModel", "[processes]") {
+
+ auto const& dummyCS =
+ corsika::RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
+ corsika::Point const origin(dummyCS, {0_m, 0_m, 0_m});
+ corsika::Vector<SpeedType::dimension_type> v(dummyCS, 0_m / second,
+ 0_m / second, 1_m / second);
+ corsika::Line line(origin, v);
+ corsika::Trajectory<corsika::Line> track(line, 10_s);
+
+ setup::Stack stack;
+ setup::Stack::ParticleType particle = stack.AddParticle(
+ std::tuple<Code, units::si::HEPEnergyType,
+ corsika::MomentumVector, corsika::Point, units::si::TimeType>{
+ Code::Electron, 100_GeV,
+ corsika::MomentumVector(dummyCS, {0_GeV, 0_GeV, -1_GeV}),
+ corsika::Point(dummyCS, {0_m, 0_m, 10_km}), 0_ns});
+
+ SECTION("interface") {
+
+ NullModel model(10_m);
+
+ model.Init();
+ [[maybe_unused]] const EProcessReturn ret =
+ model.DoContinuous(particle, track);
+ LengthType const length = model.MaxStepLength(particle, track);
+
+ CHECK((length / 10_m) == Approx(1));
+ }
+}
diff --git a/tests/modules/testObservationPlane.cpp b/tests/modules/testObservationPlane.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..ba0733da78d7ae1c6042731b0b1c88746d71363e
--- /dev/null
+++ b/tests/modules/testObservationPlane.cpp
@@ -0,0 +1,94 @@
+/*
+ * (c) Copyright 2019 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#include <catch2/catch.hpp>
+
+#include <corsika/modules/observation_plane/ObservationPlane.hpp>
+
+#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/framework/geometry/RootCoordinateSystem.hpp>
+#include <corsika/framework/geometry/Vector.hpp>
+
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+
+using namespace corsika::units::si;
+using namespace corsika::observation_plane;
+using namespace corsika;
+
+TEST_CASE("ContinuousProcess interface", "[proccesses][observation_plane]") {
+
+ auto const& rootCS = RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
+
+ /*
+ Test with downward going 1_GeV neutrino, starting at 0,1_m,10m
+
+ ObservationPlane has origin at 0,0,0
+ */
+
+ Point const start(rootCS, {0_m, 1_m, 10_m});
+ Vector<units::si::SpeedType::dimension_type> vec(rootCS, 0_m / second, 0_m / second,
+ -units::constants::c);
+ Line line(start, vec);
+ Trajectory<Line> track(line, 12_m / units::constants::c);
+
+ // setup particle stack, and add primary particle
+ setup::Stack stack;
+ stack.Clear();
+ {
+ auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
+ return sqrt((Elab - m) * (Elab + m));
+ };
+ stack.AddParticle(
+ std::tuple<Code, units::si::HEPEnergyType, corsika::MomentumVector, Point,
+ units::si::TimeType>{
+ Code::NuMu, 1_GeV,
+ corsika::MomentumVector(
+ rootCS, {0_GeV, 0_GeV, -elab2plab(1_GeV, NuMu::GetMass())}),
+ Point(rootCS, {1_m, 1_m, 10_m}), 0_ns});
+ }
+ auto particle = stack.GetNextParticle();
+
+ SECTION("horizontal plane") {
+
+ Plane const obsPlane(Point(rootCS, {0_m, 0_m, 0_m}),
+ Vector<dimensionless_d>(rootCS, {0., 0., 1.}));
+ ObservationPlane obs(obsPlane, "particles.dat", true);
+
+ obs.Init();
+ const LengthType length = obs.MaxStepLength(particle, track);
+ const EProcessReturn ret = obs.DoContinuous(particle, track);
+
+ REQUIRE(length / 10_m == Approx(1).margin(1e-4));
+ REQUIRE(ret == EProcessReturn::eParticleAbsorbed);
+
+ /*
+ SECTION("horizontal plane") {
+ REQUIRE(true); // todo: we have to check content of output file...
+
+ }
+ */
+ }
+
+ SECTION("inclined plane") {}
+
+ SECTION("transparent plane") {
+ Plane const obsPlane(Point(rootCS, {0_m, 0_m, 0_m}),
+ Vector<dimensionless_d>(rootCS, {0., 0., 1.}));
+ ObservationPlane obs(obsPlane, "particles.dat", false);
+
+ obs.Init();
+ const LengthType length = obs.MaxStepLength(particle, track);
+ const EProcessReturn ret = obs.DoContinuous(particle, track);
+
+ REQUIRE(length / 10_m == Approx(1).margin(1e-4));
+ REQUIRE(ret == EProcessReturn::eOk);
+ }
+}
diff --git a/tests/modules/testParticleCut.cpp b/tests/modules/testParticleCut.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..2375600c8bacb5908950dbeab0efdd0b46c19b3f
--- /dev/null
+++ b/tests/modules/testParticleCut.cpp
@@ -0,0 +1,107 @@
+
+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#include <corsika/process/particle_cut/ParticleCut.h>
+
+#include <corsika/environment/Environment.h>
+#include <corsika/geometry/Point.h>
+#include <corsika/geometry/RootCoordinateSystem.h>
+#include <corsika/geometry/Vector.h>
+#include <corsika/units/PhysicalUnits.h>
+#include <corsika/utl/CorsikaFenv.h>
+
+#include <corsika/setup/SetupStack.h>
+
+#include <catch2/catch.hpp>
+
+using namespace corsika;
+using namespace corsika::process::particle_cut;
+using namespace corsika::units;
+using namespace corsika::units::si;
+
+TEST_CASE("ParticleCut", "[processes]") {
+ feenableexcept(FE_INVALID);
+ using EnvType = environment::Environment<setup::IEnvironmentModel>;
+ EnvType env;
+ const geometry::CoordinateSystem& rootCS = env.GetCoordinateSystem();
+
+ // setup empty particle stack
+ setup::Stack stack;
+ stack.Clear();
+ // two energies
+ const HEPEnergyType Eabove = 1_TeV;
+ const HEPEnergyType Ebelow = 10_GeV;
+ // list of arbitrary particles
+ std::vector<particles::Code> particleList = {
+ particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus,
+ particles::Code::KMinus, particles::Code::K0Long, particles::Code::K0Short,
+ particles::Code::Electron, particles::Code::MuPlus, particles::Code::NuE,
+ particles::Code::Neutron};
+
+ SECTION("cut on particle type") {
+
+ ParticleCut cut(20_GeV);
+
+ // add primary particle to stack
+ auto particle = stack.AddParticle(
+ std::tuple<particles::Code, units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
+ particles::Code::Proton, Eabove,
+ corsika::stack::MomentumVector(rootCS, {0_GeV, 0_GeV, 0_GeV}),
+ geometry::Point(rootCS, 0_m, 0_m, 0_m), 0_ns});
+ // view on secondary particles
+ corsika::stack::SecondaryView view(particle);
+ // ref. to primary particle through the secondary view.
+ // only this way the secondary view is populated
+ auto projectile = view.GetProjectile();
+ // add secondaries, all with energies above the threshold
+ // only cut is by species
+ for (auto proType : particleList)
+ projectile.AddSecondary(std::tuple<particles::Code, units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, geometry::Point,
+ units::si::TimeType>{
+ proType, Eabove, corsika::stack::MomentumVector(rootCS, {0_GeV, 0_GeV, 0_GeV}),
+ geometry::Point(rootCS, 0_m, 0_m, 0_m), 0_ns});
+
+ cut.DoSecondaries(view);
+
+ REQUIRE(view.GetSize() == 8);
+ }
+
+ SECTION("cut low energy") {
+ ParticleCut cut(20_GeV);
+
+ // add primary particle to stack
+ auto particle = stack.AddParticle(
+ std::tuple<particles::Code, units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
+ particles::Code::Proton, Eabove,
+ corsika::stack::MomentumVector(rootCS, {0_GeV, 0_GeV, 0_GeV}),
+ geometry::Point(rootCS, 0_m, 0_m, 0_m), 0_ns});
+ // view on secondary particles
+ corsika::stack::SecondaryView view(particle);
+ // ref. to primary particle through the secondary view.
+ // only this way the secondary view is populated
+ auto projectile = view.GetProjectile();
+ // add secondaries, all with energies below the threshold
+ // only cut is by species
+ for (auto proType : particleList)
+ projectile.AddSecondary(std::tuple<particles::Code, units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, geometry::Point,
+ units::si::TimeType>{
+ proType, Ebelow, corsika::stack::MomentumVector(rootCS, {0_GeV, 0_GeV, 0_GeV}),
+ geometry::Point(rootCS, 0_m, 0_m, 0_m), 0_ns});
+
+ cut.DoSecondaries(view);
+
+ REQUIRE(view.GetSize() == 0);
+ }
+}
diff --git a/tests/modules/testPythia.cpp b/tests/modules/testPythia.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..f229b907f92968d7c26ff7bed70b4476ce56282d
--- /dev/null
+++ b/tests/modules/testPythia.cpp
@@ -0,0 +1,165 @@
+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#include <Pythia8/Pythia.h>
+#include <corsika/process/pythia/Decay.h>
+#include <corsika/process/pythia/Interaction.h>
+
+#include <corsika/random/RNGManager.h>
+
+#include <corsika/particles/ParticleProperties.h>
+
+#include <corsika/geometry/Point.h>
+#include <corsika/units/PhysicalUnits.h>
+
+#include <catch2/catch.hpp>
+
+TEST_CASE("Pythia", "[processes]") {
+
+ SECTION("linking pythia") {
+ using namespace Pythia8;
+ using std::cout;
+ using std::endl;
+
+ // Generator. Process selection. LHC initialization. Histogram.
+ Pythia pythia;
+
+ pythia.readString("Next:numberShowInfo = 0");
+ pythia.readString("Next:numberShowProcess = 0");
+ pythia.readString("Next:numberShowEvent = 0");
+
+ pythia.readString("ProcessLevel:all = off");
+
+ pythia.init();
+
+ Event& event = pythia.event;
+ event.reset();
+
+ pythia.particleData.mayDecay(321, true);
+ double pz = 100.;
+ double m = 0.49368;
+ event.append(321, 1, 0, 0, 0., 0., 100., sqrt(pz * pz + m * m), m);
+
+ if (!pythia.next())
+ cout << "decay failed!" << endl;
+ else
+ cout << "particles after decay: " << event.size() << endl;
+ event.list();
+
+ // loop over final state
+ for (int i = 0; i < pythia.event.size(); ++i)
+ if (pythia.event[i].isFinal()) {
+ cout << "particle: id=" << pythia.event[i].id() << endl;
+ }
+ }
+
+ SECTION("pythia interface") {
+ using namespace corsika;
+
+ const std::vector<particles::Code> particleList = {
+ particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus,
+ particles::Code::KMinus, particles::Code::K0Long, particles::Code::K0Short};
+
+ random::RNGManager::GetInstance().RegisterRandomStream("pythia");
+
+ process::pythia::Decay model(particleList);
+
+ model.Init();
+ }
+}
+
+#include <corsika/geometry/Point.h>
+#include <corsika/geometry/RootCoordinateSystem.h>
+#include <corsika/geometry/Vector.h>
+
+#include <corsika/units/PhysicalUnits.h>
+
+#include <corsika/particles/ParticleProperties.h>
+#include <corsika/setup/SetupStack.h>
+#include <corsika/setup/SetupTrajectory.h>
+
+#include <corsika/environment/Environment.h>
+#include <corsika/environment/HomogeneousMedium.h>
+#include <corsika/environment/NuclearComposition.h>
+
+using namespace corsika;
+using namespace corsika::units::si;
+
+TEST_CASE("pythia process") {
+
+ // setup environment, geometry
+ environment::Environment<environment::IMediumModel> env;
+
+ 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>(
+ 1_kg / (1_m * 1_m * 1_m),
+ environment::NuclearComposition(
+ 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
+
+ SECTION("pythia decay") {
+
+ setup::Stack stack;
+ const HEPEnergyType E0 = 10_GeV;
+ HEPMomentumType P0 =
+ sqrt(E0 * E0 - particles::PiPlus::GetMass() * particles::PiPlus::GetMass());
+ 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(
+ std::tuple<particles::Code, units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
+ particles::Code::PiPlus, E0, plab, pos, 0_ns});
+
+ const std::vector<particles::Code> particleList = {
+ particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus,
+ particles::Code::KMinus, particles::Code::K0Long, particles::Code::K0Short};
+
+ random::RNGManager::GetInstance().RegisterRandomStream("pythia");
+
+ corsika::stack::SecondaryView view(particle);
+ auto projectile = view.GetProjectile();
+
+ process::pythia::Decay model(particleList);
+ model.Init();
+ model.DoDecay(projectile);
+ [[maybe_unused]] const TimeType time = model.GetLifetime(particle);
+ }
+
+ SECTION("pythia interaction") {
+
+ setup::Stack stack;
+ const HEPEnergyType E0 = 100_GeV;
+ HEPMomentumType P0 =
+ sqrt(E0 * E0 - particles::PiPlus::GetMass() * particles::PiPlus::GetMass());
+ 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(
+ 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;
+ model.Init();
+ [[maybe_unused]] const process::EProcessReturn ret = model.DoInteraction(projectile);
+ [[maybe_unused]] const GrammageType length = model.GetInteractionLength(particle);
+ }
+}
diff --git a/tests/modules/testQGSJetII.cpp b/tests/modules/testQGSJetII.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..4f284e056992a9987ae4cd70b4a5f471f5caca6a
--- /dev/null
+++ b/tests/modules/testQGSJetII.cpp
@@ -0,0 +1,126 @@
+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#include <corsika/process/qgsjetII/Interaction.h>
+#include <corsika/process/qgsjetII/ParticleConversion.h>
+
+#include <corsika/random/RNGManager.h>
+
+#include <corsika/particles/ParticleProperties.h>
+
+#include <corsika/geometry/Point.h>
+#include <corsika/units/PhysicalUnits.h>
+
+#include <catch2/catch.hpp>
+
+using namespace corsika;
+using namespace corsika::process::qgsjetII;
+
+TEST_CASE("QgsjetII", "[processes]") {
+
+ SECTION("QgsjetII -> Corsika") {
+ REQUIRE(particles::PiPlus::GetCode() == process::qgsjetII::ConvertFromQgsjetII(
+ process::qgsjetII::QgsjetIICode::PiPlus));
+ }
+
+ SECTION("Corsika -> QgsjetII") {
+ REQUIRE(process::qgsjetII::ConvertToQgsjetII(particles::PiMinus::GetCode()) ==
+ process::qgsjetII::QgsjetIICode::PiMinus);
+ REQUIRE(process::qgsjetII::ConvertToQgsjetIIRaw(particles::Proton::GetCode()) == 2);
+ }
+
+ SECTION("canInteractInQgsjetII") {
+
+ REQUIRE(process::qgsjetII::CanInteract(particles::Proton::GetCode()));
+ REQUIRE(process::qgsjetII::CanInteract(particles::Code::KPlus));
+ REQUIRE(process::qgsjetII::CanInteract(particles::Nucleus::GetCode()));
+ // REQUIRE(process::qgsjetII::CanInteract(particles::Helium::GetCode()));
+
+ REQUIRE_FALSE(process::qgsjetII::CanInteract(particles::EtaC::GetCode()));
+ REQUIRE_FALSE(process::qgsjetII::CanInteract(particles::SigmaC0::GetCode()));
+ }
+
+ SECTION("cross-section type") {
+
+ REQUIRE(process::qgsjetII::GetQgsjetIIXSCode(particles::Code::Neutron) == 2);
+ REQUIRE(process::qgsjetII::GetQgsjetIIXSCode(particles::Code::K0Long) == 3);
+ REQUIRE(process::qgsjetII::GetQgsjetIIXSCode(particles::Code::Proton) == 2);
+ REQUIRE(process::qgsjetII::GetQgsjetIIXSCode(particles::Code::PiMinus) == 1);
+ }
+}
+
+#include <corsika/geometry/Point.h>
+#include <corsika/geometry/RootCoordinateSystem.h>
+#include <corsika/geometry/Vector.h>
+
+#include <corsika/units/PhysicalUnits.h>
+
+#include <corsika/particles/ParticleProperties.h>
+#include <corsika/setup/SetupStack.h>
+#include <corsika/setup/SetupTrajectory.h>
+
+#include <corsika/environment/Environment.h>
+#include <corsika/environment/HomogeneousMedium.h>
+#include <corsika/environment/NuclearComposition.h>
+#include <corsika/process/qgsjetII/qgsjet-II-04.h>
+
+using namespace corsika::units::si;
+using namespace corsika::units;
+
+TEST_CASE("QgsjetIIInterface", "[processes]") {
+
+ // setup environment, geometry
+ environment::Environment<environment::IMediumModel> env;
+ auto& universe = *(env.GetUniverse());
+
+ auto theMedium =
+ environment::Environment<environment::IMediumModel>::CreateNode<geometry::Sphere>(
+ geometry::Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m},
+ 1_km * std::numeric_limits<double>::infinity());
+
+ using MyHomogeneousModel = environment::HomogeneousMedium<environment::IMediumModel>;
+ theMedium->SetModelProperties<MyHomogeneousModel>(
+ 1_kg / (1_m * 1_m * 1_m),
+ environment::NuclearComposition(
+ std::vector<particles::Code>{particles::Code::Oxygen}, std::vector<float>{1.}));
+
+ auto const* nodePtr = theMedium.get();
+ universe.AddChild(std::move(theMedium));
+
+ const geometry::CoordinateSystem& cs = env.GetCoordinateSystem();
+
+ random::RNGManager::GetInstance().RegisterRandomStream("qgran");
+
+ SECTION("InteractionInterface") {
+
+ setup::Stack stack;
+ const HEPEnergyType E0 = 100_GeV;
+ HEPMomentumType P0 =
+ sqrt(E0 * E0 - particles::Proton::GetMass() * particles::Proton::GetMass());
+ 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(std::tuple<particles::Code, units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, geometry::Point,
+ units::si::TimeType, unsigned int, unsigned int>{
+ particles::Code::Nucleus, E0, plab, pos, 0_ns, 16, 8});
+ // 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();
+
+ Interaction model;
+ model.Init();
+ [[maybe_unused]] const process::EProcessReturn ret = model.DoInteraction(projectile);
+ [[maybe_unused]] const GrammageType length = model.GetInteractionLength(particle);
+ }
+}
diff --git a/tests/modules/testStackInspector.cpp b/tests/modules/testStackInspector.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..1371f8bfaf48bd8a582cfca0abca152510b515d1
--- /dev/null
+++ b/tests/modules/testStackInspector.cpp
@@ -0,0 +1,55 @@
+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#include <catch2/catch.hpp>
+
+#include <corsika/process/stack_inspector/StackInspector.h>
+
+#include <corsika/geometry/Point.h>
+#include <corsika/geometry/RootCoordinateSystem.h>
+#include <corsika/geometry/Vector.h>
+
+#include <corsika/units/PhysicalUnits.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]") {
+
+ auto const& rootCS =
+ geometry::RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
+ geometry::Point const origin(rootCS, {0_m, 0_m, 0_m});
+ 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);
+
+ 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, E0,
+ corsika::stack::MomentumVector(rootCS, {0_GeV, 0_GeV, -1_GeV}),
+ Point(rootCS, {0_m, 0_m, 10_km}), 0_ns});
+
+ SECTION("interface") {
+
+ StackInspector<TestCascadeStack> model(1, true, E0);
+
+ model.Init();
+ [[maybe_unused]] const process::EProcessReturn ret = model.DoStack(stack);
+ }
+}
diff --git a/tests/modules/testSwitchProcess.cpp b/tests/modules/testSwitchProcess.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..510c731ff633331b8996326d3c940c68025b4796
--- /dev/null
+++ b/tests/modules/testSwitchProcess.cpp
@@ -0,0 +1,264 @@
+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#include <corsika/process/switch_process/SwitchProcess.h>
+#include <corsika/stack/SecondaryView.h>
+#include <corsika/stack/Stack.h>
+#include <corsika/units/PhysicalUnits.h>
+
+#include <catch2/catch.hpp>
+
+#include <algorithm>
+#include <random>
+
+using namespace corsika;
+using namespace corsika::process;
+using namespace corsika::units::si;
+
+class TestStackData {
+
+public:
+ // these functions are needed for the Stack interface
+ void Init() {}
+ void Clear() { fData.clear(); }
+ unsigned int GetSize() const { return fData.size(); }
+ unsigned int GetCapacity() const { return fData.size(); }
+ void Copy(int i1, int i2) { fData[i2] = fData[i1]; }
+ void Swap(int i1, int i2) { std::swap(fData[i1], fData[i2]); }
+
+ // custom data access function
+ void SetData(unsigned int i, HEPEnergyType v) { fData[i] = v; }
+ HEPEnergyType GetData(const unsigned int i) const { return fData[i]; }
+
+ // these functions are also needed by the Stack interface
+ void IncrementSize() { fData.resize(fData.size() + 1); }
+ void DecrementSize() {
+ if (fData.size() > 0) { fData.pop_back(); }
+ }
+
+ // custom private data section
+private:
+ std::vector<HEPEnergyType> fData;
+};
+
+/**
+ * From static_cast of a StackIterator over entries in the
+ * TestStackData class you get and object of type
+ * TestParticleInterface defined here
+ *
+ * It provides Get/Set methods to read and write data to the "Data"
+ * storage of TestStackData obtained via
+ * "StackIteratorInterface::GetStackData()", given the index of the
+ * iterator "StackIteratorInterface::GetIndex()"
+ *
+ */
+template <typename StackIteratorInterface>
+class TestParticleInterface
+ : public corsika::stack::ParticleBase<StackIteratorInterface> {
+
+public:
+ using corsika::stack::ParticleBase<StackIteratorInterface>::GetStackData;
+ using corsika::stack::ParticleBase<StackIteratorInterface>::GetIndex;
+
+ /*
+ The SetParticleData methods are called for creating new entries
+ on the stack. You can specifiy various parametric versions to
+ perform this task:
+ */
+
+ // default version for particle-creation from input data
+ void SetParticleData(const std::tuple<HEPEnergyType> v) { SetEnergy(std::get<0>(v)); }
+ void SetParticleData(TestParticleInterface<StackIteratorInterface>& /*parent*/,
+ std::tuple<HEPEnergyType> v) {
+ SetEnergy(std::get<0>(v));
+ }
+
+ // here are the fundamental methods for access to TestStackData data
+ void SetEnergy(HEPEnergyType v) { GetStackData().SetData(GetIndex(), v); }
+ HEPEnergyType GetEnergy() const { return GetStackData().GetData(GetIndex()); }
+};
+
+using SimpleStack = corsika::stack::Stack<TestStackData, TestParticleInterface>;
+
+// see issue 161
+#if defined(__clang__)
+using StackTestView = corsika::stack::SecondaryView<TestStackData, TestParticleInterface>;
+#elif defined(__GNUC__) || defined(__GNUG__)
+using StackTestView = corsika::stack::MakeView<SimpleStack>::type;
+#endif
+
+auto constexpr kgMSq = 1_kg / (1_m * 1_m);
+
+template <int N>
+struct DummyProcess : InteractionProcess<DummyProcess<N>> {
+ void Init() {}
+
+ template <typename TParticle>
+ corsika::units::si::GrammageType GetInteractionLength(TParticle const&) const {
+ return N * kgMSq;
+ }
+
+ template <typename TSecondaries>
+ corsika::process::EProcessReturn DoInteraction(TSecondaries& vSec) {
+ // to figure out which process was selected in the end, we produce N
+ // secondaries for DummyProcess<N>
+
+ for (int i = 0; i < N; ++i) {
+ vSec.AddSecondary(std::tuple<HEPEnergyType>{vSec.GetEnergy() / N});
+ }
+
+ return EProcessReturn::eOk;
+ }
+};
+
+using DummyLowEnergyProcess = DummyProcess<1>;
+using DummyHighEnergyProcess = DummyProcess<2>;
+using DummyAdditionalProcess = DummyProcess<3>;
+
+TEST_CASE("SwitchProcess from InteractionProcess") {
+ DummyLowEnergyProcess low;
+ DummyHighEnergyProcess high;
+ DummyAdditionalProcess proc;
+
+ switch_process::SwitchProcess switchProcess(low, high, 1_TeV);
+ auto seq = switchProcess << proc;
+
+ SimpleStack stack;
+
+ SECTION("low energy") {
+ stack.AddParticle(std::tuple<HEPEnergyType>{0.5_TeV});
+ auto p = stack.GetNextParticle();
+
+ // low energy process returns 1 kg/m²
+ SECTION("interaction length") {
+ REQUIRE(switchProcess.GetInteractionLength(p) / kgMSq == Approx(1));
+ REQUIRE(seq.GetTotalInteractionLength(p) / kgMSq == Approx(3. / 4));
+ }
+
+ // low energy process creates 1 secondary
+ //~ SECTION("SelectInteraction") {
+ //~ typename SimpleStack::ParticleType theParticle =
+ //~ stack.GetNextParticle(); // as in corsika::Cascade
+ //~ StackTestView view(theParticle);
+ //~ auto projectile = view.GetProjectile();
+
+ //~ InverseGrammageType invLambda = 0 / kgMSq;
+ //~ switchProcess.SelectInteraction(p, projectile, 0.01 / kgMSq, invLambda);
+
+ //~ REQUIRE(view.GetSize() == 1);
+ //~ }
+ }
+
+ SECTION("high energy") {
+ stack.AddParticle(std::tuple<HEPEnergyType>{4_TeV});
+ auto p = stack.GetNextParticle();
+
+ // high energy process returns 2 kg/m²
+ SECTION("interaction length") {
+ REQUIRE(switchProcess.GetInteractionLength(p) / kgMSq == Approx(2));
+ REQUIRE(seq.GetTotalInteractionLength(p) / kgMSq == Approx(6. / 5));
+ }
+
+ // high energy process creates 2 secondaries
+ SECTION("SelectInteraction") {
+ typename SimpleStack::ParticleType theParticle =
+ stack.GetNextParticle(); // as in corsika::Cascade
+ StackTestView view(theParticle);
+ auto projectile = view.GetProjectile();
+
+ InverseGrammageType invLambda = 0 / kgMSq;
+ switchProcess.SelectInteraction(p, projectile, 0.01 / kgMSq, invLambda);
+
+ REQUIRE(view.GetSize() == 2);
+ }
+ }
+}
+
+TEST_CASE("SwitchProcess from ProcessSequence") {
+ DummyProcess<1> innerA;
+ DummyProcess<2> innerB;
+ DummyProcess<3> outer;
+ DummyProcess<4> additional;
+
+ auto seq = innerA << innerB;
+
+ switch_process::SwitchProcess switchProcess(seq, outer, 1_TeV);
+ auto completeSeq = switchProcess << additional;
+
+ SimpleStack stack;
+
+ SECTION("low energy") {
+ stack.AddParticle(std::tuple<HEPEnergyType>{0.5_TeV});
+ auto p = stack.GetNextParticle();
+
+ SECTION("interaction length") {
+ REQUIRE(switchProcess.GetInteractionLength(p) / kgMSq == Approx(2. / 3));
+ REQUIRE(completeSeq.GetTotalInteractionLength(p) / kgMSq == Approx(4. / 7));
+ }
+
+ SECTION("SelectInteraction") {
+ std::vector<int> numberOfSecondaries;
+
+ for (int i = 0; i < 1000; ++i) {
+ typename SimpleStack::ParticleType theParticle =
+ stack.GetNextParticle(); // as in corsika::Cascade
+ StackTestView view(theParticle);
+ auto projectile = view.GetProjectile();
+
+ double r = i / 1000.;
+ InverseGrammageType invLambda = r * 7. / 4 / kgMSq;
+
+ InverseGrammageType accumulator = 0 / kgMSq;
+ completeSeq.SelectInteraction(p, projectile, invLambda, accumulator);
+
+ numberOfSecondaries.push_back(view.GetSize());
+ }
+
+ auto const mean =
+ std::accumulate(numberOfSecondaries.cbegin(), numberOfSecondaries.cend(), 0.) /
+ numberOfSecondaries.size();
+ REQUIRE(mean == Approx(12. / 7.).margin(0.01));
+ }
+ }
+
+ SECTION("high energy") {
+ stack.AddParticle(std::tuple<HEPEnergyType>{3.0_TeV});
+ auto p = stack.GetNextParticle();
+
+ SECTION("interaction length") {
+ REQUIRE(switchProcess.GetInteractionLength(p) / kgMSq == Approx(3));
+ REQUIRE(completeSeq.GetTotalInteractionLength(p) / kgMSq == Approx(12. / 7.));
+ }
+
+ SECTION("SelectInteraction") {
+ std::vector<int> numberOfSecondaries;
+
+ for (int i = 0; i < 1000; ++i) {
+ typename SimpleStack::ParticleType theParticle =
+ stack.GetNextParticle(); // as in corsika::Cascade
+ StackTestView view(theParticle);
+ auto projectile = view.GetProjectile();
+
+ double r = i / 1000.;
+ InverseGrammageType invLambda = r * 7. / 12. / kgMSq;
+
+ InverseGrammageType accumulator = 0 / kgMSq;
+ completeSeq.SelectInteraction(p, projectile, invLambda, accumulator);
+
+ numberOfSecondaries.push_back(view.GetSize());
+ }
+
+ auto const mean =
+ std::accumulate(numberOfSecondaries.cbegin(), numberOfSecondaries.cend(), 0.) /
+ numberOfSecondaries.size();
+ REQUIRE(mean == Approx(24. / 7.).margin(0.01));
+ }
+ }
+}
diff --git a/tests/modules/testTrackingLine.cpp b/tests/modules/testTrackingLine.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..638e43bc7f83f8f45d1446b75dbfdb1d5092dee2
--- /dev/null
+++ b/tests/modules/testTrackingLine.cpp
@@ -0,0 +1,100 @@
+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#include <corsika/process/tracking_line/TrackingLine.h>
+#include <testTrackingLineStack.h> // test-build, and include file is obtained from CMAKE_CURRENT_SOURCE_DIR
+
+#include <corsika/environment/Environment.h>
+#include <corsika/particles/ParticleProperties.h>
+
+#include <corsika/geometry/Point.h>
+#include <corsika/geometry/Sphere.h>
+#include <corsika/geometry/Vector.h>
+
+#include <corsika/setup/SetupTrajectory.h>
+using corsika::setup::Trajectory;
+
+#include <catch2/catch.hpp>
+
+using namespace corsika;
+using namespace corsika::process;
+using namespace corsika::units;
+using namespace corsika::geometry;
+
+#include <iostream>
+using namespace std;
+using namespace corsika::units::si;
+
+TEST_CASE("TrackingLine") {
+ environment::Environment<environment::Empty> env; // dummy environment
+ auto const& cs = env.GetCoordinateSystem();
+
+ tracking_line::TrackingLine tracking;
+
+ SECTION("intersection with sphere") {
+ Point const origin(cs, {0_m, 0_m, -5_m});
+ Point const center(cs, {0_m, 0_m, 10_m});
+ Sphere const sphere(center, 1_m);
+ Vector<corsika::units::si::SpeedType::dimension_type> v(cs, 0_m / second,
+ 0_m / second, 1_m / second);
+ Line line(origin, v);
+
+ setup::Trajectory traj(line, 12345_s);
+
+ auto const opt =
+ tracking_line::TimeOfIntersection(traj, Sphere(Point(cs, {0_m, 0_m, 10_m}), 1_m));
+ REQUIRE(opt.has_value());
+
+ auto [t1, t2] = opt.value();
+ REQUIRE(t1 / 14_s == Approx(1));
+ REQUIRE(t2 / 16_s == Approx(1));
+
+ auto const optNoIntersection =
+ tracking_line::TimeOfIntersection(traj, Sphere(Point(cs, {5_m, 0_m, 10_m}), 1_m));
+ REQUIRE_FALSE(optNoIntersection.has_value());
+ }
+
+ SECTION("maximally possible propagation") {
+ auto& universe = *(env.GetUniverse());
+
+ auto const radius = 20_m;
+
+ auto theMedium = environment::Environment<environment::Empty>::CreateNode<Sphere>(
+ Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m}, radius);
+ auto const* theMediumPtr = theMedium.get();
+ universe.AddChild(std::move(theMedium));
+
+ TestTrackingLineStack stack;
+ stack.AddParticle(
+ std::tuple<particles::Code, units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
+ particles::Code::MuPlus,
+ 1_GeV,
+ {cs, {0_GeV, 0_GeV, 1_GeV}},
+ {cs, {0_m, 0_m, 0_km}},
+ 0_ns});
+ auto p = stack.GetNextParticle();
+ p.SetNode(theMediumPtr);
+
+ Point const origin(cs, {0_m, 0_m, 0_m});
+ Vector<corsika::units::si::SpeedType::dimension_type> v(cs, 0_m / second,
+ 0_m / second, 1_m / second);
+ 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)
+ .norm()
+ .magnitude() == Approx(0).margin(1e-4));
+ }
+}
diff --git a/tests/modules/testUrQMD.cpp b/tests/modules/testUrQMD.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..3339fb9d32cc84939a3d9c4b0bbb6dbc031c4012
--- /dev/null
+++ b/tests/modules/testUrQMD.cpp
@@ -0,0 +1,227 @@
+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#include <corsika/process/urqmd/UrQMD.h>
+#include <corsika/random/RNGManager.h>
+
+#include <corsika/geometry/Point.h>
+#include <corsika/geometry/RootCoordinateSystem.h>
+#include <corsika/geometry/Vector.h>
+
+#include <corsika/units/PhysicalConstants.h>
+#include <corsika/units/PhysicalUnits.h>
+
+#include <corsika/utl/CorsikaFenv.h>
+
+#include <corsika/particles/ParticleProperties.h>
+#include <corsika/setup/SetupStack.h>
+#include <corsika/setup/SetupTrajectory.h>
+
+#include <corsika/environment/Environment.h>
+#include <corsika/environment/HomogeneousMedium.h>
+#include <corsika/environment/NuclearComposition.h>
+
+#include <tuple>
+#include <utility>
+
+#include <catch2/catch.hpp>
+
+using namespace corsika;
+using namespace corsika::process::UrQMD;
+using namespace corsika::units::si;
+
+template <typename TStackView>
+auto sumCharge(TStackView const& view) {
+ int totalCharge = 0;
+
+ for (auto const& p : view) { totalCharge += particles::GetChargeNumber(p.GetPID()); }
+
+ return totalCharge;
+}
+
+template <typename TStackView>
+auto sumMomentum(TStackView const& view, geometry::CoordinateSystem const& vCS) {
+ geometry::Vector<hepenergy_d> sum{vCS, 0_eV, 0_eV, 0_eV};
+
+ for (auto const& p : view) { sum += p.GetMomentum(); }
+
+ return sum;
+}
+
+auto setupEnvironment(particles::Code vTargetCode) {
+ // setup environment, geometry
+ auto env = std::make_unique<environment::Environment<environment::IMediumModel>>();
+ auto& universe = *(env->GetUniverse());
+ const geometry::CoordinateSystem& 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>(
+ 1_kg / (1_m * 1_m * 1_m),
+ environment::NuclearComposition(std::vector<particles::Code>{vTargetCode},
+ std::vector<float>{1.}));
+
+ auto const* nodePtr = theMedium.get();
+ universe.AddChild(std::move(theMedium));
+
+ return std::make_tuple(std::move(env), &cs, nodePtr);
+}
+
+template <typename TNodeType>
+auto setupStack(int vA, int vZ, HEPEnergyType vMomentum, TNodeType* vNodePtr,
+ geometry::CoordinateSystem const& cs) {
+ auto stack = std::make_unique<setup::Stack>();
+ auto constexpr mN = corsika::units::constants::nucleonMass;
+
+ geometry::Point const origin(cs, {0_m, 0_m, 0_m});
+ corsika::stack::MomentumVector const pLab(cs, {vMomentum, 0_GeV, 0_GeV});
+
+ HEPEnergyType const E0 =
+ sqrt(units::si::detail::static_pow<2>(mN * vA) + pLab.squaredNorm());
+ 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, origin, 0_ns, vA, vZ});
+
+ particle.SetNode(vNodePtr);
+ return std::make_tuple(
+ std::move(stack),
+ std::make_unique<decltype(corsika::stack::SecondaryView(particle))>(particle));
+}
+
+template <typename TNodeType>
+auto setupStack(particles::Code vProjectileType, HEPEnergyType vMomentum,
+ TNodeType* vNodePtr, geometry::CoordinateSystem const& cs) {
+ auto stack = std::make_unique<setup::Stack>();
+
+ geometry::Point const origin(cs, {0_m, 0_m, 0_m});
+ corsika::stack::MomentumVector const pLab(cs, {vMomentum, 0_GeV, 0_GeV});
+
+ HEPEnergyType const E0 =
+ sqrt(units::si::detail::static_pow<2>(particles::GetMass(vProjectileType)) +
+ pLab.squaredNorm());
+ auto particle = stack->AddParticle(
+ std::tuple<particles::Code, units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
+ vProjectileType, E0, pLab, origin, 0_ns});
+
+ particle.SetNode(vNodePtr);
+ return std::make_tuple(
+ std::move(stack),
+ std::make_unique<decltype(corsika::stack::SecondaryView(particle))>(particle));
+}
+
+TEST_CASE("UrQMD") {
+ SECTION("conversion") {
+ REQUIRE_THROWS(process::UrQMD::ConvertFromUrQMD(106, 0));
+ REQUIRE(process::UrQMD::ConvertFromUrQMD(101, 0) == particles::Code::Pi0);
+ REQUIRE(process::UrQMD::ConvertToUrQMD(particles::Code::PiPlus) ==
+ std::make_pair<int, int>(101, 2));
+ }
+
+ feenableexcept(FE_INVALID);
+ corsika::random::RNGManager::GetInstance().RegisterRandomStream("UrQMD");
+ UrQMD urqmd;
+
+ SECTION("cross sections") {
+ auto [env, csPtr, nodePtr] = setupEnvironment(particles::Code::Unknown);
+ auto const& cs = *csPtr;
+
+ particles::Code validProjectileCodes[] = {
+ particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::Proton,
+ particles::Code::Neutron, particles::Code::KPlus, particles::Code::KMinus,
+ particles::Code::K0, particles::Code::K0Bar, particles::Code::K0Long};
+
+ for (auto code : validProjectileCodes) {
+ auto [stack, view] = setupStack(code, 100_GeV, nodePtr, cs);
+ REQUIRE(stack->GetSize() == 1);
+
+ // simple check whether the cross-section is non-vanishing
+ REQUIRE(urqmd.GetCrossSection(view->GetProjectile(), particles::Code::Proton) /
+ 1_mb >
+ 0);
+ REQUIRE(urqmd.GetCrossSection(view->GetProjectile(), particles::Code::Nitrogen) /
+ 1_mb >
+ 0);
+ REQUIRE(urqmd.GetCrossSection(view->GetProjectile(), particles::Code::Oxygen) /
+ 1_mb >
+ 0);
+ REQUIRE(urqmd.GetCrossSection(view->GetProjectile(), particles::Code::Argon) /
+ 1_mb >
+ 0);
+ }
+ }
+
+ SECTION("nucleon projectile") {
+ auto [env, csPtr, nodePtr] = setupEnvironment(particles::Code::Oxygen);
+ unsigned short constexpr A = 14, Z = 7;
+ auto [stackPtr, secViewPtr] = setupStack(A, Z, 400_GeV, nodePtr, *csPtr);
+
+ // must be assigned to variable, cannot be used as rvalue?!
+ auto projectile = secViewPtr->GetProjectile();
+ auto const projectileMomentum = projectile.GetMomentum();
+ [[maybe_unused]] process::EProcessReturn const ret = urqmd.DoInteraction(projectile);
+
+ REQUIRE(sumCharge(*secViewPtr) ==
+ Z + particles::GetChargeNumber(particles::Code::Oxygen));
+
+ auto const secMomSum =
+ sumMomentum(*secViewPtr, projectileMomentum.GetCoordinateSystem());
+ REQUIRE((secMomSum - projectileMomentum).norm() / projectileMomentum.norm() ==
+ Approx(0).margin(1e-2));
+ }
+
+ SECTION("\"special\" projectile") {
+ auto [env, csPtr, nodePtr] = setupEnvironment(particles::Code::Oxygen);
+ auto [stackPtr, secViewPtr] =
+ setupStack(particles::Code::PiPlus, 400_GeV, nodePtr, *csPtr);
+
+ // must be assigned to variable, cannot be used as rvalue?!
+ auto projectile = secViewPtr->GetProjectile();
+ auto const projectileMomentum = projectile.GetMomentum();
+
+ [[maybe_unused]] process::EProcessReturn const ret = urqmd.DoInteraction(projectile);
+
+ REQUIRE(sumCharge(*secViewPtr) ==
+ particles::GetChargeNumber(particles::Code::PiPlus) +
+ particles::GetChargeNumber(particles::Code::Oxygen));
+
+ auto const secMomSum =
+ sumMomentum(*secViewPtr, projectileMomentum.GetCoordinateSystem());
+ REQUIRE((secMomSum - projectileMomentum).norm() / projectileMomentum.norm() ==
+ Approx(0).margin(1e-2));
+ }
+
+ SECTION("K0Long projectile") {
+ auto [env, csPtr, nodePtr] = setupEnvironment(particles::Code::Oxygen);
+ auto [stackPtr, secViewPtr] =
+ setupStack(particles::Code::K0Long, 400_GeV, nodePtr, *csPtr);
+
+ // must be assigned to variable, cannot be used as rvalue?!
+ auto projectile = secViewPtr->GetProjectile();
+ auto const projectileMomentum = projectile.GetMomentum();
+
+ [[maybe_unused]] process::EProcessReturn const ret = urqmd.DoInteraction(projectile);
+
+ REQUIRE(sumCharge(*secViewPtr) ==
+ particles::GetChargeNumber(particles::Code::K0Long) +
+ particles::GetChargeNumber(particles::Code::Oxygen));
+
+ auto const secMomSum =
+ sumMomentum(*secViewPtr, projectileMomentum.GetCoordinateSystem());
+ REQUIRE((secMomSum - projectileMomentum).norm() / projectileMomentum.norm() ==
+ Approx(0).margin(1e-2));
+ }
+}