corsika/detail/media/LayeredSphericalAtmosphereBuilder.hpp

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+namespace corsika {
+
+  namespace detail {
+
+    struct NoExtraModelInner {};
+
+    template <typename M>
+    struct NoExtraModel {};
+
+    template <template <typename> typename M>
+    struct has_extra_models : std::true_type {};
+
+    template <>
+    struct has_extra_models<NoExtraModel> : std::false_type {};
+
+  } // namespace detail
+
+} // namespace corsika

corsika/detail/media/LayeredSphericalAtmosphereBuilder.inl

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/core/Logging.hpp>
+
+#include <corsika/media/FlatExponential.hpp>
+#include <corsika/media/HomogeneousMedium.hpp>
+#include <corsika/media/SlidingPlanarExponential.hpp>
+#include <corsika/media/SlidingPlanarTabular.hpp>
+
+namespace corsika {
+
+  template <typename TMediumInterface, template <typename> typename TMediumModelExtra,
+            typename... TModelArgs>
+  inline void LayeredSphericalAtmosphereBuilder<
+      TMediumInterface, TMediumModelExtra, TModelArgs...>::checkRadius(LengthType const r)
+      const {
+    if (r <= previousRadius_) {
+      throw std::runtime_error("radius must be greater than previous");
+    }
+  }
+
+  template <typename TMediumInterface, template <typename> typename TMediumModelExtra,
+            typename... TModelArgs>
+  inline void LayeredSphericalAtmosphereBuilder<
+      TMediumInterface, TMediumModelExtra,
+      TModelArgs...>::setNuclearComposition(NuclearComposition const& composition) {
+    composition_ = std::make_unique<NuclearComposition>(composition);
+  }
+
+  template <typename TMediumInterface, template <typename> typename TMediumModelExtra,
+            typename... TModelArgs>
+  inline typename LayeredSphericalAtmosphereBuilder<TMediumInterface, TMediumModelExtra,
+                                                    TModelArgs...>::volume_tree_node*
+  LayeredSphericalAtmosphereBuilder<TMediumInterface, TMediumModelExtra, TModelArgs...>::
+      addExponentialLayer(GrammageType const b, LengthType const scaleHeight,
+                          LengthType const upperBoundary) {
+
+    // outer radius
+    auto const radius = planetRadius_ + upperBoundary;
+    checkRadius(radius);
+    previousRadius_ = radius;
+
+    auto node = std::make_unique<VolumeTreeNode<TMediumInterface>>(
+        std::make_unique<Sphere>(center_, radius));
+
+    auto const rho0 = b / scaleHeight;
+
+    if constexpr (detail::has_extra_models<TMediumModelExtra>::value) {
+      // helper lambda in which the last 5 arguments to make_shared<...> are bound
+      auto lastBound = [&](auto... argPack) {
+        return std::make_shared<
+            TMediumModelExtra<SlidingPlanarExponential<TMediumInterface>>>(
+            argPack..., center_, rho0, -scaleHeight, *composition_, planetRadius_);
+      };
+
+      // now unpack the additional arguments
+      auto model = std::apply(lastBound, additionalModelArgs_);
+      node->setModelProperties(std::move(model));
+    } else {
+      node->template setModelProperties<SlidingPlanarExponential<TMediumInterface>>(
+          center_, rho0, -scaleHeight, *composition_, planetRadius_);
+    }
+
+    layers_.push(std::move(node));
+    return layers_.top().get();
+  }
+
+  template <typename TMediumInterface, template <typename> typename TMediumModelExtra,
+            typename... TModelArgs>
+  inline void LayeredSphericalAtmosphereBuilder<
+      TMediumInterface, TMediumModelExtra,
+      TModelArgs...>::addLinearLayer(GrammageType const b, LengthType const scaleHeight,
+                                     LengthType const upperBoundary) {
+    // outer radius
+    auto const radius = planetRadius_ + upperBoundary;
+    checkRadius(radius);
+    previousRadius_ = radius;
+
+    auto node = std::make_unique<VolumeTreeNode<TMediumInterface>>(
+        std::make_unique<Sphere>(center_, radius));
+
+    auto const rho0 = b / scaleHeight;
+
+    if constexpr (detail::has_extra_models<TMediumModelExtra>::value) {
+      // helper lambda in which the last 2 arguments to make_shared<...> are bound
+      auto lastBound = [&](auto... argPack) {
+        return std::make_shared<TMediumModelExtra<HomogeneousMedium<TMediumInterface>>>(
+            argPack..., rho0, *composition_);
+      };
+
+      // now unpack the additional arguments
+      auto model = std::apply(lastBound, additionalModelArgs_);
+      node->setModelProperties(std::move(model));
+    } else {
+      node->template setModelProperties<HomogeneousMedium<TMediumInterface>>(
+          rho0, *composition_);
+    }
+
+    layers_.push(std::move(node));
+  }
+
+  template <typename TMediumInterface, template <typename> typename TMediumModelExtra,
+            typename... TModelArgs>
+  inline void
+  LayeredSphericalAtmosphereBuilder<TMediumInterface, TMediumModelExtra, TModelArgs...>::
+      addTabularLayer(std::function<MassDensityType(LengthType)> const& funcRho,
+                      unsigned int const nBins, LengthType const deltaHeight,
+                      LengthType const upperBoundary) {
+
+    auto const radius = planetRadius_ + upperBoundary;
+    checkRadius(radius);
+    previousRadius_ = radius;
+
+    auto node = std::make_unique<VolumeTreeNode<TMediumInterface>>(
+        std::make_unique<Sphere>(center_, radius));
+
+    if constexpr (detail::has_extra_models<TMediumModelExtra>::value) {
+      // helper lambda in which the last 5 arguments to make_shared<...> are bound
+      auto lastBound = [&](auto... argPack) {
+        return std::make_shared<
+            TMediumModelExtra<SlidingPlanarTabular<TMediumInterface>>>(
+            argPack..., center_, funcRho, nBins, deltaHeight, *composition_,
+            planetRadius_);
+      };
+
+      // now unpack the additional arguments
+      auto model = std::apply(lastBound, additionalModelArgs_);
+      node->setModelProperties(std::move(model));
+    } else {
+      node->template setModelProperties<SlidingPlanarTabular<TMediumInterface>>(
+          center_, funcRho, nBins, deltaHeight, *composition_, planetRadius_);
+    }
+    layers_.push(std::move(node));
+  }
+
+  template <typename TMediumInterface, template <typename> typename TMediumModelExtra,
+            typename... TModelArgs>
+  inline Environment<TMediumInterface> LayeredSphericalAtmosphereBuilder<
+      TMediumInterface, TMediumModelExtra, TModelArgs...>::assemble() {
+    Environment<TMediumInterface> env;
+    assemble(env);
+    return env;
+  }
+
+  template <typename TMediumInterface, template <typename> typename TMediumModelExtra,
+            typename... TModelArgs>
+  inline void LayeredSphericalAtmosphereBuilder<
+      TMediumInterface, TMediumModelExtra,
+      TModelArgs...>::assemble(Environment<TMediumInterface>& env) {
+    auto& universe = env.getUniverse();
+    auto* outmost = universe.get();
+
+    while (!layers_.empty()) {
+      auto l = std::move(layers_.top());
+      auto* tmp = l.get();
+      outmost->addChild(std::move(l));
+      layers_.pop();
+      outmost = tmp;
+    }
+  }
+
+  template <typename TMediumInterface, template <typename> typename MExtraEnvirnoment>
+  struct make_layered_spherical_atmosphere_builder {
+    template <typename... TArgs>
+    static auto create(Point const& center, LengthType const planetRadius,
+                       TArgs... args) {
+      return LayeredSphericalAtmosphereBuilder<TMediumInterface, MExtraEnvirnoment,
+                                               TArgs...>{std::forward<TArgs>(args)...,
+                                                         center, planetRadius};
+    }
+  };
+
+} // namespace corsika

corsika/detail/media/LinearApproximationIntegrator.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/media/LinearApproximationIntegrator.hpp>
+
+namespace corsika {
+
+  template <typename TDerived>
+  inline TDerived const& LinearApproximationIntegrator<TDerived>::getImplementation()
+      const {
+    return *static_cast<TDerived const*>(this);
+  }
+
+  template <typename TDerived>
+  inline GrammageType LinearApproximationIntegrator<TDerived>::getIntegrateGrammage(
+      BaseTrajectory const& line) const {
+    LengthType const length = line.getLength();
+    auto const c0 = getImplementation().evaluateAt(line.getPosition(0));
+    auto const c1 = getImplementation().rho_.getFirstDerivative(line.getPosition(0),
+                                                                line.getDirection(0));
+    CORSIKA_LOG_INFO("length={} c0={} c1={} pos={} dir={} return={}", length, c0, c1,
+                     line.getPosition(0), line.getDirection(0),
+                     (c0 + 0.5 * c1 * length) * length);
+    return (c0 + 0.5 * c1 * length) * length;
+  }
+
+  template <typename TDerived>
+  inline LengthType LinearApproximationIntegrator<TDerived>::getArclengthFromGrammage(
+      BaseTrajectory const& line, GrammageType grammage) const {
+    auto const c0 = getImplementation().rho_(line.getPosition(0));
+    auto const c1 = getImplementation().rho_.getFirstDerivative(line.getPosition(0),
+                                                                line.getDirection(0));
+
+    return (1 - 0.5 * grammage * c1 / (c0 * c0)) * grammage / c0;
+  }
+
+  template <typename TDerived>
+  inline LengthType LinearApproximationIntegrator<TDerived>::getMaximumLength(
+      BaseTrajectory const& line, [[maybe_unused]] double relError) const {
+    [[maybe_unused]] auto const c1 = getImplementation().rho_.getSecondDerivative(
+        line.getPosition(0), line.getDirection(0));
+
+    // todo: provide a real, working implementation
+    return 1_m * std::numeric_limits<double>::infinity();
+  }
+
+} // namespace corsika

corsika/detail/media/MediumPropertyModel.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/media/IMediumPropertyModel.hpp>
+
+namespace corsika {
+
+  template <typename T>
+  template <typename... Args>
+  inline MediumPropertyModel<T>::MediumPropertyModel(Medium const medium, Args&&... args)
+      : T(std::forward<Args>(args)...)
+      , medium_(medium) {}
+
+  template <typename T>
+  inline Medium MediumPropertyModel<T>::getMedium() const {
+    return medium_;
+  }
+
+  template <typename T>
+  inline void MediumPropertyModel<T>::setMedium(Medium const medium) {
+    medium_ = medium;
+  }
+
+} // namespace corsika

corsika/detail/media/NuclearComposition.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>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+
+#include <boost/iterator/zip_iterator.hpp>
+#include <boost/iterator/transform_iterator.hpp>
+
+#include <cassert>
+#include <functional>
+#include <numeric>
+#include <random>
+#include <stdexcept>
+#include <vector>
+
+namespace corsika {
+
+  inline NuclearComposition::NuclearComposition(std::vector<Code> const& pComponents,
+                                                std::vector<double> const& pFractions)
+      : numberFractions_(pFractions)
+      , components_(pComponents)
+      , avgMassNumber_(getWeightedSum([](Code const compID) -> double {
+        if (is_nucleus(compID)) {
+          return get_nucleus_A(compID);
+        } else {
+          return get_mass(compID) / convert_SI_to_HEP(constants::u);
+        }
+      })) {
+    if (pComponents.size() != pFractions.size()) {
+      throw std::runtime_error(
+          "Cannot construct NuclearComposition from vectors of different sizes.");
+    }
+    auto const sumFractions = std::accumulate(pFractions.cbegin(), pFractions.cend(), 0.);
+
+    if (!(0.999 < sumFractions && sumFractions < 1.001)) {
+      throw std::runtime_error("element fractions do not add up to 1");
+    }
+    this->updateHash();
+  }
+
+  template <typename TFunction>
+  inline auto NuclearComposition::getWeighted(TFunction func) const {
+    using ResultQuantity = decltype(func(std::declval<Code>()));
+    auto const product = [&](auto const compID, auto const fraction) {
+      return func(compID) * fraction;
+    };
+
+    if constexpr (phys::units::is_quantity_v<ResultQuantity>) {
+      std::vector<ResultQuantity> result(components_.size(), ResultQuantity::zero());
+      std::transform(components_.cbegin(), components_.cend(), numberFractions_.cbegin(),
+                     result.begin(), product);
+      return result;
+    } else {
+      std::vector<ResultQuantity> result(components_.size(), ResultQuantity(0));
+      std::transform(components_.cbegin(), components_.cend(), numberFractions_.cbegin(),
+                     result.begin(), product);
+      return result;
+    }
+  } // namespace corsika
+
+  template <typename TFunction>
+  inline auto NuclearComposition::getWeightedSum(TFunction func) const
+      -> decltype(func(std::declval<Code>())) {
+    using ResultQuantity = decltype(func(std::declval<Code>()));
+
+    auto const prod = [&](auto const compID, auto const fraction) {
+      return func(compID) * fraction;
+    };
+
+    if constexpr (phys::units::is_quantity_v<ResultQuantity>) {
+      return std::inner_product(
+          components_.cbegin(), components_.cend(), numberFractions_.cbegin(),
+          ResultQuantity::zero(), // .zero() is defined for quantity types only
+          std::plus<ResultQuantity>(), prod);
+    } else {
+      return std::inner_product(
+          components_.cbegin(), components_.cend(), numberFractions_.cbegin(),
+          ResultQuantity(0), // in other cases we have to use a bare 0
+          std::plus<ResultQuantity>(), prod);
+    }
+  }
+
+  inline size_t NuclearComposition::getSize() const { return numberFractions_.size(); }
+
+  inline std::vector<double> const& NuclearComposition::getFractions() const {
+    return numberFractions_;
+  }
+
+  inline std::vector<Code> const& NuclearComposition::getComponents() const {
+    return components_;
+  }
+
+  inline double const NuclearComposition::getAverageMassNumber() const {
+    return avgMassNumber_;
+  }
+
+  template <class TRNG>
+  inline Code NuclearComposition::sampleTarget(std::vector<CrossSectionType> const& sigma,
+                                               TRNG&& randomStream) const {
+    if (sigma.size() != numberFractions_.size()) {
+      throw std::runtime_error("incompatible vector sigma as input");
+    }
+
+    auto zip_beg = boost::make_zip_iterator(
+        boost::make_tuple(numberFractions_.cbegin(), sigma.cbegin()));
+    auto zip_end = boost::make_zip_iterator(
+        boost::make_tuple(numberFractions_.cend(), sigma.cend()));
+    using zip_iter_type = decltype(zip_beg);
+
+    auto const mult_func = [](zip_iter_type::value_type const& zipit) -> double {
+      return zipit.get<0>() * zipit.get<1>().magnitude();
+    };
+
+    using transform_iter_type =
+        boost::transform_iterator<decltype(mult_func), zip_iter_type, double, double>;
+
+    auto trans_beg = transform_iter_type{zip_beg, mult_func};
+    auto trans_end = transform_iter_type{zip_end, mult_func};
+
+    std::discrete_distribution channelDist{trans_beg, trans_end};
+
+    auto const iChannel = channelDist(randomStream);
+    return components_[iChannel];
+  }
+
+  // Note: when this class ever modifies its internal data, the hash
+  // must be updated, too!
+  // the hash value is important to find tables, etc.
+  inline size_t NuclearComposition::getHash() const { return hash_; }
+
+  inline bool NuclearComposition::operator==(NuclearComposition const& v) const {
+    return v.hash_ == hash_;
+  }
+
+  inline void NuclearComposition::updateHash() {
+    std::vector<std::size_t> hashes;
+    for (double ifrac : this->getFractions())
+      hashes.push_back(std::hash<double>{}(ifrac));
+    for (Code icode : this->getComponents())
+      hashes.push_back(std::hash<int>{}(static_cast<int>(icode)));
+    std::size_t h = std::hash<double>{}(this->getAverageMassNumber());
+    for (std::size_t ih : hashes) h = h ^ (ih << 1);
+    hash_ = h;
+  }
+
+} // namespace corsika

corsika/detail/media/ShowerAxis.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.
+ */
+
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/core/Logging.hpp>
+
+#include <string>
+
+namespace corsika {
+
+  template <typename TEnvModel>
+  inline ShowerAxis::ShowerAxis(Point const& pStart, Vector<length_d> const& length,
+                                Environment<TEnvModel> const& env, bool const doThrow,
+                                int const steps)
+      : pointStart_(pStart)
+      , length_(length)
+      , throw_(doThrow)
+      , max_length_(length_.getNorm())
+      , steplength_(max_length_ / steps)
+      , axis_normalized_(length / max_length_)
+      , X_(steps + 1) {
+    auto const* const universe = env.getUniverse().get();
+
+    auto rho = [pStart, length, universe, doThrow](double x) {
+      auto const p = pStart + length * x;
+      auto const* node = universe->getContainingNode(p);
+      if (!node->hasModelProperties()) {
+        CORSIKA_LOG_CRITICAL(
+            "Unable to construct ShowerAxis. ShowerAxis includes volume "
+            "with no model properties at point {}.",
+            p);
+        if (doThrow) throw std::runtime_error("Unable to construct ShowerAxis.");
+      }
+      return node->getModelProperties().getMassDensity(p).magnitude();
+    };
+
+    double error;
+    int k = 0;
+    X_[0] = GrammageType::zero();
+    auto sum = GrammageType::zero();
+
+    for (int i = 1; i <= steps; ++i) {
+      auto const x_prev = (i - 1.) / steps;
+      auto const d_prev = max_length_ * x_prev;
+      auto const x = double(i) / steps;
+      auto const r = boost::math::quadrature::gauss_kronrod<double, 15>::integrate(
+          rho, x_prev, x, 15, 1e-9, &error);
+      auto const result =
+          MassDensityType(phys::units::detail::magnitude_tag, r) * max_length_;
+
+      sum += result;
+      X_[i] = sum;
+
+      for (; sum > k * X_binning_; ++k) {
+        d_.emplace_back(d_prev + k * X_binning_ * steplength_ / result);
+      }
+    }
+
+    assert(std::is_sorted(X_.cbegin(), X_.cend()));
+    assert(std::is_sorted(d_.cbegin(), d_.cend()));
+  }
+
+  inline GrammageType ShowerAxis::getX(LengthType l) const {
+    double const fractionalBin = l / steplength_;
+    int const lower = fractionalBin; // indices of nearest X support points
+    double const fraction = fractionalBin - lower;
+    unsigned int const upper = lower + 1;
+
+    if (fractionalBin < 0) {
+      CORSIKA_LOG_TRACE("cannot extrapolate to points behind point of injection l={} m",
+                        l / 1_m);
+      if (throw_) {
+        throw std::runtime_error(
+            "cannot extrapolate to points behind point of injection");
+      }
+      return getMinimumX();
+    }
+
+    if (upper >= X_.size()) {
+      CORSIKA_LOG_TRACE(
+          "shower axis too short, cannot extrapolate (l / max_length_ = {} )",
+          l / max_length_);
+      if (throw_) {
+        const std::string err = fmt::format(
+            "shower axis too short, cannot extrapolate (l / max_length_ = {} )",
+            l / max_length_);
+        throw std::runtime_error(err.c_str());
+      }
+      return getMaximumX();
+    }
+    CORSIKA_LOG_TRACE("showerAxis::X frac={}, fractionalBin={}, lower={}, upper={}",
+                      fraction, fractionalBin, lower, upper);
+
+    assert(0 <= fraction && fraction <= 1.);
+
+    CORSIKA_LOG_TRACE("ShowerAxis::getX l={} m, lower={}, fraction={}, upper={}", l / 1_m,
+                      lower, fraction, upper);
+
+    // linear interpolation between getX[lower] and X[upper]
+    return X_[upper] * fraction + X_[lower] * (1 - fraction);
+  }
+
+  inline LengthType ShowerAxis::getSteplength() const { return steplength_; }
+
+  inline GrammageType ShowerAxis::getMaximumX() const { return *X_.rbegin(); }
+
+  inline GrammageType ShowerAxis::getMinimumX() const { return GrammageType::zero(); }
+
+  inline GrammageType ShowerAxis::getProjectedX(Point const& p) const {
+    auto const projectedLength = (p - pointStart_).dot(axis_normalized_);
+    return getX(projectedLength);
+  }
+
+  inline DirectionVector const& ShowerAxis::getDirection() const {
+    return axis_normalized_;
+  }
+
+  inline Point const& ShowerAxis::getStart() const { return pointStart_; }
+
+} // namespace corsika

corsika/detail/media/SlidingPlanarExponential.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
+
+namespace corsika {
+
+  template <typename TDerived>
+  inline SlidingPlanarExponential<TDerived>::SlidingPlanarExponential(
+      Point const& p0, MassDensityType const rho0, LengthType const lambda,
+      NuclearComposition const& nuclComp, LengthType const referenceHeight)
+      : BaseExponential<SlidingPlanarExponential<TDerived>>(p0, referenceHeight, rho0,
+                                                            lambda)
+      , nuclComp_(nuclComp) {}
+
+  template <typename TDerived>
+  inline MassDensityType SlidingPlanarExponential<TDerived>::getMassDensity(
+      Point const& point) const {
+    auto const heightFull =
+        (point - BaseExponential<SlidingPlanarExponential<TDerived>>::getAnchorPoint())
+            .getNorm();
+    return BaseExponential<SlidingPlanarExponential<TDerived>>::getMassDensity(
+        heightFull);
+  }
+
+  template <typename TDerived>
+  inline NuclearComposition const&
+  SlidingPlanarExponential<TDerived>::getNuclearComposition() const {
+    return nuclComp_;
+  }
+
+  template <typename TDerived>
+  inline GrammageType SlidingPlanarExponential<TDerived>::getIntegratedGrammage(
+      BaseTrajectory const& traj) const {
+    auto const axis =
+        (traj.getPosition(0) -
+         BaseExponential<SlidingPlanarExponential<TDerived>>::getAnchorPoint())
+            .normalized();
+    return BaseExponential<SlidingPlanarExponential<TDerived>>::getIntegratedGrammage(
+        traj, axis);
+  }
+
+  template <typename TDerived>
+  inline LengthType SlidingPlanarExponential<TDerived>::getArclengthFromGrammage(
+      BaseTrajectory const& traj, GrammageType const grammage) const {
+    auto const axis =
+        (traj.getPosition(0) -
+         BaseExponential<SlidingPlanarExponential<TDerived>>::getAnchorPoint())
+            .normalized();
+    return BaseExponential<SlidingPlanarExponential<TDerived>>::getArclengthFromGrammage(
+        traj, grammage, axis);
+  }
+
+} // namespace corsika

corsika/detail/media/SlidingPlanarTabular.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
+
+namespace corsika {
+
+  template <typename T>
+  inline SlidingPlanarTabular<T>::SlidingPlanarTabular(
+      Point const& p0, std::function<MassDensityType(LengthType)> const& rho,
+      unsigned int const nBins, LengthType const deltaHeight,
+      NuclearComposition const& nuclComp, LengthType referenceHeight)
+      : BaseTabular<SlidingPlanarTabular<T>>(p0, referenceHeight, rho, nBins, deltaHeight)
+      , nuclComp_(nuclComp) {}
+
+  template <typename T>
+  inline MassDensityType SlidingPlanarTabular<T>::getMassDensity(
+      Point const& point) const {
+    auto const heightFull =
+        (point - BaseTabular<SlidingPlanarTabular<T>>::getAnchorPoint()).getNorm();
+    return BaseTabular<SlidingPlanarTabular<T>>::getMassDensity(heightFull);
+  }
+
+  template <typename T>
+  inline NuclearComposition const& SlidingPlanarTabular<T>::getNuclearComposition()
+      const {
+    return nuclComp_;
+  }
+
+  template <typename T>
+  inline GrammageType SlidingPlanarTabular<T>::getIntegratedGrammage(
+      BaseTrajectory const& traj) const {
+    return BaseTabular<SlidingPlanarTabular<T>>::getIntegratedGrammage(traj);
+  }
+
+  template <typename T>
+  inline LengthType SlidingPlanarTabular<T>::getArclengthFromGrammage(
+      BaseTrajectory const& traj, GrammageType const grammage) const {
+    return BaseTabular<SlidingPlanarTabular<T>>::getArclengthFromGrammage(traj, grammage);
+  }
+
+} // namespace corsika

corsika/detail/media/UniformRefractiveIndex.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/media/IRefractiveIndexModel.hpp>
+
+namespace corsika {
+
+  template <typename T>
+  template <typename... Args>
+  inline UniformRefractiveIndex<T>::UniformRefractiveIndex(double const n, Args&&... args)
+      : T(std::forward<Args>(args)...)
+      , n_(n) {}
+
+  template <typename T>
+  inline double UniformRefractiveIndex<T>::getRefractiveIndex(Point const&) const {
+    return n_;
+  }
+
+  template <typename T>
+  inline void UniformRefractiveIndex<T>::setRefractiveIndex(double const n) {
+    n_ = n;
+  }
+
+} // namespace corsika

corsika/detail/media/Universe.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/media/Environment.hpp>
+
+namespace corsika {
+
+  inline Universe::Universe(CoordinateSystemPtr const& pCS)
+      : corsika::Sphere(Point{pCS, 0 * meter, 0 * meter, 0 * meter},
+                        meter * std::numeric_limits<double>::infinity()) {}
+
+  inline bool Universe::contains(corsika::Point const&) const { return true; }
+
+} // namespace corsika

corsika/detail/media/VolumeTreeNode.inl

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/geometry/IVolume.hpp>
+#include <corsika/media/IMediumModel.hpp>
+
+namespace corsika {
+
+  template <typename IModelProperties>
+  inline bool VolumeTreeNode<IModelProperties>::contains(Point const& p) const {
+    return geoVolume_->contains(p);
+  }
+
+  template <typename IModelProperties>
+  inline VolumeTreeNode<IModelProperties> const*
+  VolumeTreeNode<IModelProperties>::excludes(Point const& p) const {
+    auto exclContainsIter =
+        std::find_if(excludedNodes_.cbegin(), excludedNodes_.cend(),
+                     [&](auto const& s) { return bool(s->contains(p)); });
+
+    return exclContainsIter != excludedNodes_.cend() ? *exclContainsIter : nullptr;
+  }
+
+  /** returns a pointer to the sub-VolumeTreeNode which is "responsible" for the given
+   * \class Point \p p, or nullptr iff \p p is not contained in this volume.
+   */
+  template <typename IModelProperties>
+  inline VolumeTreeNode<IModelProperties> const*
+  VolumeTreeNode<IModelProperties>::getContainingNode(Point const& p) const {
+    if (!contains(p)) { return nullptr; }
+
+    if (auto const childContainsIter =
+            std::find_if(childNodes_.cbegin(), childNodes_.cend(),
+                         [&](auto const& s) { return bool(s->contains(p)); });
+        childContainsIter == childNodes_.cend()) // not contained in any of the children
+    {
+      if (auto const exclContainsIter = excludes(p)) // contained in any excluded nodes
+      {
+        return exclContainsIter->getContainingNode(p);
+      } else {
+        return this;
+      }
+    } else {
+      return (*childContainsIter)->getContainingNode(p);
+    }
+  }
+
+  template <typename IModelProperties>
+  inline VolumeTreeNode<IModelProperties>*
+  VolumeTreeNode<IModelProperties>::getContainingNode(Point const& p) {
+    // see Scott Meyers, Effective C++ 3rd ed., Item 3
+    return const_cast<VTN_type*>(std::as_const(*this).getContainingNode(p));
+  }
+
+  template <typename IModelProperties>
+  inline void VolumeTreeNode<IModelProperties>::addChildToContainingNode(Point const& p,
+                                                                         VTNUPtr pChild) {
+    VolumeTreeNode<IModelProperties>* node = getContainingNode(p);
+    if (!node) {
+      CORSIKA_LOG_ERROR("Adding child at {} failed!. No containing node", p);
+      throw std::runtime_error("Failed adding child node. No parent at chosen location");
+    }
+    node->addChild(std::move(pChild));
+  }
+
+  template <typename IModelProperties>
+  template <typename TCallable, bool preorder>
+  inline void VolumeTreeNode<IModelProperties>::walk(TCallable func) const {
+    if constexpr (preorder) { func(*this); }
+
+    std::for_each(childNodes_.begin(), childNodes_.end(),
+                  [&](auto const& v) { v->walk(func); });
+
+    if constexpr (!preorder) { func(*this); };
+  }
+
+  template <typename IModelProperties>
+  inline void VolumeTreeNode<IModelProperties>::addChild(
+      typename VolumeTreeNode<IModelProperties>::VTNUPtr pChild) {
+    pChild->parentNode_ = this;
+    childNodes_.push_back(std::move(pChild));
+    // It is a bad idea to return an iterator to the inserted element
+    // because it might get invalidated when the vector needs to grow
+    // later and the caller won't notice.
+  }
+
+  template <typename IModelProperties>
+  inline void VolumeTreeNode<IModelProperties>::excludeOverlapWith(
+      typename VolumeTreeNode<IModelProperties>::VTNUPtr const& pNode) {
+    excludedNodes_.push_back(pNode.get());
+  }
+
+} // namespace corsika

corsika/detail/modules/HadronicElasticModel.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/HadronicElasticModel.hpp>
+
+#include <corsika/media/Environment.hpp>
+#include <corsika/media/NuclearComposition.hpp>
+#include <corsika/framework/geometry/FourVector.hpp>
+#include <corsika/framework/random/ExponentialDistribution.hpp>
+#include <corsika/framework/utility/COMBoost.hpp>
+
+#include <iomanip>
+#include <iostream>
+
+namespace corsika {
+
+  inline HadronicElasticInteraction::HadronicElasticInteraction(CrossSectionType x,
+                                                                CrossSectionType y)
+      : parX_(x)
+      , parY_(y) {}
+
+  template <typename TParticle>
+  inline GrammageType HadronicElasticInteraction::getInteractionLength(
+      TParticle const& p) {
+    if (p.getPID() == Code::Proton) {
+      auto const* currentNode = p.getNode();
+      auto const& mediumComposition =
+          currentNode->getModelProperties().getNuclearComposition();
+
+      auto const& components = mediumComposition.getComponents();
+      auto const& fractions = mediumComposition.getFractions();
+
+      auto const projectileMomentum = p.getMomentum();
+      auto const projectileMomentumSquaredNorm = projectileMomentum.getSquaredNorm();
+      auto const projectileEnergy = p.getEnergy();
+
+      auto const avgCrossSection = [&]() {
+        CrossSectionType avgCrossSection = 0_b;
+
+        for (size_t i = 0; i < fractions.size(); ++i) {
+          auto const targetMass = get_mass(components[i]);
+          auto const s = static_pow<2>(projectileEnergy + targetMass) -
+                         projectileMomentumSquaredNorm;
+          avgCrossSection += getCrossSection(s) * fractions[i];
+        }
+
+        CORSIKA_LOG_DEBUG("avgCrossSection: {} mb", avgCrossSection / 1_mb);
+
+        return avgCrossSection;
+      }();
+
+      auto const avgTargetMassNumber = mediumComposition.getAverageMassNumber();
+
+      GrammageType const interactionLength =
+          avgTargetMassNumber * constants::u / avgCrossSection;
+
+      return interactionLength;
+    } else {
+      return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
+    }
+  }
+
+  template <typename TParticle>
+  inline ProcessReturn HadronicElasticInteraction::doInteraction(TParticle& p) {
+    if (p.getPID() != Code::Proton) { return ProcessReturn::Ok; }
+
+    const auto* currentNode = p.getNode();
+    const auto& composition = currentNode->getModelProperties().getNuclearComposition();
+    const auto& components = composition.getComponents();
+
+    std::vector<CrossSectionType> cross_section_of_components(
+        composition.getComponents().size());
+
+    auto const projectileMomentum = p.getMomentum();
+    auto const projectileMomentumSquaredNorm = projectileMomentum.getSquaredNorm();
+    auto const projectileEnergy = p.getEnergy();
+
+    for (size_t i = 0; i < components.size(); ++i) {
+      auto const targetMass = get_mass(components[i]);
+      auto const s =
+          static_pow<2>(projectileEnergy + targetMass) - projectileMomentumSquaredNorm;
+      cross_section_of_components[i] = CrossSection(s);
+    }
+
+    const auto targetCode = composition.SampleTarget(cross_section_of_components, RNG_);
+
+    auto const targetMass = get_mass(targetCode);
+
+    std::uniform_real_distribution phiDist(0., 2 * M_PI);
+
+    FourVector const projectileLab(projectileEnergy, projectileMomentum);
+    FourVector const targetLab(
+        targetMass,
+        MomentumVector(projectileMomentum.getCoordinateSystem(), {0_eV, 0_eV, 0_eV}));
+    COMBoost const boost(projectileLab, targetMass);
+
+    auto const projectileCoM = boost.toCoM(projectileLab);
+    auto const targetCoM = boost.toCoM(targetLab);
+
+    auto const pProjectileCoMSqNorm =
+        projectileCoM.getSpaceLikeComponents().getSquaredNorm();
+    auto const pProjectileCoMNorm = sqrt(pProjectileCoMSqNorm);
+
+    auto const eProjectileCoM = projectileCoM.getTimeLikeComponent();
+    auto const eTargetCoM = targetCoM.getTimeLikeComponent();
+
+    auto const sqrtS = eProjectileCoM + eTargetCoM;
+    auto const s = static_pow<2>(sqrtS);
+
+    auto const B = this->B(s);
+    CORSIKA_LOG_DEBUG(B);
+
+    ExponentialDistribution tDist(1 / B);
+    auto const absT = [&]() {
+      decltype(tDist(RNG_)) absT;
+      auto const maxT = 4 * pProjectileCoMSqNorm;
+
+      do {
+        // |t| cannot become arbitrarily large, max. given by GER eq. (4.16), so we just
+        // throw again until we have an acceptable value. Note that the formula holds in
+        // any frame despite of what is stated in the book.
+        absT = tDist(RNG_);
+      } while (absT >= maxT);
+
+      return absT;
+    }();
+
+    CORSIKA_LOG_DEBUG(
+        "HadronicElasticInteraction: s = {}"
+        " GeV²; absT = {} "
+        " GeV² (max./GeV² = {})",
+        s * constants::invGeVsq, absT * constants::invGeVsq,
+        4 * constants::invGeVsq * projectileMomentumSquaredNorm);
+
+    auto const theta = 2 * asin(sqrt(absT / (4 * pProjectileCoMSqNorm)));
+    auto const phi = phiDist(RNG_);
+
+    auto const projectileScatteredLab =
+        boost.fromCoM(FourVector<HEPEnergyType, MomentumVector>(
+            eProjectileCoM, MomentumVector(projectileMomentum.getCoordinateSystem(),
+                                           {pProjectileCoMNorm * sin(theta) * cos(phi),
+                                            pProjectileCoMNorm * sin(theta) * sin(phi),
+                                            pProjectileCoMNorm * cos(theta)})));
+
+    p.setMomentum(projectileScatteredLab.getSpaceLikeComponents());
+    p.setEnergy(
+        sqrt(projectileScatteredLab.getSpaceLikeComponents().getSquaredNorm() +
+             static_pow<2>(get_mass(
+                 p.getPID())))); // Don't use energy from boost. It can be smaller than
+                                 // the momentum due to limited numerical accuracy.
+
+    return ProcessReturn::Ok;
+  }
+
+  inline HadronicElasticInteraction::inveV2 HadronicElasticInteraction::B(eV2 s) const {
+    auto constexpr b_p = 2.3;
+    auto const result =
+        (2 * b_p + 2 * b_p + 4 * pow(s * constants::invGeVsq, gfEpsilon) - 4.2) *
+        constants::invGeVsq;
+    CORSIKA_LOG_DEBUG("B({}) = {}  GeV¯²", s, result / constants::invGeVsq);
+
+    return result;
+  }
+
+  inline CrossSectionType HadronicElasticInteraction::getCrossSection(
+      SquaredHEPEnergyType s) const {
+    // assuming every target behaves like a proton, parX_ and parY_ are universal
+    CrossSectionType const sigmaTotal = parX_ * pow(s * constants::invGeVsq, gfEpsilon) +
+                                        parY_ * pow(s * constants::invGeVsq, -gfEta);
+
+    // according to Schuler & Sjöstrand, PRD 49, 2257 (1994)
+    // (we ignore rho because rho^2 is just ~2 %)
+    auto const sigmaElastic =
+        static_pow<2>(sigmaTotal) /
+        (16 * constants::pi * convert_HEP_to_SI<CrossSectionType::dimension_type>(B(s)));
+
+    CORSIKA_LOG_DEBUG("HEM sigmaTot = {} mb", sigmaTotal / 1_mb);
+    CORSIKA_LOG_DEBUG("HEM sigmaElastic = {} mb", sigmaElastic / 1_mb);
+    return sigmaElastic;
+  }
+
+} // namespace corsika

corsika/detail/modules/LongitudinalProfile.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.
+ */
+
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <corsika/framework/core/Logging.hpp>
+#include <corsika/framework/core/Step.hpp>
+#include <corsika/modules/LongitudinalProfile.hpp>
+
+#include <cmath>
+#include <iomanip>
+#include <limits>
+#include <utility>
+
+namespace corsika {
+
+  template <typename TOutput>
+  template <typename... TArgs>
+  inline LongitudinalProfile<TOutput>::LongitudinalProfile(TArgs&&... args)
+      : TOutput(std::forward<TArgs>(args)...) {}
+
+  template <typename TOutput>
+  template <typename TParticle>
+  inline ProcessReturn LongitudinalProfile<TOutput>::doContinuous(
+      Step<TParticle> const& step, bool const) {
+
+    auto const pid = step.getParticlePre().getPID();
+    this->write(step.getPositionPre(), step.getPositionPost(), pid,
+                step.getParticlePre().getWeight()); // weight hardcoded so far
+    return ProcessReturn::Ok;
+  }
+
+  template <typename TOutput>
+  inline YAML::Node LongitudinalProfile<TOutput>::getConfig() const {
+    YAML::Node node;
+    node["type"] = "LongitudinalProfile";
+
+    return node;
+  }
+
+} // namespace corsika

corsika/detail/modules/ObservationPlane.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.
+ */
+
+namespace corsika {
+
+  template <typename TTracking, typename TOutput>
+  template <typename... TArgs>
+  ObservationPlane<TTracking, TOutput>::ObservationPlane(Plane const& obsPlane,
+                                                         DirectionVector const& x_axis,
+                                                         bool const deleteOnHit,
+                                                         TArgs&&... args)
+      : TOutput(std::forward<TArgs>(args)...)
+      , plane_(obsPlane)
+      , xAxis_(x_axis.normalized())
+      , yAxis_(obsPlane.getNormal().cross(xAxis_))
+      , deleteOnHit_(deleteOnHit) {}
+
+  template <typename TTracking, typename TOutput>
+  template <typename TParticle>
+  inline ProcessReturn ObservationPlane<TTracking, TOutput>::doContinuous(
+      Step<TParticle>& step, bool const stepLimit) {
+    /*
+       The current step did not yet reach the ObservationPlane, do nothing now and wait:
+     */
+    if (!stepLimit) {
+      // @todo this is actually needed to fix small instabilities of the leap-frog
+      // tracking: Note, this is NOT a general solution and should be clearly revised with
+      // a more robust tracking. #ifdef DEBUG
+      if (deleteOnHit_) {
+        // since this is basically a bug, it cannot be tested LCOV_EXCL_START
+        LengthType const check =
+            (step.getPositionPost() - plane_.getCenter()).dot(plane_.getNormal());
+        if (check < 0_m) {
+          CORSIKA_LOG_WARN("PARTICLE AVOIDED OBSERVATIONPLANE {}", check);
+          CORSIKA_LOG_WARN("Temporary fix: write and remove particle.");
+        } else
+          return ProcessReturn::Ok;
+        // LCOV_EXCL_STOP
+      } else
+        // #endif
+        return ProcessReturn::Ok;
+    }
+
+    HEPEnergyType const kineticEnergy = step.getEkinPost();
+    Point const pointOfIntersection = step.getPositionPost();
+    Vector const displacement = pointOfIntersection - plane_.getCenter();
+    DirectionVector const direction = step.getDirectionPost();
+
+    // add our particles to the output file stream
+    double const weight = step.getParticlePre().getWeight();
+    this->write(step.getParticlePre().getPID(), kineticEnergy, displacement.dot(xAxis_),
+                displacement.dot(yAxis_), 0_m, direction.dot(xAxis_),
+                direction.dot(yAxis_), direction.dot(plane_.getNormal()),
+                step.getTimePost(), weight);
+
+    CORSIKA_LOG_TRACE("Particle detected absorbed={}", deleteOnHit_);
+
+    if (deleteOnHit_) {
+      return ProcessReturn::ParticleAbsorbed;
+    } else {
+      return ProcessReturn::Ok;
+    }
+  } // namespace corsika
+
+  template <typename TTracking, typename TOutput>
+  template <typename TParticle, typename TTrajectory>
+  inline LengthType ObservationPlane<TTracking, TOutput>::getMaxStepLength(
+      TParticle const& particle, TTrajectory const& trajectory) {
+
+    CORSIKA_LOG_TRACE("getMaxStepLength, particle={}, pos={}, dir={}, plane={}",
+                      particle.asString(), particle.getPosition(),
+                      particle.getDirection(), plane_.asString());
+
+    auto const intersection = TTracking::intersect(particle, plane_);
+
+    TimeType const timeOfIntersection = intersection.getEntry();
+    CORSIKA_LOG_TRACE("timeOfIntersection={}", timeOfIntersection);
+    if (timeOfIntersection <= TimeType::zero()) {
+      return std::numeric_limits<double>::infinity() * 1_m;
+    }
+    if (timeOfIntersection > trajectory.getDuration()) {
+      return std::numeric_limits<double>::infinity() * 1_m;
+    }
+    double const fractionOfIntersection = timeOfIntersection / trajectory.getDuration();
+    CORSIKA_LOG_TRACE("ObservationPlane: getMaxStepLength dist={} m, pos={}",
+                      trajectory.getLength(fractionOfIntersection) / 1_m,
+                      trajectory.getPosition(fractionOfIntersection));
+    return trajectory.getLength(fractionOfIntersection);
+  }
+
+  template <typename TTracking, typename TOutput>
+  inline YAML::Node ObservationPlane<TTracking, TOutput>::getConfig() const {
+    using namespace units::si;
+
+    // construct the top-level node
+    YAML::Node node;
+
+    // basic info
+    node["type"] = "ObservationPlane";
+    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{plane_.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);
+
+    // the normal vector of the plane
+    auto const normal{plane_.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());
+
+    // the x-axis vector
+    auto const xAxis_coords{xAxis_.getComponents(xAxis_.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());
+
+    // the y-axis vector
+    auto const yAxis_coords{yAxis_.getComponents(yAxis_.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());
+
+    node["delete_on_hit"] = deleteOnHit_;
+
+    return node;
+  }
+
+} // namespace corsika

corsika/detail/modules/ObservationVolume.inl

+/*
+ * (c) Copyright 2023 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.
+ */
+
+namespace corsika {
+
+  template <typename TTracking, typename TVolume, typename TOutput>
+  ObservationVolume<TTracking, TVolume, TOutput>::ObservationVolume(TVolume vol)
+      : vol_(vol)
+      , energy_(0_GeV)
+      , count_(0) {}
+
+  template <typename TTracking, typename TVolume, typename TOutput>
+  template <typename TParticle>
+  inline ProcessReturn ObservationVolume<TTracking, TVolume, TOutput>::doContinuous(
+      Step<TParticle>& step, bool const stepLimit) {
+    /*
+       The current step did not yet reach the ObservationVolume, do nothing now and
+       wait:
+     */
+    if (!stepLimit) { return ProcessReturn::Ok; }
+
+    HEPEnergyType const kineticEnergy = step.getEkinPost();
+    Point const pointOfIntersection = step.getPositionPost();
+    DirectionVector const dirction = step.getDirectionPost();
+    double const weight = step.getParticlePre().getWeight();
+
+    // add particles to the output file stream
+    auto cs = vol_.getCoordinateSystem();
+    this->write(step.getParticlePre().getPID(), kineticEnergy,
+                pointOfIntersection.getX(cs), pointOfIntersection.getY(cs),
+                pointOfIntersection.getZ(cs), dirction.getX(cs), dirction.getY(cs),
+                dirction.getZ(cs), step.getTimePost(), weight);
+
+    // always absorb particles
+    count_++;
+    energy_ += kineticEnergy + get_mass(step.getParticlePre().getPID());
+    return ProcessReturn::ParticleAbsorbed;
+  }
+
+  template <typename TTracking, typename TVolume, typename TOutput>
+  template <typename TParticle, typename TTrajectory>
+  inline LengthType ObservationVolume<TTracking, TVolume, TOutput>::getMaxStepLength(
+      TParticle const& particle, TTrajectory const& trajectory) {
+
+    CORSIKA_LOG_TRACE("getMaxStepLength, particle={}, pos={}, dir={}, Box={}",
+                      particle.asString(), particle.getPosition(),
+                      particle.getDirection(), vol_.asString());
+
+    auto const intersection = TTracking::intersect(particle, vol_);
+
+    TimeType const timeOfEntry = intersection.getEntry();
+    TimeType const timeOfExit = intersection.getExit();
+    CORSIKA_LOG_TRACE("timeOfEntry={}, timeOfExit={}", timeOfEntry, timeOfExit);
+    if (timeOfEntry < TimeType::zero()) {
+      if (timeOfExit < TimeType::zero()) {
+        // opposite direction
+        return std::numeric_limits<double>::infinity() * 1_m;
+      } else {
+        // inside box: not zero but 1 pm, to allow short lived particles to decay
+        return 1e-12 * 1_m;
+      }
+    }
+    if (timeOfEntry > trajectory.getDuration()) {
+      // can not reach
+      return std::numeric_limits<double>::infinity() * 1_m;
+    }
+    double const fractionOfIntersection = timeOfEntry / trajectory.getDuration();
+    CORSIKA_LOG_TRACE("ObservationVolume: getMaxStepLength dist={} m, pos={}",
+                      trajectory.getLength(fractionOfIntersection) / 1_m,
+                      trajectory.getPosition(fractionOfIntersection));
+    return trajectory.getLength(fractionOfIntersection);
+  }
+
+  template <typename TTracking, typename TVolume, typename TOutput>
+  inline void ObservationVolume<TTracking, TVolume, TOutput>::showResults() const {
+    CORSIKA_LOG_INFO(
+        " ******************************\n"
+        " ObservationVolume: \n"
+        " energy at Box (GeV)     :  {}\n"
+        " no. of particles at Box :  {}\n"
+        " ******************************",
+        energy_ / 1_GeV, count_);
+  }
+
+  template <typename TTracking, typename TVolume, typename TOutput>
+  inline YAML::Node ObservationVolume<TTracking, TVolume, TOutput>::getConfig() const {
+    using namespace units::si;
+    auto cs = vol_.getCoordinateSystem();
+    auto center = vol_.getCenter();
+
+    // construct the top-level node
+    YAML::Node node;
+
+    // basic info
+    node["type"] = "ObservationVolume";
+    node["units"]["length"] = "m";
+    node["units"]["energy"] = "GeV";
+    node["units"]["time"] = "s";
+
+    // save each component in its native coordinate system
+    auto const root_cs = get_root_CoordinateSystem();
+    node["center"].push_back(center.getX(root_cs) / 1_m);
+    node["center"].push_back(center.getY(root_cs) / 1_m);
+    node["center"].push_back(center.getZ(root_cs) / 1_m);
+
+    // the x-axis vector
+    DirectionVector const x_axis = DirectionVector{cs, {1, 0, 0}};
+    node["x-axis"].push_back(x_axis.getX(root_cs).magnitude());
+    node["x-axis"].push_back(x_axis.getY(root_cs).magnitude());
+    node["x-axis"].push_back(x_axis.getZ(root_cs).magnitude());
+
+    // the y-axis vector
+    DirectionVector const y_axis = DirectionVector{cs, {0, 1, 0}};
+    node["y-axis"].push_back(y_axis.getX(root_cs).magnitude());
+    node["y-axis"].push_back(y_axis.getY(root_cs).magnitude());
+    node["y-axis"].push_back(y_axis.getZ(root_cs).magnitude());
+
+    // the x-axis vector
+    DirectionVector const z_axis = DirectionVector{cs, {0, 0, 1}};
+    node["z-axis"].push_back(z_axis.getX(root_cs).magnitude());
+    node["z-axis"].push_back(z_axis.getY(root_cs).magnitude());
+    node["z-axis"].push_back(z_axis.getZ(root_cs).magnitude());
+
+    return node;
+  }
+
+  template <typename TTracking, typename TVolume, typename TOutput>
+  inline void ObservationVolume<TTracking, TVolume, TOutput>::reset() {
+    energy_ = 0_GeV;
+    count_ = 0;
+  }
+
+} // namespace corsika

corsika/detail/modules/OnShellCheck.inl

+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/geometry/FourVector.hpp>
+#include <corsika/modules/OnShellCheck.hpp>
+#include <corsika/framework/core/Logging.hpp>
+
+namespace corsika {
+
+  inline OnShellCheck::OnShellCheck(double const vMassTolerance,
+                                    double const vEnergyTolerance, bool const vError)
+      : mass_tolerance_(vMassTolerance)
+      , energy_tolerance_(vEnergyTolerance)
+      , throw_error_(vError) {
+    CORSIKA_LOGGER_DEBUG(logger_, "mass tolerance is set to {:3.2f}%",
+                         mass_tolerance_ * 100);
+    CORSIKA_LOGGER_DEBUG(logger_, "energy tolerance is set to {:3.2f}%",
+                         energy_tolerance_ * 100);
+  }
+
+  inline OnShellCheck::~OnShellCheck() {
+    logger_->info(
+        " summary \n"
+        " particles shifted: {} \n"
+        " average energy shift (%): {} \n"
+        " max. energy shift (%): {} ",
+        int(count_), (count_ ? average_shift_ / count_ * 100 : 0), max_shift_ * 100.);
+  }
+
+  template <typename TView>
+  inline void OnShellCheck::doSecondaries(TView& vS) {
+    for (auto& p : vS) {
+      auto const pid = p.getPID();
+      if (is_nucleus(pid)) continue;
+      auto const e_original = p.getEnergy();
+      auto const p_original = p.getMomentum();
+      auto const Plab = FourVector(e_original, p_original);
+      auto const m_kinetic = Plab.getNorm();
+      auto const m_corsika = get_mass(pid);
+      auto const m_err_abs = abs(m_kinetic - m_corsika);
+      if (m_err_abs >= mass_tolerance_ * m_corsika) {
+        const HEPEnergyType e_shifted =
+            sqrt(p_original.getSquaredNorm() + m_corsika * m_corsika);
+        auto const e_shift_relative = (e_shifted / e_original - 1);
+        count_ = count_ + 1;
+        average_shift_ += abs(e_shift_relative);
+        if (abs(e_shift_relative) > max_shift_) max_shift_ = abs(e_shift_relative);
+        CORSIKA_LOGGER_TRACE(
+            logger_,
+            "shift particle mass for {} \n"
+            "{:>45} {:7.5f} \n"
+            "{:>45} {:7.5f} \n"
+            "{:>45} {:7.5f} \n"
+            "{:>45} {:7.5f} \n",
+            pid, "corsika mass (GeV):", m_corsika / 1_GeV,
+            "kinetic mass (GeV): ", m_kinetic / 1_GeV,
+            "m_kin-m_cor (GeV): ", m_err_abs / 1_GeV,
+            "mass tolerance (GeV): ", (m_corsika * mass_tolerance_) / 1_GeV);
+        /*
+          For now we warn if the necessary shift is larger than 1%.
+          we could promote this to an error.
+        */
+        if (abs(e_shift_relative) > energy_tolerance_) {
+          logger_->warn("warning! shifted particle energy by {} %",
+                        e_shift_relative * 100);
+          if (throw_error_) {
+            throw std::runtime_error(
+                "OnShellCheck: error! shifted energy by large amount!");
+          }
+        }
+
+        // reset energy
+        p.setEnergy(e_shifted);
+      } else {
+        CORSIKA_LOGGER_DEBUG(logger_, "particle mass for {} OK", pid);
+      }
+    }
+  } // namespace corsika
+
+} // namespace corsika

corsika/detail/modules/ParticleCut.inl

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/core/Logging.hpp>
+
+namespace corsika {
+
+  template <typename TOutput>
+  template <typename... TArgs>
+  inline ParticleCut<TOutput>::ParticleCut(HEPEnergyType const eEleCut,
+                                           HEPEnergyType const ePhoCut,
+                                           HEPEnergyType const eHadCut,
+                                           HEPEnergyType const eMuCut,
+                                           HEPEnergyType const eTauCut, bool const inv,
+                                           TArgs&&... outputArgs)
+      : TOutput(std::forward<TArgs>(outputArgs)...)
+      , cut_electrons_(eEleCut)
+      , cut_photons_(ePhoCut)
+      , cut_hadrons_(eHadCut)
+      , cut_muons_(eMuCut)
+      , cut_tau_(eTauCut)
+      , doCutInv_(inv) {
+    for (auto const p : get_all_particles()) {
+      if (is_hadron(p)) // nuclei are also hadrons
+        set_kinetic_energy_propagation_threshold(p, eHadCut);
+      else if (is_muon(p))
+        set_kinetic_energy_propagation_threshold(p, eMuCut);
+      else if (p == Code::TauMinus || p == Code::TauPlus)
+        set_kinetic_energy_propagation_threshold(p, eTauCut);
+      else if (p == Code::Electron || p == Code::Positron)
+        set_kinetic_energy_propagation_threshold(p, eEleCut);
+      else if (p == Code::Photon)
+        set_kinetic_energy_propagation_threshold(p, ePhoCut);
+    }
+    set_kinetic_energy_propagation_threshold(Code::Nucleus, eHadCut);
+    CORSIKA_LOG_DEBUG(
+        "setting kinetic energy thresholds: electrons = {} GeV, photons = {} GeV, "
+        "hadrons = {} GeV, "
+        "muons = {} GeV",
+        "tau = {} GeV", eEleCut / 1_GeV, ePhoCut / 1_GeV, eHadCut / 1_GeV, eMuCut / 1_GeV,
+        eTauCut / 1_GeV);
+  }
+
+  template <typename TOutput>
+  template <typename... TArgs>
+  inline ParticleCut<TOutput>::ParticleCut(HEPEnergyType const eCut, bool const inv,
+                                           TArgs&&... outputArgs)
+      : TOutput(std::forward<TArgs>(outputArgs)...)
+      , doCutInv_(inv) {
+    for (auto p : get_all_particles()) {
+      set_kinetic_energy_propagation_threshold(p, eCut);
+    }
+    set_kinetic_energy_propagation_threshold(Code::Nucleus, eCut);
+    CORSIKA_LOG_DEBUG("setting kinetic energy threshold {} GeV", eCut / 1_GeV);
+  }
+
+  template <typename TOutput>
+  template <typename... TArgs>
+  inline ParticleCut<TOutput>::ParticleCut(
+      std::unordered_map<Code const, HEPEnergyType const> const& eCuts, bool const inv,
+      TArgs&&... args)
+      : TOutput(std::forward<TArgs>(args)...)
+      , doCutInv_(inv) {
+    for (auto const& cut : eCuts) {
+      set_kinetic_energy_propagation_threshold(cut.first, cut.second);
+    }
+    CORSIKA_LOG_DEBUG("setting threshold particles individually");
+  }
+
+  template <typename TOutput>
+  inline bool ParticleCut<TOutput>::isBelowEnergyCut(
+      Code const pid, HEPEnergyType const energyLab) const {
+    // nuclei
+    if (is_nucleus(pid)) {
+      // calculate energy per nucleon
+      auto const ElabNuc = energyLab / get_nucleus_A(pid);
+      return (ElabNuc < get_kinetic_energy_propagation_threshold(pid));
+    } else {
+      return (energyLab < get_kinetic_energy_propagation_threshold(pid));
+    }
+  }
+
+  template <typename TOutput>
+  inline bool ParticleCut<TOutput>::checkCutParticle(Code const pid,
+                                                     HEPEnergyType const kine_energy,
+                                                     TimeType const timePost) const {
+
+    HEPEnergyType const energy = kine_energy + get_mass(pid);
+    CORSIKA_LOG_DEBUG(
+        "ParticleCut: checking {} ({}), E_kin= {} GeV, E={} GeV, m={} "
+        "GeV",
+        pid, get_PDG(pid), kine_energy / 1_GeV, energy / 1_GeV, get_mass(pid) / 1_GeV);
+    if (doCutInv_ && is_neutrino(pid)) {
+      CORSIKA_LOG_DEBUG("removing inv. particle...");
+      return true;
+    } else if (isBelowEnergyCut(pid, kine_energy)) {
+      CORSIKA_LOG_DEBUG("removing low en. particle...");
+      return true;
+    } else if (timePost > 10_ms) {
+      CORSIKA_LOG_DEBUG("removing OLD particle...");
+      return true;
+    } else {
+      for (auto const& cut : cuts_) {
+        if (pid == cut.first && kine_energy < cut.second) { return true; }
+      }
+    }
+    return false; // this particle will not be removed/cut
+  }
+
+  template <typename TOutput>
+  template <typename TStackView>
+  inline void ParticleCut<TOutput>::doSecondaries(TStackView& vS) {
+    HEPEnergyType energy_event = 0_GeV; // per event counting for printout
+    auto particle = vS.begin();
+    while (particle != vS.end()) {
+      Code pid = particle.getPID();
+      HEPEnergyType Ekin = particle.getKineticEnergy();
+      if (checkCutParticle(pid, Ekin, particle.getTime())) {
+        this->write(particle.getPosition(), pid, particle.getWeight() * Ekin);
+        particle.erase();
+      }
+      ++particle; // next entry in SecondaryView
+    }
+    CORSIKA_LOG_DEBUG("Event cut: {} GeV", energy_event / 1_GeV);
+  }
+
+  template <typename TOutput>
+  template <typename TParticle>
+  inline ProcessReturn ParticleCut<TOutput>::doContinuous(Step<TParticle>& step,
+                                                          bool const) {
+    if (checkCutParticle(step.getParticlePre().getPID(), step.getEkinPost(),
+                         step.getTimePost())) {
+      this->write(
+          step.getPositionPost(), step.getParticlePre().getPID(),
+          step.getParticlePre().getWeight() *
+              step.getEkinPost()); // ToDO: should the cut happen at the start of the
+                                   // track? For now, I set it to happen at the start
+      CORSIKA_LOG_TRACE("removing during continuous");
+      // signal to upstream code that this particle was deleted
+      return ProcessReturn::ParticleAbsorbed;
+    }
+    return ProcessReturn::Ok;
+  }
+
+  template <typename TOutput>
+  inline void ParticleCut<TOutput>::printThresholds() const {
+
+    CORSIKA_LOG_DEBUG("kinetic energy threshold for electrons is {} GeV",
+                      cut_electrons_ / 1_GeV);
+    CORSIKA_LOG_DEBUG("kinetic energy threshold for photons is {} GeV",
+                      cut_photons_ / 1_GeV);
+    CORSIKA_LOG_DEBUG("kinetic energy threshold for muons is {} GeV", cut_muons_ / 1_GeV);
+    CORSIKA_LOG_DEBUG("kinetic energy threshold for tau is {} GeV", cut_tau_ / 1_GeV);
+    CORSIKA_LOG_DEBUG("kinetic energy threshold for hadrons is {} GeV",
+                      cut_hadrons_ / 1_GeV);
+
+    for (auto const& cut : cuts_) {
+      CORSIKA_LOG_DEBUG("kinetic energy threshold for particle {} is {} GeV", cut.first,
+                        cut.second / 1_GeV);
+    }
+  }
+
+  template <typename TOutput>
+  inline YAML::Node ParticleCut<TOutput>::getConfig() const {
+
+    YAML::Node node;
+    node["type"] = "ParticleCut";
+    node["units"]["energy"] = "GeV";
+    node["cut_electrons"] = cut_electrons_ / 1_GeV;
+    node["cut_photons"] = cut_photons_ / 1_GeV;
+    node["cut_muons"] = cut_muons_ / 1_GeV;
+    node["cut_hadrons"] = cut_hadrons_ / 1_GeV;
+    node["cut_tau"] = cut_tau_ / 1_GeV;
+    node["cut_invisibles"] = doCutInv_;
+    for (auto const& cut : cuts_) {
+      node[fmt::format("cut_{}", cut.first)] = cut.second / 1_GeV;
+    }
+    return node;
+  }
+
+} // namespace corsika

corsika/detail/modules/StackInspector.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/StackInspector.hpp>
+
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/geometry/RootCoordinateSystem.hpp>
+
+#include <chrono>
+#include <iomanip>
+#include <iostream>
+#include <limits>
+
+namespace corsika {
+
+  template <typename TStack>
+  inline StackInspector<TStack>::StackInspector(int const vNStep, bool const vReportStack,
+                                                HEPEnergyType const vE0)
+      : StackProcess<StackInspector<TStack>>(vNStep)
+      , ReportStack_(vReportStack)
+      , E0_(vE0)
+      , StartTime_(std::chrono::system_clock::now())
+      , energyPostInit_(HEPEnergyType::zero())
+      , timePostInit_(std::chrono::system_clock::now()) {}
+
+  template <typename TStack>
+  inline StackInspector<TStack>::~StackInspector() {}
+
+  template <typename TStack>
+  inline void StackInspector<TStack>::doStack(TStack const& vS) {
+
+    [[maybe_unused]] int i = 0;
+    HEPEnergyType Etot = 0_GeV;
+
+    for (auto const& iterP : vS) {
+      HEPEnergyType E = iterP.getEnergy();
+      Etot += E;
+      if (ReportStack_) {
+        CoordinateSystemPtr const& rootCS = get_root_CoordinateSystem(); // for printout
+        auto pos = iterP.getPosition().getCoordinates(rootCS);
+        CORSIKA_LOG_INFO(
+            "StackInspector: i= {:5d}"
+            ", id= {:^15}"
+            " E={:15e} GeV, "
+            " pos= {}"
+            " node = {}",
+            (i++), iterP.getPID(), (E / 1_GeV), pos, fmt::ptr(iterP.getNode()));
+      }
+    }
+
+    if ((E0_ - Etot) < dE_threshold_) return;
+
+    // limit number of printouts
+    if (PrintoutCounter_ < MaxNumberOfPrintouts_) {
+      std::chrono::system_clock::time_point const now = std::chrono::system_clock::now();
+      std::chrono::duration<double> const elapsed_seconds = now - StartTime_; // seconds
+
+      // Select reference times and energies using either the true start
+      // or the delayed start to avoid counting overhead in ETA
+      bool const usePostVals = (energyPostInit_ != HEPEnergyType::zero());
+      auto const dE = (usePostVals ? energyPostInit_ : E0_) - Etot;
+      std::chrono::duration<double> const usedSeconds = now - timePostInit_;
+      double const dT = usedSeconds.count();
+
+      double const progress = (E0_ - Etot) / E0_;
+      // for printout
+      std::time_t const now_time = std::chrono::system_clock::to_time_t(now);
+
+      double const ETA_seconds = (1.0 - progress) * dT * (E0_ / dE);
+      std::chrono::system_clock::time_point const ETA =
+          now + std::chrono::seconds((int)ETA_seconds);
+
+      std::time_t const ETA_time = std::chrono::system_clock::to_time_t(ETA);
+
+      int const yday0 = std::localtime(&now_time)->tm_yday;
+      int const yday1 = std::localtime(&ETA_time)->tm_yday;
+      int const dyday = yday1 - yday0;
+
+      std::stringstream ETA_string;
+      ETA_string << std::put_time(std::localtime(&ETA_time), "%T %Z");
+
+      CORSIKA_LOG_DEBUG(
+          "StackInspector: "
+          " time={}"
+          ", running={:.1f} seconds"
+          " ( {:.1f}%)"
+          ", nStep={}"
+          ", stackSize={}"
+          ", Estack={:.1f} GeV"
+          ", ETA={}{}",
+          std::put_time(std::localtime(&now_time), "%T %Z"), elapsed_seconds.count(),
+          (progress * 100), getStep(), vS.getSize(), Etot / 1_GeV,
+          (dyday == 0 ? "" : fmt::format("+{}d ", dyday)), ETA_string.str());
+
+      PrintoutCounter_++;
+
+      if (PrintoutCounter_ == MaxNumberOfPrintouts_) {
+        CORSIKA_LOG_DEBUG("StackInspector reached allowed maximum of {} lines printout",
+                          MaxNumberOfPrintouts_);
+      }
+
+      // Change reference time once the shower has begin (avoid counting overhead time)
+      if (progress > 0.02 && energyPostInit_ == HEPEnergyType::zero()) {
+        energyPostInit_ = Etot;
+        timePostInit_ = std::chrono::system_clock::now();
+      }
+    }
+  }
+
+} // namespace corsika

corsika/detail/modules/TimeCut.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/TimeCut.hpp>
+
+namespace corsika {
+
+  inline TimeCut::TimeCut(const TimeType time)
+      : time_(time) {}
+
+  template <typename Particle>
+  inline ProcessReturn TimeCut::doContinuous(Step<Particle> const& step, bool const) {
+    CORSIKA_LOG_TRACE("TimeCut::doContinuous");
+    if (step.getTimePost() >= time_) {
+      CORSIKA_LOG_TRACE("stopping continuous process");
+      return ProcessReturn::ParticleAbsorbed;
+    }
+    return ProcessReturn::Ok;
+  }
+
+} // namespace corsika

corsika/detail/modules/TrackWriter.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>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+
+#include <limits>
+
+namespace corsika {
+
+  template <typename TOutput>
+  inline TrackWriter<TOutput>::TrackWriter() {}
+
+  template <typename TOutput>
+  template <typename TParticle>
+  inline ProcessReturn TrackWriter<TOutput>::doContinuous(Step<TParticle> const& step,
+                                                          bool const) {
+
+    auto const start = step.getPositionPre().getCoordinates();
+    auto const end = step.getPositionPost().getCoordinates();
+
+    // write the track to the file
+    TOutput::write(step.getParticlePre().getPID(), step.getEkinPre(),
+                   step.getParticlePre().getWeight(), start, step.getTimePre(), end,
+                   step.getEkinPost(), step.getTimePost());
+
+    return ProcessReturn::Ok;
+  }
+
+  template <typename TOutput>
+  template <typename TParticle, typename TTrack>
+  inline LengthType TrackWriter<TOutput>::getMaxStepLength(TParticle const&,
+                                                           TTrack const&) {
+    return meter * std::numeric_limits<double>::infinity();
+  }
+
+  template <typename TOutput>
+  YAML::Node TrackWriter<TOutput>::getConfig() const {
+    using namespace units::si;
+
+    YAML::Node node;
+
+    // add default units for values
+    node["type"] = "TrackWriter";
+    node["units"] = "GeV | m | ns";
+
+    return node;
+  }
+} // namespace corsika