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 GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
+ * 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
@@ -29,15 +28,17 @@ namespace corsika {
           particle.getMomentum() / particle.getEnergy() * constants::c;
 
       Point const initialPosition = particle.getPosition();
+
       CORSIKA_LOG_DEBUG(
-          "TrackingB pid: {}"
-          " , E = {} GeV",
-          particle.getPID(), particle.getEnergy() / 1_GeV);
-      CORSIKA_LOG_DEBUG("TrackingB pos: {}", initialPosition.getCoordinates());
-      CORSIKA_LOG_DEBUG("TrackingB   E: {} GeV", particle.getEnergy() / 1_GeV);
-      CORSIKA_LOG_DEBUG("TrackingB   p: {} GeV",
-                        particle.getMomentum().getComponents() / 1_GeV);
-      CORSIKA_LOG_DEBUG("TrackingB   v: {} ", initialVelocity.getComponents());
+          "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();
@@ -50,19 +51,18 @@ namespace corsika {
       }
 
       // charge of the particle, and magnetic field
-      const int chargeNumber = particle.getChargeNumber();
+      auto const charge = particle.getCharge();
       auto magneticfield =
           volumeNode->getModelProperties().getMagneticField(initialPosition);
       auto const magnitudeB = magneticfield.getNorm();
-      CORSIKA_LOG_DEBUG("field={} uT, chargeNumber={}, magnitudeB={} uT",
-                        magneticfield.getComponents() / 1_uT, chargeNumber,
-                        magnitudeB / 1_T);
-      bool const no_deflection = chargeNumber == 0 || magnitudeB == 0_T;
+      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; }
+      //{ [[maybe_unused]] auto& initialTrackNextVolume_dum = initialTrackNextVolume; }
       auto const initialTrackLength = initialTrack.getLength(1);
 
       CORSIKA_LOG_DEBUG("initialTrack(0)={}, initialTrack(1)={}, initialTrackLength={}",
@@ -75,61 +75,62 @@ namespace corsika {
         return std::make_tuple(initialTrack, initialTrackNextVolume);
       }
 
-      HEPMomentumType const pAlongB_delta =
+      HEPMomentumType const p_perp =
           (particle.getMomentum() -
            particle.getMomentum().getParallelProjectionOnto(magneticfield))
               .getNorm();
 
-      if (pAlongB_delta == 0_GeV) {
-        // particle travel along, parallel to magnetic field. Rg is
-        // "0", but for purpose of step limit we return infinity here.
-        CORSIKA_LOG_TRACE("pAlongB_delta is 0_GeV --> parallel");
+      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 =
-          (pAlongB_delta * 1_V / (constants::c * abs(chargeNumber) * magnitudeB * 1_eV));
+      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.01;
+      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"
-      auto const initialMomentum = particle.getMomentum();
-      auto const absMomentum = initialMomentum.getNorm();
-      DirectionVector const direction = initialVelocity.normalized();
+      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, initialTrackLength * firstFraction_);
+          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 + direction * firstHalveSteplength;
-      auto const k =
-          chargeNumber * (constants::c * 1_eV / 1_V) / particle.getMomentum().getNorm();
-      auto const new_direction =
-          direction + direction.cross(magneticfield) * firstHalveSteplength * 2 * k;
+      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, direction.dot(new_direction) / new_direction_norm)) * 180 /
-              M_PI);
+          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.setMomentum(new_direction * absMomentum);
+      particle.setDirection(new_direction);
       auto const [finalTrack, finalTrackNextVolume] =
           tracking_line::Tracking::getTrack(particle);
-      particle.setPosition(initialPosition); // this is not nice...
-      particle.setMomentum(initialMomentum); // this is not nice...
+      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;

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 GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
+ * 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/Vector.hpp>
 #include <corsika/framework/geometry/StraightTrajectory.hpp>
-#include <corsika/framework/geometry/Intersections.hpp>
-#include <corsika/framework/core/Logging.hpp>
+#include <corsika/framework/geometry/Vector.hpp>
 #include <corsika/modules/tracking/Intersect.hpp>
 
+#include <cmath>
 #include <type_traits>
 #include <utility>
 
@@ -28,16 +28,15 @@ namespace corsika::tracking_line {
     VelocityVector const initialVelocity =
         particle.getMomentum() / particle.getEnergy() * constants::c;
 
-    auto const initialPosition = particle.getPosition();
+    auto const& initialPosition = particle.getPosition();
     CORSIKA_LOG_DEBUG(
-        "Tracking pid: {}"
-        " , E = {} GeV",
-        particle.getPID(), particle.getEnergy() / 1_GeV);
-    CORSIKA_LOG_DEBUG("Tracking pos: {}", initialPosition.getCoordinates());
-    CORSIKA_LOG_DEBUG("Tracking   E: {} GeV", particle.getEnergy() / 1_GeV);
-    CORSIKA_LOG_DEBUG("Tracking   p: {} GeV",
-                      particle.getMomentum().getComponents() / 1_GeV);
-    CORSIKA_LOG_DEBUG("Tracking   v: {} ", initialVelocity.getComponents());
+        "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
@@ -51,7 +50,8 @@ namespace corsika::tracking_line {
   template <typename TParticle>
   inline Intersections Tracking::intersect(TParticle const& particle,
                                            Sphere const& sphere) {
-    auto const delta = particle.getPosition() - sphere.getCenter();
+    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();
@@ -63,19 +63,76 @@ namespace corsika::tracking_line {
     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) {
-    Sphere const* sphere = dynamic_cast<Sphere const*>(&volumeNode.getVolume());
-    if (sphere) { return Tracking::intersect<TParticle>(particle, *sphere); }
-    throw std::runtime_error(
-        "The Volume type provided is not supported in Intersect(particle, node)");
+    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>
@@ -89,9 +146,16 @@ namespace corsika::tracking_line {
     CORSIKA_LOG_TRACE("n_dot_v={}, delta={}, momentum={}", n_dot_v, delta,
                       particle.getMomentum());
 
-    return Intersections(n_dot_v.magnitude() == 0
-                             ? std::numeric_limits<TimeType::value_type>::infinity() * 1_s
-                             : n.dot(delta) / n_dot_v);
+    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 GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
+ * 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>
@@ -29,17 +33,28 @@
 
 namespace corsika::urqmd {
 
-  inline UrQMD::UrQMD(boost::filesystem::path const& xs_file) {
+  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 CrossSectionType UrQMD::getTabulatedCrossSection(Code projectileCode,
-                                                          Code targetCode,
-                                                          HEPEnergyType labEnergy) const {
+  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(projectileCode);
+    auto const kinEnergy = labEnergy - get_mass(projectileId);
 
     assert(kinEnergy >= HEPEnergyType::zero());
 
@@ -53,7 +68,7 @@ namespace corsika::urqmd {
     w[2 - 1] = w[2 - 1] - 2 * w[3 - 1];
 
     int projectileIndex;
-    switch (projectileCode) {
+    switch (projectileId) {
       case Code::Proton:
         projectileIndex = 0;
         break;
@@ -87,13 +102,14 @@ namespace corsika::urqmd {
         projectileIndex = 8;
         break;
       default: { // LCOV_EXCL_START since this can never happen due to canInteract
-        CORSIKA_LOG_WARN("UrQMD cross-section not tabulated for {}", projectileCode);
-        return CrossSectionType::zero(); // LCOV_EXCL_STOP
+        CORSIKA_LOG_WARN("UrQMD cross-section not tabulated for {}", projectileId);
+        return CrossSectionType::zero();
+        // LCOV_EXCL_STOP
       }
     }
 
     int targetIndex;
-    switch (targetCode) {
+    switch (targetId) {
       case Code::Nitrogen:
         targetIndex = 0;
         break;
@@ -105,7 +121,7 @@ namespace corsika::urqmd {
         break;
       default:
         std::stringstream ss;
-        ss << "UrQMD cross-section not tabluated for target " << targetCode;
+        ss << "UrQMD cross-section not tabluated for target " << targetId;
         throw std::runtime_error(ss.str().data());
     }
 
@@ -117,54 +133,17 @@ namespace corsika::urqmd {
 
     CORSIKA_LOG_TRACE(
         "UrQMD::GetTabulatedCrossSection proj={}, targ={}, E={} GeV, sigma={}",
-        get_name(projectileCode), get_name(targetCode), labEnergy / 1_GeV, result);
+        get_name(projectileId), get_name(targetId), labEnergy / 1_GeV, result);
 
     return result;
   }
 
-  inline CrossSectionType UrQMD::getCrossSection(Code vProjectileCode, Code vTargetCode,
-                                                 HEPEnergyType vLabEnergy,
-                                                 int vAProjectile = 1) {
-
-    // the following is a translation of ptsigtot() into C++
-    if (vProjectileCode != Code::Nucleus &&
-        !is_nucleus(vTargetCode)) { // both particles are "special"
-      auto const mProj = get_mass(vProjectileCode);
-      auto const mTar = get_mass(vTargetCode);
-      double sqrtS =
-          sqrt(static_pow<2>(mProj) + static_pow<2>(mTar) + 2 * vLabEnergy * mTar) *
-          (1 / 1_GeV);
-
-      // we must set some UrQMD globals first...
-      auto const [ityp, iso3] = convertToUrQMD(vProjectileCode);
-      ::urqmd::inputs_.spityp[0] = ityp;
-      ::urqmd::inputs_.spiso3[0] = iso3;
-
-      auto const [itypTar, iso3Tar] = convertToUrQMD(vTargetCode);
-      ::urqmd::inputs_.spityp[1] = itypTar;
-      ::urqmd::inputs_.spiso3[1] = iso3Tar;
-
-      int one = 1;
-      int two = 2;
-      return ::urqmd::sigtot_(one, two, sqrtS) * 1_mb;
-    } else {
-      int const Ap = vAProjectile;
-      int const At = is_nucleus(vTargetCode) ? get_nucleus_A(vTargetCode) : 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 TParticle> // need template here, as this is called both with
-                                // SetupParticle as well as SetupProjectile
-  inline CrossSectionType UrQMD::getCrossSection(TParticle const& projectile,
-                                                 Code targetCode) const {
-    auto const projectileCode = projectile.getPID();
+  inline CrossSectionType UrQMD::getCrossSection(Code const projectileId,
+                                                 Code const targetId,
+                                                 FourMomentum const& projectileP4,
+                                                 FourMomentum const& targetP4) const {
 
-    if (projectileCode == Code::Nucleus) {
+    if (!isValid(projectileId, targetId)) {
       /*
        * unfortunately unavoidable at the moment until we have tools to get the actual
        * inealstic cross-section from UrQMD
@@ -172,72 +151,61 @@ namespace corsika::urqmd {
       return CrossSectionType::zero();
     }
 
-    return getTabulatedCrossSection(projectileCode, targetCode, projectile.getEnergy());
-  }
-
-  inline bool UrQMD::canInteract(Code vCode) const {
-    // 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
+    // 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));
 
-    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};
+    bool const tabulated = true;
+    if (tabulated) { return getTabulatedCrossSection(projectileId, targetId, Elab); }
 
-    return std::find(std::cbegin(validProjectileCodes), std::cend(validProjectileCodes),
-                     vCode) != std::cend(validProjectileCodes);
-  }
+    // the following is a translation of ptsigtot() into C++
+    if (!is_nucleus(projectileId) &&
+        !is_nucleus(targetId)) { // both particles are "special"
 
-  template <typename TParticle>
-  inline GrammageType UrQMD::getInteractionLength(TParticle const& vParticle) const {
+      double sqrtS = sqrt(sqrtS2) / 1_GeV;
 
-    if (!canInteract(vParticle.getPID())) {
-      // we could do the canInteract check in getCrossSection, too but if
-      // we do it here we have the advantage of avoiding the loop
-      return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
-    }
+      // 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& mediumComposition =
-        vParticle.getNode()->getModelProperties().getNuclearComposition();
-    using namespace std::placeholders;
+      auto const [itypTar, iso3Tar] = corsika::urqmd::convertToUrQMD(targetId);
+      ::urqmd::inputs_.spityp[1] = itypTar;
+      ::urqmd::inputs_.spiso3[1] = iso3Tar;
 
-    CrossSectionType const weightedProdCrossSection = mediumComposition.getWeightedSum(
-        std::bind(&UrQMD::getCrossSection<decltype(vParticle)>, this, vParticle, _1));
+      int one = 1;
+      int two = 2;
+      return ::urqmd::sigtot_(one, two, sqrtS) * 1_mb;
+    }
 
-    return mediumComposition.getAverageMassNumber() * constants::u /
-           weightedProdCrossSection;
+    // 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) {
-
-    auto projectile = view.getProjectile();
-
-    auto projectileCode = projectile.getPID();
-    auto const projectileEnergyLab = projectile.getEnergy();
-    auto const& projectileMomentumLab = projectile.getMomentum();
-    auto const& projectilePosition = projectile.getPosition();
-    auto const projectileTime = projectile.getTime();
-
-    // sample target particle
-    auto const& mediumComposition =
-        projectile.getNode()->getModelProperties().getNuclearComposition();
-    auto const componentCrossSections = std::invoke([&]() {
-      auto const& components = mediumComposition.getComponents();
-      std::vector<CrossSectionType> crossSections;
-      crossSections.reserve(components.size());
-
-      for (auto const c : components) {
-        crossSections.push_back(getCrossSection(projectile, c));
-      }
-
-      return crossSections;
-    });
+  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");
+    }
 
-    auto const targetCode = mediumComposition.sampleTarget(componentCrossSections, RNG_);
-    auto const targetA = get_nucleus_A(targetCode);
-    auto const targetZ = get_nucleus_Z(targetCode);
+    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
@@ -245,14 +213,13 @@ namespace corsika::urqmd {
     ::urqmd::sys_.nsteps = 1;
 
     // initialization regarding projectile
-    if (Code::Nucleus == projectileCode) {
+    if (is_nucleus(projectileId)) {
       // is this everything?
       ::urqmd::inputs_.prspflg = 0;
 
-      ::urqmd::sys_.Ap = projectile.getNuclearA();
-      ::urqmd::sys_.Zp = projectile.getNuclearZ();
-      ::urqmd::rsys_.ebeam = (projectileEnergyLab - projectile.getMass()) * (1 / 1_GeV) /
-                             projectile.getNuclearA();
+      ::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) +
@@ -266,22 +233,19 @@ namespace corsika::urqmd {
           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 = (projectileEnergyLab - projectile.getMass()) * (1 / 1_GeV);
-
-      if (projectileCode == Code::K0Long) {
-        projectileCode = booleanDist_(RNG_) ? Code::K0 : Code::K0Bar;
-      } else if (projectileCode == Code::K0Short) {
-        throw std::runtime_error("K0Short should not interact");
-      }
+      ::urqmd::rsys_.ebeam = (Elab - get_mass(projectileId)) / 1_GeV;
 
-      auto const [ityp, iso3] = convertToUrQMD(projectileCode);
+      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;
     }
 
-    // initilazation regarding target
-    if (is_nucleus(targetCode)) {
+    // initialization regarding target
+    if (is_nucleus(targetId)) {
       ::urqmd::sys_.Zt = targetZ;
       ::urqmd::sys_.At = targetA;
       ::urqmd::inputs_.trspflg = 0; // nucleus as target
@@ -289,120 +253,57 @@ namespace corsika::urqmd {
       ::urqmd::cascinit_(::urqmd::sys_.Zt, ::urqmd::sys_.At, id);
     } else {
       ::urqmd::inputs_.trspflg = 1; // special particle as target
-      auto const [ityp, iso3] = convertToUrQMD(targetCode);
+      auto const [ityp, iso3] = corsika::urqmd::convertToUrQMD(targetId);
       ::urqmd::inputs_.spityp[1] = ityp;
       ::urqmd::inputs_.spiso3[1] = iso3;
     }
 
-    int iflb = 0; // flag for retrying interaction in case of empty event, 0 means retry
+    int iflb =
+        iflb_; // flag for retrying interaction in case of empty event, 0 means retry
     ::urqmd::urqmd_(iflb);
 
     // now retrieve secondaries from UrQMD
-    auto const& originalCS = projectileMomentumLab.getCoordinateSystem();
-    CoordinateSystemPtr const& zAxisFrame =
-        make_rotationToZ(originalCS, projectileMomentumLab);
+    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 = convertFromUrQMD(::urqmd::isys_.ityp[i], ::urqmd::isys_.iso3[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
-      auto momentum =
-          Vector(zAxisFrame,
-                 QuantityVector<dimensionless_d>{
-                     ::urqmd::coor_.px[i], ::urqmd::coor_.py[i], ::urqmd::coor_.pz[i]} *
-                     1_GeV);
-
-      auto const energy = sqrt(momentum.getSquaredNorm() + square(get_mass(code)));
+      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());
 
-      projectile.addSecondary(
-          std::make_tuple(code, energy, momentum, projectilePosition, projectileTime));
+      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 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);
-  }
-
-  inline std::pair<int, int> convertToUrQMD(Code code) {
-    static const std::map<int, std::pair<int, int>> mapPDGToUrQMD{
-        // data mostly from github.com/afedynitch/ParticleDataTool
-        {22, {100, 0}},      // gamma
-        {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 void UrQMD::readXSFile(boost::filesystem::path const& filename) {
+  inline void UrQMD::readXSFile(boost::filesystem::path const filename) {
     boost::filesystem::ifstream file(filename, std::ios::in);
 
     if (!file.is_open()) {
-      throw std::runtime_error(
-          filename.native() +
-          " could not be opened."); // LCOV_EXCL_LINE since this is pointless to test
+      // 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;
@@ -410,7 +311,7 @@ namespace corsika::urqmd {
     std::getline(file, line);
     std::stringstream ss(line);
 
-    char dummy;
+    char dummy; // this is '#'
     int nTargets, nProjectiles, nSupports;
     ss >> dummy >> nTargets >> nProjectiles >> nSupports;
 
@@ -431,6 +332,7 @@ namespace corsika::urqmd {
         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/ObservationPlaneWriterParquet.inl

-/*
- * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#pragma once
-
-namespace corsika {
-
-  ObservationPlaneWriterParquet::ObservationPlaneWriterParquet()
-      : output_() {}
-
-  void ObservationPlaneWriterParquet::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("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);
-
-    // and build the streamer
-    output_.buildStreamer();
-
-    setInit(true);
-  }
-
-  void ObservationPlaneWriterParquet::endOfShower() { ++shower_; }
-
-  void ObservationPlaneWriterParquet::endOfLibrary() { output_.closeStreamer(); }
-
-  void ObservationPlaneWriterParquet::write(Code const& pid, HEPEnergyType const& energy,
-                                            LengthType const& x, LengthType const& y) {
-    if (!isInit()) {
-      std::runtime_error(
-          "ObservationPlaneWriterParquet not initialized. Either 1) add the "
-          "corresponding module to "
-          "the OutputManager, or 2) declare the module to write no output using "
-          "NoOutput.");
-    }
-
-    // write the next row - we must write `shower_` first.
-    *(output_.getWriter()) << shower_ << static_cast<int>(get_PDG(pid))
-                           << static_cast<float>(energy / 1_GeV)
-                           << static_cast<float>(x / 1_m) << static_cast<float>(y / 1_m)
-                           << parquet::EndRow;
-  }
-
-} // namespace corsika

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 GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
+ * 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 {
 
-  TrackWriterParquet::TrackWriterParquet()
-      : output_() {}
+  inline TrackWriterParquet::TrackWriterParquet()
+      : output_()
+      , showerId_(0) {}
 
-  void TrackWriterParquet::startOfLibrary(boost::filesystem::path const& directory) {
+  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("energy", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
+    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);
@@ -29,48 +35,58 @@ namespace corsika {
                      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();
-
-    setInit(true);
+    showerId_ = 0;
   }
 
-  void TrackWriterParquet::endOfShower() { ++shower_; }
+  inline void TrackWriterParquet::startOfShower(unsigned int const showerId) {
+    showerId_ = showerId;
+  }
 
-  void TrackWriterParquet::endOfLibrary() { output_.closeStreamer(); }
+  inline void TrackWriterParquet::endOfShower(unsigned int const) {}
 
-  void TrackWriterParquet::write(Code const& pid, HEPEnergyType const& energy,
-                                 QuantityVector<length_d> const& start,
-                                 QuantityVector<length_d> const& end) {
+  inline void TrackWriterParquet::endOfLibrary() { output_.closeStreamer(); }
 
-    if (!isInit()) {
-      std::runtime_error(
-          "TrackWriterParquet not initialized. Either 1) add the corresponding module to "
-          "the OutputManager, or 2) declare the module to write no output using "
-          "NoOutput.");
-    }
+  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())
-        << shower_
-        << static_cast<int>(get_PDG(pid))
-        << static_cast<float>(energy / 1_GeV)
-        << static_cast<float>(start[0] / 1_m)
-        << static_cast<float>(start[1] / 1_m)
-        << static_cast<float>(start[2] / 1_m)
-        << static_cast<float>(end[0] / 1_m)
-        << static_cast<float>(end[1] / 1_m)
-        << static_cast<float>(end[2] / 1_m)
-        << parquet::EndRow;
+   *(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
+} // namespace corsika
\ No newline at end of file

corsika/detail/output/BaseOutput.inl

-/*
- * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-#pragma once
-
-namespace corsika {
-
-  BaseOutput::BaseOutput()
-      : shower_(0) {}
-
-} // namespace corsika

corsika/detail/output/DummyOutputManager.inl

 /*
  * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu
  *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
+ * 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 {
 
-  DummyOutputManager::DummyOutputManager() {}
+  inline DummyOutputManager::DummyOutputManager() {}
 
-  DummyOutputManager::~DummyOutputManager() {}
+  inline DummyOutputManager::~DummyOutputManager() {}
 
   template <typename TOutput>
-  void DummyOutputManager::add(std::string const& name, TOutput& output) {}
+  inline void DummyOutputManager::add(std::string const&, TOutput&) {}
 
-  void DummyOutputManager::startOfLibrary() {}
+  inline void DummyOutputManager::startOfLibrary() {}
 
-  void DummyOutputManager::startOfShower() {}
+  inline void DummyOutputManager::startOfShower() {}
 
-  void DummyOutputManager::endOfShower() {}
+  inline void DummyOutputManager::endOfShower() {}
 
-  void DummyOutputManager::endOfLibrary() {}
+  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 GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
+ * 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
 
@@ -15,46 +14,126 @@
 #include <ctime>
 #include <sstream>
 
+#include <boost/filesystem.hpp>
+
 #include <fmt/core.h>
 #include <fmt/chrono.h>
 
 namespace corsika {
 
-  void OutputManager::writeNode(YAML::Node const& node,
-                                boost::filesystem::path const& path) const {
+  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() {
 
-    // construct a YAML emitter for this config file
-    YAML::Emitter out;
+    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();
 
-    // and write the node to the output
-    out << node;
+      // Ensure we use "./" if the parent path is empty
+      std::string parent_path = parent.empty() ? "./" : parent.string() + "/";
 
-    // open the output file - this is <output name>.yaml
-    std::ofstream file(path.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());
 
-    // dump the YAML to the file
-    file << out.c_str() << std::endl;
+      if (returnCode) {
+        CORSIKA_LOG_ERROR("Compression returned with error code {}", returnCode);
+      } else {
+        // remove the original directory
+        boost::filesystem::remove_all(root_);
+      }
+    }
   }
 
-  void OutputManager::writeTopLevelConfig() const {
+  inline int OutputManager::getEventId() const { return count_; }
+
+  inline YAML::Node OutputManager::getConfig() const {
 
     YAML::Node config;
 
     // some basic info
-    config["name"] = name_;               // the simulation name
+    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
-
-    // write the node to a file
-    writeNode(config, root_ / ("config.yaml"));
+    config["args"] = cmnd_line_args_;     // the command line parameters
+    return config;
   }
 
-  void OutputManager::writeTopLevelSummary() const {
+  inline YAML::Node OutputManager::getSummary() const {
 
-    YAML::Node config;
+    YAML::Node summary;
 
     // the total number of showers contained in the library
-    config["showers"] = count_;
+    summary["showers"] = count_;
+
+    summary["seed"] = seed_;
 
     // this next section handles writing some time and duration information
 
@@ -75,95 +154,31 @@ namespace corsika {
 
     auto end_time{std::chrono::system_clock::now()};
 
-    // now let's construct an estimate of the runtime
-    auto runtime{end_time - start_time};
+    // 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
-    config["start time"] = timeToString(start_time);
-    config["end time"] = timeToString(end_time);
-    config["runtime"] = fmt::format("{:%H:%M:%S}", runtime);
-
-    // write the node to a file
-    writeNode(config, root_ / ("summary.yaml"));
-  }
-
-  void OutputManager::initOutput(std::string const& name) const {
-    // construct the path to this directory
-    auto const path{root_ / name};
-
-    // create the directory for this process.
-    boost::filesystem::create_directory(path);
+    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);
 
-    // get the config for this output
-    auto config = outputs_.at(name).get().getConfig();
+    std::vector<std::string> output_dirs;
+    for (auto const& outs : outputs_) { output_dirs.push_back(outs.first); }
+    summary["output_dirs"] = output_dirs;
 
-    // and assign the name for this output
-    config["name"] = name;
-
-    // write the config for this output to the file
-    writeNode(config, path / "config.yaml");
+    return summary;
   }
 
-  OutputManager::OutputManager(
-      std::string const& name,
-      boost::filesystem::path const& dir = boost::filesystem::current_path())
-      : name_(name)
-      , root_(dir / name) {
-
-    // check if this directory already exists
-    if (boost::filesystem::exists(root_)) {
-      logger->warn(
-          "Output directory '{}' already exists! This is currenty not supported.",
-          root_.string());
-      throw std::runtime_error("Output directory already exists.");
-    }
-
-    // construct the directory for this library
-    boost::filesystem::create_directory(root_);
-
-    // write the top level config file
-    writeTopLevelConfig();
-  }
-
-  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.
-      logger->warn(
-          "OutputManager was destroyed before endOfShower() called."
-          " The last shower in this libray may be incomplete.");
-      endOfShower();
-    }
-
-    // write the top level summary file (summary.yaml)
-    writeTopLevelSummary();
-
-    // 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(); }
-  }
+  inline void OutputManager::writeSummary() const {
 
-  template <typename TOutput>
-  void OutputManager::add(std::string const& name, TOutput& output) {
-
-    // check if that name is already in the map
-    if (outputs_.count(name) > 0) {
-      logger->warn("'{}' is already registered. All outputs must have unique names.",
-                   name);
-      return;
-    }
-
-    // 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)));
-
-    // and initialize this output
-    initOutput(name);
+    // write the node to a file
+    writeYAML(getSummary(), root_ / ("summary.yaml"));
   }
 
-  void OutputManager::startOfLibrary() {
+  inline void OutputManager::startOfLibrary() {
 
     // this is only valid when we haven't started a library
     // or have already finished a library
@@ -175,42 +190,49 @@ namespace corsika {
     // we now forward this signal to all of our outputs
     for (auto& [name, output] : outputs_) {
 
-      // construct the path to this output subdirectory
-      auto const path{root_ / name};
-
       // and start the library
-      output.get().startOfLibrary(path);
+      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
   }
 
-  void OutputManager::startOfShower() {
+  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(); }
-
-    // increment our shower count
-    ++count_;
+    for (auto& [name, output] : outputs_) { output.get().startOfShower(count_); }
 
     // and transition to the in progress state
     state_ = OutputState::ShowerInProgress;
   }
 
-  void OutputManager::endOfShower() {
+  inline void OutputManager::endOfShower() {
 
-    for (auto& [name, output] : outputs_) { output.get().endOfShower(); }
+    for (auto& [name, output] : outputs_) { output.get().endOfShower(count_); }
 
     // switch back to the initialized state
     state_ = OutputState::LibraryReady;
+
+    // increment our shower count
+    ++count_;
   }
 
-  void OutputManager::endOfLibrary() {
+  inline void OutputManager::endOfLibrary() {
 
     // we can only call endOfLibrary when we have already started
     if (state_ == OutputState::NoInit) {
@@ -219,12 +241,11 @@ namespace corsika {
 
     // write the summary for each output and forward the endOfLibrary call()
     for (auto& [name, output] : outputs_) {
-
-      // we get the summary for each output as a YAML node
-      auto summary{outputs_.at(name).get().getSummary()};
-
-      // write the summary for this output to the file
-      writeNode(summary, root_ / name / "summary.yaml");
+      // 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();
@@ -233,5 +254,4 @@ namespace corsika {
     // 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 GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
+ * 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 {
 
-  ParquetStreamer::ParquetStreamer() {}
+  inline ParquetStreamer::ParquetStreamer()
+      : isInit_(false) {}
 
-  void ParquetStreamer::initStreamer(std::string const& filepath) {
+  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));
@@ -22,20 +22,20 @@ namespace corsika {
 
     // add run and event tags to the file
     addField("shower", parquet::Repetition::REQUIRED, parquet::Type::INT32,
-             parquet::ConvertedType::INT_32);
+             parquet::ConvertedType::UINT_32);
   }
 
   template <typename... TArgs>
-  void ParquetStreamer::addField(TArgs&&... args) {
+  inline void ParquetStreamer::addField(TArgs&&... args) {
     fields_.push_back(parquet::schema::PrimitiveNode::Make(args...));
   }
 
-  void ParquetStreamer::enableCompression(int const /*level*/) {
-    // builder_.compression(parquet::Compression::ZSTD);
-    // builder_.compression_level(level);
+  inline void ParquetStreamer::enableCompression(int const level) {
+    builder_.compression(parquet::Compression::LZ4);
+    builder_.compression_level(level);
   }
 
-  void ParquetStreamer::buildStreamer() {
+  inline void ParquetStreamer::buildStreamer() {
 
     // build the top level schema
     schema_ = std::static_pointer_cast<parquet::schema::GroupNode>(
@@ -45,13 +45,26 @@ namespace corsika {
     // 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;
   }
 
-  void ParquetStreamer::closeStreamer() {
+  inline void ParquetStreamer::closeStreamer() {
     writer_.reset();
     [[maybe_unused]] auto status = outfile_->Close();
+    isInit_ = false;
   }
 
-  std::shared_ptr<parquet::StreamWriter> ParquetStreamer::getWriter() { return writer_; }
+  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/SetupEnvironment.inl

-/*
- * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#pragma once
-
-#include <corsika/framework/geometry/Point.hpp>
-#include <corsika/framework/geometry/CoordinateSystem.hpp>
-
-#include <limits>

corsika/detail/setup/SetupStack.hpp

 /*
  * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
  *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
+ * 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/NuclearStackExtension.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/setup/SetupEnvironment.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>;
 
-    // ------------------------------------------
-    // add geometry node tracking data to stack:
-
-    // the GeometryNode stack needs to know the type of geometry-nodes from the
-    // environment:
-    template <typename TStackIter>
-    using SetupGeometryDataInterface =
-        typename node::MakeGeometryDataInterface<TStackIter, setup::Environment>::type;
-
-    // combine particle data stack with geometry information for tracking
-    template <typename TStackIter>
-    using StackWithGeometryInterface =
-        CombinedParticleInterface<nuclear_stack::ParticleDataStack::pi_type,
-                                  SetupGeometryDataInterface, TStackIter>;
-
-    using StackWithGeometry = CombinedStack<
-        typename nuclear_stack::ParticleDataStack::stack_implementation_type,
-        node::GeometryData<setup::Environment>, StackWithGeometryInterface>;
-
-    // ------------------------------------------
-    // Add [optional] history data to stack, too:
-
-    // combine dummy stack with geometry information for tracking
-    template <typename TStackIter>
-    using StackWithHistoryInterface =
-        CombinedParticleInterface<StackWithGeometry::pi_type,
-                                  history::HistoryEventDataInterface, TStackIter>;
-
-    using StackWithHistory =
-        CombinedStack<typename StackWithGeometry::stack_implementation_type,
-                      history::HistoryEventData, StackWithHistoryInterface>;
+    public:
+      using StackWithHistory =
+          CombinedStack<typename StackWithWeight::stack_data_type,
+                        history::HistoryEventData, StackWithHistoryInterface,
+                        history::HistorySecondaryProducer>;
+    };
 
   } // namespace setup::detail