corsika/detail/modules/tracking/TrackingLeapFrogStraight.inl

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/geometry/Line.hpp>
+#include <corsika/framework/geometry/StraightTrajectory.hpp>
+#include <corsika/framework/geometry/Vector.hpp>
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+
+#include <corsika/modules/tracking/TrackingStraight.hpp>
+
+#include <type_traits>
+#include <utility>
+
+namespace corsika {
+
+  namespace tracking_leapfrog_straight {
+
+    template <typename Particle>
+    inline auto Tracking::getTrack(Particle& particle) {
+      VelocityVector initialVelocity =
+          particle.getMomentum() / particle.getEnergy() * constants::c;
+
+      Point const initialPosition = particle.getPosition();
+
+      CORSIKA_LOG_DEBUG(
+          "TrackingLeapFrogStraight pid: {}"
+          " , E = {} GeV \n"
+          "\tTracking pos: {} \n"
+          "\tTracking   p: {} GeV \n"
+          "\tTracking   v: {} ",
+          particle.getPID(), particle.getEnergy() / 1_GeV,
+          initialPosition.getCoordinates(),
+          particle.getMomentum().getComponents() / 1_GeV,
+          initialVelocity.getComponents());
+
+      typedef decltype(particle.getNode()) node_type;
+      node_type const volumeNode = particle.getNode();
+
+      // check if particle is moving at all
+      auto const absVelocity = initialVelocity.getNorm();
+      if (absVelocity * 1_s == 0_m) {
+        return std::make_tuple(
+            StraightTrajectory(Line(initialPosition, initialVelocity), 0_s), volumeNode);
+      }
+
+      // charge of the particle, and magnetic field
+      auto const charge = particle.getCharge();
+      auto magneticfield =
+          volumeNode->getModelProperties().getMagneticField(initialPosition);
+      auto const magnitudeB = magneticfield.getNorm();
+      CORSIKA_LOG_DEBUG("field={} uT, charge={}, magnitudeB={} uT",
+                        magneticfield.getComponents() / 1_uT, charge, magnitudeB / 1_T);
+      bool const no_deflection = (charge == 0 * constants::e) || magnitudeB == 0_T;
+
+      // check, where the first halve-step direction has geometric intersections
+      auto const [initialTrack, initialTrackNextVolume] =
+          tracking_line::Tracking::getTrack(particle);
+      //{ [[maybe_unused]] auto& initialTrackNextVolume_dum = initialTrackNextVolume; }
+      auto const initialTrackLength = initialTrack.getLength(1);
+
+      CORSIKA_LOG_DEBUG("initialTrack(0)={}, initialTrack(1)={}, initialTrackLength={}",
+                        initialTrack.getPosition(0).getCoordinates(),
+                        initialTrack.getPosition(1).getCoordinates(), initialTrackLength);
+
+      // if particle is non-deflectable, we are done:
+      if (no_deflection) {
+        CORSIKA_LOG_DEBUG("no deflection. tracking finished");
+        return std::make_tuple(initialTrack, initialTrackNextVolume);
+      }
+
+      HEPMomentumType const p_perp =
+          (particle.getMomentum() -
+           particle.getMomentum().getParallelProjectionOnto(magneticfield))
+              .getNorm();
+
+      if (p_perp == 0_GeV) {
+        // particle travels along, parallel to magnetic field. Rg is
+        // "0", return straight track here.
+        CORSIKA_LOG_TRACE("p_perp is 0_GeV --> parallel");
+        return std::make_tuple(initialTrack, initialTrackNextVolume);
+      }
+
+      LengthType const gyroradius = (convert_HEP_to_SI<MassType::dimension_type>(p_perp) *
+                                     constants::c / (abs(charge) * magnitudeB));
+
+      // we need to limit maximum step-length since we absolutely
+      // need to follow strongly curved trajectories segment-wise,
+      // at least if we don't employ concepts as "Helix
+      // Trajectories" or similar
+      double const maxRadians = 0.001;
+      LengthType const steplimit = 2 * cos(maxRadians) * sin(maxRadians) * gyroradius;
+      CORSIKA_LOG_DEBUG("gyroradius {}, Steplimit: {}", gyroradius, steplimit);
+
+      // calculate first halve step for "steplimit"
+      DirectionVector const initialDirection = particle.getDirection();
+
+      // avoid any intersections within first halve steplength
+      // aim for intersection in second halve step
+      LengthType const firstHalveSteplength =
+          std::min(steplimit / 2, initialTrackLength * firstFraction_);
+
+      CORSIKA_LOG_DEBUG("first halve step length {}, steplimit={}, initialTrackLength={}",
+                        firstHalveSteplength, steplimit, initialTrackLength);
+      // perform the first halve-step
+      Point const position_mid =
+          initialPosition + initialDirection * firstHalveSteplength;
+      auto const k = charge / (constants::c * convert_HEP_to_SI<MassType::dimension_type>(
+                                                  particle.getMomentum().getNorm()));
+      DirectionVector const new_direction =
+          initialDirection +
+          initialDirection.cross(magneticfield) * firstHalveSteplength * 2 * k;
+      auto const new_direction_norm = new_direction.getNorm(); // by design this is >1
+      CORSIKA_LOG_DEBUG(
+          "position_mid={}, new_direction={}, (new_direction_norm)={}, deflection={}",
+          position_mid.getCoordinates(), new_direction.getComponents(),
+          new_direction_norm,
+          acos(std::min(1.0, initialDirection.dot(new_direction) / new_direction_norm)) *
+              180 / M_PI);
+
+      // check, where the second halve-step direction has geometric intersections
+      particle.setPosition(position_mid);
+      particle.setDirection(new_direction);
+      auto const [finalTrack, finalTrackNextVolume] =
+          tracking_line::Tracking::getTrack(particle);
+      particle.setPosition(initialPosition);   // this is not nice...
+      particle.setDirection(initialDirection); // this is not nice...
+
+      LengthType const finalTrackLength = finalTrack.getLength(1);
+      LengthType const secondLeapFrogLength = firstHalveSteplength * new_direction_norm;
+
+      // check if volume transition is obvious, OR
+      // for numerical reasons, particles slighly bend "away" from a
+      // volume boundary have a very hard time to cross the border,
+      // thus, if secondLeapFrogLength is just slighly shorter (1e-4m) than
+      // finalTrackLength we better just [extend the
+      // secondLeapFrogLength slightly and] force the volume
+      // crossing:
+      bool const switch_volume = finalTrackLength - 0.0001_m <= secondLeapFrogLength;
+      LengthType const secondHalveStepLength =
+          std::min(secondLeapFrogLength, finalTrackLength);
+
+      CORSIKA_LOG_DEBUG(
+          "finalTrack(0)={}, finalTrack(1)={}, finalTrackLength={}, "
+          "secondLeapFrogLength={}, secondHalveStepLength={}, "
+          "secondLeapFrogLength-finalTrackLength={}, "
+          "secondHalveStepLength-finalTrackLength={}, "
+          "nextVol={}, transition={}",
+          finalTrack.getPosition(0).getCoordinates(),
+          finalTrack.getPosition(1).getCoordinates(), finalTrackLength,
+          secondLeapFrogLength, secondHalveStepLength,
+          secondLeapFrogLength - finalTrackLength,
+          secondHalveStepLength - finalTrackLength, fmt::ptr(finalTrackNextVolume),
+          switch_volume);
+
+      // perform the second halve-step
+      auto const new_direction_normalized = new_direction.normalized();
+      Point const finalPosition =
+          position_mid + new_direction_normalized * secondHalveStepLength;
+
+      LengthType const totalStep = firstHalveSteplength + secondHalveStepLength;
+      auto const delta_pos = finalPosition - initialPosition;
+      auto const distance = delta_pos.getNorm();
+
+      return std::make_tuple(
+          StraightTrajectory(
+              Line(initialPosition,
+                   (distance == 0_m ? initialVelocity
+                                    : delta_pos.normalized() * absVelocity)),
+              distance / absVelocity,  // straight distance
+              totalStep / absVelocity, // bend distance
+              initialVelocity,
+              new_direction_normalized * absVelocity), // trajectory
+          (switch_volume ? finalTrackNextVolume : volumeNode));
+    }
+
+  } // namespace tracking_leapfrog_straight
+
+} // namespace corsika

corsika/detail/modules/tracking/TrackingStraight.inl

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/core/Logging.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/geometry/Intersections.hpp>
+#include <corsika/framework/geometry/Line.hpp>
+#include <corsika/framework/geometry/Plane.hpp>
+#include <corsika/framework/geometry/Sphere.hpp>
+#include <corsika/framework/geometry/StraightTrajectory.hpp>
+#include <corsika/framework/geometry/Vector.hpp>
+#include <corsika/modules/tracking/Intersect.hpp>
+
+#include <cmath>
+#include <type_traits>
+#include <utility>
+
+namespace corsika::tracking_line {
+
+  template <typename TParticle>
+  inline auto Tracking::getTrack(TParticle const& particle) {
+    VelocityVector const initialVelocity =
+        particle.getMomentum() / particle.getEnergy() * constants::c;
+
+    auto const& initialPosition = particle.getPosition();
+    CORSIKA_LOG_DEBUG(
+        "TrackingStraight pid: {}"
+        " , E = {} GeV \n"
+        "\tTracking pos: {} \n"
+        "\tTracking   p: {} GeV \n"
+        "\tTracking   v: {}",
+        particle.getPID(), particle.getEnergy() / 1_GeV, initialPosition.getCoordinates(),
+        particle.getMomentum().getComponents() / 1_GeV, initialVelocity.getComponents());
+
+    // traverse the environment volume tree and find next
+    // intersection
+    auto [minTime, minNode] = nextIntersect(particle);
+
+    return std::make_tuple(StraightTrajectory(Line(initialPosition, initialVelocity),
+                                              minTime), // trajectory
+                           minNode);                    // next volume node
+  }
+
+  template <typename TParticle>
+  inline Intersections Tracking::intersect(TParticle const& particle,
+                                           Sphere const& sphere) {
+    auto const& position = particle.getPosition();
+    auto const delta = position - sphere.getCenter();
+    auto const velocity = particle.getMomentum() / particle.getEnergy() * constants::c;
+    auto const vSqNorm = velocity.getSquaredNorm();
+    auto const R = sphere.getRadius();
+
+    auto const vDotDelta = velocity.dot(delta);
+    auto const discriminant =
+        vDotDelta * vDotDelta - vSqNorm * (delta.getSquaredNorm() - R * R);
+
+    if (discriminant.magnitude() > 0) {
+      auto const sqDisc = sqrt(discriminant);
+      auto const invDenom = 1 / vSqNorm;
+
+      CORSIKA_LOG_TRACE("numericallyInside={}", sphere.contains(position));
+      return Intersections((-vDotDelta - sqDisc) * invDenom,
+                           (-vDotDelta + sqDisc) * invDenom);
+    }
+    return Intersections();
+  }
+
+  template <typename TParticle>
+  inline Intersections Tracking::intersect(TParticle const& particle, Box const& box) {
+    Point const& position = particle.getPosition();
+    VelocityVector const velocity =
+        particle.getMomentum() / particle.getEnergy() * constants::c;
+    CoordinateSystemPtr const& cs = box.getCoordinateSystem();
+    LengthType x0 = position.getX(cs);
+    LengthType y0 = position.getY(cs);
+    LengthType z0 = position.getZ(cs);
+    SpeedType vx = velocity.getX(cs);
+    SpeedType vy = velocity.getY(cs);
+    SpeedType vz = velocity.getZ(cs);
+    CORSIKA_LOG_TRACE(
+        "particle in box coordinate: position: ({:.3f}, {:.3f}, "
+        "{:.3f}) m, veolocity: ({:.3f}, {:.3f}, {:.3f}) m/ns",
+        x0 / 1_m, y0 / 1_m, z0 / 1_m, vx / (1_m / 1_ns), vy / (1_m / 1_ns),
+        vz / (1_m / 1_ns));
+
+    auto get_intersect_min_max = [](LengthType x0, SpeedType v0, LengthType dx) {
+      auto t1 = (dx - x0) / v0;
+      auto t2 = (-dx - x0) / v0;
+      if (t1 > t2)
+        return std::make_pair(t1, t2);
+      else
+        return std::make_pair(t2, t1);
+    };
+
+    auto [tx_max, tx_min] = get_intersect_min_max(x0, vx, box.getX());
+    auto [ty_max, ty_min] = get_intersect_min_max(y0, vy, box.getY());
+    auto [tz_max, tz_min] = get_intersect_min_max(z0, vz, box.getZ());
+
+    TimeType t_exit = std::min(std::min(tx_max, ty_max), tz_max);
+    TimeType t_enter = std::max(std::max(tx_min, ty_min), tz_min);
+
+    CORSIKA_LOG_DEBUG("t_enter: {} ns, t_exit: {} ns", t_enter / 1_ns, t_exit / 1_ns);
+    if ((t_exit > t_enter)) {
+      if (t_enter < 0_s && t_exit > 0_s)
+        CORSIKA_LOG_DEBUG("numericallyInside={}", box.contains(position));
+      else if (t_enter < 0_s && t_exit < 0_s)
+        CORSIKA_LOG_DEBUG("oppisite direction");
+      return Intersections(std::move(t_enter), std::move(t_exit));
+    } else
+      return Intersections();
+  }
+
+  template <typename TParticle, typename TBaseNodeType>
+  inline Intersections Tracking::intersect(TParticle const& particle,
+                                           TBaseNodeType const& volumeNode) {
+    if (Sphere const* sphere = dynamic_cast<Sphere const*>(&volumeNode.getVolume());
+        sphere) {
+      return Tracking::intersect<TParticle>(particle, *sphere);
+    } else if (Box const* box = dynamic_cast<Box const*>(&volumeNode.getVolume()); box) {
+      return Tracking::intersect<TParticle>(particle, *box);
+    } else if (SeparationPlane const* sepPlane =
+                   dynamic_cast<SeparationPlane const*>(&volumeNode.getVolume());
+               sepPlane) {
+      return Tracking::intersect<TParticle>(particle, *sepPlane);
+    } else {
+      throw std::runtime_error(
+          "The Volume type provided is not supported in "
+          "TrackingStraight::intersect(particle, node)");
+    }
+  }
+
+  template <typename TParticle>
+  inline Intersections Tracking::intersect(TParticle const& particle,
+                                           Plane const& plane) {
+    auto const delta = plane.getCenter() - particle.getPosition();
+    auto const velocity = particle.getMomentum() / particle.getEnergy() * constants::c;
+    auto const n = plane.getNormal();
+    auto const n_dot_v = n.dot(velocity);
+
+    CORSIKA_LOG_TRACE("n_dot_v={}, delta={}, momentum={}", n_dot_v, delta,
+                      particle.getMomentum());
+
+    if (n_dot_v.magnitude() == 0)
+      return Intersections();
+    else
+      return Intersections(n.dot(delta) / n_dot_v);
+  }
+
+  template <typename TParticle>
+  inline Intersections Tracking::intersect(TParticle const& particle,
+                                           SeparationPlane const& sepPlane) {
+    return intersect(particle, sepPlane.getPlane());
+  }
+
+} // namespace corsika::tracking_line

corsika/detail/modules/urqmd/ParticleConversion.inl

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/modules/urqmd/ParticleConversion.hpp>
+
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+
+#include <urqmd.hpp>
+
+namespace corsika::urqmd {
+
+  inline bool canInteract(Code const vCode) {
+    // According to the manual, UrQMD can use all mesons, baryons and nucleons
+    // which are modeled also as input particles. I think it is safer to accept
+    // only the usual long-lived species as input.
+    // TODO: Charmed mesons should be added to the list, too
+
+    static std::array constexpr validProjectileCodes{
+        Code::Proton,  Code::AntiProton, Code::Neutron, Code::AntiNeutron, Code::PiPlus,
+        Code::PiMinus, Code::KPlus,      Code::KMinus,  Code::K0Short,     Code::K0Long};
+
+    return std::find(std::cbegin(validProjectileCodes), std::cend(validProjectileCodes),
+                     vCode) != std::cend(validProjectileCodes);
+  }
+
+  inline std::pair<int, int> convertToUrQMD(Code const code) {
+    static const std::map<int, std::pair<int, int>> mapPDGToUrQMD{
+        // data mostly from github.com/afedynitch/ParticleDataTool
+        {22, {100, 0}},      // photon
+        {111, {101, 0}},     // pi0
+        {211, {101, 2}},     // pi+
+        {-211, {101, -2}},   // pi-
+        {321, {106, 1}},     // K+
+        {-321, {-106, -1}},  // K-
+        {311, {106, -1}},    // K0
+        {-311, {-106, 1}},   // K0bar
+        {2212, {1, 1}},      // p
+        {2112, {1, -1}},     // n
+        {-2212, {-1, -1}},   // pbar
+        {-2112, {-1, 1}},    // nbar
+        {221, {102, 0}},     // eta
+        {213, {104, 2}},     // rho+
+        {-213, {104, -2}},   // rho-
+        {113, {104, 0}},     // rho0
+        {323, {108, 2}},     // K*+
+        {-323, {108, -2}},   // K*-
+        {313, {108, 0}},     // K*0
+        {-313, {-108, 0}},   // K*0-bar
+        {223, {103, 0}},     // omega
+        {333, {109, 0}},     // phi
+        {3222, {40, 2}},     // Sigma+
+        {3212, {40, 0}},     // Sigma0
+        {3112, {40, -2}},    // Sigma-
+        {3322, {49, 0}},     // Xi0
+        {3312, {49, -1}},    // Xi-
+        {3122, {27, 0}},     // Lambda0
+        {2224, {17, 4}},     // Delta++
+        {2214, {17, 2}},     // Delta+
+        {2114, {17, 0}},     // Delta0
+        {1114, {17, -2}},    // Delta-
+        {3224, {41, 2}},     // Sigma*+
+        {3214, {41, 0}},     // Sigma*0
+        {3114, {41, -2}},    // Sigma*-
+        {3324, {50, 0}},     // Xi*0
+        {3314, {50, -1}},    // Xi*-
+        {3334, {55, 0}},     // Omega-
+        {411, {133, 2}},     // D+
+        {-411, {133, -2}},   // D-
+        {421, {133, 0}},     // D0
+        {-421, {-133, 0}},   // D0-bar
+        {441, {107, 0}},     // etaC
+        {431, {138, 1}},     // Ds+
+        {-431, {138, -1}},   // Ds-
+        {433, {139, 1}},     // Ds*+
+        {-433, {139, -1}},   // Ds*-
+        {413, {134, 1}},     // D*+
+        {-413, {134, -1}},   // D*-
+        {10421, {134, 0}},   // D*0
+        {-10421, {-134, 0}}, // D*0-bar
+        {443, {135, 0}},     // jpsi
+    };
+
+    return mapPDGToUrQMD.at(static_cast<int>(get_PDG(code)));
+  }
+
+  inline Code convertFromUrQMD(int vItyp, int vIso3) {
+    int const pdgInt =
+        ::urqmd::pdgid_(vItyp, vIso3); // use the conversion function provided by UrQMD
+    if (pdgInt == 0) {                 // ::urqmd::pdgid_ returns 0 on error
+      throw std::runtime_error("UrQMD pdgid() returned 0");
+    }
+    auto const pdg = static_cast<PDGCode>(pdgInt);
+    return convert_from_PDG(pdg);
+  }
+
+} // namespace corsika::urqmd

corsika/detail/modules/urqmd/UrQMD.inl

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/modules/Random.hpp>
+#include <corsika/modules/urqmd/UrQMD.hpp>
+#include <corsika/modules/urqmd/ParticleConversion.hpp>
+
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/core/EnergyMomentumOperations.hpp>
+#include <corsika/framework/geometry/QuantityVector.hpp>
+#include <corsika/framework/geometry/Vector.hpp>
+#include <corsika/framework/utility/COMBoost.hpp>
+
+#include <boost/filesystem.hpp>
+#include <boost/filesystem/fstream.hpp>
+#include <boost/multi_array.hpp>
+
+#include <algorithm>
+#include <cassert>
+#include <functional>
+#include <iostream>
+#include <fstream>
+#include <sstream>
+
+#include <urqmd.hpp>
+
+namespace corsika::urqmd {
+
+  inline UrQMD::UrQMD(boost::filesystem::path xs_file, int const retryFlag)
+      : iflb_(retryFlag) {
+    readXSFile(xs_file);
+    corsika::connect_random_stream(RNG_, ::urqmd::set_rng_function);
+    ::urqmd::iniurqmdc8_();
+  }
+
+  inline bool UrQMD::isValid(Code const projectileId, Code const targetId) const {
+
+    if (!is_hadron(projectileId) || !corsika::urqmd::canInteract(projectileId)) {
+      return false;
+    }
+    if (!is_nucleus(targetId)) { return false; }
+    return true;
+  }
+
+  inline CrossSectionType UrQMD::getTabulatedCrossSection(
+      Code const projectileId, Code const targetId, HEPEnergyType const labEnergy) const {
+
+    // translated to C++ from CORSIKA 7 subroutine cxtot_u
+
+    auto const kinEnergy = labEnergy - get_mass(projectileId);
+
+    assert(kinEnergy >= HEPEnergyType::zero());
+
+    double const logKinEnergy = std::log10(kinEnergy * (1 / 1_GeV));
+    double const ye = std::max(10 * logKinEnergy + 10.5, 1.);
+    int const je = std::min(int(ye), int(xs_interp_support_table_.shape()[2] - 2));
+    std::array<double, 3> w;
+    w[2 - 1] = ye - je;
+    w[3 - 1] = w[2 - 1] * (w[2 - 1] - 1.) * .5;
+    w[1 - 1] = 1 - w[2 - 1] + w[3 - 1];
+    w[2 - 1] = w[2 - 1] - 2 * w[3 - 1];
+
+    int projectileIndex;
+    switch (projectileId) {
+      case Code::Proton:
+        projectileIndex = 0;
+        break;
+      case Code::AntiProton:
+        projectileIndex = 1;
+        break;
+      case Code::Neutron:
+        projectileIndex = 2;
+        break;
+      case Code::AntiNeutron:
+        projectileIndex = 3;
+        break;
+      case Code::PiPlus:
+        projectileIndex = 4;
+        break;
+      case Code::PiMinus:
+        projectileIndex = 5;
+        break;
+      case Code::KPlus:
+        projectileIndex = 6;
+        break;
+      case Code::KMinus:
+        projectileIndex = 7;
+        break;
+      case Code::K0Short:
+      case Code::K0Long:
+      /* since K0Short and K0Long are treated the same, we can also add K0 and K0Bar
+       * to the list. This is a deviation from CORSIKA 7. */
+      case Code::K0:
+      case Code::K0Bar:
+        projectileIndex = 8;
+        break;
+      default: { // LCOV_EXCL_START since this can never happen due to canInteract
+        CORSIKA_LOG_WARN("UrQMD cross-section not tabulated for {}", projectileId);
+        return CrossSectionType::zero();
+        // LCOV_EXCL_STOP
+      }
+    }
+
+    int targetIndex;
+    switch (targetId) {
+      case Code::Nitrogen:
+        targetIndex = 0;
+        break;
+      case Code::Oxygen:
+        targetIndex = 1;
+        break;
+      case Code::Argon:
+        targetIndex = 2;
+        break;
+      default:
+        std::stringstream ss;
+        ss << "UrQMD cross-section not tabluated for target " << targetId;
+        throw std::runtime_error(ss.str().data());
+    }
+
+    auto result = CrossSectionType::zero();
+    for (int i = 0; i < 3; ++i) {
+      result +=
+          xs_interp_support_table_[projectileIndex][targetIndex][je + i - 1 - 1] * w[i];
+    }
+
+    CORSIKA_LOG_TRACE(
+        "UrQMD::GetTabulatedCrossSection proj={}, targ={}, E={} GeV, sigma={}",
+        get_name(projectileId), get_name(targetId), labEnergy / 1_GeV, result);
+
+    return result;
+  }
+
+  inline CrossSectionType UrQMD::getCrossSection(Code const projectileId,
+                                                 Code const targetId,
+                                                 FourMomentum const& projectileP4,
+                                                 FourMomentum const& targetP4) const {
+
+    if (!isValid(projectileId, targetId)) {
+      /*
+       * unfortunately unavoidable at the moment until we have tools to get the actual
+       * inealstic cross-section from UrQMD
+       */
+      return CrossSectionType::zero();
+    }
+
+    // define projectile, in lab frame
+    auto const sqrtS2 = (projectileP4 + targetP4).getNormSqr();
+    HEPEnergyType const Elab = (sqrtS2 - static_pow<2>(get_mass(projectileId)) -
+                                static_pow<2>(get_mass(targetId))) /
+                               (2 * get_mass(targetId));
+
+    bool const tabulated = true;
+    if (tabulated) { return getTabulatedCrossSection(projectileId, targetId, Elab); }
+
+    // the following is a translation of ptsigtot() into C++
+    if (!is_nucleus(projectileId) &&
+        !is_nucleus(targetId)) { // both particles are "special"
+
+      double sqrtS = sqrt(sqrtS2) / 1_GeV;
+
+      // we must set some UrQMD globals first...
+      auto const [ityp, iso3] = corsika::urqmd::convertToUrQMD(projectileId);
+      ::urqmd::inputs_.spityp[0] = ityp;
+      ::urqmd::inputs_.spiso3[0] = iso3;
+
+      auto const [itypTar, iso3Tar] = corsika::urqmd::convertToUrQMD(targetId);
+      ::urqmd::inputs_.spityp[1] = itypTar;
+      ::urqmd::inputs_.spiso3[1] = iso3Tar;
+
+      int one = 1;
+      int two = 2;
+      return ::urqmd::sigtot_(one, two, sqrtS) * 1_mb;
+    }
+
+    // at least one of them is a nucleus
+    int const Ap = is_nucleus(projectileId) ? get_nucleus_A(projectileId) : 1;
+    int const At = is_nucleus(targetId) ? get_nucleus_A(targetId) : 1;
+    double const maxImpact = ::urqmd::nucrad_(Ap) + ::urqmd::nucrad_(At) +
+                             2 * ::urqmd::options_.CTParam[30 - 1];
+    return 10_mb * M_PI * static_pow<2>(maxImpact);
+    // is a constant cross-section really reasonable?
+  }
+
+  template <typename TView>
+  inline void UrQMD::doInteraction(TView& view, Code const projectileId,
+                                   Code const targetId, FourMomentum const& projectileP4,
+                                   FourMomentum const& targetP4) {
+
+    // define projectile, in lab frame
+    auto const sqrtS2 = (projectileP4 + targetP4).getNormSqr();
+    HEPEnergyType const Elab = (sqrtS2 - static_pow<2>(get_mass(projectileId)) -
+                                static_pow<2>(get_mass(targetId))) /
+                               (2 * get_mass(targetId));
+
+    if (!isValid(projectileId, targetId)) {
+      throw std::runtime_error("invalid target,projectile,energy combination");
+    }
+
+    size_t const targetA = get_nucleus_A(targetId);
+    size_t const targetZ = get_nucleus_Z(targetId);
+
+    ::urqmd::inputs_.nevents = 1;
+    ::urqmd::sys_.eos = 0; // could be configurable in principle
+    ::urqmd::inputs_.outsteps = 1;
+    ::urqmd::sys_.nsteps = 1;
+
+    // initialization regarding projectile
+    if (is_nucleus(projectileId)) {
+      // is this everything?
+      ::urqmd::inputs_.prspflg = 0;
+
+      ::urqmd::sys_.Ap = get_nucleus_A(projectileId);
+      ::urqmd::sys_.Zp = get_nucleus_Z(projectileId);
+      ::urqmd::rsys_.ebeam = (Elab - get_mass(projectileId)) / 1_GeV / ::urqmd::sys_.Ap;
+
+      ::urqmd::rsys_.bdist = ::urqmd::nucrad_(targetA) +
+                             ::urqmd::nucrad_(::urqmd::sys_.Ap) +
+                             2 * ::urqmd::options_.CTParam[30 - 1];
+
+      int const id = 1;
+      ::urqmd::cascinit_(::urqmd::sys_.Zp, ::urqmd::sys_.Ap, id);
+    } else {
+      ::urqmd::inputs_.prspflg = 1;
+      ::urqmd::sys_.Ap =
+          1; // even for non-baryons this has to be set, see vanilla UrQMD.f
+      ::urqmd::rsys_.bdist = ::urqmd::nucrad_(targetA) + ::urqmd::nucrad_(1) +
+                             2 * ::urqmd::options_.CTParam[30 - 1];
+      ::urqmd::rsys_.ebeam = (Elab - get_mass(projectileId)) / 1_GeV;
+
+      auto const [ityp, iso3] = corsika::urqmd::convertToUrQMD(
+          (projectileId == Code::K0Long || projectileId == Code::K0Short)
+              ? (booleanDist_(RNG_) ? Code::K0 : Code::K0Bar)
+              : projectileId);
+      // todo: conversion of K_long/short into strong eigenstates;
+      ::urqmd::inputs_.spityp[0] = ityp;
+      ::urqmd::inputs_.spiso3[0] = iso3;
+    }
+
+    // initialization regarding target
+    if (is_nucleus(targetId)) {
+      ::urqmd::sys_.Zt = targetZ;
+      ::urqmd::sys_.At = targetA;
+      ::urqmd::inputs_.trspflg = 0; // nucleus as target
+      int const id = 2;
+      ::urqmd::cascinit_(::urqmd::sys_.Zt, ::urqmd::sys_.At, id);
+    } else {
+      ::urqmd::inputs_.trspflg = 1; // special particle as target
+      auto const [ityp, iso3] = corsika::urqmd::convertToUrQMD(targetId);
+      ::urqmd::inputs_.spityp[1] = ityp;
+      ::urqmd::inputs_.spiso3[1] = iso3;
+    }
+
+    int iflb =
+        iflb_; // flag for retrying interaction in case of empty event, 0 means retry
+    ::urqmd::urqmd_(iflb);
+
+    // now retrieve secondaries from UrQMD
+    COMBoost const boost(projectileP4, targetP4);
+    auto const& originalCS = boost.getOriginalCS();
+    auto const& csPrime = boost.getRotatedCS();
+
+    for (int i = 0; i < ::urqmd::sys_.npart; ++i) {
+      auto code = corsika::urqmd::convertFromUrQMD(::urqmd::isys_.ityp[i],
+                                                   ::urqmd::isys_.iso3[i]);
+      if (code == Code::K0 || code == Code::K0Bar) {
+        code = booleanDist_(RNG_) ? Code::K0Short : Code::K0Long;
+      }
+
+      // "coor_.p0[i] * 1_GeV" is likely off-shell as UrQMD doesn't preserve masses well
+      MomentumVector momentum{csPrime,
+                              {::urqmd::coor_.px[i] * 1_GeV, ::urqmd::coor_.py[i] * 1_GeV,
+                               ::urqmd::coor_.pz[i] * 1_GeV}};
+
+      momentum.rebase(originalCS); // transform back into standard lab frame
+      CORSIKA_LOG_DEBUG(" {} {} {} ", i, code, momentum.getComponents());
+
+      HEPEnergyType const mass = get_mass(code);
+      HEPEnergyType const Ekin = calculate_kinetic_energy(momentum.getNorm(), mass);
+      if (Ekin <= 0_GeV) {
+        if (Ekin < 0_GeV) {
+          CORSIKA_LOG_WARN("Negative kinetic energy {} {}. Skipping.", code, Ekin);
+        }
+        view.addSecondary(
+            std::make_tuple(code, 0_eV, DirectionVector{originalCS, {0, 0, 0}}));
+      } else {
+        view.addSecondary(std::make_tuple(code, Ekin, momentum.normalized()));
+      }
+    }
+    CORSIKA_LOG_DEBUG("UrQMD generated {} secondaries!", ::urqmd::sys_.npart);
+  }
+
+  inline void UrQMD::readXSFile(boost::filesystem::path const filename) {
+    boost::filesystem::ifstream file(filename, std::ios::in);
+
+    if (!file.is_open()) {
+      // LCOV_EXCL_START since this is pointless to test
+      throw std::runtime_error(filename.native() + " could not be opened.");
+      // LCOV_EXCL_STOP
+    }
+
+    std::string line;
+
+    std::getline(file, line);
+    std::stringstream ss(line);
+
+    char dummy; // this is '#'
+    int nTargets, nProjectiles, nSupports;
+    ss >> dummy >> nTargets >> nProjectiles >> nSupports;
+
+    decltype(xs_interp_support_table_)::extent_gen extents;
+    xs_interp_support_table_.resize(extents[nProjectiles][nTargets][nSupports]);
+
+    for (int i = 0; i < nTargets; ++i) {
+      for (int j = 0; j < nProjectiles; ++j) {
+        for (int k = 0; k < nSupports; ++k) {
+          std::getline(file, line);
+          std::stringstream s(line);
+          double energy, sigma;
+          s >> energy >> sigma;
+          xs_interp_support_table_[j][i][k] = sigma * 1_mb;
+        }
+
+        std::getline(file, line);
+        std::getline(file, line);
+      }
+    }
+    file.close();
+  }
+
+} // namespace corsika::urqmd

corsika/detail/modules/writers/EnergyLossWriter.inl

+/*
+ * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/utility/FindXmax.hpp>
+
+#include <corsika/media/ShowerAxis.hpp>
+
+#include <exception>
+
+namespace corsika {
+
+  template <typename TOutput>
+  inline EnergyLossWriter<TOutput>::EnergyLossWriter(ShowerAxis const& axis,
+                                                     GrammageType dX,
+                                                     GrammageType dX_threshold)
+      : EnergyLossWriter<TOutput>{axis,
+                                  static_cast<unsigned int>(axis.getMaximumX() / dX) + 1,
+                                  dX, dX_threshold} {}
+
+  template <typename TOutput>
+  inline EnergyLossWriter<TOutput>::EnergyLossWriter(ShowerAxis const& axis,
+                                                     unsigned int const nBins,
+                                                     GrammageType dX,
+                                                     GrammageType dX_threshold)
+      : TOutput(dEdX_output::ProfileIndexNames)
+      , showerAxis_(axis)
+      , dX_(dX)
+      , nBins_(nBins)
+      , dX_threshold_(dX_threshold) {}
+
+  template <typename TOutput>
+  inline void EnergyLossWriter<TOutput>::startOfLibrary(
+      boost::filesystem::path const& directory) {
+    TOutput::startOfLibrary(directory);
+  }
+
+  template <typename TOutput>
+  inline void EnergyLossWriter<TOutput>::startOfShower(unsigned int const showerId) {
+    TOutput::startOfShower(showerId);
+    // reset profile
+    profile_.clear();
+    profile_.resize(nBins_);
+  }
+
+  template <typename TOutput>
+  inline void EnergyLossWriter<TOutput>::endOfLibrary() {
+    TOutput::endOfLibrary();
+  }
+
+  template <typename TOutput>
+  inline void EnergyLossWriter<TOutput>::endOfShower(unsigned int const showerId) {
+
+    int iRow{0};
+    for (dEdX_output::Profile const& row : profile_) {
+      // here: write to underlying writer (e.g. parquet)
+      TOutput::write(showerId, iRow * dX_, row);
+      iRow++;
+    }
+
+    TOutput::endOfShower(showerId);
+  }
+
+  template <typename TOutput>
+  inline void EnergyLossWriter<TOutput>::write(Point const& p0, Point const& p1,
+                                               Code const PID, HEPEnergyType const dE) {
+
+    GrammageType grammageStart = showerAxis_.getProjectedX(p0);
+    GrammageType grammageEnd = showerAxis_.getProjectedX(p1);
+
+    if (grammageStart > grammageEnd) { // particle going upstream
+      std::swap(grammageStart, grammageEnd);
+    }
+
+    GrammageType const deltaX = grammageEnd - grammageStart;
+
+    CORSIKA_LOGGER_TRACE(
+        TOutput::getLogger(),
+        "dE={} GeV, grammageStart={} g/cm2, End={}g /cm2, deltaX={} g/cm2", dE / 1_GeV,
+        grammageStart / 1_g * square(1_cm), grammageEnd / 1_g * square(1_cm),
+        deltaX / 1_g * square(1_cm));
+
+    if (deltaX < dX_threshold_) {
+      CORSIKA_LOGGER_TRACE(TOutput::getLogger(), "Point-like dE");
+      this->write(p0, PID, dE);
+      return;
+    }
+
+    // only register the range that is covered by the profile
+    int const maxBin = int(profile_.size() - 1);
+    int binStart = grammageStart / dX_;
+    if (binStart < 0) binStart = 0;
+    if (binStart > maxBin) binStart = maxBin;
+    int binEnd = grammageEnd / dX_;
+    if (binEnd < 0) binEnd = 0;
+    if (binEnd > maxBin) binEnd = maxBin;
+
+    CORSIKA_LOGGER_TRACE(TOutput::getLogger(), "maxBin={}, binStart={}, binEnd={}",
+                         maxBin, binStart, binEnd);
+
+    auto energyCount = HEPEnergyType::zero();
+
+    auto const factor = dE / deltaX; // [ energy / grammage ]
+    auto fill = [&](int const bin, GrammageType const weight) {
+      auto const increment = factor * weight;
+      CORSIKA_LOGGER_TRACE(TOutput::getLogger(),
+                           "filling bin={} with weight {} : dE={} GeV ", bin, weight,
+                           increment / 1_GeV);
+      profile_[bin][static_cast<int>(dEdX_output::ProfileIndex::Total)] += increment;
+      energyCount += increment;
+    };
+
+    // fill longitudinal profile
+    if (binStart == binEnd) {
+      fill(binStart, deltaX);
+    } else {
+      fill(binStart, ((1 + binStart) * dX_ - grammageStart));
+      fill(binEnd, (grammageEnd - binEnd * dX_));
+      for (int bin = binStart + 1; bin < binEnd; ++bin) { fill(bin, dX_); }
+    }
+
+    CORSIKA_LOGGER_TRACE(TOutput::getLogger(), "total energy added to histogram: {} GeV ",
+                         energyCount / 1_GeV);
+  }
+
+  template <typename TOutput>
+  inline void EnergyLossWriter<TOutput>::write(Point const& point, Code const,
+                                               HEPEnergyType const dE) {
+    GrammageType grammage = showerAxis_.getProjectedX(point);
+    int const maxBin = int(profile_.size() - 1);
+    int bin = grammage / dX_;
+    if (bin < 0) bin = 0;
+    if (bin > maxBin) bin = maxBin;
+
+    CORSIKA_LOGGER_TRACE(TOutput::getLogger(), "add local energy loss bin={} dE={} GeV ",
+                         bin, dE / 1_GeV);
+
+    profile_[bin][static_cast<int>(dEdX_output::ProfileIndex::Total)] += dE;
+  }
+
+  template <typename TOutput>
+  inline void EnergyLossWriter<TOutput>::write(GrammageType const Xstart,
+                                               GrammageType const Xend, Code const,
+                                               HEPEnergyType const dE) {
+    double const bstart = Xstart / dX_;
+    double const bend = Xend / dX_;
+
+    if (abs(bstart - floor(bstart + 0.5)) > 1e-2 ||
+        abs(bend - floor(bend + 0.5)) > 1e-2 || abs(bend - bstart - 1) > 1e-2) {
+      CORSIKA_LOGGER_ERROR(
+          TOutput::getLogger(),
+          "CascadeEquation (CONEX) and Corsika8 dX grammage binning are not the same! "
+          "Xstart={} Xend={} dX={} g/cm2",
+          Xstart / 1_g * square(1_cm), Xend / 1_g * square(1_cm),
+          dX_ / 1_g * square(1_cm));
+      throw std::runtime_error(
+          "CONEX and Corsika8 dX grammage binning are not the same!");
+    }
+
+    size_t const bin = size_t((bend + bstart) / 2);
+    CORSIKA_LOGGER_TRACE(TOutput::getLogger(),
+                         "add binned energy loss {} {} bin={} dE={} GeV ", bstart, bend,
+                         bin, dE / 1_GeV);
+    if (bin >= profile_.size()) {
+      CORSIKA_LOGGER_WARN(TOutput::getLogger(),
+                          "Grammage bin {} outside of profile {}. skipping.", bin,
+                          profile_.size());
+      return;
+    }
+    profile_[bin][static_cast<int>(dEdX_output::ProfileIndex::Total)] += dE;
+  }
+
+  template <typename TOutput>
+  inline HEPEnergyType EnergyLossWriter<TOutput>::getEnergyLost() const {
+    HEPEnergyType tot = HEPEnergyType::zero();
+    for (dEdX_output::Profile const& row : profile_)
+      tot += row.at(static_cast<int>(dEdX_output::ProfileIndex::Total));
+    return tot;
+  }
+
+  template <typename TOutput>
+  inline YAML::Node EnergyLossWriter<TOutput>::getConfig() const {
+
+    YAML::Node node;
+
+    node["type"] = "EnergyLoss";
+    node["units"]["energy"] = "GeV";
+    node["units"]["grammage"] = "g/cm^2";
+    node["bin-size"] = dX_ / (1_g / square(1_cm));
+    node["nbins"] = nBins_;
+    node["grammage_threshold"] = dX_threshold_ / (1_g / square(1_cm));
+
+    return node;
+  }
+
+  template <typename TOutput>
+  inline YAML::Node EnergyLossWriter<TOutput>::getSummary() const {
+
+    // determined Xmax and dEdXmax from quadratic interpolation
+    double maximum = 0;
+    size_t iMaximum = 0;
+    for (size_t i = 0; i < profile_.size() - 3; ++i) {
+      double value =
+          (profile_[i + 0].at(static_cast<int>(dEdX_output::ProfileIndex::Total)) +
+           profile_[i + 1].at(static_cast<int>(dEdX_output::ProfileIndex::Total)) +
+           profile_[i + 2].at(static_cast<int>(dEdX_output::ProfileIndex::Total))) /
+          1_GeV;
+      if (value > maximum) {
+        maximum = value;
+        iMaximum = i;
+      }
+    }
+
+    double const dX = dX_ / 1_g * square(1_cm);
+
+    auto [Xmax, dEdXmax] = FindXmax::interpolateProfile(
+        dX * (0.5 + iMaximum), dX * (1.5 + iMaximum), dX * (2.5 + iMaximum),
+        profile_[iMaximum + 0].at(static_cast<int>(dEdX_output::ProfileIndex::Total)) /
+            1_GeV,
+        profile_[iMaximum + 1].at(static_cast<int>(dEdX_output::ProfileIndex::Total)) /
+            1_GeV,
+        profile_[iMaximum + 2].at(static_cast<int>(dEdX_output::ProfileIndex::Total)) /
+            1_GeV);
+
+    YAML::Node summary;
+    summary["sum_dEdX"] = getEnergyLost() / 1_GeV;
+    summary["Xmax"] = Xmax;
+    summary["dEdXmax"] = dEdXmax;
+    return summary;
+  }
+
+} // namespace corsika
\ No newline at end of file

corsika/detail/modules/writers/EnergyLossWriterParquet.inl

+/*
+ * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/utility/FindXmax.hpp>
+
+#include <corsika/media/ShowerAxis.hpp>
+
+#include <exception>
+
+namespace corsika {
+
+  template <size_t NColumns>
+  inline EnergyLossWriterParquet<NColumns>::EnergyLossWriterParquet(
+      std::array<const char*, NColumns> const& colNames)
+      : columns_(colNames) {}
+
+  template <size_t NColumns>
+  inline void EnergyLossWriterParquet<NColumns>::startOfLibrary(
+      boost::filesystem::path const& directory) {
+
+    // setup the streamer
+    output_.initStreamer((directory / "dEdX.parquet").string());
+
+    // enable compression with the default level
+    output_.enableCompression();
+
+    // build the schema
+    output_.addField("X", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
+                     parquet::ConvertedType::NONE);
+    for (auto const& col : columns_) {
+      output_.addField(col, parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
+                       parquet::ConvertedType::NONE);
+    }
+
+    // and build the streamer
+    output_.buildStreamer();
+  }
+
+  template <size_t NColumns>
+  inline void EnergyLossWriterParquet<NColumns>::write(
+      unsigned int const showerId, GrammageType const grammage,
+      std::array<HEPEnergyType, NColumns> const& data) {
+
+    // and write the data into the column
+    *(output_.getWriter()) << showerId
+                           << static_cast<float>(grammage / 1_g * square(1_cm));
+    for (HEPEnergyType const& dedx : data) {
+      *(output_.getWriter()) << static_cast<float>(dedx / 1_GeV);
+    }
+    *(output_.getWriter()) << parquet::EndRow;
+  }
+
+  template <size_t NColumns>
+  inline void EnergyLossWriterParquet<NColumns>::startOfShower(unsigned int const) {}
+
+  template <size_t NColumns>
+  inline void EnergyLossWriterParquet<NColumns>::endOfShower(unsigned int const) {}
+
+  template <size_t NColumns>
+  inline void EnergyLossWriterParquet<NColumns>::endOfLibrary() {
+    output_.closeStreamer();
+  }
+
+} // namespace corsika

corsika/detail/modules/writers/InteractionWriter.inl

+/*
+ * (c) Copyright 2024 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/core/Logging.hpp>
+
+namespace corsika {
+
+  template <typename TTracking, typename TOutput>
+  inline InteractionWriter<TTracking, TOutput>::InteractionWriter(
+      ShowerAxis const& axis, ObservationPlane<TTracking, TOutput> const& obsPlane)
+      : obsPlane_(obsPlane)
+      , showerAxis_(axis)
+      , interactionCounter_(0)
+      , showerId_(0) {}
+
+  template <typename TTracking, typename TOutput>
+  template <typename TStackView>
+  inline void InteractionWriter<TTracking, TOutput>::doSecondaries(TStackView& vS) {
+    // Dumps the stack to the output parquet stream
+    // The primary and secondaries are all written along with the slant depth
+
+    if (interactionCounter_++) { return; } // Only run on the first interaction
+
+    auto primary = vS.getProjectile();
+    auto dX = showerAxis_.getProjectedX(primary.getPosition());
+    CORSIKA_LOG_INFO("First interaction at dX {}", dX);
+    CORSIKA_LOG_INFO("Primary: {}, E_kin {}", primary.getPID(),
+                     primary.getKineticEnergy());
+
+    // get the location of the primary w.r.t. observation plane
+    Vector const displacement = primary.getPosition() - obsPlane_.getPlane().getCenter();
+
+    auto const x = displacement.dot(obsPlane_.getXAxis());
+    auto const y = displacement.dot(obsPlane_.getYAxis());
+    auto const z = displacement.dot(obsPlane_.getPlane().getNormal());
+    auto const nx = primary.getDirection().dot(obsPlane_.getXAxis());
+    auto const ny = primary.getDirection().dot(obsPlane_.getYAxis());
+    auto const nz = primary.getDirection().dot(obsPlane_.getPlane().getNormal());
+
+    auto const px = primary.getMomentum().dot(obsPlane_.getXAxis());
+    auto const py = primary.getMomentum().dot(obsPlane_.getYAxis());
+    auto const pz = primary.getMomentum().dot(obsPlane_.getPlane().getNormal());
+
+    summary_["shower_" + std::to_string(showerId_)]["pdg"] =
+        static_cast<int>(get_PDG(primary.getPID()));
+    summary_["shower_" + std::to_string(showerId_)]["name"] =
+        static_cast<std::string>(get_name(primary.getPID()));
+    summary_["shower_" + std::to_string(showerId_)]["total_energy"] =
+        (primary.getKineticEnergy() + get_mass(primary.getPID())) / 1_GeV;
+    summary_["shower_" + std::to_string(showerId_)]["kinetic_energy"] =
+        primary.getKineticEnergy() / 1_GeV;
+    summary_["shower_" + std::to_string(showerId_)]["x"] = x / 1_m;
+    summary_["shower_" + std::to_string(showerId_)]["y"] = y / 1_m;
+    summary_["shower_" + std::to_string(showerId_)]["z"] = z / 1_m;
+    summary_["shower_" + std::to_string(showerId_)]["nx"] = static_cast<double>(nx);
+    summary_["shower_" + std::to_string(showerId_)]["ny"] = static_cast<double>(ny);
+    summary_["shower_" + std::to_string(showerId_)]["nz"] = static_cast<double>(nz);
+    summary_["shower_" + std::to_string(showerId_)]["px"] =
+        static_cast<double>(px / 1_GeV);
+    summary_["shower_" + std::to_string(showerId_)]["py"] =
+        static_cast<double>(py / 1_GeV);
+    summary_["shower_" + std::to_string(showerId_)]["pz"] =
+        static_cast<double>(pz / 1_GeV);
+    summary_["shower_" + std::to_string(showerId_)]["time"] = primary.getTime() / 1_s;
+    summary_["shower_" + std::to_string(showerId_)]["slant_depth"] =
+        dX / (1_g / 1_cm / 1_cm);
+
+    uint nSecondaries = 0;
+
+    // Loop through secondaries
+    auto particle = vS.begin();
+    while (particle != vS.end()) {
+
+      // Get the momentum of the secondary, write w.r.t. observation plane
+      auto const p_2nd = particle.getMomentum();
+
+      *(output_.getWriter())
+          << showerId_ << static_cast<int>(get_PDG(particle.getPID()))
+          << static_cast<float>(p_2nd.dot(obsPlane_.getXAxis()) / 1_GeV)
+          << static_cast<float>(p_2nd.dot(obsPlane_.getYAxis()) / 1_GeV)
+          << static_cast<float>(p_2nd.dot(obsPlane_.getPlane().getNormal()) / 1_GeV)
+          << parquet::EndRow;
+
+      CORSIKA_LOG_INFO(" 2ndary: {}, E_kin {}", particle.getPID(),
+                       particle.getKineticEnergy());
+      ++particle;
+      ++nSecondaries;
+    }
+
+    summary_["shower_" + std::to_string(showerId_)]["n_secondaries"] = nSecondaries;
+  }
+
+  template <typename TTracking, typename TOutput>
+  inline void InteractionWriter<TTracking, TOutput>::startOfLibrary(
+      boost::filesystem::path const& directory) {
+    output_.initStreamer((directory / ("interactions.parquet")).string());
+
+    // enable compression with the default level
+    output_.enableCompression();
+
+    output_.addField("pdg", parquet::Repetition::REQUIRED, parquet::Type::INT32,
+                     parquet::ConvertedType::INT_32);
+    output_.addField("px", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
+                     parquet::ConvertedType::NONE);
+    output_.addField("py", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
+                     parquet::ConvertedType::NONE);
+    output_.addField("pz", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
+                     parquet::ConvertedType::NONE);
+
+    output_.buildStreamer();
+
+    showerId_ = 0;
+    interactionCounter_ = 0;
+    summary_ = YAML::Node();
+  }
+
+  template <typename TTracking, typename TOutput>
+  inline void InteractionWriter<TTracking, TOutput>::startOfShower(
+      unsigned int const showerId) {
+    showerId_ = showerId;
+    interactionCounter_ = 0;
+  }
+
+  template <typename TTracking, typename TOutput>
+  inline void InteractionWriter<TTracking, TOutput>::endOfShower(unsigned int const) {}
+
+  template <typename TTracking, typename TOutput>
+  inline void InteractionWriter<TTracking, TOutput>::endOfLibrary() {
+    output_.closeStreamer();
+  }
+
+  template <typename TTracking, typename TOutput>
+  inline YAML::Node InteractionWriter<TTracking, TOutput>::getConfig() const {
+    YAML::Node node;
+    node["type"] = "Interactions";
+    node["units"]["energy"] = "GeV";
+    node["units"]["length"] = "m";
+    node["units"]["time"] = "ns";
+    node["units"]["grammage"] = "g/cm^2";
+    return node;
+  }
+
+  template <typename TTracking, typename TOutput>
+  inline YAML::Node InteractionWriter<TTracking, TOutput>::getSummary() const {
+    return summary_;
+  }
+
+} // namespace corsika

corsika/detail/modules/writers/LongitudinalProfileWriterParquet.inl

+/*
+ * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/utility/FindXmax.hpp>
+
+#include <corsika/media/ShowerAxis.hpp>
+
+#include <string>
+#include <exception>
+
+namespace corsika {
+
+  template <size_t NColumns>
+  inline LongitudinalProfileWriterParquet<NColumns>::LongitudinalProfileWriterParquet(
+      std::array<const char*, NColumns> const& columns)
+      : columns_(columns) {}
+
+  template <size_t NColumns>
+  inline void LongitudinalProfileWriterParquet<NColumns>::startOfLibrary(
+      boost::filesystem::path const& directory) {
+    // setup the streamer
+    output_.initStreamer((directory / "profile.parquet").string());
+
+    // enable compression with the default level
+    output_.enableCompression();
+
+    // build the schema
+    output_.addField("X", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
+                     parquet::ConvertedType::NONE);
+    for (auto const& col : columns_) {
+      output_.addField(col, parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
+                       parquet::ConvertedType::NONE);
+    }
+
+    // and build the streamer
+    output_.buildStreamer();
+  }
+
+  template <size_t NColumns>
+  inline void LongitudinalProfileWriterParquet<NColumns>::startOfShower(
+      unsigned int const) {}
+
+  template <size_t NColumns>
+  inline void LongitudinalProfileWriterParquet<NColumns>::endOfShower(
+      unsigned int const) {}
+
+  template <size_t NColumns>
+  inline void LongitudinalProfileWriterParquet<NColumns>::endOfLibrary() {
+    output_.closeStreamer();
+  }
+
+  template <size_t NColumns>
+  inline void LongitudinalProfileWriterParquet<NColumns>::write(
+      unsigned int const showerId, GrammageType const grammage,
+      std::array<double, NColumns> const& data) {
+
+    // and write the data into the column
+    *(output_.getWriter()) << showerId
+                           << static_cast<float>(grammage / 1_g * square(1_cm));
+    for (double const weight : data) {
+      *(output_.getWriter()) << static_cast<float>(weight);
+    }
+    *(output_.getWriter()) << parquet::EndRow;
+  }
+
+} // namespace corsika

corsika/detail/modules/writers/LongitudinalWriter.inl

+/*
+ * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/utility/FindXmax.hpp>
+
+#include <corsika/media/ShowerAxis.hpp>
+
+#include <exception>
+#include <algorithm>
+#include <iostream>
+
+namespace corsika {
+
+  template <typename TOutput>
+  inline LongitudinalWriter<TOutput>::LongitudinalWriter(ShowerAxis const& axis,
+                                                         GrammageType dX)
+      : LongitudinalWriter<TOutput>{
+            axis, static_cast<unsigned int>(axis.getMaximumX() / dX) + 1, dX} {}
+
+  template <typename TOutput>
+  inline LongitudinalWriter<TOutput>::LongitudinalWriter(ShowerAxis const& axis,
+                                                         size_t nbins, GrammageType dX)
+      : TOutput(number_profile::ProfileIndexNames)
+      , showerAxis_(axis)
+      , dX_(dX)
+      , nBins_(nbins)
+      , profile_{nbins} {}
+
+  template <typename TOutput>
+  inline void LongitudinalWriter<TOutput>::startOfLibrary(
+      boost::filesystem::path const& directory) {
+    TOutput::startOfLibrary(directory);
+  }
+
+  template <typename TOutput>
+  inline void LongitudinalWriter<TOutput>::startOfShower(unsigned int const showerId) {
+    profile_.clear();
+    for (size_t i = 0; i < nBins_; ++i) { profile_.emplace_back(); }
+    TOutput::startOfShower(showerId);
+  }
+
+  template <typename TOutput>
+  inline void LongitudinalWriter<TOutput>::endOfLibrary() {
+    TOutput::endOfLibrary();
+  }
+
+  template <typename TOutput>
+  inline void LongitudinalWriter<TOutput>::endOfShower(unsigned int const showerId) {
+
+    int iRow{0};
+    for (number_profile::ProfileData const& row : profile_) {
+      // here: write to underlying writer (e.g. parquet)
+      TOutput::write(showerId, iRow * dX_, row);
+      iRow++;
+    }
+    TOutput::endOfShower(showerId);
+  }
+
+  template <typename TOutput>
+  inline void LongitudinalWriter<TOutput>::write(Point const& p0, Point const& p1,
+                                                 Code const pid, double const weight) {
+    GrammageType const grammageStart = showerAxis_.getProjectedX(p0);
+    GrammageType const grammageEnd = showerAxis_.getProjectedX(p1);
+
+    // Avoid over counting in first bin when backscattered particle goes beyond the
+    // injection point.
+    if (grammageStart == grammageEnd) { return; }
+
+    // Note: particle may go also "upward", thus, grammageEnd<grammageStart
+    size_t const binStart = std::ceil(grammageStart / dX_);
+    size_t const binEnd = std::floor(grammageEnd / dX_);
+
+    CORSIKA_LOGGER_TRACE(TOutput::getLogger(),
+                         "grammageStart={} End={} binStart={}, end={}",
+                         grammageStart / 1_g * square(1_cm),
+                         grammageEnd / 1_g * square(1_cm), binStart, binEnd);
+
+    for (size_t bin = binStart; bin <= std::min(binEnd, profile_.size() - 1); ++bin) {
+      if (pid == Code::Photon) {
+        profile_.at(bin)[static_cast<int>(number_profile::ProfileIndex::Photon)] +=
+            weight;
+      } else if (pid == Code::Positron) {
+        profile_.at(bin)[static_cast<int>(number_profile::ProfileIndex::Positron)] +=
+            weight;
+      } else if (pid == Code::Electron) {
+        profile_.at(bin)[static_cast<int>(number_profile::ProfileIndex::Electron)] +=
+            weight;
+      } else if (pid == Code::MuPlus) {
+        profile_.at(bin)[static_cast<int>(number_profile::ProfileIndex::MuPlus)] +=
+            weight;
+      } else if (pid == Code::MuMinus) {
+        profile_.at(bin)[static_cast<int>(number_profile::ProfileIndex::MuMinus)] +=
+            weight;
+      } else if (is_hadron(pid)) {
+        profile_.at(bin)[static_cast<int>(number_profile::ProfileIndex::Hadron)] +=
+            weight;
+      }
+      if (is_charged(pid)) {
+        profile_[bin][static_cast<int>(number_profile::ProfileIndex::Charged)] += weight;
+      }
+    }
+  }
+
+  template <typename TOutput>
+  inline void LongitudinalWriter<TOutput>::write(GrammageType const Xstart,
+                                                 GrammageType const Xend, Code const pid,
+                                                 double const weight) {
+    double const bstart = Xstart / dX_;
+    double const bend = Xend / dX_;
+
+    if (abs(bstart - floor(bstart + 0.5)) > 1e-2 ||
+        abs(bend - floor(bend + 0.5)) > 1e-2 || abs(bend - bstart - 1) > 1e-2) {
+      CORSIKA_LOGGER_ERROR(TOutput::getLogger(),
+                           "CONEX and Corsika8 dX grammage binning are not the same! "
+                           "Xstart={} Xend={} dX={}",
+                           Xstart / 1_g * square(1_cm), Xend / 1_g * square(1_cm),
+                           dX_ / 1_g * square(1_cm));
+      throw std::runtime_error(
+          "CONEX and Corsika8 dX grammage binning are not the same!");
+    }
+
+    size_t const bin = size_t((bend + bstart) / 2);
+    if (bin >= profile_.size()) {
+      CORSIKA_LOGGER_WARN(TOutput::getLogger(),
+                          "Grammage bin {} outside of profile {}. skipping.", bin,
+                          profile_.size());
+      return;
+    }
+
+    if (pid == Code::Photon) {
+      profile_.at(bin)[static_cast<int>(number_profile::ProfileIndex::Photon)] += weight;
+    } else if (pid == Code::Positron) {
+      profile_.at(bin)[static_cast<int>(number_profile::ProfileIndex::Positron)] +=
+          weight;
+    } else if (pid == Code::Electron) {
+      profile_.at(bin)[static_cast<int>(number_profile::ProfileIndex::Electron)] +=
+          weight;
+    } else if (pid == Code::MuPlus) {
+      profile_.at(bin)[static_cast<int>(number_profile::ProfileIndex::MuPlus)] += weight;
+    } else if (pid == Code::MuMinus) {
+      profile_.at(bin)[static_cast<int>(number_profile::ProfileIndex::MuMinus)] += weight;
+    } else if (is_hadron(pid)) {
+      profile_.at(bin)[static_cast<int>(number_profile::ProfileIndex::Hadron)] += weight;
+    }
+    if (is_charged(pid)) {
+      profile_[bin][static_cast<int>(number_profile::ProfileIndex::Charged)] += weight;
+    }
+  }
+
+  template <typename TOutput>
+  inline YAML::Node LongitudinalWriter<TOutput>::getConfig() const {
+
+    YAML::Node node;
+
+    node["type"] = "LongitudinalProfile";
+    node["units"]["grammage"] = "g/cm^2";
+    node["bin-size"] = dX_ / (1_g / square(1_cm));
+    node["nbins"] = nBins_;
+
+    return node;
+  }
+
+  template <typename TOutput>
+  inline YAML::Node LongitudinalWriter<TOutput>::getSummary() const {
+    YAML::Node summary;
+    return summary;
+  }
+
+} // namespace corsika
\ No newline at end of file

corsika/detail/modules/writers/ParticleWriterParquet.inl

+/*
+ * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/core/Logging.hpp>
+
+namespace corsika {
+
+  inline ParticleWriterParquet::ParticleWriterParquet(bool const printZ)
+      : output_()
+      , showerId_(0)
+      , printZ_(printZ) {}
+
+  inline void ParticleWriterParquet::startOfLibrary(
+      boost::filesystem::path const& directory) {
+
+    // setup the streamer
+    output_.initStreamer((directory / "particles.parquet").string());
+
+    // enable compression with the default level
+    output_.enableCompression();
+
+    // build the schema
+    output_.addField("pdg", parquet::Repetition::REQUIRED, parquet::Type::INT32,
+                     parquet::ConvertedType::INT_32);
+    output_.addField("kinetic_energy", parquet::Repetition::REQUIRED,
+                     parquet::Type::FLOAT, parquet::ConvertedType::NONE);
+    output_.addField("x", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
+                     parquet::ConvertedType::NONE);
+    output_.addField("y", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
+                     parquet::ConvertedType::NONE);
+    if (printZ_) {
+      output_.addField("z", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
+                       parquet::ConvertedType::NONE);
+    }
+    output_.addField("nx", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
+                     parquet::ConvertedType::NONE);
+    output_.addField("ny", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
+                     parquet::ConvertedType::NONE);
+    output_.addField("nz", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
+                     parquet::ConvertedType::NONE);
+    output_.addField("time", parquet::Repetition::REQUIRED, parquet::Type::DOUBLE,
+                     parquet::ConvertedType::NONE);
+    output_.addField("weight", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
+                     parquet::ConvertedType::NONE);
+
+    // and build the streamer
+    output_.buildStreamer();
+
+    showerId_ = 0;
+  }
+
+  inline void ParticleWriterParquet::startOfShower(unsigned int const showerId) {
+    showerId_ = showerId;
+
+    countHadrons_ = 0;
+    countOthers_ = 0;
+    countEM_ = 0;
+    countMuons_ = 0;
+
+    kineticEnergyHadrons_ = 0_eV;
+    kineticEnergyMuons_ = 0_eV;
+    kineticEnergyEM_ = 0_eV;
+    kineticEnergyOthers_ = 0_eV;
+
+    totalEnergyHadrons_ = 0_eV;
+    totalEnergyMuons_ = 0_eV;
+    totalEnergyEM_ = 0_eV;
+    totalEnergyOthers_ = 0_eV;
+  }
+
+  inline void ParticleWriterParquet::endOfShower(unsigned int const) {
+    summary_["shower_" + std::to_string(showerId_)]["hadron"]["count"] = countHadrons_;
+    summary_["shower_" + std::to_string(showerId_)]["hadron"]["kinetic_energy"] =
+        kineticEnergyHadrons_ / 1_GeV;
+    summary_["shower_" + std::to_string(showerId_)]["hadron"]["total_energy"] =
+        totalEnergyHadrons_ / 1_GeV;
+
+    summary_["shower_" + std::to_string(showerId_)]["muon"]["count"] = countMuons_;
+    summary_["shower_" + std::to_string(showerId_)]["muon"]["kinetic_energy"] =
+        kineticEnergyMuons_ / 1_GeV;
+    summary_["shower_" + std::to_string(showerId_)]["muon"]["total_energy"] =
+        totalEnergyMuons_ / 1_GeV;
+
+    summary_["shower_" + std::to_string(showerId_)]["em"]["count"] = countEM_;
+    summary_["shower_" + std::to_string(showerId_)]["em"]["kinetic_energy"] =
+        kineticEnergyEM_ / 1_GeV;
+    summary_["shower_" + std::to_string(showerId_)]["em"]["total_energy"] =
+        totalEnergyEM_ / 1_GeV;
+
+    summary_["shower_" + std::to_string(showerId_)]["other"]["count"] = countOthers_;
+    summary_["shower_" + std::to_string(showerId_)]["other"]["kinetic_energy"] =
+        kineticEnergyOthers_ / 1_GeV;
+    summary_["shower_" + std::to_string(showerId_)]["other"]["total_energy"] =
+        totalEnergyOthers_ / 1_GeV;
+  }
+
+  inline void ParticleWriterParquet::endOfLibrary() { output_.closeStreamer(); }
+
+  inline HEPEnergyType ParticleWriterParquet::getEnergyGround() const {
+    return totalEnergyHadrons_ + totalEnergyMuons_ + totalEnergyEM_ + totalEnergyOthers_;
+  }
+
+  inline void ParticleWriterParquet::write(Code const pid,
+                                           HEPEnergyType const kineticEnergy,
+                                           LengthType const x, LengthType const y,
+                                           LengthType const z, double const nx,
+                                           double const ny, double const nz,
+                                           TimeType const t, double const weight) {
+
+    // write the next row - we must write `shower_` first.
+    *(output_.getWriter()) << showerId_ << static_cast<int>(get_PDG(pid))
+                           << static_cast<float>(kineticEnergy / 1_GeV)
+                           << static_cast<float>(x / 1_m) << static_cast<float>(y / 1_m);
+    if (printZ_) { *(output_.getWriter()) << static_cast<float>(z / 1_m); }
+
+    *(output_.getWriter()) << static_cast<float>(nx) << static_cast<float>(ny)
+                           << static_cast<float>(nz) << static_cast<double>(t / 1_s)
+                           << static_cast<float>(weight) << parquet::EndRow;
+
+    if (is_hadron(pid)) {
+      ++countHadrons_;
+      kineticEnergyHadrons_ += kineticEnergy;
+      totalEnergyHadrons_ += kineticEnergy + get_mass(pid);
+    } else if (is_muon(pid)) {
+      ++countMuons_;
+      kineticEnergyMuons_ += kineticEnergy;
+      totalEnergyMuons_ += kineticEnergy + get_mass(pid);
+    } else if (is_em(pid)) {
+      ++countEM_;
+      kineticEnergyEM_ += kineticEnergy;
+      totalEnergyEM_ += kineticEnergy + get_mass(pid);
+    } else {
+      ++countOthers_;
+      kineticEnergyOthers_ += kineticEnergy;
+      totalEnergyOthers_ += kineticEnergy + get_mass(pid);
+    }
+  }
+
+  /**
+   * Return collected library-level summary for output.
+   */
+  inline YAML::Node ParticleWriterParquet::getSummary() const { return summary_; }
+
+} // namespace corsika

corsika/detail/modules/writers/PrimaryWriter.inl

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/core/ParticleProperties.hpp>
+
+namespace corsika {
+
+  template <typename TTracking, typename TOutput>
+  inline PrimaryWriter<TTracking, TOutput>::PrimaryWriter(
+      ObservationPlane<TTracking, TOutput> const& obsPlane)
+      : obsPlane_(obsPlane)
+      , writtenThisShower_(false)
+      , showerId_(0) {}
+
+  template <typename TTracking, typename TOutput>
+  inline PrimaryWriter<TTracking, TOutput>::~PrimaryWriter() {}
+
+  template <typename TTracking, typename TOutput>
+  inline void PrimaryWriter<TTracking, TOutput>::recordPrimary(
+      std::tuple<Code, HEPEnergyType, DirectionVector, Point, TimeType> const&
+          primaryProps) {
+    auto const [pID, kineticEnergy, direction, position, time] = primaryProps;
+
+    Vector const displacement = position - obsPlane_.getPlane().getCenter();
+
+    auto const x = displacement.dot(obsPlane_.getXAxis());
+    auto const y = displacement.dot(obsPlane_.getYAxis());
+    auto const z = displacement.dot(obsPlane_.getPlane().getNormal());
+    auto const nx = direction.dot(obsPlane_.getXAxis());
+    auto const ny = direction.dot(obsPlane_.getYAxis());
+    auto const nz = direction.dot(obsPlane_.getPlane().getNormal());
+
+    output_["shower_" + std::to_string(showerId_)]["pdg"] =
+        static_cast<int>(get_PDG(pID));
+    output_["shower_" + std::to_string(showerId_)]["name"] =
+        static_cast<std::string>(get_name(pID));
+    output_["shower_" + std::to_string(showerId_)]["total_energy"] =
+        (kineticEnergy + get_mass(pID)) / 1_GeV;
+    output_["shower_" + std::to_string(showerId_)]["kinetic_energy"] =
+        kineticEnergy / 1_GeV;
+    output_["shower_" + std::to_string(showerId_)]["x"] = x / 1_m;
+    output_["shower_" + std::to_string(showerId_)]["y"] = y / 1_m;
+    output_["shower_" + std::to_string(showerId_)]["z"] = z / 1_m;
+    output_["shower_" + std::to_string(showerId_)]["nx"] = static_cast<double>(nx);
+    output_["shower_" + std::to_string(showerId_)]["ny"] = static_cast<double>(ny);
+    output_["shower_" + std::to_string(showerId_)]["nz"] = static_cast<double>(nz);
+    output_["shower_" + std::to_string(showerId_)]["time"] = time / 1_s;
+  }
+
+  template <typename TTracking, typename TOutput>
+  inline void PrimaryWriter<TTracking, TOutput>::startOfLibrary([
+      [maybe_unused]] boost::filesystem::path const& directory) {
+    output_ = YAML::Node();
+    showerId_ = 0;
+  }
+
+  template <typename TTracking, typename TOutput>
+  inline void PrimaryWriter<TTracking, TOutput>::endOfShower([
+      [maybe_unused]] unsigned int const showerId) {
+    ++showerId_;
+  }
+
+  template <typename TTracking, typename TOutput>
+  inline void PrimaryWriter<TTracking, TOutput>::endOfLibrary() {}
+
+  template <typename TTracking, typename TOutput>
+  inline YAML::Node PrimaryWriter<TTracking, TOutput>::getConfig() const {
+    YAML::Node node;
+    node["type"] = "PrimaryParticle";
+
+    node["units"]["length"] = "m"; // add default units for values
+    node["units"]["energy"] = "GeV";
+    node["units"]["time"] = "s";
+
+    // the center of the plane
+    auto const center{obsPlane_.getPlane().getCenter()};
+
+    // save each component in its native coordinate system
+    auto const center_coords{center.getCoordinates(center.getCoordinateSystem())};
+    node["plane"]["center"].push_back(center_coords.getX() / 1_m);
+    node["plane"]["center"].push_back(center_coords.getY() / 1_m);
+    node["plane"]["center"].push_back(center_coords.getZ() / 1_m);
+
+    auto const normal{obsPlane_.getPlane().getNormal().getComponents()};
+    node["plane"]["normal"].push_back(normal.getX().magnitude());
+    node["plane"]["normal"].push_back(normal.getY().magnitude());
+    node["plane"]["normal"].push_back(normal.getZ().magnitude());
+
+    auto const xAxis_coords{
+        obsPlane_.getXAxis().getComponents(obsPlane_.getXAxis().getCoordinateSystem())};
+    node["x-axis"].push_back(xAxis_coords.getX().magnitude());
+    node["x-axis"].push_back(xAxis_coords.getY().magnitude());
+    node["x-axis"].push_back(xAxis_coords.getZ().magnitude());
+
+    auto const yAxis_coords{
+        obsPlane_.getYAxis().getComponents(obsPlane_.getYAxis().getCoordinateSystem())};
+    node["y-axis"].push_back(yAxis_coords.getX().magnitude());
+    node["y-axis"].push_back(yAxis_coords.getY().magnitude());
+    node["y-axis"].push_back(yAxis_coords.getZ().magnitude());
+
+    return node;
+  }
+
+  template <typename TTracking, typename TOutput>
+  inline YAML::Node PrimaryWriter<TTracking, TOutput>::getSummary() const {
+    return output_;
+  }
+
+} // namespace corsika

corsika/detail/modules/writers/TrackWriterParquet.inl

+/*
+ * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+namespace corsika {
+
+  inline TrackWriterParquet::TrackWriterParquet()
+      : output_()
+      , showerId_(0) {}
+
+  inline void TrackWriterParquet::startOfLibrary(
+      boost::filesystem::path const& directory) {
+
+    // setup the streamer
+    output_.initStreamer((directory / "tracks.parquet").string());
+
+    // enable compression with the default level
+    output_.enableCompression();
+
+    // build the schema
+    output_.addField("pdg", parquet::Repetition::REQUIRED, parquet::Type::INT32,
+                     parquet::ConvertedType::INT_32);
+    output_.addField("start_energy", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
+                     parquet::ConvertedType::NONE);
+    output_.addField("weight", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
+                     parquet::ConvertedType::NONE);
+    output_.addField("start_x", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
+                     parquet::ConvertedType::NONE);
+    output_.addField("start_y", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
+                     parquet::ConvertedType::NONE);
+    output_.addField("start_z", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
+                     parquet::ConvertedType::NONE);
+    output_.addField("start_t", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
+                     parquet::ConvertedType::NONE);
+    output_.addField("end_x", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
+                     parquet::ConvertedType::NONE);
+    output_.addField("end_y", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
+                     parquet::ConvertedType::NONE);
+    output_.addField("end_z", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
+                     parquet::ConvertedType::NONE);
+    output_.addField("end_energy", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
+                     parquet::ConvertedType::NONE);
+    output_.addField("end_t", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
+                     parquet::ConvertedType::NONE);
+
+    // and build the streamer
+    output_.buildStreamer();
+    showerId_ = 0;
+  }
+
+  inline void TrackWriterParquet::startOfShower(unsigned int const showerId) {
+    showerId_ = showerId;
+  }
+
+  inline void TrackWriterParquet::endOfShower(unsigned int const) {}
+
+  inline void TrackWriterParquet::endOfLibrary() { output_.closeStreamer(); }
+
+  inline void TrackWriterParquet::write(Code const pid, HEPEnergyType const KinenergyPre,
+                                        double const weight,
+                                        QuantityVector<length_d> const& start,
+                                        TimeType const t_start,
+                                        QuantityVector<length_d> const& end,
+                                        HEPEnergyType const KinenergyPost,
+                                        TimeType const t_end) {
+
+    // write the next row - we must write `shower_` first.
+    // clang-format off
+   *(output_.getWriter())
+       << showerId_
+       << static_cast<int>(get_PDG(pid))
+       << static_cast<float>(KinenergyPre / 1_GeV)
+       << static_cast<float>(weight)
+       << static_cast<float>(start[0] / 1_m)
+       << static_cast<float>(start[1] / 1_m)
+       << static_cast<float>(start[2] / 1_m)
+       << static_cast<float>(t_start / 1_ns)
+       << static_cast<float>(end[0] / 1_m)
+       << static_cast<float>(end[1] / 1_m)
+       << static_cast<float>(end[2] / 1_m)
+       << static_cast<float>(KinenergyPost / 1_GeV)
+       << static_cast<float>(t_end / 1_ns)
+       << parquet::EndRow;
+    // clang-format on
+  }
+
+} // namespace corsika
\ No newline at end of file

corsika/detail/output/DummyOutputManager.inl

+/*
+ * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+#pragma once
+
+namespace corsika {
+
+  inline DummyOutputManager::DummyOutputManager() {}
+
+  inline DummyOutputManager::~DummyOutputManager() {}
+
+  template <typename TOutput>
+  inline void DummyOutputManager::add(std::string const&, TOutput&) {}
+
+  inline void DummyOutputManager::startOfLibrary() {}
+
+  inline void DummyOutputManager::startOfShower() {}
+
+  inline void DummyOutputManager::endOfShower() {}
+
+  inline void DummyOutputManager::endOfLibrary() {}
+
+} // namespace corsika

corsika/detail/output/OutputManager.inl

+/*
+ * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+#pragma once
+
+#include <algorithm>
+#include <fstream>
+#include <functional>
+
+#include <iomanip>
+#include <ctime>
+#include <sstream>
+
+#include <boost/filesystem.hpp>
+
+#include <fmt/core.h>
+#include <fmt/chrono.h>
+
+namespace corsika {
+
+  inline OutputManager::OutputManager(std::string const& dir_path, const long& vseed = 0,
+                                      std::string const& vargs = "",
+                                      bool useCompression = false)
+      : root_(dir_path)
+      , cmnd_line_args_(vargs)
+      , useCompression_(useCompression)
+      , count_(0)
+      , seed_(vseed) {
+
+    // check if this directory already exists
+    if (boost::filesystem::exists(root_)) {
+      CORSIKA_LOGGER_ERROR(logger_,
+                           "Output directory '{}' already exists! Do not overwrite!.",
+                           root_.string());
+      throw std::runtime_error("Output directory already exists.");
+    }
+
+    // construct the directory for this library
+    boost::filesystem::create_directories(root_);
+
+    CORSIKA_LOGGER_INFO(logger_, fmt::format("Output library: \"{}\"", root_.string()));
+    writeYAML(getConfig(), root_ / ("config.yaml"));
+  }
+
+  template <typename TOutput>
+  inline void OutputManager::add(std::string const& name, TOutput& output) {
+
+    if (state_ != OutputState::NoInit) {
+      // if "add" is called after the ouptput has started, this is an ERROR.
+      CORSIKA_LOGGER_ERROR(
+          logger_, "Cannot add more outputs to OutputManager after output was started.");
+      throw std::runtime_error(
+          "Cannot add more outputs to OutputManager after output was started.");
+    }
+
+    // check if that name is already in the map
+    if (outputs_.count(name) > 0) {
+      CORSIKA_LOGGER_ERROR(
+          logger_, "'{}' is already registered. All outputs must have unique names.",
+          name);
+      throw std::runtime_error("Output already exists. Do not overwrite!");
+    }
+
+    // if we get here, the name is not already in the map
+    // so we create the output and register it into the map
+    outputs_.insert(std::make_pair(name, std::ref(output)));
+
+    // create the directory for this process.
+    boost::filesystem::create_directory(root_ / name);
+  }
+
+  inline OutputManager::~OutputManager() {
+
+    if (state_ == OutputState::ShowerInProgress) {
+      // if this the destructor is called before the shower has been explicitly
+      // ended, print a warning and end the shower before continuing.
+      CORSIKA_LOGGER_WARN(logger_,
+                          "OutputManager was destroyed before endOfShower() called."
+                          " The last shower in this libray may be incomplete.");
+      endOfShower();
+    }
+
+    // write the top level summary file (summary.yaml)
+    writeSummary();
+
+    // if we are being destructed but EndOfLibrary() has not been called,
+    // make sure that we gracefully close all the outputs. This is a supported
+    // method of operation so we don't issue a warning here
+    if (state_ == OutputState::LibraryReady) { endOfLibrary(); }
+
+    if (useCompression_) {
+      auto const parent = root_.parent_path();
+      auto const final_part = root_.filename();
+
+      // Ensure we use "./" if the parent path is empty
+      std::string parent_path = parent.empty() ? "./" : parent.string() + "/";
+
+      std::string const cmd = "tar cf " + parent_path + final_part.string() + ".tar -C " +
+                              parent_path + " " + final_part.string();
+      CORSIKA_LOG_INFO("Compressing output directory using: {}", cmd.c_str());
+      int const returnCode = std::system(cmd.c_str());
+
+      if (returnCode) {
+        CORSIKA_LOG_ERROR("Compression returned with error code {}", returnCode);
+      } else {
+        // remove the original directory
+        boost::filesystem::remove_all(root_);
+      }
+    }
+  }
+
+  inline int OutputManager::getEventId() const { return count_; }
+
+  inline YAML::Node OutputManager::getConfig() const {
+
+    YAML::Node config;
+
+    // some basic info
+    config["path"] = root_.string();      // the simulation name
+    config["creator"] = "CORSIKA8";       // a tag to identify C8 libraries
+    config["version"] = "8.0.0-prealpha"; // the current version
+    config["args"] = cmnd_line_args_;     // the command line parameters
+    return config;
+  }
+
+  inline YAML::Node OutputManager::getSummary() const {
+
+    YAML::Node summary;
+
+    // the total number of showers contained in the library
+    summary["showers"] = count_;
+
+    summary["seed"] = seed_;
+
+    // this next section handles writing some time and duration information
+
+    // create a quick lambda function to convert a time-instance to a string
+    auto timeToString = [&](auto const time) -> std::string {
+      // ISO 8601 time format
+      auto format{"%FT%T%z"};
+
+      // convert the clock to a time_t
+      auto time_tc{std::chrono::system_clock::to_time_t(time)};
+
+      // create the string and push the time onto it
+      std::ostringstream oss;
+      oss << std::put_time(std::localtime(&time_tc), format);
+
+      return oss.str();
+    };
+
+    auto end_time{std::chrono::system_clock::now()};
+
+    // calculate the number of days
+    auto const start_t = std::chrono::system_clock::to_time_t(start_time);
+    auto const end_t = std::chrono::system_clock::to_time_t(end_time);
+    int const durationDays = std::difftime(end_t, start_t) / (60 * 60 * 24);
+
+    // add the time and duration info
+    summary["start time"] = timeToString(start_time);
+    summary["end time"] = timeToString(end_time);
+    summary["runtime"] = (durationDays ? fmt::format("+{}d ", durationDays) : "") +
+                         fmt::format("{:%H:%M:%S}", end_time - start_time);
+
+    std::vector<std::string> output_dirs;
+    for (auto const& outs : outputs_) { output_dirs.push_back(outs.first); }
+    summary["output_dirs"] = output_dirs;
+
+    return summary;
+  }
+
+  inline void OutputManager::writeSummary() const {
+
+    // write the node to a file
+    writeYAML(getSummary(), root_ / ("summary.yaml"));
+  }
+
+  inline void OutputManager::startOfLibrary() {
+
+    // this is only valid when we haven't started a library
+    // or have already finished a library
+    if (!(state_ == OutputState::NoInit || state_ == OutputState::LibraryFinished)) {
+
+      throw std::runtime_error("startOfLibrary() called in invalid state.");
+    }
+
+    // we now forward this signal to all of our outputs
+    for (auto& [name, output] : outputs_) {
+
+      // and start the library
+      output.get().startOfLibrary(root_ / name);
+
+      // get the config from this output
+      auto config = output.get().getConfig();
+
+      // add the name keyword
+      config["name"] = name;
+
+      // write the output configuration to config.yaml in the output directory
+      writeYAML(config, root_ / name / ("config.yaml"));
+    }
+
+    // we have now started running
+    state_ = OutputState::LibraryReady;
+    count_ = 0; // event counter
+  }
+
+  inline void OutputManager::startOfShower() {
+
+    // if this is called and we still in the "no init" state, then
+    // this is the first shower in the library so make sure we start it
+    if (state_ == OutputState::NoInit) { startOfLibrary(); }
+
+    // now start the event for all the outputs
+    for (auto& [name, output] : outputs_) { output.get().startOfShower(count_); }
+
+    // and transition to the in progress state
+    state_ = OutputState::ShowerInProgress;
+  }
+
+  inline void OutputManager::endOfShower() {
+
+    for (auto& [name, output] : outputs_) { output.get().endOfShower(count_); }
+
+    // switch back to the initialized state
+    state_ = OutputState::LibraryReady;
+
+    // increment our shower count
+    ++count_;
+  }
+
+  inline void OutputManager::endOfLibrary() {
+
+    // we can only call endOfLibrary when we have already started
+    if (state_ == OutputState::NoInit) {
+      throw std::runtime_error("endOfLibrary() called in invalid state.");
+    }
+
+    // write the summary for each output and forward the endOfLibrary call()
+    for (auto& [name, output] : outputs_) {
+      // save eventual YAML summary
+      YAML::Node const summary = output.get().getSummary();
+      if (!summary.IsNull()) {
+        writeYAML(output.get().getSummary(), root_ / name / ("summary.yaml"));
+      }
+
+      // and forward the end of library call
+      output.get().endOfLibrary();
+    }
+
+    // and the library has finished
+    state_ = OutputState::LibraryFinished;
+  }
+} // namespace corsika

corsika/detail/output/ParquetStreamer.inl

+/*
+ * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+namespace corsika {
+
+  inline ParquetStreamer::ParquetStreamer()
+      : isInit_(false) {}
+
+  inline void ParquetStreamer::initStreamer(std::string const& filepath) {
+
+    // open the file and connect it to our pointer
+    PARQUET_ASSIGN_OR_THROW(outfile_, arrow::io::FileOutputStream::Open(filepath));
+
+    // the default builder settings
+    builder_.created_by("CORSIKA8");
+
+    // add run and event tags to the file
+    addField("shower", parquet::Repetition::REQUIRED, parquet::Type::INT32,
+             parquet::ConvertedType::UINT_32);
+  }
+
+  template <typename... TArgs>
+  inline void ParquetStreamer::addField(TArgs&&... args) {
+    fields_.push_back(parquet::schema::PrimitiveNode::Make(args...));
+  }
+
+  inline void ParquetStreamer::enableCompression(int const level) {
+    builder_.compression(parquet::Compression::LZ4);
+    builder_.compression_level(level);
+  }
+
+  inline void ParquetStreamer::buildStreamer() {
+
+    // build the top level schema
+    schema_ = std::static_pointer_cast<parquet::schema::GroupNode>(
+        parquet::schema::GroupNode::Make("schema", parquet::Repetition::REQUIRED,
+                                         fields_));
+
+    // and build the writer
+    writer_ = std::make_shared<parquet::StreamWriter>(
+        parquet::ParquetFileWriter::Open(outfile_, schema_, builder_.build()));
+
+    //  only now this object is ready to stream
+    isInit_ = true;
+  }
+
+  inline void ParquetStreamer::closeStreamer() {
+    writer_.reset();
+    [[maybe_unused]] auto status = outfile_->Close();
+    isInit_ = false;
+  }
+
+  inline std::shared_ptr<parquet::StreamWriter> ParquetStreamer::getWriter() {
+    if (!isInit()) {
+      throw std::runtime_error(
+          "ParquetStreamer not initialized. Either 1) add the "
+          "corresponding module to "
+          "the OutputManager, or 2) declare the module to write no output using "
+          "NoOutput.");
+    }
+    return writer_;
+  }
+
+} // namespace corsika

corsika/detail/output/YAMLStreamer.inl

+/*
+ * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <boost/filesystem/fstream.hpp>
+
+namespace corsika {
+
+  inline void YAMLStreamer::writeYAML(YAML::Node const& node,
+                                      boost::filesystem::path const& path) const {
+
+    // construct a YAML emitter for this config file
+    YAML::Emitter out;
+
+    // and write the node to the output
+    out << node;
+
+    // open the output file - this is <output name>.yaml
+    boost::filesystem::ofstream file(path);
+
+    // dump the YAML to the file
+    file << out.c_str() << std::endl;
+  }
+
+} // namespace corsika

corsika/detail/setup/SetupStack.hpp

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/stack/CombinedStack.hpp>
+#include <corsika/stack/GeometryNodeStackExtension.hpp>
+#include <corsika/stack/VectorStack.hpp>
+#include <corsika/stack/WeightStackExtension.hpp>
+#include <corsika/stack/history/HistorySecondaryProducer.hpp>
+#include <corsika/stack/history/HistoryStackExtension.hpp>
+#include <corsika/media/Environment.hpp>
+#include <corsika/media/IMagneticFieldModel.hpp>
+#include <corsika/media/IMediumModel.hpp>
+#include <corsika/media/IMediumPropertyModel.hpp>
+
+namespace corsika {
+
+  namespace setup::detail {
+    template <typename TEnvironment>
+    class StackGenerator {
+    private:
+      using env_type = TEnvironment;
+
+      // ------------------------------------------
+      // add geometry node data to stack. This is fundamentally needed
+      // for robust tracking through multiple volumes.
+
+      // the GeometryNode stack needs to know the type of geometry-nodes from the
+      // environment:
+
+      template <typename TStackIter>
+      using SetupGeometryDataInterface =
+          typename node::MakeGeometryDataInterface<TStackIter, env_type>::type;
+
+      // combine particle data stack with geometry information for tracking
+      template <typename TStackIter>
+      using StackWithGeometryInterface =
+          CombinedParticleInterface<VectorStack::pi_type, SetupGeometryDataInterface,
+                                    TStackIter>;
+
+      using StackWithGeometry =
+          CombinedStack<typename VectorStack::stack_data_type,
+                        node::GeometryData<env_type>, StackWithGeometryInterface,
+                        DefaultSecondaryProducer>;
+
+      template <class T>
+      using StackWithGeometry_PI_type = typename StackWithGeometry::template pi_type<T>;
+
+      // ------------------------------------------
+      // add weight data to stack. This is fundamentally needed
+      // for thinning.
+
+      // the "pure" weight stack (interface)
+      template <typename TStackIter>
+      using SetupWeightDataInterface =
+          typename weights::MakeWeightDataInterface<TStackIter>::type;
+
+      // combine geometry-node-vector data stack with weight information for tracking
+      template <typename TStackIter>
+      using StackWithWeightInterface =
+          CombinedParticleInterface<StackWithGeometry_PI_type, SetupWeightDataInterface,
+                                    TStackIter>;
+
+    public:
+      // the combined stack data: particle + geometry + weight
+      using StackWithWeight =
+          CombinedStack<typename StackWithGeometry::stack_data_type, weights::WeightData,
+                        StackWithWeightInterface, DefaultSecondaryProducer>;
+
+    private:
+      template <typename T>
+      using StackWithWeight_PI_type = typename StackWithWeight::template pi_type<T>;
+
+      // ------------------------------------------
+      // Add [OPTIONAL] history data to stack, too.
+      // This keeps the entire lineage of particles in memory.
+
+      template <typename TStackIter>
+      using StackWithHistoryInterface =
+          CombinedParticleInterface<StackWithWeight_PI_type,
+                                    history::HistoryEventDataInterface, TStackIter>;
+
+    public:
+      using StackWithHistory =
+          CombinedStack<typename StackWithWeight::stack_data_type,
+                        history::HistoryEventData, StackWithHistoryInterface,
+                        history::HistorySecondaryProducer>;
+    };
+
+  } // namespace setup::detail
+
+} // namespace corsika

corsika/detail/stack/VectorStack.inl

+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/stack/Stack.hpp>
+
+#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/framework/geometry/RootCoordinateSystem.hpp>
+#include <corsika/framework/geometry/Vector.hpp>
+
+#include <string>
+#include <tuple>
+#include <vector>
+
+namespace corsika {
+
+  template <typename StackIteratorInterface>
+  inline void ParticleInterface<StackIteratorInterface>::setParticleData(
+      particle_data_type const& v) {
+    this->setPID(std::get<0>(v));
+    this->setKineticEnergy(std::get<1>(v));
+    this->setDirection(std::get<2>(v));
+    this->setPosition(std::get<3>(v));
+    this->setTime(std::get<4>(v));
+  }
+
+  template <typename StackIteratorInterface>
+  inline void ParticleInterface<StackIteratorInterface>::setParticleData(
+      ParticleInterface<StackIteratorInterface> const& parent,
+      secondary_data_type const& v) {
+    this->setPID(std::get<0>(v));
+    this->setKineticEnergy(std::get<1>(v));
+    this->setDirection(std::get<2>(v));
+    this->setPosition(parent.getPosition()); // position
+    this->setTime(parent.getTime());         // parent time
+  }
+
+  template <typename StackIteratorInterface>
+  inline void ParticleInterface<StackIteratorInterface>::setParticleData(
+      ParticleInterface<StackIteratorInterface> const& parent,
+      secondary_extended_data_type const& v) {
+    this->setPID(std::get<0>(v));
+    this->setKineticEnergy(std::get<1>(v));
+    this->setDirection(std::get<2>(v));
+    this->setPosition(parent.getPosition() + std::get<3>(v)); // + position
+    this->setTime(parent.getTime() + std::get<4>(v));         // + parent time
+  }
+
+  template <typename StackIteratorInterface>
+  inline std::string ParticleInterface<StackIteratorInterface>::asString() const {
+    return fmt::format("particle: i={}, PID={}, Ekin={}GeV", super_type::getIndex(),
+                       get_name(this->getPID()), this->getKineticEnergy() / 1_GeV);
+  }
+
+  inline void VectorStackImpl::clear() {
+    dataPID_.clear();
+    dataEkin_.clear();
+    direction_.clear();
+    position_.clear();
+    time_.clear();
+  }
+
+  inline void VectorStackImpl::copy(size_t const i1, size_t const i2) {
+    // index range check
+    if (i1 >= getSize() || i2 >= getSize()) {
+      std::ostringstream err;
+      err << "VectorStackImpl: trying to access data beyond size of stack !";
+      throw std::runtime_error(err.str());
+    }
+    dataPID_[i2] = dataPID_[i1];
+    dataEkin_[i2] = dataEkin_[i1];
+    direction_[i2] = direction_[i1];
+    position_[i2] = position_[i1];
+    time_[i2] = time_[i1];
+  }
+
+  inline void VectorStackImpl::swap(size_t const i1, size_t const i2) {
+    // index range check
+    if (i1 >= getSize() || i2 >= getSize()) {
+      std::ostringstream err;
+      err << "VectorStackImpl: trying to access data beyond size of stack !";
+      throw std::runtime_error(err.str());
+    }
+    std::swap(dataPID_[i2], dataPID_[i1]);
+    std::swap(dataEkin_[i2], dataEkin_[i1]);
+    std::swap(direction_[i2], direction_[i1]);
+    std::swap(position_[i2], position_[i1]);
+    std::swap(time_[i2], time_[i1]);
+  }
+
+  inline void VectorStackImpl::incrementSize() {
+    dataPID_.push_back(Code::Unknown);
+    dataEkin_.push_back(0 * electronvolt);
+
+    CoordinateSystemPtr const& dummyCS = get_root_CoordinateSystem();
+
+    direction_.push_back(DirectionVector(dummyCS, {0, 0, 0}));
+
+    position_.push_back(Point(dummyCS, {0 * meter, 0 * meter, 0 * meter}));
+    time_.push_back(0 * second);
+  }
+
+  inline void VectorStackImpl::decrementSize() {
+    if (dataEkin_.size() > 0) {
+      dataPID_.pop_back();
+      dataEkin_.pop_back();
+      direction_.pop_back();
+      position_.pop_back();
+      time_.pop_back();
+    }
+  }
+
+} // namespace corsika

corsika/detail/stack/WeightStackExtension.inl

+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/core/Logging.hpp>
+#include <corsika/framework/stack/Stack.hpp>
+
+#include <tuple>
+#include <utility>
+#include <vector>
+
+namespace corsika::weights {
+
+  // default version for particle-creation from input data
+  template <typename TParentStack>
+  inline void WeightDataInterface<TParentStack>::setParticleData(
+      std::tuple<double> const v) {
+    setWeight(std::get<0>(v));
+  }
+
+  template <typename TParentStack>
+  inline void WeightDataInterface<TParentStack>::setParticleData(
+      WeightDataInterface const& parent, std::tuple<double> const) {
+    setWeight(parent.getWeight()); // copy Weight from parent particle!
+  }
+
+  template <typename TParentStack>
+  inline void WeightDataInterface<TParentStack>::setParticleData() {
+    setWeight(1);
+  } // default weight
+
+  template <typename TParentStack>
+  inline void WeightDataInterface<TParentStack>::setParticleData(
+      WeightDataInterface const& parent) {
+    setWeight(parent.getWeight()); // copy Weight from parent particle!
+  }
+
+  template <typename TParentStack>
+  inline std::string WeightDataInterface<TParentStack>::asString() const {
+    return fmt::format("weight={}", getWeight());
+  }
+
+  template <typename TParentStack>
+  inline void WeightDataInterface<TParentStack>::setWeight(double const v) {
+    super_type::getStackData().setWeight(super_type::getIndex(), v);
+  }
+
+  template <typename TParentStack>
+  inline double WeightDataInterface<TParentStack>::getWeight() const {
+    return super_type::getStackData().getWeight(super_type::getIndex());
+  }
+
+  // definition of stack-data object to store geometry information
+
+  // these functions are needed for the Stack interface
+  inline void WeightData::clear() { weight_vector_.clear(); }
+
+  inline unsigned int WeightData::getSize() const { return weight_vector_.size(); }
+
+  inline unsigned int WeightData::getCapacity() const { return weight_vector_.size(); }
+
+  inline void WeightData::copy(int const i1, int const i2) {
+    weight_vector_[i2] = weight_vector_[i1];
+  }
+
+  inline void WeightData::swap(int const i1, int const i2) {
+    std::swap(weight_vector_[i1], weight_vector_[i2]);
+  }
+
+  // custom data access function
+  inline void WeightData::setWeight(int const i, double const v) {
+    weight_vector_[i] = v;
+  }
+
+  inline double WeightData::getWeight(int const i) const { return weight_vector_[i]; }
+
+  // these functions are also needed by the Stack interface
+  inline void WeightData::incrementSize() {
+    weight_vector_.push_back(1);
+  } // default weight
+
+  inline void WeightData::decrementSize() {
+    if (weight_vector_.size() > 0) { weight_vector_.pop_back(); }
+  }
+
+} // namespace corsika::weights

corsika/detail/stack/history/HistoryObservationPlane.hpp

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <boost/histogram.hpp>
+
+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
+
+} // namespace corsika::history