diff --git a/corsika/detail/media/LayeredSphericalAtmosphereBuilder.hpp b/corsika/detail/media/LayeredSphericalAtmosphereBuilder.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..0bf88ff6743b9fc90271d509730858563d9d58d3
--- /dev/null
+++ b/corsika/detail/media/LayeredSphericalAtmosphereBuilder.hpp
@@ -0,0 +1,20 @@
+#pragma once
+
+namespace corsika {
+
+ namespace detail {
+
+ struct NoExtraModelInner {};
+
+ template <typename M>
+ struct NoExtraModel {};
+
+ template <template <typename> typename M>
+ struct has_extra_models : std::true_type {};
+
+ template <>
+ struct has_extra_models<NoExtraModel> : std::false_type {};
+
+ } // namespace detail
+
+}
diff --git a/corsika/detail/modules/qgsjetII/QGSJetIIStack.inl b/corsika/detail/modules/qgsjetII/QGSJetIIStack.inl
new file mode 100644
index 0000000000000000000000000000000000000000..8405888ccf2717f196ec3d4736fe50713970af6e
--- /dev/null
+++ b/corsika/detail/modules/qgsjetII/QGSJetIIStack.inl
@@ -0,0 +1,105 @@
+#pragma once
+
+namespace corsika::qgsjetII {
+
+ void QGSJetIIStackData::clear() {
+ qgarr12_.nsp = 0;
+ qgarr13_.nsf = 0;
+ qgarr55_.nwt = 0;
+ }
+ unsigned int QGSJetIIStackData::getSize() const { return qgarr12_.nsp; }
+ unsigned int QGSJetIIStackData::getCapacity() const { return nptmax; }
+
+ void QGSJetIIStackData::setId(const unsigned int i, const int v) {
+ qgarr14_.ich[i] = v;
+ }
+ void QGSJetIIStackData::setEnergy(const unsigned int i, const HEPEnergyType v) {
+ qgarr14_.esp[i][0] = v / 1_GeV;
+ }
+
+ void QGSJetIIStackData::setMomentum(const unsigned int i, const MomentumVector& v) {
+ auto tmp = v.getComponents();
+ qgarr14_.esp[i][2] = tmp[0] / 1_GeV;
+ qgarr14_.esp[i][3] = tmp[1] / 1_GeV;
+ qgarr14_.esp[i][1] = tmp[2] / 1_GeV;
+ }
+
+ int QGSJetIIStackData::getId(const unsigned int i) const { return qgarr14_.ich[i]; }
+ HEPEnergyType QGSJetIIStackData::getEnergy(const int i) const {
+ return qgarr14_.esp[i][0] * 1_GeV;
+ }
+ MomentumVector QGSJetIIStackData::getMomentum(
+ const unsigned int i, const CoordinateSystemPtr& CS) const {
+ QuantityVector<hepmomentum_d> components = {qgarr14_.esp[i][2] * 1_GeV,
+ qgarr14_.esp[i][3] * 1_GeV,
+ qgarr14_.esp[i][1] * 1_GeV};
+ return MomentumVector(CS, components);
+ }
+
+ void QGSJetIIStackData::copy(const unsigned int i1, const unsigned int i2) {
+ qgarr14_.ich[i2] = qgarr14_.ich[i1];
+ for (unsigned int i = 0; i < 4; ++i) qgarr14_.esp[i2][i] = qgarr14_.esp[i1][i];
+ }
+
+ void QGSJetIIStackData::swap(const unsigned int i1, const unsigned int i2) {
+ std::swap(qgarr14_.ich[i1], qgarr14_.ich[i2]);
+ for (unsigned int i = 0; i < 4; ++i)
+ std::swap(qgarr14_.esp[i1][i], qgarr14_.esp[i2][i]);
+ }
+
+ void QGSJetIIStackData::incrementSize() { qgarr12_.nsp++; }
+ void QGSJetIIStackData::decrementSize() {
+ if (qgarr12_.nsp > 0) { qgarr12_.nsp--; }
+ }
+
+ template <typename StackIteratorInterface>
+ void ParticleInterface<StackIteratorInterface>::setParticleData(
+ const int vID, const HEPEnergyType vE, const MomentumVector& vP,
+ const HEPMassType) {
+ setPID(vID);
+ setEnergy(vE);
+ setMomentum(vP);
+ }
+
+ template <typename StackIteratorInterface>
+ void ParticleInterface<StackIteratorInterface>::setParticleData(
+ ParticleInterface<StackIteratorInterface>& /*parent*/, const int vID,
+ const HEPEnergyType vE, const MomentumVector& vP, const HEPMassType) {
+ setPID(vID);
+ setEnergy(vE);
+ setMomentum(vP);
+ }
+
+ template <typename StackIteratorInterface>
+ void ParticleInterface<StackIteratorInterface>::setEnergy(const HEPEnergyType v) {
+ getStackData().setEnergy(getIndex(), v);
+ }
+
+ template <typename StackIteratorInterface>
+ HEPEnergyType ParticleInterface<StackIteratorInterface>::getEnergy() const {
+ return getStackData().getEnergy(getIndex());
+ }
+
+ template <typename StackIteratorInterface>
+ void ParticleInterface<StackIteratorInterface>::setPID(const int v) {
+ getStackData().setId(getIndex(), v);
+ }
+
+ template <typename StackIteratorInterface>
+ corsika::qgsjetII::QgsjetIICode ParticleInterface<StackIteratorInterface>::getPID()
+ const {
+ return static_cast<corsika::qgsjetII::QgsjetIICode>(getStackData().getId(getIndex()));
+ }
+
+ template <typename StackIteratorInterface>
+ MomentumVector ParticleInterface<StackIteratorInterface>::getMomentum(
+ const CoordinateSystemPtr& CS) const {
+ return getStackData().getMomentum(getIndex(), CS);
+ }
+
+ template <typename StackIteratorInterface>
+ void ParticleInterface<StackIteratorInterface>::setMomentum(const MomentumVector& v) {
+ getStackData().setMomentum(getIndex(), v);
+ }
+
+} // namespace corsika::qgsjetII
diff --git a/corsika/framework/core/PhysicalGeometry.hpp b/corsika/framework/core/PhysicalGeometry.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..fb1645786db242dc0625c2324e2574dedb57eccf
--- /dev/null
+++ b/corsika/framework/core/PhysicalGeometry.hpp
@@ -0,0 +1,26 @@
+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * 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/core/PhysicalUnits.hpp>
+#include <corsika/framework/geometry/Vector.hpp>
+
+/**
+ * \file PhysicalUnits.hpp
+ *
+ * Import and extend the phys::units package. The SI units are also imported into the
+ * `\namespace corsika`, since they are used everywhere as integral part of the framework.
+ */
+
+namespace corsika {
+
+ typedef Vector<hepmomentum_d> MomentumVector;
+
+
+} // namespace corsika
diff --git a/corsika/framework/geometry/PhysicalGeometry.hpp b/corsika/framework/geometry/PhysicalGeometry.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..fb1645786db242dc0625c2324e2574dedb57eccf
--- /dev/null
+++ b/corsika/framework/geometry/PhysicalGeometry.hpp
@@ -0,0 +1,26 @@
+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * 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/core/PhysicalUnits.hpp>
+#include <corsika/framework/geometry/Vector.hpp>
+
+/**
+ * \file PhysicalUnits.hpp
+ *
+ * Import and extend the phys::units package. The SI units are also imported into the
+ * `\namespace corsika`, since they are used everywhere as integral part of the framework.
+ */
+
+namespace corsika {
+
+ typedef Vector<hepmomentum_d> MomentumVector;
+
+
+} // namespace corsika
diff --git a/corsika/stack/history/CMakeLists.txt b/corsika/stack/history/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..2a47b86189bef946269e8f2815a51b2cfcf171f6
--- /dev/null
+++ b/corsika/stack/history/CMakeLists.txt
@@ -0,0 +1,69 @@
+set (
+ HISTORY_HEADERS
+ EventType.hpp
+ Event.hpp
+ HistorySecondaryProducer.hpp
+ HistoryObservationPlane.hpp
+ HistoryStackExtension.hpp
+ SecondaryParticle.hpp
+ )
+
+set (
+ HISTORY_NAMESPACE
+ corsika/stack/history
+ )
+
+if (WITH_HISTORY)
+ set (
+ HISTORY_SOURCES
+ HistoryObservationPlane.cpp
+ )
+ add_library (CORSIKAhistory STATIC ${HISTORY_SOURCES})
+ set (IS_INTERFACE "")
+else (WITH_HISTORY)
+ add_library (CORSIKAhistory INTERFACE)
+ set (IS_INTERFACE "INTERFACE")
+endif (WITH_HISTORY)
+
+CORSIKA_COPY_HEADERS_TO_NAMESPACE (CORSIKAhistory ${HISTORY_NAMESPACE} ${HISTORY_HEADERS})
+
+target_link_libraries (
+ CORSIKAhistory
+ ${IS_INTERFACE}
+ CORSIKAparticles
+ CORSIKAgeometry
+ CORSIKAprocesssequence # for HistoryObservationPlane
+ CORSIKAsetup # for HistoryObservationPlane
+ CORSIKAlogging
+ SuperStupidStack
+ NuclearStackExtension # for testHistoryView.cc
+ C8::ext::boost # for HistoryObservationPlane
+ )
+
+target_include_directories (
+ CORSIKAhistory
+ INTERFACE
+ $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
+ $<INSTALL_INTERFACE:include/include>
+ )
+
+install (
+ FILES ${HISTORY_HEADERS}
+ DESTINATION include/${HISTORY_NAMESPACE}
+ )
+
+# ----------------
+# code unit testing
+CORSIKA_ADD_TEST(testHistoryStack)
+target_link_libraries (
+ testHistoryStack
+ CORSIKAhistory
+ CORSIKAtesting
+ DummyStack
+ )
+CORSIKA_ADD_TEST(testHistoryView)
+target_link_libraries (
+ testHistoryView
+ CORSIKAhistory
+ CORSIKAtesting
+ )
diff --git a/corsika/stack/history/Event.hpp b/corsika/stack/history/Event.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..c00f5511b3ce91f748cdcb372d345f6d4da8ff13
--- /dev/null
+++ b/corsika/stack/history/Event.hpp
@@ -0,0 +1,68 @@
+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * 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/logging/Logging.hpp>
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <corsika/stack/history/EventType.hpp>
+#include <corsika/stack/history/SecondaryParticle.hpp>
+
+#include <memory>
+#include <optional>
+#include <vector>
+
+namespace corsika::history {
+
+ class Event;
+ using EventPtr = std::shared_ptr<history::Event>;
+
+ class Event {
+
+ size_t projectile_index_ = 0; //!< index of projectile on stack
+ std::vector<SecondaryParticle> secondaries_;
+ EventPtr parent_event_;
+
+ EventType type_ = EventType::Uninitialized;
+
+ std::optional<Code> targetCode_;
+
+ public:
+ Event() = default;
+
+ void setParentEvent(EventPtr const& evt) { parent_event_ = evt; }
+
+ bool hasParentEvent() const { return bool(parent_event_); }
+ EventPtr& parentEvent() { return parent_event_; }
+ EventPtr const& parentEvent() const { return parent_event_; }
+
+ void setProjectileIndex(size_t i) { projectile_index_ = i; }
+ size_t projectileIndex() const { return projectile_index_; }
+
+ size_t addSecondary(HEPEnergyType energy,
+ Vector<hepmomentum_d> const& momentum,
+ Code pid) {
+ secondaries_.emplace_back(energy, momentum, pid);
+ return secondaries_.size() - 1;
+ }
+
+ std::vector<SecondaryParticle> const& secondaries() const { return secondaries_; }
+
+ void setTargetCode(const Code t) { targetCode_ = t; }
+
+ std::string as_string() const {
+ return fmt::format("hasParent={}, projIndex={}, Nsec={}", hasParentEvent(),
+ projectile_index_, secondaries_.size());
+ }
+
+ EventType eventType() const { return type_; }
+
+ void setEventType(EventType t) { type_ = t; }
+ };
+
+} // namespace corsika::history
diff --git a/corsika/stack/history/EventType.hpp b/corsika/stack/history/EventType.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..4781d5c64c085684a8ea139abda2aa856b56cd66
--- /dev/null
+++ b/corsika/stack/history/EventType.hpp
@@ -0,0 +1,13 @@
+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * 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
+
+namespace corsika::history {
+ enum class EventType { Uninitialized, Interaction, Decay };
+}
diff --git a/corsika/stack/history/HistoryObservationPlane.cpp b/corsika/stack/history/HistoryObservationPlane.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..b1425f18be3ff20fe68f8ffbaf2f6ff9244d6e78
--- /dev/null
+++ b/corsika/stack/history/HistoryObservationPlane.cpp
@@ -0,0 +1,84 @@
+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * 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/logging/Logging.h>
+#include <corsika/stack/history/HistoryObservationPlane.hpp>
+
+#include <boost/histogram/ostream.hpp>
+
+#include <iomanip>
+#include <iostream>
+
+using namespace corsika::units::si;
+using namespace corsika::history;
+using namespace corsika;
+
+HistoryObservationPlane::HistoryObservationPlane(setup::Stack const& stack,
+ geometry::Plane const& obsPlane,
+ bool deleteOnHit)
+ : stack_{stack}
+ , plane_{obsPlane}
+ , deleteOnHit_{deleteOnHit} {}
+
+corsika::process::EProcessReturn HistoryObservationPlane::DoContinuous(
+ setup::Stack::ParticleType const& particle, setup::Trajectory const& trajectory) {
+ TimeType const timeOfIntersection =
+ (plane_.GetCenter() - trajectory.GetR0()).dot(plane_.GetNormal()) /
+ trajectory.GetV0().dot(plane_.GetNormal());
+
+ if (timeOfIntersection < TimeType::zero()) { return process::EProcessReturn::eOk; }
+
+ if (plane_.IsAbove(trajectory.GetR0()) == plane_.IsAbove(trajectory.GetPosition(1))) {
+ return process::EProcessReturn::eOk;
+ }
+
+ C8LOG_DEBUG(fmt::format("HistoryObservationPlane: Particle detected: pid={}",
+ particle.GetPID()));
+
+ auto const pid = particle.GetPID();
+ if (particles::IsMuon(pid)) { fillHistoryHistogram(particle); }
+
+ if (deleteOnHit_) {
+ return process::EProcessReturn::eParticleAbsorbed;
+ } else {
+ return process::EProcessReturn::eOk;
+ }
+}
+
+LengthType HistoryObservationPlane::MaxStepLength(setup::Stack::ParticleType const&,
+ setup::Trajectory const& trajectory) {
+ 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;
+}
+
+void HistoryObservationPlane::fillHistoryHistogram(
+ setup::Stack::ParticleType const& muon) {
+ double const muon_energy = muon.GetEnergy() / 1_GeV;
+
+ int genctr{0};
+ Event const* event = muon.GetEvent().get();
+ while (event) {
+ auto const projectile = stack_.cfirst() + event->projectileIndex();
+ if (event->eventType() == EventType::Interaction) {
+ genctr++;
+ double const projEnergy = projectile.GetEnergy() / 1_GeV;
+ int const pdg = static_cast<int>(particles::GetPDG(projectile.GetPID()));
+
+ histogram_(muon_energy, projEnergy, pdg);
+ }
+ event = event->parentEvent().get(); // projectile.GetEvent().get();
+ }
+}
diff --git a/corsika/stack/history/HistoryObservationPlane.hpp b/corsika/stack/history/HistoryObservationPlane.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..fb67ca7f9b8897cb6d6f2b11c6ce02101a5ea653
--- /dev/null
+++ b/corsika/stack/history/HistoryObservationPlane.hpp
@@ -0,0 +1,60 @@
+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * 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/geometry/Plane.h>
+#include <corsika/process/ContinuousProcess.h>
+#include <corsika/setup/SetupStack.h>
+#include <corsika/setup/SetupTrajectory.h>
+#include <corsika/units/PhysicalUnits.h>
+
+#include <boost/histogram.hpp>
+
+#include <functional>
+
+namespace corsika::history {
+ namespace detail {
+ inline auto hist_factory() {
+ namespace bh = boost::histogram;
+ namespace bha = bh::axis;
+ auto h = bh::make_histogram(
+ bha::regular<float, bha::transform::log>{11 * 5, 1e0, 1e11, "muon energy"},
+ bha::regular<float, bha::transform::log>{11 * 5, 1e0, 1e11,
+ "projectile energy"},
+ bha::category<int, bh::use_default, bha::option::growth_t>{
+ {211, -211, 2212, -2212}, "projectile PDG"});
+ return h;
+ }
+ } // namespace detail
+
+ class HistoryObservationPlane
+ : public corsika::process::ContinuousProcess<HistoryObservationPlane> {
+ public:
+ HistoryObservationPlane(setup::Stack const&, geometry::Plane const&, bool = true);
+
+ corsika::units::si::LengthType MaxStepLength(
+ corsika::setup::Stack::ParticleType const&,
+ corsika::setup::Trajectory const& vTrajectory);
+
+ corsika::process::EProcessReturn DoContinuous(
+ corsika::setup::Stack::ParticleType const& vParticle,
+ corsika::setup::Trajectory const& vTrajectory);
+
+ auto const& histogram() const { return histogram_; }
+
+ private:
+ void fillHistoryHistogram(setup::Stack::ParticleType const&);
+
+ setup::Stack const& stack_;
+ geometry::Plane const plane_;
+ bool const deleteOnHit_;
+
+ decltype(detail::hist_factory()) histogram_ = detail::hist_factory();
+ };
+} // namespace corsika::history
diff --git a/corsika/stack/history/HistorySecondaryProducer.hpp b/corsika/stack/history/HistorySecondaryProducer.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..f7317a22b0c290203880af3110b07a63fb4b2379
--- /dev/null
+++ b/corsika/stack/history/HistorySecondaryProducer.hpp
@@ -0,0 +1,59 @@
+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * 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/stack/SecondaryView.hpp>
+#include <corsika/stack/history/Event.hpp>
+
+#include <corsika/framework/logging/Logging.hpp>
+
+#include <boost/type_index.hpp>
+
+#include <memory>
+#include <type_traits>
+#include <utility>
+
+namespace corsika::history {
+
+ //! mix-in class for SecondaryView that fills secondaries into an \class Event
+ template <class T1, template <class> class T2>
+ class HistorySecondaryProducer {
+ public:
+ EventPtr event_;
+ static bool constexpr has_event{true};
+
+ public:
+ template <typename Particle>
+ HistorySecondaryProducer(Particle const& p)
+ : event_{std::make_shared<Event>()} {
+ CORSIKA_LOG_TRACE("HistorySecondaryProducer::HistorySecondaryProducer");
+ event_->setProjectileIndex(p.getIndex());
+ event_->setParentEvent(p.getEvent());
+ }
+
+ /**
+ * Method is called after a new Secondary has been created on the
+ * SecondaryView. Extra logic can be introduced here.
+ *
+ * The input Particle is the new secondary that was produced and
+ * is of course a reference into the SecondaryView itself.
+ */
+ template <typename Particle>
+ auto new_secondary(Particle& sec) {
+ CORSIKA_LOG_TRACE("HistorySecondaryProducer::new_secondary(sec)");
+
+ // store particles at production time in Event here
+ auto const sec_index =
+ event_->addSecondary(sec.getEnergy(), sec.getMomentum(), sec.getPID());
+ sec.setParentEventIndex(sec_index);
+ sec.setEvent(event_);
+ }
+ };
+
+} // namespace corsika::history
diff --git a/corsika/stack/history/HistoryStackExtension.hpp b/corsika/stack/history/HistoryStackExtension.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..2d79af4c4fa12a009a019e1fe27e576e130a184c
--- /dev/null
+++ b/corsika/stack/history/HistoryStackExtension.hpp
@@ -0,0 +1,139 @@
+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * 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/logging/Logging.hpp>
+#include <corsika/framework/stack/Stack.hpp>
+
+#include <memory>
+#include <utility>
+#include <vector>
+
+namespace corsika::history {
+
+ /**
+ * @class HistoryData
+ *
+ * definition of stack-data object to store history information this
+ * is vector with shared_ptr<TEvent>, where TEvent is a free
+ * template parameter for customization.
+ */
+ template <typename TEvent>
+ class HistoryData {
+ using EventPtr =
+ std::shared_ptr<TEvent>; //!< Pointer to the event where this particle was created
+ using ParentEventIndex = int; //!< index to TEvent::secondaries_
+ using DataType = std::pair<EventPtr, ParentEventIndex>;
+
+ public:
+ // these functions are needed for the Stack interface
+ void clear() { historyData_.clear(); }
+ unsigned int getSize() const { return historyData_.size(); }
+ unsigned int getCapacity() const { return historyData_.size(); }
+ void copy(const int i1, const int i2) { historyData_[i2] = historyData_[i1]; }
+ void swap(const int i1, const int i2) {
+ std::swap(historyData_[i1], historyData_[i2]);
+ }
+
+ // custom data access function
+ void setEvent(const int i, EventPtr v) { historyData_[i].first = std::move(v); }
+ const EventPtr& getEvent(const int i) const { return historyData_[i].first; }
+ EventPtr& getEvent(const int i) { return historyData_[i].first; }
+
+ void setParentEventIndex(const int i, ParentEventIndex v) {
+ historyData_[i].second = std::move(v);
+ }
+ ParentEventIndex getParentEventIndex(const int i) const {
+ return historyData_[i].second;
+ }
+
+ // these functions are also needed by the Stack interface
+ void incrementSize() { historyData_.push_back(DataType{}); }
+ void decrementSize() {
+ if (historyData_.size() > 0) { historyData_.pop_back(); }
+ }
+
+ // custom private data section
+ private:
+ std::vector<DataType> historyData_;
+ };
+
+ /**
+ * @class HistoryDataInterface
+ *
+ * corresponding defintion of a stack-readout object, the iteractor
+ * dereference operator will deliver access to these function
+ // defintion of a stack-readout object, the iteractor dereference
+ // operator will deliver access to these function
+ */
+ template <typename T, typename TEvent>
+ class HistoryDataInterface : public T {
+ protected:
+ using T::getStack;
+ using T::getStackData;
+
+ public:
+ using T::getIndex;
+
+ public:
+ // create a new particle from scratch
+ void setParticleData() {
+ CORSIKA_LOG_TRACE("HistoyDatatInterface::setParticleData()");
+ getStackData().setParentEventIndex(getIndex(), -1);
+ }
+
+ // create a new particle as secondary of a parent
+ void setParticleData(HistoryDataInterface& /*parent*/) {
+ CORSIKA_LOG_TRACE("HistoyDatatInterface::setParticleData(parnt)");
+ setParticleData();
+ }
+
+ void setEvent(const std::shared_ptr<TEvent>& v) {
+ getStackData().setEvent(getIndex(), v);
+ }
+
+ void setParentEventIndex(int index) {
+ getStackData().setParentEventIndex(getIndex(), index);
+ }
+
+ std::shared_ptr<TEvent> getEvent() const {
+ return getStackData().getEvent(getIndex());
+ }
+
+ int getParentEventIndex() const {
+ return getStackData().getParentEventIndex(getIndex());
+ }
+
+ std::string asString() const {
+ return fmt::format("i_parent={}, [evt: {}]", getParentEventIndex(),
+ (bool(getEvent()) ? getEvent()->as_string() : "n/a"));
+ }
+ };
+
+ template <typename T, typename TEvent>
+ struct MakeHistoryDataInterface {
+ typedef HistoryDataInterface<T, TEvent> type;
+ };
+
+} // namespace corsika::history
+
+// for user-friendlyness we create the HistoryDataInterface type
+// with the histoy::Event data content right here:
+
+#include <corsika/stack/history/Event.hpp>
+
+namespace corsika::history {
+
+ template <typename TStackIter>
+ using HistoryEventDataInterface =
+ typename history::MakeHistoryDataInterface<TStackIter, history::Event>::type;
+
+ using HistoryEventData = history::HistoryData<history::Event>;
+
+} // namespace corsika::history
diff --git a/corsika/stack/history/SecondaryParticle.hpp b/corsika/stack/history/SecondaryParticle.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..4a36abf37f50eb5ce61f5f7e5cd8f16720ea7893
--- /dev/null
+++ b/corsika/stack/history/SecondaryParticle.hpp
@@ -0,0 +1,35 @@
+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * 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/Vector.hpp>
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+
+#include <vector>
+
+namespace corsika::history {
+
+ /**
+ * This class stores the non-common properties of secondaries in an event. All
+ * other (common) properties are available via the event itself or its projectile.
+ */
+ struct SecondaryParticle {
+ HEPEnergyType const energy_;
+ Vector<hepmomentum_d> const momentum_;
+ Code const pid_;
+
+ public:
+ SecondaryParticle(HEPEnergyType energy, Vector<hepmomentum_d> momentum, Code pid)
+ : energy_{energy}
+ , momentum_{momentum}
+ , pid_{pid} {}
+ };
+
+} // namespace corsika::history
diff --git a/examples/em_shower.cpp b/examples/em_shower.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..594a95845c9a2e8d8640aeda4c7847ed591d9420
--- /dev/null
+++ b/examples/em_shower.cpp
@@ -0,0 +1,185 @@
+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#include <corsika/cascade/Cascade.h>
+#include <corsika/environment/Environment.h>
+#include <corsika/environment/LayeredSphericalAtmosphereBuilder.h>
+#include <corsika/environment/NuclearComposition.h>
+#include <corsika/environment/ShowerAxis.h>
+#include <corsika/geometry/Plane.h>
+#include <corsika/geometry/Sphere.h>
+#include <corsika/process/ProcessSequence.h>
+#include <corsika/process/StackProcess.h>
+#include <corsika/process/longitudinal_profile/LongitudinalProfile.h>
+#include <corsika/process/observation_plane/ObservationPlane.h>
+#include <corsika/process/particle_cut/ParticleCut.h>
+#include <corsika/process/proposal/ContinuousProcess.h>
+#include <corsika/process/proposal/Interaction.h>
+#include <corsika/process/track_writer/TrackWriter.h>
+#include <corsika/process/tracking_line/TrackingLine.h>
+#include <corsika/random/RNGManager.h>
+#include <corsika/setup/SetupStack.h>
+#include <corsika/setup/SetupTrajectory.h>
+#include <corsika/units/PhysicalUnits.h>
+#include <corsika/utl/CorsikaFenv.h>
+#include <corsika/process/interaction_counter/InteractionCounter.hpp>
+
+#include <corsika/logging/Logging.h>
+
+#include <iomanip>
+#include <iostream>
+#include <limits>
+#include <string>
+#include <typeinfo>
+
+using namespace corsika;
+using namespace corsika::process;
+using namespace corsika::units;
+using namespace corsika::particles;
+using namespace corsika::random;
+using namespace corsika::geometry;
+using namespace corsika::environment;
+
+using namespace std;
+using namespace corsika::units::si;
+
+void registerRandomStreams() {
+ random::RNGManager::GetInstance().RegisterRandomStream("cascade");
+ random::RNGManager::GetInstance().RegisterRandomStream("proposal");
+ random::RNGManager::GetInstance().SeedAll();
+}
+
+template <typename T>
+using MyExtraEnv = environment::MediumPropertyModel<environment::UniformMagneticField<T>>;
+
+int main(int argc, char** argv) {
+
+ logging::SetLevel(logging::level::info);
+
+ if (argc != 2) {
+ std::cerr << "usage: em_shower <energy/GeV>" << std::endl;
+ return 1;
+ }
+ feenableexcept(FE_INVALID);
+ // initialize random number sequence(s)
+ registerRandomStreams();
+
+ // setup environment, geometry
+ using EnvType = setup::Environment;
+ EnvType env;
+ const CoordinateSystem& rootCS = env.GetCoordinateSystem();
+ Point const center{rootCS, 0_m, 0_m, 0_m};
+ auto builder = environment::make_layered_spherical_atmosphere_builder<
+ setup::EnvironmentInterface,
+ MyExtraEnv>::create(center, units::constants::EarthRadius::Mean,
+ environment::Medium::AirDry1Atm,
+ geometry::Vector{rootCS, 0_T, 0_T, 1_T});
+ builder.setNuclearComposition(
+ {{particles::Code::Nitrogen, particles::Code::Oxygen},
+ {0.7847f, 1.f - 0.7847f}}); // values taken from AIRES manual, Ar removed for now
+
+ builder.addExponentialLayer(1222.6562_g / (1_cm * 1_cm), 994186.38_cm, 4_km);
+ builder.addExponentialLayer(1144.9069_g / (1_cm * 1_cm), 878153.55_cm, 10_km);
+ builder.addExponentialLayer(1305.5948_g / (1_cm * 1_cm), 636143.04_cm, 40_km);
+ builder.addExponentialLayer(540.1778_g / (1_cm * 1_cm), 772170.16_cm, 100_km);
+ builder.addLinearLayer(1e9_cm, 112.8_km);
+ builder.assemble(env);
+
+ // setup particle stack, and add primary particle
+ setup::Stack stack;
+ stack.Clear();
+ const Code beamCode = Code::Electron;
+ auto const mass = particles::GetMass(beamCode);
+ const HEPEnergyType E0 = 1_GeV * std::stof(std::string(argv[1]));
+ double theta = 0.;
+ auto const thetaRad = theta / 180. * M_PI;
+
+ auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
+ return sqrt((Elab - m) * (Elab + m));
+ };
+ HEPMomentumType P0 = elab2plab(E0, mass);
+ auto momentumComponents = [](double thetaRad, HEPMomentumType ptot) {
+ return std::make_tuple(ptot * sin(thetaRad), 0_eV, -ptot * cos(thetaRad));
+ };
+
+ auto const [px, py, pz] = momentumComponents(thetaRad, P0);
+ auto plab = corsika::stack::MomentumVector(rootCS, {px, py, pz});
+ cout << "input particle: " << beamCode << endl;
+ cout << "input angles: theta=" << theta << endl;
+ cout << "input momentum: " << plab.GetComponents() / 1_GeV << ", norm = " << plab.norm()
+ << endl;
+
+ auto const observationHeight = 1.4_km + builder.getEarthRadius();
+ auto const injectionHeight = 112.75_km + builder.getEarthRadius();
+ auto const t = -observationHeight * cos(thetaRad) +
+ sqrt(-static_pow<2>(sin(thetaRad) * observationHeight) +
+ static_pow<2>(injectionHeight));
+ Point const showerCore{rootCS, 0_m, 0_m, observationHeight};
+ Point const injectionPos =
+ showerCore +
+ Vector<dimensionless_d>{rootCS, {-sin(thetaRad), 0, cos(thetaRad)}} * t;
+
+ std::cout << "point of injection: " << injectionPos.GetCoordinates() << std::endl;
+
+ stack.AddParticle(
+ std::tuple<particles::Code, units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
+ beamCode, E0, plab, injectionPos, 0_ns});
+
+ std::cout << "shower axis length: " << (showerCore - injectionPos).norm() * 1.02
+ << std::endl;
+
+ environment::ShowerAxis const showerAxis{injectionPos,
+ (showerCore - injectionPos) * 1.02, env};
+
+ // setup processes, decays and interactions
+
+ // PROPOSAL processs proposal{...};
+ process::particle_cut::ParticleCut cut(10_GeV, false, true);
+ process::proposal::Interaction proposal(env, cut.GetECut());
+ process::proposal::ContinuousProcess em_continuous(env, cut.GetECut());
+ process::interaction_counter::InteractionCounter proposalCounted(proposal);
+
+ process::track_writer::TrackWriter trackWriter("tracks.dat");
+
+ // long. profile; columns for gamma, e+, e- still need to be added
+ process::longitudinal_profile::LongitudinalProfile longprof{showerAxis};
+
+ Plane const obsPlane(showerCore, Vector<dimensionless_d>(rootCS, {0., 0., 1.}));
+ process::observation_plane::ObservationPlane observationLevel(obsPlane,
+ "particles.dat");
+
+ auto sequence = process::sequence(proposalCounted, em_continuous, longprof, cut,
+ observationLevel, trackWriter);
+ // define air shower object, run simulation
+ tracking_line::TrackingLine tracking;
+ cascade::Cascade EAS(env, tracking, sequence, stack);
+
+ // to fix the point of first interaction, uncomment the following two lines:
+ // EAS.SetNodes();
+ // EAS.forceInteraction();
+
+ EAS.Run();
+
+ cut.ShowResults();
+ em_continuous.showResults();
+ observationLevel.ShowResults();
+ const HEPEnergyType Efinal = cut.GetCutEnergy() + cut.GetInvEnergy() +
+ cut.GetEmEnergy() + em_continuous.energyLost() +
+ observationLevel.GetEnergyGround();
+ cout << "total cut energy (GeV): " << Efinal / 1_GeV << endl
+ << "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl;
+ observationLevel.Reset();
+ cut.Reset();
+ em_continuous.reset();
+
+ auto const hists = proposalCounted.GetHistogram();
+ hists.saveLab("inthist_lab_emShower.npz");
+ hists.saveCMS("inthist_cms_emShower.npz");
+ longprof.save("longprof_emShower.txt");
+}
diff --git a/examples/hybrid_MC.cpp b/examples/hybrid_MC.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..c6647209bb13d4843020c2bce4d0354486f97753
--- /dev/null
+++ b/examples/hybrid_MC.cpp
@@ -0,0 +1,281 @@
+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * 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.
+ */
+
+/* clang-format off */
+// InteractionCounter used boost/histogram, which
+// fails if boost/type_traits have been included before. Thus, we have
+// to include it first...
+#include <corsika/process/interaction_counter/InteractionCounter.hpp>
+/* clang-format on */
+#include <corsika/cascade/Cascade.h>
+#include <corsika/environment/Environment.h>
+#include <corsika/environment/FlatExponential.h>
+#include <corsika/environment/LayeredSphericalAtmosphereBuilder.h>
+#include <corsika/environment/NuclearComposition.h>
+#include <corsika/environment/ShowerAxis.h>
+#include <corsika/geometry/Plane.h>
+#include <corsika/geometry/Sphere.h>
+#include <corsika/logging/Logging.h>
+#include <corsika/process/ProcessSequence.h>
+#include <corsika/process/SwitchProcessSequence.h>
+#include <corsika/process/StackProcess.h>
+#include <corsika/process/conex_source_cut/CONEXSourceCut.h>
+#include <corsika/process/energy_loss/EnergyLoss.h>
+#include <corsika/process/longitudinal_profile/LongitudinalProfile.h>
+#include <corsika/process/observation_plane/ObservationPlane.h>
+#include <corsika/process/on_shell_check/OnShellCheck.h>
+#include <corsika/process/particle_cut/ParticleCut.h>
+#include <corsika/process/pythia/Decay.h>
+#include <corsika/process/sibyll/Decay.h>
+#include <corsika/process/sibyll/Interaction.h>
+#include <corsika/process/sibyll/NuclearInteraction.h>
+#include <corsika/process/tracking_line/TrackingLine.h>
+#include <corsika/process/urqmd/UrQMD.h>
+#include <corsika/random/RNGManager.h>
+#include <corsika/setup/SetupStack.h>
+#include <corsika/setup/SetupTrajectory.h>
+#include <corsika/units/PhysicalUnits.h>
+#include <corsika/utl/CorsikaFenv.h>
+
+#include <iomanip>
+#include <iostream>
+#include <limits>
+#include <string>
+
+using namespace corsika;
+using namespace corsika::process;
+using namespace corsika::units;
+using namespace corsika::particles;
+using namespace corsika::random;
+using namespace corsika::geometry;
+using namespace corsika::environment;
+
+using namespace std;
+using namespace corsika::units::si;
+
+void registerRandomStreams(const int seed) {
+ random::RNGManager::GetInstance().RegisterRandomStream("cascade");
+ random::RNGManager::GetInstance().RegisterRandomStream("qgsjet");
+ random::RNGManager::GetInstance().RegisterRandomStream("sibyll");
+ random::RNGManager::GetInstance().RegisterRandomStream("pythia");
+ random::RNGManager::GetInstance().RegisterRandomStream("urqmd");
+ random::RNGManager::GetInstance().RegisterRandomStream("proposal");
+
+ if (seed == 0)
+ random::RNGManager::GetInstance().SeedAll();
+ else
+ random::RNGManager::GetInstance().SeedAll(seed);
+}
+
+template <typename T>
+using MyExtraEnv = environment::MediumPropertyModel<environment::UniformMagneticField<T>>;
+
+int main(int argc, char** argv) {
+
+ logging::SetLevel(logging::level::info);
+
+ C8LOG_INFO("hybrid_MC");
+
+ if (argc < 4) {
+ std::cerr << "usage: hybrid_MC <A> <Z> <energy/GeV> [seed]" << std::endl;
+ std::cerr << " if no seed is given, a random seed is chosen" << std::endl;
+ return 1;
+ }
+ feenableexcept(FE_INVALID);
+
+ int seed = 0;
+ if (argc > 4) seed = std::stoi(std::string(argv[4]));
+ // initialize random number sequence(s)
+ registerRandomStreams(seed);
+
+ // setup environment, geometry
+ using EnvType = setup::Environment;
+ EnvType env;
+ const CoordinateSystem& rootCS = env.GetCoordinateSystem();
+ Point const center{rootCS, 0_m, 0_m, 0_m};
+ auto builder = environment::make_layered_spherical_atmosphere_builder<
+ setup::EnvironmentInterface,
+ MyExtraEnv>::create(center, units::constants::EarthRadius::Mean,
+ environment::Medium::AirDry1Atm,
+ geometry::Vector{rootCS, 0_T, 0_T, 1_T});
+ builder.setNuclearComposition(
+ {{particles::Code::Nitrogen, particles::Code::Oxygen},
+ {0.7847f, 1.f - 0.7847f}}); // values taken from AIRES manual, Ar removed for now
+
+ builder.addExponentialLayer(1222.6562_g / (1_cm * 1_cm), 994186.38_cm, 4_km);
+ builder.addExponentialLayer(1144.9069_g / (1_cm * 1_cm), 878153.55_cm, 10_km);
+ builder.addExponentialLayer(1305.5948_g / (1_cm * 1_cm), 636143.04_cm, 40_km);
+ builder.addExponentialLayer(540.1778_g / (1_cm * 1_cm), 772170.16_cm, 100_km);
+ builder.addLinearLayer(1e9_cm, 112.8_km);
+ builder.assemble(env);
+
+ // setup particle stack, and add primary particle
+ setup::Stack stack;
+ stack.Clear();
+ const Code beamCode = Code::Nucleus;
+ unsigned short const A = std::stoi(std::string(argv[1]));
+ unsigned short Z = std::stoi(std::string(argv[2]));
+ auto const mass = particles::GetNucleusMass(A, Z);
+ const HEPEnergyType E0 = 1_GeV * std::stof(std::string(argv[3]));
+ double theta = 0.;
+ auto const thetaRad = theta / 180. * M_PI;
+
+ auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
+ return sqrt((Elab - m) * (Elab + m));
+ };
+ HEPMomentumType P0 = elab2plab(E0, mass);
+ auto momentumComponents = [](double thetaRad, HEPMomentumType ptot) {
+ return std::make_tuple(ptot * sin(thetaRad), 0_eV, -ptot * cos(thetaRad));
+ };
+
+ auto const [px, py, pz] = momentumComponents(thetaRad, P0);
+ auto plab = corsika::stack::MomentumVector(rootCS, {px, py, pz});
+ cout << "input particle: " << beamCode << endl;
+ cout << "input angles: theta=" << theta << endl;
+ cout << "input momentum: " << plab.GetComponents() / 1_GeV << ", norm = " << plab.norm()
+ << endl;
+
+ auto const observationHeight = 0_km + builder.getEarthRadius();
+ auto const injectionHeight = 112.75_km + builder.getEarthRadius();
+ auto const t = -observationHeight * cos(thetaRad) +
+ sqrt(-units::static_pow<2>(sin(thetaRad) * observationHeight) +
+ units::static_pow<2>(injectionHeight));
+ Point const showerCore{rootCS, 0_m, 0_m, observationHeight};
+ Point const injectionPos =
+ showerCore +
+ Vector<dimensionless_d>{rootCS, {-sin(thetaRad), 0, cos(thetaRad)}} * t;
+
+ std::cout << "point of injection: " << injectionPos.GetCoordinates() << std::endl;
+
+ if (A != 1) {
+ stack.AddParticle(std::tuple<particles::Code, units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, geometry::Point,
+ units::si::TimeType, unsigned short, unsigned short>{
+ beamCode, E0, plab, injectionPos, 0_ns, A, Z});
+
+ } else {
+ stack.AddParticle(
+ std::tuple<particles::Code, units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
+ particles::Code::Proton, E0, plab, injectionPos, 0_ns});
+ }
+
+ std::cout << "shower axis length: " << (showerCore - injectionPos).norm() * 1.02
+ << std::endl;
+
+ environment::ShowerAxis const showerAxis{injectionPos,
+ (showerCore - injectionPos) * 1.02, env};
+
+ // setup processes, decays and interactions
+
+ process::sibyll::Interaction sibyll;
+ process::interaction_counter::InteractionCounter sibyllCounted(sibyll);
+
+ process::sibyll::NuclearInteraction sibyllNuc(sibyll, env);
+ process::interaction_counter::InteractionCounter sibyllNucCounted(sibyllNuc);
+
+ process::pythia::Decay decayPythia;
+
+ // use sibyll decay routine for decays of particles unknown to pythia
+ process::sibyll::Decay decaySibyll{{
+ Code::N1440Plus,
+ Code::N1440MinusBar,
+ Code::N1440_0,
+ Code::N1440_0Bar,
+ Code::N1710Plus,
+ Code::N1710MinusBar,
+ Code::N1710_0,
+ Code::N1710_0Bar,
+
+ Code::Pi1300Plus,
+ Code::Pi1300Minus,
+ Code::Pi1300_0,
+
+ Code::KStar0_1430_0,
+ Code::KStar0_1430_0Bar,
+ Code::KStar0_1430_Plus,
+ Code::KStar0_1430_MinusBar,
+ }};
+
+ decaySibyll.PrintDecayConfig();
+
+ process::particle_cut::ParticleCut cut{60_GeV, false, true};
+ process::energy_loss::EnergyLoss eLoss{showerAxis, cut.GetECut()};
+
+ corsika::process::conex_source_cut::CONEXSourceCut conex(
+ center, showerAxis, t, injectionHeight, E0,
+ particles::GetPDG(particles::Code::Proton));
+
+ process::on_shell_check::OnShellCheck reset_particle_mass(1.e-3, 1.e-1, false);
+
+ process::longitudinal_profile::LongitudinalProfile longprof{showerAxis};
+
+ Plane const obsPlane(showerCore, Vector<dimensionless_d>(rootCS, {0., 0., 1.}));
+ process::observation_plane::ObservationPlane observationLevel(obsPlane,
+ "particles.dat");
+
+ process::UrQMD::UrQMD urqmd;
+ process::interaction_counter::InteractionCounter urqmdCounted{urqmd};
+
+ // assemble all processes into an ordered process list
+ struct EnergySwitch {
+ HEPEnergyType cutE_;
+ EnergySwitch(HEPEnergyType cutE)
+ : cutE_(cutE) {}
+ process::SwitchResult operator()(const setup::Stack::ParticleType& p) {
+ if (p.GetEnergy() < cutE_)
+ return process::SwitchResult::First;
+ else
+ return process::SwitchResult::Second;
+ }
+ };
+ auto hadronSequence =
+ process::select(urqmdCounted, process::sequence(sibyllNucCounted, sibyllCounted),
+ EnergySwitch(55_GeV));
+ auto decaySequence = process::sequence(decayPythia, decaySibyll);
+ auto sequence = process::sequence(hadronSequence, reset_particle_mass, decaySequence,
+ eLoss, cut, conex, longprof, observationLevel);
+
+ // define air shower object, run simulation
+ tracking_line::TrackingLine tracking;
+ cascade::Cascade EAS(env, tracking, sequence, stack);
+
+ // to fix the point of first interaction, uncomment the following two lines:
+ // EAS.SetNodes();
+ // EAS.forceInteraction();
+
+ EAS.Run();
+
+ cut.ShowResults();
+ eLoss.showResults();
+ observationLevel.ShowResults();
+ const HEPEnergyType Efinal = cut.GetCutEnergy() + cut.GetInvEnergy() +
+ cut.GetEmEnergy() + eLoss.energyLost() +
+ observationLevel.GetEnergyGround();
+ cout << "total cut energy (GeV): " << Efinal / 1_GeV << endl
+ << "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl;
+ observationLevel.Reset();
+ cut.Reset();
+ eLoss.reset();
+
+ conex.SolveCE();
+
+ auto const hists = sibyllCounted.GetHistogram() + sibyllNucCounted.GetHistogram() +
+ urqmdCounted.GetHistogram();
+
+ hists.saveLab("inthist_lab.txt");
+ hists.saveCMS("inthist_cms.txt");
+
+ hists.saveLab("inthist_lab.txt");
+ hists.saveCMS("inthist_cms.txt");
+
+ longprof.save("longprof.txt");
+
+ std::ofstream finish("finished");
+ finish << "run completed without error" << std::endl;
+}
diff --git a/examples/particle_list_example.cpp b/examples/particle_list_example.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..1280398ab6cafd8aabcae8785be75432124f93d3
--- /dev/null
+++ b/examples/particle_list_example.cpp
@@ -0,0 +1,59 @@
+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * 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/framework/core/ParticleProperties.hpp>
+#include <corsika/modules/QGSJetII.hpp>
+#include <corsika/modules/Sibyll.hpp>
+#include <corsika/setup/SetupEnvironment.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+
+#include <iomanip>
+#include <string>
+
+using namespace corsika;
+using namespace std;
+
+//
+// The example main program for a particle list
+//
+int main() {
+
+ logging::set_level(spdlog::level::info);
+
+ logging::info(
+ "------------------------------------------\n"
+ "particles in CORSIKA\n"
+ "------------------------------------------\n"
+ "Name | "
+ "PDG-id | "
+ "SIBYLL-id | "
+ "QGSJETII-id| "
+ "PDG-mass (GeV) | "
+ "SIBYLL-mass (GeV)|\n"
+ "{:-}",
+ "", 104);
+ for (auto p : get_all_particles()) {
+ if (!is_nucleus(p)) {
+ corsika::sibyll::SibyllCode sib_id = corsika::sibyll::convertToSibyll(p);
+ auto const sib_mass = (sib_id != corsika::sibyll::SibyllCode::Unknown
+ ? to_string(corsika::sibyll::getSibyllMass(p) / 1_GeV)
+ : "--");
+ auto const qgs_id = corsika::qgsjetII::convertToQgsjetII(p);
+ logging::info("{:20} | {:10} | {:10} | {:10.5} | {:18.5}", p,
+ static_cast<int>(get_PDG(p)),
+ (sib_id != corsika::sibyll::SibyllCode::Unknown
+ ? to_string(static_cast<int>(sib_id))
+ : "--"),
+ (qgs_id != corsika::qgsjetII::QgsjetIICode::Unknown
+ ? to_string(static_cast<int>(qgs_id))
+ : "--"),
+ get_mass(p) / 1_GeV, sib_mass);
+ }
+ }
+ logging::info("{:-104}", "");
+}
diff --git a/tests/common/SetupEnvironment.hpp b/tests/common/SetupEnvironment.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..69d1f6233cb390cb074347bcdd7a3a1cbbdcb61e
--- /dev/null
+++ b/tests/common/SetupEnvironment.hpp
@@ -0,0 +1,45 @@
+#include <corsika/media/HomogeneousMedium.hpp>
+#include <corsika/media/InhomogeneousMedium.hpp>
+#include <corsika/media/MediumPropertyModel.hpp>
+#include <corsika/media/UniformMagneticField.hpp>
+
+#include <tuple>
+#include <memory>
+
+/**
+ * \function setup_environment
+ *
+ * standard environment for unit testing.
+ *
+ */
+namespace corsika::setup::testing {
+
+ inline std::tuple<std::unique_ptr<setup::Environment>, CoordinateSystem const*,
+ setup::Environment::BaseNodeType const*>
+ setup_environment(Code vTargetCode) {
+ auto env = std::make_unique<setup::Environment>();
+ auto& universe = *(env->getUniverse());
+ CoordinateSystemPtr const& cs = env->getCoordinateSystem();
+
+ /**
+ * our world is a sphere at 0,0,0 with R=infty
+ */
+ auto world = setup::Environment::createNode<Sphere>(
+ Point{cs, 0_m, 0_m, 0_m}, 1_km * std::numeric_limits<double>::infinity());
+
+ /**
+ * construct suited environment medium model:
+ */
+ using MyHomogeneousModel = MediumPropertyModel<
+ UniformMagneticField<HomogeneousMedium<setup::EnvironmentInterface>>>;
+
+ world->setModelProperties<MyHomogeneousModel>(
+ Medium::AirDry1Atm, Vector(cs, 0_T, 0_T, 1_T), 1_kg / (1_m * 1_m * 1_m),
+ NuclearComposition(std::vector<Code>{vTargetCode}, std::vector<float>{1.}));
+
+ setup::Environment::BaseNodeType const* nodePtr = world.get();
+ universe.addChild(std::move(world));
+
+ return std::make_tuple(std::move(env), &cs, nodePtr);
+ }
+} // namespace corsika::setup::testing
diff --git a/tests/common/SetupTestEnvironment.hpp b/tests/common/SetupTestEnvironment.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..50d5561dc38021bfeed55f3f2ae0a69613069785
--- /dev/null
+++ b/tests/common/SetupTestEnvironment.hpp
@@ -0,0 +1,45 @@
+#pragma once
+
+#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/framework/geometry/CoordinateSystem.hpp>
+
+#include <corsika/setup/SetupEnvironment.hpp>
+#include <corsika/media/UniformMagneticField.hpp>
+#include <corsika/media/MediumPropertyModel.hpp>
+#include <corsika/media/HomogeneousMedium.hpp>
+
+#include <limits>
+
+namespace corsika::setup::testing {
+
+ inline std::tuple<std::unique_ptr<setup::Environment>, CoordinateSystemPtr const*,
+ setup::Environment::BaseNodeType const*>
+ setup_environment(Code vTargetCode) {
+
+ auto env = std::make_unique<setup::Environment>();
+ auto& universe = *(env->getUniverse());
+ CoordinateSystemPtr const& cs = env->getCoordinateSystem();
+
+ /**
+ * our world is a sphere at 0,0,0 with R=infty
+ */
+ auto world = setup::Environment::createNode<Sphere>(
+ Point{cs, 0_m, 0_m, 0_m}, 1_km * std::numeric_limits<double>::infinity());
+
+ /**
+ * construct suited environment medium model:
+ */
+ using MyHomogeneousModel = MediumPropertyModel<
+ UniformMagneticField<HomogeneousMedium<setup::EnvironmentInterface>>>;
+
+ world->setModelProperties<MyHomogeneousModel>(
+ Medium::AirDry1Atm, Vector(cs, 0_T, 0_T, 1_T), 1_kg / (1_m * 1_m * 1_m),
+ NuclearComposition(std::vector<Code>{vTargetCode}, std::vector<float>{1.}));
+
+ setup::Environment::BaseNodeType const* nodePtr = world.get();
+ universe.addChild(std::move(world));
+
+ return std::make_tuple(std::move(env), &cs, nodePtr);
+ }
+
+} // namespace corsika::setup::testing
diff --git a/tests/common/SetupTestStack.hpp b/tests/common/SetupTestStack.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..35103012b2865aeab9f8af8a4a2bbc4ec43881a9
--- /dev/null
+++ b/tests/common/SetupTestStack.hpp
@@ -0,0 +1,56 @@
+#pragma once
+
+#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/framework/geometry/RootCoordinateSystem.hpp>
+#include <corsika/framework/geometry/Vector.hpp>
+
+/**
+ * \file SetupTestStack
+ *
+ * standard stack setup for unit tests.
+ **/
+
+namespace corsika::setup::testing {
+
+ /**
+ * \function setup_stack
+ *
+ * standard stack setup for unit tests.
+ *
+ *
+ *
+ *
+ * \return a tuple with element 0 being a Stack object filled with
+ * one particle, and element 1 the StackView on it.
+ **/
+
+ inline std::tuple<std::unique_ptr<setup::Stack>, std::unique_ptr<setup::StackView>>
+ setup_stack(Code vProjectileType, int vA, int vZ, HEPEnergyType vMomentum,
+ setup::Environment::BaseNodeType* const vNodePtr,
+ CoordinateSystemPtr const& cs) {
+
+ auto stack = std::make_unique<setup::Stack>();
+
+ Point const origin(cs, {0_m, 0_m, 0_m});
+ MomentumVector const pLab(cs, {vMomentum, 0_GeV, 0_GeV});
+
+ if (vProjectileType == Code::Nucleus) {
+ auto constexpr mN = constants::nucleonMass;
+ HEPEnergyType const E0 = sqrt(static_pow<2>(mN * vA) + pLab.getSquaredNorm());
+ auto particle = stack->addParticle(
+ std::make_tuple(Code::Nucleus, E0, pLab, origin, 0_ns, vA, vZ));
+ particle.setNode(vNodePtr);
+ return std::make_tuple(std::move(stack),
+ std::make_unique<setup::StackView>(particle));
+ } else { // not a nucleus
+ HEPEnergyType const E0 =
+ sqrt(static_pow<2>(get_mass(vProjectileType)) + pLab.getSquaredNorm());
+ auto particle =
+ stack->addParticle(std::make_tuple(vProjectileType, E0, pLab, origin, 0_ns));
+ particle.setNode(vNodePtr);
+ return std::make_tuple(std::move(stack),
+ std::make_unique<setup::StackView>(particle));
+ }
+ }
+
+} // namespace corsika::setup::testing
diff --git a/tests/framework/CMakeCache.txt b/tests/framework/CMakeCache.txt
new file mode 100644
index 0000000000000000000000000000000000000000..792df53c17d42a1b2a9f008c4c4450be089dbfec
--- /dev/null
+++ b/tests/framework/CMakeCache.txt
@@ -0,0 +1,342 @@
+# This is the CMakeCache file.
+# For build in directory: /home/rulrich/work/corsika/corsika/tests/framework
+# It was generated by CMake: /usr/bin/cmake
+# You can edit this file to change values found and used by cmake.
+# If you do not want to change any of the values, simply exit the editor.
+# If you do want to change a value, simply edit, save, and exit the editor.
+# The syntax for the file is as follows:
+# KEY:TYPE=VALUE
+# KEY is the name of a variable in the cache.
+# TYPE is a hint to GUIs for the type of VALUE, DO NOT EDIT TYPE!.
+# VALUE is the current value for the KEY.
+
+########################
+# EXTERNAL cache entries
+########################
+
+//Path to a program.
+CMAKE_AR:FILEPATH=/usr/bin/ar
+
+//For backwards compatibility, what version of CMake commands and
+// syntax should this version of CMake try to support.
+CMAKE_BACKWARDS_COMPATIBILITY:STRING=2.4
+
+//Choose the type of build, options are: None(CMAKE_CXX_FLAGS or
+// CMAKE_C_FLAGS used) Debug Release RelWithDebInfo MinSizeRel.
+CMAKE_BUILD_TYPE:STRING=
+
+//Enable/Disable color output during build.
+CMAKE_COLOR_MAKEFILE:BOOL=ON
+
+//CXX compiler
+CMAKE_CXX_COMPILER:FILEPATH=/usr/bin/c++
+
+//A wrapper around 'ar' adding the appropriate '--plugin' option
+// for the GCC compiler
+CMAKE_CXX_COMPILER_AR:FILEPATH=/usr/bin/gcc-ar-7
+
+//A wrapper around 'ranlib' adding the appropriate '--plugin' option
+// for the GCC compiler
+CMAKE_CXX_COMPILER_RANLIB:FILEPATH=/usr/bin/gcc-ranlib-7
+
+//Flags used by the compiler during all build types.
+CMAKE_CXX_FLAGS:STRING=
+
+//Flags used by the compiler during debug builds.
+CMAKE_CXX_FLAGS_DEBUG:STRING=-g
+
+//Flags used by the compiler during release builds for minimum
+// size.
+CMAKE_CXX_FLAGS_MINSIZEREL:STRING=-Os -DNDEBUG
+
+//Flags used by the compiler during release builds.
+CMAKE_CXX_FLAGS_RELEASE:STRING=-O3 -DNDEBUG
+
+//Flags used by the compiler during release builds with debug info.
+CMAKE_CXX_FLAGS_RELWITHDEBINFO:STRING=-O2 -g -DNDEBUG
+
+//C compiler
+CMAKE_C_COMPILER:FILEPATH=/usr/bin/cc
+
+//A wrapper around 'ar' adding the appropriate '--plugin' option
+// for the GCC compiler
+CMAKE_C_COMPILER_AR:FILEPATH=/usr/bin/gcc-ar-7
+
+//A wrapper around 'ranlib' adding the appropriate '--plugin' option
+// for the GCC compiler
+CMAKE_C_COMPILER_RANLIB:FILEPATH=/usr/bin/gcc-ranlib-7
+
+//Flags used by the compiler during all build types.
+CMAKE_C_FLAGS:STRING=
+
+//Flags used by the compiler during debug builds.
+CMAKE_C_FLAGS_DEBUG:STRING=-g
+
+//Flags used by the compiler during release builds for minimum
+// size.
+CMAKE_C_FLAGS_MINSIZEREL:STRING=-Os -DNDEBUG
+
+//Flags used by the compiler during release builds.
+CMAKE_C_FLAGS_RELEASE:STRING=-O3 -DNDEBUG
+
+//Flags used by the compiler during release builds with debug info.
+CMAKE_C_FLAGS_RELWITHDEBINFO:STRING=-O2 -g -DNDEBUG
+
+//Flags used by the linker.
+CMAKE_EXE_LINKER_FLAGS:STRING=
+
+//Flags used by the linker during debug builds.
+CMAKE_EXE_LINKER_FLAGS_DEBUG:STRING=
+
+//Flags used by the linker during release minsize builds.
+CMAKE_EXE_LINKER_FLAGS_MINSIZEREL:STRING=
+
+//Flags used by the linker during release builds.
+CMAKE_EXE_LINKER_FLAGS_RELEASE:STRING=
+
+//Flags used by the linker during Release with Debug Info builds.
+CMAKE_EXE_LINKER_FLAGS_RELWITHDEBINFO:STRING=
+
+//Enable/Disable output of compile commands during generation.
+CMAKE_EXPORT_COMPILE_COMMANDS:BOOL=OFF
+
+//Install path prefix, prepended onto install directories.
+CMAKE_INSTALL_PREFIX:PATH=/usr/local
+
+//Path to a program.
+CMAKE_LINKER:FILEPATH=/usr/bin/ld
+
+//Path to a program.
+CMAKE_MAKE_PROGRAM:FILEPATH=/usr/bin/make
+
+//Flags used by the linker during the creation of modules.
+CMAKE_MODULE_LINKER_FLAGS:STRING=
+
+//Flags used by the linker during debug builds.
+CMAKE_MODULE_LINKER_FLAGS_DEBUG:STRING=
+
+//Flags used by the linker during release minsize builds.
+CMAKE_MODULE_LINKER_FLAGS_MINSIZEREL:STRING=
+
+//Flags used by the linker during release builds.
+CMAKE_MODULE_LINKER_FLAGS_RELEASE:STRING=
+
+//Flags used by the linker during Release with Debug Info builds.
+CMAKE_MODULE_LINKER_FLAGS_RELWITHDEBINFO:STRING=
+
+//Path to a program.
+CMAKE_NM:FILEPATH=/usr/bin/nm
+
+//Path to a program.
+CMAKE_OBJCOPY:FILEPATH=/usr/bin/objcopy
+
+//Path to a program.
+CMAKE_OBJDUMP:FILEPATH=/usr/bin/objdump
+
+//Value Computed by CMake
+CMAKE_PROJECT_NAME:STATIC=Project
+
+//Path to a program.
+CMAKE_RANLIB:FILEPATH=/usr/bin/ranlib
+
+//Flags used by the linker during the creation of dll's.
+CMAKE_SHARED_LINKER_FLAGS:STRING=
+
+//Flags used by the linker during debug builds.
+CMAKE_SHARED_LINKER_FLAGS_DEBUG:STRING=
+
+//Flags used by the linker during release minsize builds.
+CMAKE_SHARED_LINKER_FLAGS_MINSIZEREL:STRING=
+
+//Flags used by the linker during release builds.
+CMAKE_SHARED_LINKER_FLAGS_RELEASE:STRING=
+
+//Flags used by the linker during Release with Debug Info builds.
+CMAKE_SHARED_LINKER_FLAGS_RELWITHDEBINFO:STRING=
+
+//If set, runtime paths are not added when installing shared libraries,
+// but are added when building.
+CMAKE_SKIP_INSTALL_RPATH:BOOL=NO
+
+//If set, runtime paths are not added when using shared libraries.
+CMAKE_SKIP_RPATH:BOOL=NO
+
+//Flags used by the linker during the creation of static libraries.
+CMAKE_STATIC_LINKER_FLAGS:STRING=
+
+//Flags used by the linker during debug builds.
+CMAKE_STATIC_LINKER_FLAGS_DEBUG:STRING=
+
+//Flags used by the linker during release minsize builds.
+CMAKE_STATIC_LINKER_FLAGS_MINSIZEREL:STRING=
+
+//Flags used by the linker during release builds.
+CMAKE_STATIC_LINKER_FLAGS_RELEASE:STRING=
+
+//Flags used by the linker during Release with Debug Info builds.
+CMAKE_STATIC_LINKER_FLAGS_RELWITHDEBINFO:STRING=
+
+//Path to a program.
+CMAKE_STRIP:FILEPATH=/usr/bin/strip
+
+//If this value is on, makefiles will be generated without the
+// .SILENT directive, and all commands will be echoed to the console
+// during the make. This is useful for debugging only. With Visual
+// Studio IDE projects all commands are done without /nologo.
+CMAKE_VERBOSE_MAKEFILE:BOOL=FALSE
+
+//Single output directory for building all executables.
+EXECUTABLE_OUTPUT_PATH:PATH=
+
+//Single output directory for building all libraries.
+LIBRARY_OUTPUT_PATH:PATH=
+
+//Value Computed by CMake
+Project_BINARY_DIR:STATIC=/home/rulrich/work/corsika/corsika/tests/framework
+
+//Value Computed by CMake
+Project_SOURCE_DIR:STATIC=/home/rulrich/work/corsika/corsika/tests
+
+
+########################
+# INTERNAL cache entries
+########################
+
+//ADVANCED property for variable: CMAKE_AR
+CMAKE_AR-ADVANCED:INTERNAL=1
+//This is the directory where this CMakeCache.txt was created
+CMAKE_CACHEFILE_DIR:INTERNAL=/home/rulrich/work/corsika/corsika/tests/framework
+//Major version of cmake used to create the current loaded cache
+CMAKE_CACHE_MAJOR_VERSION:INTERNAL=3
+//Minor version of cmake used to create the current loaded cache
+CMAKE_CACHE_MINOR_VERSION:INTERNAL=10
+//Patch version of cmake used to create the current loaded cache
+CMAKE_CACHE_PATCH_VERSION:INTERNAL=2
+//ADVANCED property for variable: CMAKE_COLOR_MAKEFILE
+CMAKE_COLOR_MAKEFILE-ADVANCED:INTERNAL=1
+//Path to CMake executable.
+CMAKE_COMMAND:INTERNAL=/usr/bin/cmake
+//Path to cpack program executable.
+CMAKE_CPACK_COMMAND:INTERNAL=/usr/bin/cpack
+//Path to ctest program executable.
+CMAKE_CTEST_COMMAND:INTERNAL=/usr/bin/ctest
+//ADVANCED property for variable: CMAKE_CXX_COMPILER
+CMAKE_CXX_COMPILER-ADVANCED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_CXX_COMPILER_AR
+CMAKE_CXX_COMPILER_AR-ADVANCED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_CXX_COMPILER_RANLIB
+CMAKE_CXX_COMPILER_RANLIB-ADVANCED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_CXX_FLAGS
+CMAKE_CXX_FLAGS-ADVANCED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_CXX_FLAGS_DEBUG
+CMAKE_CXX_FLAGS_DEBUG-ADVANCED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_CXX_FLAGS_MINSIZEREL
+CMAKE_CXX_FLAGS_MINSIZEREL-ADVANCED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_CXX_FLAGS_RELEASE
+CMAKE_CXX_FLAGS_RELEASE-ADVANCED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_CXX_FLAGS_RELWITHDEBINFO
+CMAKE_CXX_FLAGS_RELWITHDEBINFO-ADVANCED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_C_COMPILER
+CMAKE_C_COMPILER-ADVANCED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_C_COMPILER_AR
+CMAKE_C_COMPILER_AR-ADVANCED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_C_COMPILER_RANLIB
+CMAKE_C_COMPILER_RANLIB-ADVANCED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_C_FLAGS
+CMAKE_C_FLAGS-ADVANCED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_C_FLAGS_DEBUG
+CMAKE_C_FLAGS_DEBUG-ADVANCED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_C_FLAGS_MINSIZEREL
+CMAKE_C_FLAGS_MINSIZEREL-ADVANCED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_C_FLAGS_RELEASE
+CMAKE_C_FLAGS_RELEASE-ADVANCED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_C_FLAGS_RELWITHDEBINFO
+CMAKE_C_FLAGS_RELWITHDEBINFO-ADVANCED:INTERNAL=1
+//Path to cache edit program executable.
+CMAKE_EDIT_COMMAND:INTERNAL=/usr/bin/ccmake
+//Executable file format
+CMAKE_EXECUTABLE_FORMAT:INTERNAL=ELF
+//ADVANCED property for variable: CMAKE_EXE_LINKER_FLAGS
+CMAKE_EXE_LINKER_FLAGS-ADVANCED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_EXE_LINKER_FLAGS_DEBUG
+CMAKE_EXE_LINKER_FLAGS_DEBUG-ADVANCED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_EXE_LINKER_FLAGS_MINSIZEREL
+CMAKE_EXE_LINKER_FLAGS_MINSIZEREL-ADVANCED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_EXE_LINKER_FLAGS_RELEASE
+CMAKE_EXE_LINKER_FLAGS_RELEASE-ADVANCED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_EXE_LINKER_FLAGS_RELWITHDEBINFO
+CMAKE_EXE_LINKER_FLAGS_RELWITHDEBINFO-ADVANCED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_EXPORT_COMPILE_COMMANDS
+CMAKE_EXPORT_COMPILE_COMMANDS-ADVANCED:INTERNAL=1
+//Name of external makefile project generator.
+CMAKE_EXTRA_GENERATOR:INTERNAL=
+//Name of generator.
+CMAKE_GENERATOR:INTERNAL=Unix Makefiles
+//Name of generator platform.
+CMAKE_GENERATOR_PLATFORM:INTERNAL=
+//Name of generator toolset.
+CMAKE_GENERATOR_TOOLSET:INTERNAL=
+//Source directory with the top level CMakeLists.txt file for this
+// project
+CMAKE_HOME_DIRECTORY:INTERNAL=/home/rulrich/work/corsika/corsika/tests
+//Install .so files without execute permission.
+CMAKE_INSTALL_SO_NO_EXE:INTERNAL=1
+//ADVANCED property for variable: CMAKE_LINKER
+CMAKE_LINKER-ADVANCED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_MAKE_PROGRAM
+CMAKE_MAKE_PROGRAM-ADVANCED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_MODULE_LINKER_FLAGS
+CMAKE_MODULE_LINKER_FLAGS-ADVANCED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_MODULE_LINKER_FLAGS_DEBUG
+CMAKE_MODULE_LINKER_FLAGS_DEBUG-ADVANCED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_MODULE_LINKER_FLAGS_MINSIZEREL
+CMAKE_MODULE_LINKER_FLAGS_MINSIZEREL-ADVANCED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_MODULE_LINKER_FLAGS_RELEASE
+CMAKE_MODULE_LINKER_FLAGS_RELEASE-ADVANCED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_MODULE_LINKER_FLAGS_RELWITHDEBINFO
+CMAKE_MODULE_LINKER_FLAGS_RELWITHDEBINFO-ADVANCED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_NM
+CMAKE_NM-ADVANCED:INTERNAL=1
+//number of local generators
+CMAKE_NUMBER_OF_MAKEFILES:INTERNAL=3
+//ADVANCED property for variable: CMAKE_OBJCOPY
+CMAKE_OBJCOPY-ADVANCED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_OBJDUMP
+CMAKE_OBJDUMP-ADVANCED:INTERNAL=1
+//Platform information initialized
+CMAKE_PLATFORM_INFO_INITIALIZED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_RANLIB
+CMAKE_RANLIB-ADVANCED:INTERNAL=1
+//Path to CMake installation.
+CMAKE_ROOT:INTERNAL=/usr/share/cmake-3.10
+//ADVANCED property for variable: CMAKE_SHARED_LINKER_FLAGS
+CMAKE_SHARED_LINKER_FLAGS-ADVANCED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_SHARED_LINKER_FLAGS_DEBUG
+CMAKE_SHARED_LINKER_FLAGS_DEBUG-ADVANCED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_SHARED_LINKER_FLAGS_MINSIZEREL
+CMAKE_SHARED_LINKER_FLAGS_MINSIZEREL-ADVANCED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_SHARED_LINKER_FLAGS_RELEASE
+CMAKE_SHARED_LINKER_FLAGS_RELEASE-ADVANCED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_SHARED_LINKER_FLAGS_RELWITHDEBINFO
+CMAKE_SHARED_LINKER_FLAGS_RELWITHDEBINFO-ADVANCED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_SKIP_INSTALL_RPATH
+CMAKE_SKIP_INSTALL_RPATH-ADVANCED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_SKIP_RPATH
+CMAKE_SKIP_RPATH-ADVANCED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_STATIC_LINKER_FLAGS
+CMAKE_STATIC_LINKER_FLAGS-ADVANCED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_STATIC_LINKER_FLAGS_DEBUG
+CMAKE_STATIC_LINKER_FLAGS_DEBUG-ADVANCED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_STATIC_LINKER_FLAGS_MINSIZEREL
+CMAKE_STATIC_LINKER_FLAGS_MINSIZEREL-ADVANCED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_STATIC_LINKER_FLAGS_RELEASE
+CMAKE_STATIC_LINKER_FLAGS_RELEASE-ADVANCED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_STATIC_LINKER_FLAGS_RELWITHDEBINFO
+CMAKE_STATIC_LINKER_FLAGS_RELWITHDEBINFO-ADVANCED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_STRIP
+CMAKE_STRIP-ADVANCED:INTERNAL=1
+//uname command
+CMAKE_UNAME:INTERNAL=/bin/uname
+//ADVANCED property for variable: CMAKE_VERBOSE_MAKEFILE
+CMAKE_VERBOSE_MAKEFILE-ADVANCED:INTERNAL=1
+
diff --git a/tests/framework/testSwitchProcessSequence.cpp b/tests/framework/testSwitchProcessSequence.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..a778f2efb5dd587b3fc2341dbb7cb5e205c671ab
--- /dev/null
+++ b/tests/framework/testSwitchProcessSequence.cpp
@@ -0,0 +1,247 @@
+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * 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 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>> {
+
+ 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") {
+ CHECK(switchProcess.GetInteractionLength(p) / kgMSq == Approx(1));
+ CHECK(seq.GetInteractionLength(p) / kgMSq == Approx(3. / 4));
+ }
+ }
+
+ SECTION("high energy") {
+ stack.AddParticle(std::tuple<HEPEnergyType>{4_TeV});
+ auto p = stack.GetNextParticle();
+
+ // high energy process returns 2 kg/m²
+ SECTION("interaction length") {
+ CHECK(switchProcess.GetInteractionLength(p) / kgMSq == Approx(2));
+ CHECK(seq.GetInteractionLength(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);
+
+ CHECK(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") {
+ CHECK(switchProcess.GetInteractionLength(p) / kgMSq == Approx(2. / 3));
+ CHECK(completeSeq.GetInteractionLength(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();
+ CHECK(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") {
+ CHECK(switchProcess.GetInteractionLength(p) / kgMSq == Approx(3));
+ CHECK(completeSeq.GetInteractionLength(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();
+ CHECK(mean == Approx(24. / 7.).margin(0.01));
+ }
+ }
+}
diff --git a/tests/media/testMagneticField.cpp b/tests/media/testMagneticField.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..da0a1a5d80e8702d0533c4d132db95c16a984d4a
--- /dev/null
+++ b/tests/media/testMagneticField.cpp
@@ -0,0 +1,87 @@
+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * 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/framework/core/ParticleProperties.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/geometry/RootCoordinateSystem.hpp>
+#include <corsika/media/HomogeneousMedium.hpp>
+#include <corsika/media/IMediumModel.hpp>
+#include <corsika/media/UniformMagneticField.hpp>
+#include <corsika/media/IMagneticFieldModel.hpp>
+#include <corsika/media/VolumeTreeNode.hpp>
+
+#include <catch2/catch.hpp>
+
+using namespace corsika;
+
+
+TEST_CASE("UniformMagneticField w/ Homogeneous Medium") {
+
+ CoordinateSystemPtr const& gCS = get_root_CoordinateSystem();
+ Point const gOrigin(gCS, {0_m, 0_m, 0_m});
+
+
+ // setup our interface types
+ using IModelInterface = IMagneticFieldModel<IMediumModel>;
+ using AtmModel = UniformMagneticField<HomogeneousMedium<IModelInterface>>;
+
+ // the composition we use for the homogenous medium
+ NuclearComposition const protonComposition(std::vector<Code>{Code::Proton},
+ std::vector<float>{1.f});
+
+ // create a magnetic field vector
+ Vector B0(gCS, 0_T, 0_T, 0_T);
+
+ // the constant density
+ const auto density{19.2_g / cube(1_cm)};
+
+ // create our atmospheric model
+ AtmModel medium(B0, density, protonComposition);
+
+ // and test at several locations
+ CHECK(B0.getComponents(gCS) ==
+ medium.getMagneticField(Point(gCS, -10_m, 4_m, 35_km)).getComponents(gCS));
+ CHECK(
+ B0.getComponents(gCS) ==
+ medium.getMagneticField(Point(gCS, 1000_km, -1000_km, 1000_km)).getComponents(gCS));
+ CHECK(B0.getComponents(gCS) ==
+ medium.getMagneticField(Point(gCS, 0_m, 0_m, 0_m)).getComponents(gCS));
+
+ // create a new magnetic field vector
+ Vector B1(gCS, 23_T, 57_T, -4_T);
+
+ // and update this atmospheric model
+ medium.setMagneticField(B1);
+
+ // and test at several locations
+ CHECK(B1.getComponents(gCS) ==
+ medium.getMagneticField(Point(gCS, -10_m, 4_m, 35_km)).getComponents(gCS));
+ CHECK(
+ B1.getComponents(gCS) ==
+ medium.getMagneticField(Point(gCS, 1000_km, -1000_km, 1000_km)).getComponents(gCS));
+ CHECK(B1.getComponents(gCS) ==
+ medium.getMagneticField(Point(gCS, 0_m, 0_m, 0_m)).getComponents(gCS));
+
+ // check the density and nuclear composition
+ CHECK(density == medium.getMassDensity(Point(gCS, 0_m, 0_m, 0_m)));
+ medium.getNuclearComposition();
+
+ // create a line of length 1 m
+ Line const line(gOrigin, Vector<SpeedType::dimension_type>(
+ gCS, {1_m / second, 0_m / second, 0_m / second}));
+
+ // the end time of our line
+ auto const tEnd = 1_s;
+
+ // and the associated trajectory
+ Trajectory<Line> const trajectory(line, tEnd);
+
+ // and check the integrated grammage
+ CHECK((medium.getIntegratedGrammage(trajectory, 3_m) / (density * 3_m)) == Approx(1));
+ CHECK((medium.getArclengthFromGrammage(trajectory, density * 5_m) / 5_m) == Approx(1));
+}
diff --git a/tests/modules/testCONEXSourceCut.cpp b/tests/modules/testCONEXSourceCut.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..378cd9259488e8ded32678929e987ba507387e07
--- /dev/null
+++ b/tests/modules/testCONEXSourceCut.cpp
@@ -0,0 +1,150 @@
+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * 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/setup/SetupEnvironment.h>
+
+#include <corsika/environment/Environment.h>
+#include <corsika/environment/LayeredSphericalAtmosphereBuilder.h>
+#include <corsika/environment/MediumPropertyModel.h>
+#include <corsika/environment/UniformMagneticField.h>
+
+#include <corsika/geometry/Point.h>
+#include <corsika/geometry/RootCoordinateSystem.h>
+#include <corsika/geometry/Vector.h>
+
+#include <corsika/particles/ParticleProperties.h>
+
+#include <corsika/process/conex_source_cut/CONEXSourceCut.h>
+#include <corsika/process/sibyll/Interaction.h>
+#include <corsika/process/sibyll/NuclearInteraction.h>
+
+#include <corsika/random/RNGManager.h>
+
+#include <corsika/units/PhysicalUnits.h>
+#include <corsika/utl/CorsikaFenv.h>
+
+#include <catch2/catch.hpp>
+
+using namespace corsika;
+using namespace corsika::environment;
+using namespace corsika::geometry;
+using namespace corsika::units::si;
+
+const std::string refDataDir = std::string(REFDATADIR); // from cmake
+
+template <typename T>
+using MExtraEnvirnoment =
+ environment::MediumPropertyModel<environment::UniformMagneticField<T>>;
+
+TEST_CASE("CONEXSourceCut") {
+ random::RNGManager::GetInstance().RegisterRandomStream("cascade");
+ random::RNGManager::GetInstance().RegisterRandomStream("sibyll");
+
+ feenableexcept(FE_INVALID);
+
+ // setup environment, geometry
+ setup::Environment env;
+ const CoordinateSystem& rootCS = env.GetCoordinateSystem();
+ Point const center{rootCS, 0_m, 0_m, 0_m};
+
+ auto builder = environment::make_layered_spherical_atmosphere_builder<
+ setup::EnvironmentInterface,
+ MExtraEnvirnoment>::create(center, conex::earthRadius,
+ environment::Medium::AirDry1Atm,
+ geometry::Vector{rootCS, 0_T, 50_mT, 0_T});
+
+ builder.setNuclearComposition(
+ {{particles::Code::Nitrogen, particles::Code::Oxygen},
+ {0.7847f, 1.f - 0.7847f}}); // values taken from AIRES manual, Ar removed for now
+
+ builder.addExponentialLayer(1222.6562_g / (1_cm * 1_cm), 994186.38_cm, 4_km);
+ builder.addExponentialLayer(1144.9069_g / (1_cm * 1_cm), 878153.55_cm, 10_km);
+ builder.addExponentialLayer(1305.5948_g / (1_cm * 1_cm), 636143.04_cm, 40_km);
+ builder.addExponentialLayer(540.1778_g / (1_cm * 1_cm), 772170.16_cm, 100_km);
+ builder.addLinearLayer(1e9_cm, 112.8_km);
+
+ builder.assemble(env);
+
+ const HEPEnergyType E0 = 1_PeV;
+ double thetaDeg = 60.;
+ auto const thetaRad = thetaDeg / 180. * M_PI;
+
+ auto const observationHeight = 1.4_km + conex::earthRadius;
+ auto const injectionHeight = 112.75_km + conex::earthRadius;
+ auto const t = -observationHeight * cos(thetaRad) +
+ sqrt(-units::static_pow<2>(sin(thetaRad) * observationHeight) +
+ units::static_pow<2>(injectionHeight));
+ Point const showerCore{rootCS, 0_m, 0_m, observationHeight};
+ Point const injectionPos =
+ showerCore +
+ Vector<dimensionless_d>{rootCS, {-sin(thetaRad), 0, cos(thetaRad)}} * t;
+
+ environment::ShowerAxis const showerAxis{injectionPos,
+ (showerCore - injectionPos) * 1.02, env};
+
+ // need to initialize Sibyll, done in constructor:
+ process::sibyll::Interaction sibyll;
+ [[maybe_unused]] process::sibyll::NuclearInteraction sibyllNuc(sibyll, env);
+
+ corsika::process::conex_source_cut::CONEXSourceCut conex(
+ center, showerAxis, t, injectionHeight, E0,
+ particles::GetPDG(particles::Code::Proton));
+
+ HEPEnergyType const Eem{1_PeV};
+ auto const momentum = showerAxis.GetDirection() * Eem;
+
+ auto const emPosition = showerCore + showerAxis.GetDirection() * (-20_km);
+
+ std::cout << "position injection: "
+ << injectionPos.GetCoordinates(conex.GetObserverCS()) << " "
+ << injectionPos.GetCoordinates(rootCS) << std::endl;
+ std::cout << "position core: " << showerCore.GetCoordinates(conex.GetObserverCS())
+ << " " << showerCore.GetCoordinates(rootCS) << std::endl;
+ std::cout << "position EM: " << emPosition.GetCoordinates(conex.GetObserverCS()) << " "
+ << emPosition.GetCoordinates(rootCS) << std::endl;
+
+ conex.addParticle(particles::Code::Proton, Eem, 0_eV, emPosition, momentum.normalized(),
+ 0_s);
+ // supperimpose a photon
+ auto const momentumPhoton = showerAxis.GetDirection() * 1_TeV;
+ conex.addParticle(particles::Code::Gamma, 1_TeV, 0_eV, emPosition,
+ momentumPhoton.normalized(), 0_s);
+ conex.SolveCE();
+}
+
+#include <algorithm>
+#include <iterator>
+#include <string>
+#include <fstream>
+
+TEST_CASE("ConexOutput", "[output validation]") {
+
+ auto file = GENERATE(as<std::string>{}, "conex_fit", "conex_output");
+
+ SECTION(std::string("check saved data, ") + file + ".txt") {
+
+ // compare to binary reference data
+ std::ifstream file1(file + ".txt");
+ std::ifstream file1ref(refDataDir + "/" + file + "_REF.txt");
+
+ std::istreambuf_iterator<char> begin1(file1);
+ std::istreambuf_iterator<char> begin1ref(file1ref);
+
+ std::istreambuf_iterator<char> end;
+
+ while (begin1 != end && begin1ref != end) {
+ CHECK(*begin1 == *begin1ref);
+ ++begin1;
+ ++begin1ref;
+ }
+ CHECK(begin1 == end);
+ CHECK(begin1ref == end);
+ file1.close();
+ file1ref.close();
+ }
+}
diff --git a/tests/modules/testExecTime.cpp b/tests/modules/testExecTime.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..e2fc710e7a201096184edac399befe01a8d04356
--- /dev/null
+++ b/tests/modules/testExecTime.cpp
@@ -0,0 +1,168 @@
+/*
+ * (c) Copyright 2019 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * 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.
+ */
+
+#define CATCH_CONFIG_MAIN // This tells Catch to provide a main() - only do this in one
+ // cpp file
+#include <catch2/catch.hpp>
+
+#include <corsika/process/analytic_processors/ExecTime.h>
+
+#include <corsika/process/example_processors/DummyBoundaryCrossingProcess.h>
+#include <corsika/process/example_processors/DummyContinuousProcess.h>
+#include <corsika/process/example_processors/DummyDecayProcess.h>
+#include <corsika/process/example_processors/DummyInteractionProcess.h>
+#include <corsika/process/example_processors/DummySecondariesProcess.h>
+
+#include <cmath>
+#include <random>
+#include <vector>
+
+using namespace corsika::process;
+using namespace corsika::process::analytic_processors;
+using namespace corsika::process::example_processors;
+
+TEST_CASE("Timing process", "[proccesses][analytic_processors ExecTime]") {
+
+ int tmp = 0;
+
+ SECTION("BoundaryCrossing") {
+ ExecTime<DummyBoundaryCrossingProcess<10>> execTime;
+ auto start = std::chrono::steady_clock::now();
+ CHECK(execTime.DoBoundaryCrossing(tmp, 0, 0) == EProcessReturn::eOk);
+ auto end = std::chrono::steady_clock::now();
+ CHECK(std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count() ==
+ Approx(10).margin(5));
+
+ for (int i = 0; i < 100; i++) execTime.DoBoundaryCrossing(tmp, 0, 0);
+
+ CHECK(execTime.mean() == Approx(10 * 1000).margin(2 * 1000));
+
+ CHECK(execTime.sumTime() == Approx(10 * 100 * 1000).margin((10 * 100) * 1000));
+
+ CHECK(fabs(execTime.var()) < 100000);
+ }
+
+ SECTION("Continuous") {
+ ExecTime<DummyContinuousProcess<50>> execTime;
+ auto start = std::chrono::steady_clock::now();
+ CHECK(execTime.DoContinuous(tmp, tmp) == EProcessReturn::eOk);
+ auto end = std::chrono::steady_clock::now();
+ CHECK(std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count() ==
+ Approx(50).margin(5));
+
+ for (int i = 0; i < 100; i++) execTime.DoContinuous(tmp, tmp);
+
+ CHECK(execTime.mean() == Approx(50 * 1000).margin(2 * 1000));
+
+ CHECK(execTime.sumTime() == Approx(50 * 100 * 1000).margin((10 * 100) * 1000));
+
+ CHECK(fabs(execTime.var()) < 100000);
+ }
+
+ SECTION("Decay") {
+ ExecTime<DummyDecayProcess<10>> execTime;
+ auto start = std::chrono::steady_clock::now();
+ CHECK(execTime.DoDecay(tmp) == EProcessReturn::eOk);
+ auto end = std::chrono::steady_clock::now();
+ CHECK(std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count() ==
+ Approx(10).margin(5));
+
+ for (int i = 0; i < 100; i++) execTime.DoDecay(tmp);
+
+ CHECK(execTime.mean() == Approx(10 * 1000).margin(2 * 1000));
+
+ CHECK(execTime.sumTime() == Approx(10 * 100 * 100).margin((10 * 100) * 1000));
+
+ CHECK(fabs(execTime.var()) < 100000);
+ }
+
+ SECTION("Interaction") {
+ ExecTime<DummyInteractionProcess<10>> execTime;
+ auto start = std::chrono::steady_clock::now();
+ CHECK(execTime.DoInteraction(tmp) == EProcessReturn::eOk);
+ auto end = std::chrono::steady_clock::now();
+ CHECK(std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count() ==
+ Approx(10).margin(5));
+
+ for (int i = 0; i < 100; i++) execTime.DoInteraction(tmp);
+
+ CHECK(execTime.mean() == Approx(10 * 1000).margin(2 * 1000));
+
+ CHECK(execTime.sumTime() == Approx(10 * 100 * 1000).margin((10 * 100) * 1000));
+
+ CHECK(fabs(execTime.var()) < 100000);
+ }
+
+ SECTION("Secondaries") {
+ ExecTime<DummySecondariesProcess<10>> execTime;
+ auto start = std::chrono::steady_clock::now();
+ CHECK(execTime.DoSecondaries(tmp) == EProcessReturn::eOk);
+ auto end = std::chrono::steady_clock::now();
+ CHECK(std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count() ==
+ Approx(10).margin(5));
+
+ for (int i = 0; i < 100; i++) execTime.DoSecondaries(tmp);
+
+ CHECK(execTime.mean() == Approx(10 * 1000).margin(2 * 1000));
+
+ CHECK(execTime.sumTime() == Approx(10 * 100 * 1000).margin((10 * 100) * 1000));
+
+ CHECK(fabs(execTime.var()) < 100000);
+ }
+
+ SECTION("TestMeanAlgo") {
+ std::default_random_engine generator;
+ std::normal_distribution<double> distribution(10000.0, 200.0);
+
+ double fElapsedSum = 0;
+ double fMean = 0;
+ double fMean2 = 0;
+ long long fMin = std::numeric_limits<long long>::max();
+ long long fMax = std::numeric_limits<long long>::min();
+ int fN = 0;
+
+ std::vector<double> elems;
+
+ for (int i = 0; i < 1000; i++) {
+ double timeDiv = distribution(generator);
+
+ elems.push_back(timeDiv);
+
+ fElapsedSum += timeDiv;
+ fN = fN + 1;
+
+ if (fMax < timeDiv) fMax = timeDiv;
+
+ if (timeDiv < fMin) fMin = timeDiv;
+
+ double delta = timeDiv - fMean;
+ fMean += delta / static_cast<double>(fN);
+
+ double delta2 = timeDiv - fMean;
+
+ fMean2 += delta * delta2;
+ }
+
+ CHECK(fN == 1000);
+
+ double mean = 0;
+ std::for_each(elems.begin(), elems.end(), [&](double i) { mean += i; });
+ mean = mean / fN;
+
+ double var = 0;
+ std::for_each(elems.begin(), elems.end(),
+ [&](double i) { var += (mean - i) * (mean - i); });
+ var = var / fN;
+
+ CHECK(mean == Approx(10000.0).margin(10));
+ CHECK(var == Approx(200.0 * 200).margin(2000));
+
+ CHECK(fMean2 / fN == Approx(200 * 200).margin(2000)); // Varianz
+ CHECK(fMean == Approx(10000).margin(10));
+ }
+}
diff --git a/tests/modules/testInteractionCounter.cpp b/tests/modules/testInteractionCounter.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..5417f3e3ec90ea7955e705e6fb836efd22198f8f
--- /dev/null
+++ b/tests/modules/testInteractionCounter.cpp
@@ -0,0 +1,135 @@
+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * 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/interaction_counter/InteractionCounter.hpp>
+
+#include <corsika/environment/Environment.h>
+#include <corsika/environment/HomogeneousMedium.h>
+#include <corsika/environment/NuclearComposition.h>
+#include <corsika/geometry/Point.h>
+#include <corsika/geometry/RootCoordinateSystem.h>
+#include <corsika/geometry/Vector.h>
+#include <corsika/units/PhysicalUnits.h>
+
+#include <corsika/setup/SetupStack.h>
+
+#include <catch2/catch.hpp>
+
+#include <numeric>
+
+using namespace corsika;
+using namespace corsika::process::interaction_counter;
+using namespace corsika::units;
+using namespace corsika::units::si;
+
+const std::string refDataDir = std::string(REFDATADIR); // from cmake
+
+struct DummyProcess {
+ template <typename TParticle>
+ GrammageType GetInteractionLength([[maybe_unused]] TParticle const& particle) {
+ return 100_g / 1_cm / 1_cm;
+ }
+
+ template <typename TParticle>
+ auto DoInteraction([[maybe_unused]] TParticle& projectile) {
+ return nullptr;
+ }
+};
+
+TEST_CASE("InteractionCounter", "[process]") {
+
+ logging::SetLevel(logging::level::debug);
+
+ DummyProcess d;
+ InteractionCounter countedProcess(d);
+
+ SECTION("GetInteractionLength") {
+ REQUIRE(countedProcess.GetInteractionLength(nullptr) == 100_g / 1_cm / 1_cm);
+ }
+
+ auto [env, csPtr, nodePtr] = setup::testing::setupEnvironment(particles::Code::Oxygen);
+ [[maybe_unused]] auto& env_dummy = env;
+
+ SECTION("DoInteraction nucleus") {
+ unsigned short constexpr A = 14, Z = 7;
+ auto [stackPtr, secViewPtr] = setup::testing::setupStack(particles::Code::Nucleus, A,
+ Z, 105_TeV, nodePtr, *csPtr);
+ REQUIRE(stackPtr->getEntries() == 1);
+ REQUIRE(secViewPtr->getEntries() == 0);
+
+ auto const ret = countedProcess.DoInteraction(*secViewPtr);
+ REQUIRE(ret == nullptr);
+
+ auto const& h = countedProcess.GetHistogram().labHist();
+ REQUIRE(h.at(h.axis(0).index(1'000'070'140), h.axis(1).index(1.05e14)) == 1);
+ REQUIRE(std::accumulate(h.cbegin(), h.cend(), 0) == 1);
+
+ auto const& h2 = countedProcess.GetHistogram().CMSHist();
+ REQUIRE(h2.at(h2.axis(0).index(1'000'070'140), h2.axis(1).index(1.6e12)) == 1);
+ REQUIRE(std::accumulate(h2.cbegin(), h2.cend(), 0) == 1);
+
+ countedProcess.GetHistogram().saveLab("testInteractionCounter_file1.npz",
+ utl::SaveMode::overwrite);
+ countedProcess.GetHistogram().saveCMS("testInteractionCounter_file2.npz",
+ utl::SaveMode::overwrite);
+ }
+
+ SECTION("DoInteraction Lambda") {
+ auto constexpr code = particles::Code::Lambda0;
+ auto [stackPtr, secViewPtr] =
+ setup::testing::setupStack(code, 0, 0, 105_TeV, nodePtr, *csPtr);
+ REQUIRE(stackPtr->getEntries() == 1);
+ REQUIRE(secViewPtr->getEntries() == 0);
+
+ auto const ret = countedProcess.DoInteraction(*secViewPtr);
+ REQUIRE(ret == nullptr);
+
+ auto const& h = countedProcess.GetHistogram().labHist();
+ REQUIRE(h.at(h.axis(0).index(3122), h.axis(1).index(1.05e14)) == 1);
+ REQUIRE(std::accumulate(h.cbegin(), h.cend(), 0) == 1);
+
+ auto const& h2 = countedProcess.GetHistogram().CMSHist();
+ REQUIRE(h2.at(h2.axis(0).index(3122), h2.axis(1).index(1.6e12)) == 1);
+ REQUIRE(std::accumulate(h2.cbegin(), h2.cend(), 0) == 1);
+ }
+}
+
+#include <algorithm>
+#include <iterator>
+#include <string>
+#include <fstream>
+
+TEST_CASE("InteractionCounterOutput", "[output validation]") {
+
+ auto file = GENERATE(as<std::string>{}, "testInteractionCounter_file1",
+ "testInteractionCounter_file2");
+
+ SECTION(std::string("check saved data, ") + file + ".npz") {
+
+ std::cout << file + ".npz vs " << refDataDir + "/" + file + "_REF.npz" << std::endl;
+
+ // compare to binary reference data
+ std::ifstream file1(file + ".npz");
+ std::ifstream file1ref(refDataDir + "/" + file + "_REF.npz");
+
+ std::istreambuf_iterator<char> begin1(file1);
+ std::istreambuf_iterator<char> begin1ref(file1ref);
+
+ std::istreambuf_iterator<char> end;
+
+ while (begin1 != end && begin1ref != end) {
+ CHECK(*begin1 == *begin1ref);
+ ++begin1;
+ ++begin1ref;
+ }
+ CHECK(begin1 == end);
+ CHECK(begin1ref == end);
+ file1.close();
+ file1ref.close();
+ }
+}
diff --git a/tests/modules/testOnShellCheck.cpp b/tests/modules/testOnShellCheck.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..77c64e178da2726fabcd7c3792af60d270417c25
--- /dev/null
+++ b/tests/modules/testOnShellCheck.cpp
@@ -0,0 +1,88 @@
+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * 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/on_shell_check/OnShellCheck.h>
+
+#include <corsika/environment/Environment.h>
+#include <corsika/geometry/FourVector.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::on_shell_check;
+using namespace corsika::units;
+using namespace corsika::units::si;
+
+TEST_CASE("OnShellCheck", "[processes]") {
+ feenableexcept(FE_INVALID);
+ using EnvType = setup::Environment;
+ EnvType env;
+ const geometry::CoordinateSystem& rootCS = env.GetCoordinateSystem();
+
+ // setup empty particle stack
+ setup::Stack stack;
+ stack.Clear();
+ // two energies
+ const HEPEnergyType E = 10_GeV;
+ // list of arbitrary particles
+ std::array const particleList{particles::Code::PiPlus, particles::Code::PiMinus,
+ particles::Code::Helium, particles::Code::Gamma};
+
+ std::array const mass_shifts{1.1, 1.001, 1.0, 1.0};
+
+ SECTION("check particle masses") {
+
+ OnShellCheck check(1.e-2, 0.01, false);
+
+ check.Init();
+
+ // 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, E,
+ 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
+ setup::StackView 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
+ int count = -1;
+ for (auto proType : particleList) {
+ count++;
+ const auto pz = sqrt((E - particles::GetMass(proType) * mass_shifts[count]) *
+ (E + particles::GetMass(proType) * mass_shifts[count]));
+ auto const momentum = corsika::stack::MomentumVector(rootCS, {0_GeV, 0_GeV, pz});
+ projectile.AddSecondary(std::tuple<particles::Code, units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, geometry::Point,
+ units::si::TimeType>{
+ proType, E, momentum, geometry::Point(rootCS, 0_m, 0_m, 0_m), 0_ns});
+ }
+ check.DoSecondaries(view);
+ int i = -1;
+ for (auto& p : view) {
+ i++;
+ auto const Plab = corsika::geometry::FourVector(p.GetEnergy(), p.GetMomentum());
+ auto const m_kinetic = Plab.GetNorm();
+ if (i == 0)
+ CHECK(m_kinetic / particles::PiPlus::GetMass() == Approx(1));
+ else if (i == 1)
+ CHECK_FALSE(m_kinetic / particles::PiMinus::GetMass() == Approx(1));
+ }
+ }
+}
diff --git a/tests/stack/testHistoryStack.cpp b/tests/stack/testHistoryStack.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..e636208300aafb24883817fa9d619898636c9389
--- /dev/null
+++ b/tests/stack/testHistoryStack.cpp
@@ -0,0 +1,64 @@
+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * 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/framework/stack/CombinedStack.hpp>
+#include <corsika/stack/DummyStack.hpp>
+#include <corsika/stack/history/HistoryStackExtension.hpp>
+
+#include <catch2/catch.hpp>
+
+using namespace corsika;
+
+// this is our dummy environment, it only knows its trivial BaseNodeType
+class DummyEvent {
+private:
+ size_t parent_;
+ std::vector<int> secondaries_;
+
+public:
+ DummyEvent() {}
+ DummyEvent(const size_t parent) { parent_ = parent; }
+
+ size_t getParentIndex() { return parent_; }
+ void addSecondary(const int particle) { secondaries_.push_back(particle); }
+ int multiplicity() const { return secondaries_.size(); }
+};
+
+// the GeometryNode stack needs to know the type of geometry-nodes from the DummyEnv:
+template <typename TStackIter>
+using DummyHistoryDataInterface =
+ typename history::MakeHistoryDataInterface<TStackIter, DummyEvent>::type;
+
+// combine dummy stack with geometry information for tracking
+template <typename TStackIter>
+using StackWithHistoryInterface =
+ CombinedParticleInterface<dummy_stack::DummyStack::pi_type, DummyHistoryDataInterface,
+ TStackIter>;
+
+using TestStack =
+ CombinedStack<typename dummy_stack::DummyStack::stack_implementation_type,
+ history::HistoryData<DummyEvent>, StackWithHistoryInterface>;
+
+using EvtPtr = std::shared_ptr<DummyEvent>;
+
+TEST_CASE("HistoryStackExtension", "[stack]") {
+
+ [[maybe_unused]] CoordinateSystemPtr const& dummyCS = get_root_CoordinateSystem();
+
+ const dummy_stack::NoData noData;
+ TestStack s;
+
+ auto p = s.addParticle(std::tuple<dummy_stack::NoData>{noData});
+
+ SECTION("add lone particle") {
+ CHECK(s.getEntries() == 1);
+
+ EvtPtr evt = p.getEvent();
+ CHECK(evt == nullptr);
+ }
+}
diff --git a/tests/stack/testHistoryView.cpp b/tests/stack/testHistoryView.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..c702f4c8ae8adf111c06821b1d0178b6f3355b74
--- /dev/null
+++ b/tests/stack/testHistoryView.cpp
@@ -0,0 +1,216 @@
+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * 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/stack/history/Event.hpp>
+#include <corsika/stack/history/HistorySecondaryProducer.hpp>
+#include <corsika/stack/history/HistoryStackExtension.hpp>
+
+#include <corsika/framework/stack/CombinedStack.hpp>
+#include <corsika/stack/DummyStack.hpp>
+#include <corsika/stack/NuclearStackExtension.hpp>
+
+#include <corsika/framework/logging/Logging.hpp>
+
+#include <catch2/catch.hpp>
+
+using namespace corsika;
+
+/**
+ Need to replicate setup::SetupStack in a maximally simplified
+ way, but with real particle data
+ */
+
+// combine dummy stack with geometry information for tracking
+template <typename TStackIter>
+using StackWithHistoryInterface = CombinedParticleInterface<
+ nuclear_stack::ParticleDataStack::pi_type,
+ history::HistoryEventDataInterface, TStackIter>;
+
+using TestStack = CombinedStack<
+ typename nuclear_stack::ParticleDataStack::stack_implementation_type,
+ history::HistoryEventData, StackWithHistoryInterface>;
+
+/*
+ See Issue 161
+
+ unfortunately clang does not support this in the same way (yet) as
+ gcc, so we have to distinguish here. If clang cataches up, we
+ could remove the clang branch here and also in
+ corsika::Cascade. The gcc code is much more generic and
+ universal. If we could do the gcc version, we won't had to define
+ StackView globally, we could do it with MakeView whereever it is
+ actually needed. Keep an eye on this!
+ */
+#if defined(__clang__)
+using TheTestStackView = SecondaryView<typename TestStack::StackImpl,
+ StackWithHistoryInterface,
+ history::HistorySecondaryProducer>;
+#elif defined(__GNUC__) || defined(__GNUG__)
+using TheTestStackView =
+ MakeView<TestStack, history::HistorySecondaryProducer>::type;
+#endif
+
+using TestStackView = TheTestStackView;
+
+template <typename Event>
+int count_generations(Event const* event) {
+ int genCounter = 0;
+ while (event) {
+ event = event->parentEvent().get();
+ genCounter++;
+ }
+
+ return genCounter;
+}
+
+TEST_CASE("HistoryStackExtensionView", "[stack]") {
+
+ CoordinateSystemPtr const& dummyCS = get_root_CoordinateSystem();
+
+ // in this test we only use one singel stack !
+ TestStack stack;
+
+ // add primary particle
+ auto p0 = stack.addParticle(
+ std::make_tuple(Code::Electron, 1.5_GeV,
+ MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
+ Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s));
+
+ CHECK(stack.getEntries() == 1);
+ corsika::history::EventPtr evt = p0.getEvent();
+ CHECK(evt == nullptr);
+ CHECK(count_generations(evt.get()) == 0);
+
+ SECTION("interface test, view") {
+
+ // add secondaries, 1st generation
+ TestStackView hview0(p0);
+
+ auto const ev0 = p0.getEvent();
+ CHECK(ev0 == nullptr);
+
+ CORSIKA_LOG_DEBUG("loop VIEW");
+
+ // add 5 secondaries
+ for (int i = 0; i < 5; ++i) {
+ auto sec = hview0.addSecondary(
+ std::make_tuple(Code::Electron, 1.5_GeV,
+ MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
+ Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s));
+
+ CHECK(sec.getParentEventIndex() == i);
+ CHECK(sec.getEvent() != nullptr);
+ CHECK(sec.getEvent()->parentEvent() == nullptr);
+ CHECK(count_generations(sec.getEvent().get()) == 1);
+ }
+
+ // read 1st genertion particle particle
+ auto p1 = stack.getNextParticle();
+ CHECK(count_generations(p1.getEvent().get()) == 1);
+
+ TestStackView hview1(p1);
+
+ auto const ev1 = p1.getEvent();
+
+ // add second generation of secondaries
+ // add 10 secondaries
+ for (int i = 0; i < 10; ++i) {
+ auto sec = hview1.addSecondary(
+ std::make_tuple(Code::Electron, 1.5_GeV,
+ MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
+ Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s));
+
+ CHECK(sec.getParentEventIndex() == i);
+ CHECK(sec.getEvent()->parentEvent() == ev1);
+ CHECK(sec.getEvent()->parentEvent()->parentEvent() == ev0);
+
+ CHECK(count_generations(sec.getEvent().get()) == 2);
+
+ const auto org_projectile = stack.at(sec.getEvent()->projectileIndex());
+ CHECK(org_projectile.getEvent() == sec.getEvent()->parentEvent());
+ }
+
+ // read 2nd genertion particle particle
+ auto p2 = stack.getNextParticle();
+
+ TestStackView hview2(p2);
+
+ auto const ev2 = p2.getEvent();
+
+ // add third generation of secondaries
+ // add 15 secondaries
+ for (int i = 0; i < 15; ++i) {
+ CORSIKA_LOG_TRACE("loop, view: " + std::to_string(i));
+
+ auto sec = hview2.addSecondary(
+ std::make_tuple(Code::Electron, 1.5_GeV,
+ MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
+ Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s));
+ CORSIKA_LOG_TRACE("loop, ---- ");
+
+ CHECK(sec.getParentEventIndex() == i);
+ CHECK(sec.getEvent()->parentEvent() == ev2);
+ CHECK(sec.getEvent()->parentEvent()->parentEvent() == ev1);
+ CHECK(sec.getEvent()->parentEvent()->parentEvent()->parentEvent() == ev0);
+
+ CHECK(count_generations(sec.getEvent().get()) == 3);
+ }
+ }
+
+ SECTION("also test projectile access") {
+
+ CORSIKA_LOG_TRACE("projectile test");
+
+ // add secondaries, 1st generation
+ TestStackView hview0(p0);
+ auto proj0 = hview0.getProjectile();
+ auto const ev0 = p0.getEvent();
+ CHECK(ev0 == nullptr);
+
+ CORSIKA_LOG_TRACE("loop");
+
+ // add 5 secondaries
+ for (int i = 0; i < 5; ++i) {
+ CORSIKA_LOG_TRACE("loop " + std::to_string(i));
+ auto sec = proj0.addSecondary(
+ std::make_tuple(Code::Electron, 1.5_GeV,
+ MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
+ Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s));
+
+ CHECK(sec.getParentEventIndex() == i);
+ CHECK(sec.getEvent() != nullptr);
+ CHECK(sec.getEvent()->parentEvent() == nullptr);
+ }
+ CHECK(stack.getEntries() == 6);
+
+ // read 1st genertion particle particle
+ auto p1 = stack.getNextParticle();
+
+ TestStackView hview1(p1);
+ auto proj1 = hview1.getProjectile();
+ auto const ev1 = p1.getEvent();
+
+ // add second generation of secondaries
+ // add 10 secondaries
+ for (unsigned int i = 0; i < 10; ++i) {
+ auto sec = proj1.addSecondary(
+ std::make_tuple(Code::Electron, 1.5_GeV,
+ MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
+ Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s));
+
+ CHECK(sec.getParentEventIndex() == int(i));
+ CHECK(sec.getEvent()->parentEvent() == ev1);
+ CHECK(sec.getEvent()->secondaries().size() == i + 1);
+ CHECK(sec.getEvent()->parentEvent()->parentEvent() == ev0);
+
+ const auto org_projectile = stack.at(sec.getEvent()->projectileIndex());
+ CHECK(org_projectile.getEvent() == sec.getEvent()->parentEvent());
+ }
+ CHECK(stack.getEntries() == 16);
+ }
+}