corsika/detail/media/SlidingPlanarExponential.inl

 /*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
  *
- * See file AUTHORS for a list of contributors.
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
+ * 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/SlidingPlanarExponential.hpp>
-
 namespace corsika {
 
-  template <class T>
-  units::si::MassDensityType SlidingPlanarExponential<T>::GetMassDensity(
-      Point const& p) const {
-    auto const height = (p - Base::fP0).norm() - referenceHeight_;
-    return Base::fRho0 * exp(Base::fInvLambda * height);
+  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 <class T>
-  units::si::GrammageType SlidingPlanarExponential<T>::IntegratedGrammage(
-      Trajectory<Line> const& line, units::si::LengthType l) const {
-    auto const axis = (line.GetR0() - Base::fP0).normalized();
-    return Base::IntegratedGrammage(line, l, axis);
+  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 <class T>
-  units::si::LengthType SlidingPlanarExponential<T>::ArclengthFromGrammage(
-      Trajectory<Line> const& line, units::si::GrammageType grammage) const {
-    auto const axis = (line.GetR0() - Base::fP0).normalized();
-    return Base::ArclengthFromGrammage(line, grammage, 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
\ No newline at end of file
+} // 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
  *
- * See file AUTHORS for a list of contributors.
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
+ * 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
@@ -14,12 +11,10 @@
 
 namespace corsika {
 
-  Universe::Universe(corsika::CoordinateSystem const& pCS)
-      : corsika::Sphere(
-            corsika::Point{pCS, 0 * corsika::units::si::meter,
-                           0 * corsika::units::si::meter, 0 * corsika::units::si::meter},
-            corsika::units::si::meter * std::numeric_limits<double>::infinity()) {}
+  inline Universe::Universe(CoordinateSystemPtr const& pCS)
+      : corsika::Sphere(Point{pCS, 0 * meter, 0 * meter, 0 * meter},
+                        meter * std::numeric_limits<double>::infinity()) {}
 
-  bool Universe::Contains(corsika::Point const&) const { return true; }
+  inline bool Universe::contains(corsika::Point const&) const { return true; }
 
-} // namespace corsika
\ No newline at end of file
+} // namespace corsika

corsika/detail/media/VolumeTreeNode.inl

 /*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
  *
- * See file AUTHORS for a list of contributors.
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
+ * 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/Volume.hpp>
+#include <corsika/framework/geometry/IVolume.hpp>
 #include <corsika/media/IMediumModel.hpp>
 
 namespace corsika {
 
-  //! convenience function equivalent to Volume::Contains
   template <typename IModelProperties>
-  bool VolumeTreeNode<IModelProperties>::Contains(corsika::Point const& p) const {
-    return fGeoVolume->Contains(p);
+  inline bool VolumeTreeNode<IModelProperties>::contains(Point const& p) const {
+    return geoVolume_->contains(p);
   }
 
   template <typename IModelProperties>
   inline VolumeTreeNode<IModelProperties> const*
-  VolumeTreeNode<IModelProperties>::Excludes(corsika::Point const& p) const {
+  VolumeTreeNode<IModelProperties>::excludes(Point const& p) const {
     auto exclContainsIter =
-        std::find_if(fExcludedNodes.cbegin(), fExcludedNodes.cend(),
-                     [&](auto const& s) { return bool(s->Contains(p)); });
+        std::find_if(excludedNodes_.cbegin(), excludedNodes_.cend(),
+                     [&](auto const& s) { return bool(s->contains(p)); });
 
-    return exclContainsIter != fExcludedNodes.cend() ? *exclContainsIter : nullptr;
+    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>
-  VolumeTreeNode<IModelProperties> const*
-  VolumeTreeNode<IModelProperties>::GetContainingNode(corsika::Point const& p) const {
-    if (!Contains(p)) { return nullptr; }
+  inline VolumeTreeNode<IModelProperties> const*
+  VolumeTreeNode<IModelProperties>::getContainingNode(Point const& p) const {
+    if (!contains(p)) { return nullptr; }
 
     if (auto const childContainsIter =
-            std::find_if(fChildNodes.cbegin(), fChildNodes.cend(),
-                         [&](auto const& s) { return bool(s->Contains(p)); });
-        childContainsIter == fChildNodes.cend()) // not contained in any of the children
+            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
+      if (auto const exclContainsIter = excludes(p)) // contained in any excluded nodes
       {
-        return exclContainsIter->GetContainingNode(p);
+        return exclContainsIter->getContainingNode(p);
       } else {
         return this;
       }
     } else {
-      return (*childContainsIter)->GetContainingNode(p);
+      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>
-  void VolumeTreeNode<IModelProperties>::walk(TCallable func) {
+  inline void VolumeTreeNode<IModelProperties>::walk(TCallable func) const {
     if constexpr (preorder) { func(*this); }
 
-    std::for_each(fChildNodes.begin(), fChildNodes.end(),
-                  [&](auto& v) { v->walk(func); });
+    std::for_each(childNodes_.begin(), childNodes_.end(),
+                  [&](auto const& v) { v->walk(func); });
 
     if constexpr (!preorder) { func(*this); };
   }
 
   template <typename IModelProperties>
-  void VolumeTreeNode<IModelProperties>::AddChild(typename VolumeTreeNode<IModelProperties>::VTNUPtr pChild) {
-    pChild->fParentNode = this;
-    fChildNodes.push_back(std::move(pChild));
+  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>
-  void VolumeTreeNode<IModelProperties>::ExcludeOverlapWith(typename VolumeTreeNode<IModelProperties>::VTNUPtr const& pNode) {
-    fExcludedNodes.push_back(pNode.get());
-  }
-
-  template <typename IModelProperties>
-  template <class MediumType, typename... Args>
-  auto VolumeTreeNode<IModelProperties>::CreateMedium(Args&&... args) {
-    static_assert(std::is_base_of_v<IMediumModel, MediumType>,
-                  "unusable type provided, needs to be derived from \"IMediumModel\"");
-
-    return std::make_shared<MediumType>(std::forward<Args>(args)...);
+  inline void VolumeTreeNode<IModelProperties>::excludeOverlapWith(
+      typename VolumeTreeNode<IModelProperties>::VTNUPtr const& pNode) {
+    excludedNodes_.push_back(pNode.get());
   }
 
-} // namespace corsika
\ No newline at end of file
+} // 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

corsika/detail/modules/conex/CONEXhybrid.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/Logging.hpp>
+#include <corsika/media/CORSIKA7Atmospheres.hpp>
+#include <corsika/modules/conex/CONEXhybrid.hpp>
+#include <corsika/modules/conex/CONEXrandom.hpp>
+#include <corsika/modules/Random.hpp>
+#include <corsika/modules/conex/CONEX_f.hpp>
+#include <corsika/framework/random/RNGManager.hpp>
+#include <corsika/framework/core/PhysicalConstants.hpp>
+#include <corsika/framework/core/Logging.hpp>
+
+#include <conexConfig.h>
+
+#include <algorithm>
+#include <fstream>
+#include <numeric>
+#include <utility>
+
+namespace corsika {
+
+  template <typename TOutputE, typename TOutputN>
+  inline CONEXhybrid<TOutputE, TOutputN>::CONEXhybrid(
+      Point const& center, ShowerAxis const& showerAxis, LengthType groundDist,
+      LengthType injectionHeight, HEPEnergyType primaryEnergy, PDGCode primaryPDG,
+      TOutputE& args1, TOutputN& args2)
+      : SubWriter<TOutputE>(args1)
+      , SubWriter<TOutputN>(args2)
+      , center_{center}
+      , showerAxis_{showerAxis}
+      , groundDist_{groundDist}
+      , injectionHeight_{injectionHeight}
+      , primaryEnergy_{primaryEnergy}
+      , primaryPDG_{primaryPDG}
+      , showerCore_{showerAxis_.getStart() + showerAxis_.getDirection() * groundDist_}
+      , conexObservationCS_{std::invoke([&]() {
+        auto const& c8cs = center.getCoordinateSystem();
+        auto const translation = showerCore_ - center;
+        auto const intermediateCS =
+            make_translation(c8cs, translation.getComponents(c8cs));
+        auto const transformCS = make_rotationToZ(intermediateCS, translation);
+
+        CORSIKA_LOG_DEBUG("translation C8/CONEX obs: ", translation.getComponents());
+
+        /*
+        auto const transform = CoordinateSystem::getTransformation(
+            intermediateCS2, c8cs); // either this way or vice versa... TODO: test this!
+        */
+        return transformCS;
+      })}
+      , x_sf_{std::invoke([&]() {
+        Vector<length_d> const a{conexObservationCS_, 0._m, 0._m, 1._m};
+        auto b = a.cross(showerAxis_.getDirection());
+        auto const lengthB = b.getNorm();
+        if (lengthB < 1e-10_m) {
+          b = Vector<length_d>{conexObservationCS_, 1_m, 0_m, 0_m};
+        }
+
+        return b.normalized();
+      })}
+      , y_sf_{showerAxis_.getDirection().cross(x_sf_)} {
+
+    CORSIKA_LOG_DEBUG("x_sf (conexObservationCS): {}",
+                      x_sf_.getComponents(conexObservationCS_));
+    CORSIKA_LOG_DEBUG("x_sf (C8): {}", x_sf_.getComponents(center.getCoordinateSystem()));
+
+    CORSIKA_LOG_DEBUG("y_sf (conexObservationCS): {}",
+                      y_sf_.getComponents(conexObservationCS_));
+
+    CORSIKA_LOG_DEBUG("y_sf (C8): {}", y_sf_.getComponents(center.getCoordinateSystem()));
+
+    CORSIKA_LOG_DEBUG("showerAxisDirection (conexObservationCS): {}",
+                      showerAxis_.getDirection().getComponents(conexObservationCS_));
+    CORSIKA_LOG_DEBUG(
+        "showerAxisDirection (C8): {}",
+        showerAxis_.getDirection().getComponents(center.getCoordinateSystem()));
+
+    CORSIKA_LOG_DEBUG("showerCore (conexObservationCS): {}",
+                      showerCore_.getCoordinates(conexObservationCS_));
+    CORSIKA_LOG_DEBUG("showerCore (C8): {}",
+                      showerCore_.getCoordinates(center.getCoordinateSystem()));
+
+    auto const& components = ::corsika::standardAirComposition.getComponents();
+    auto const& fractions = ::corsika::standardAirComposition.getFractions();
+    if (::corsika::standardAirComposition.getSize() != 3) {
+      throw std::runtime_error{"CONEXhybrid only usable with standard 3-component air"};
+    }
+
+    std::transform(components.cbegin(), components.cend(), ::conex::cxair_.aira.begin(),
+                   get_nucleus_A);
+    std::transform(components.cbegin(), components.cend(), ::conex::cxair_.airz.begin(),
+                   get_nucleus_Z);
+    std::copy(fractions.cbegin(), fractions.cend(), ::conex::cxair_.airw.begin());
+
+    ::conex::cxair_.airava =
+        std::inner_product(::conex::cxair_.airw.cbegin(), ::conex::cxair_.airw.cend(),
+                           ::conex::cxair_.aira.cbegin(), 0.);
+    ::conex::cxair_.airavz =
+        std::inner_product(::conex::cxair_.airw.cbegin(), ::conex::cxair_.airw.cend(),
+                           ::conex::cxair_.airz.cbegin(), 0.);
+
+    // this is the CONEX default but actually unused there
+    ::conex::cxair_.airi = {82.0e-09, 95.0e-09, 188.e-09};
+
+    int randomSeeds[3] = {1234, 0, 0}; // SEEDS ARE NOT USED. All random numbers are
+                                       // obtained from the CORSIKA 8 stream "conex"
+    corsika::connect_random_stream("conex", ::conex::set_rng_function);
+    int heModel = eSibyll23;
+
+    int nShower = 1; // large to avoid final stats.
+    int maxDetail = 0;
+#ifdef CONEX_EXTENSIONS
+    int particleListMode = 0;
+#endif
+
+    std::string configPath = CONEX_CONFIG_PATH;
+    ::conex::initconex_(nShower, randomSeeds, heModel, maxDetail,
+#ifdef CONEX_EXTENSIONS
+                        particleListMode,
+#endif
+                        configPath.c_str(), configPath.size());
+  }
+
+  template <typename TOutputE, typename TOutputN>
+  inline void CONEXhybrid<TOutputE, TOutputN>::initCascadeEquations() {
+
+    // set phi, theta
+    Vector<length_d> ez{conexObservationCS_, {0._m, 0._m, -1_m}};
+    auto const c = showerAxis_.getDirection().dot(ez) / 1_m;
+    double theta = std::acos(c) * 180 / M_PI;
+
+    auto const showerAxisConex =
+        showerAxis_.getDirection().getComponents(conexObservationCS_);
+    double phi = std::atan2(-showerAxisConex.getY().magnitude(),
+                            showerAxisConex.getX().magnitude()) *
+                 180 / M_PI;
+
+    CORSIKA_LOG_DEBUG(
+        "theta (deg) = {}"
+        "; phi (deg) = {}",
+        theta, phi);
+
+    int ipart = static_cast<int>(primaryPDG_);
+
+    double dimpact = 0.; // valid only if shower core is fixed on the observation plane;
+                         // for skimming showers an offset is needed like in CONEX
+
+    // SEEDS ARE NOT USED. All random numbers are obtained from
+    // the CORSIKA 8 stream "conex" and "epos"!
+    std::array<int, 3> ioseed{1, 1, 1};
+
+    double eprima = primaryEnergy_ / 1_GeV;
+    double xminp = injectionHeight_ / 1_m;
+
+    ::conex::conexrun_(ipart, eprima, theta, phi, xminp, dimpact, ioseed.data());
+  }
+
+  template <typename TOutputE, typename TOutputN>
+  template <typename TStackView>
+  inline void CONEXhybrid<TOutputE, TOutputN>::doSecondaries(TStackView& vS) {
+    auto p = vS.begin();
+    while (p != vS.end()) {
+      Code const pid = p.getPID();
+      if (addParticle(pid, p.getEnergy(), p.getMass(), p.getPosition(),
+                      p.getMomentum().normalized(), p.getTime())) {
+        p.erase();
+      }
+      ++p;
+    }
+  }
+
+  template <typename TOutputE, typename TOutputN>
+  inline bool CONEXhybrid<TOutputE, TOutputN>::addParticle(
+      Code pid, HEPEnergyType energy, HEPEnergyType mass, Point const& position,
+      DirectionVector const& direction, TimeType t, double weight) {
+
+    auto const it = std::find_if(egs_em_codes_.cbegin(), egs_em_codes_.cend(),
+                                 [=](auto const& p) { return pid == p.first; });
+    if (it == egs_em_codes_.cend()) { return false; }
+
+    // EM particle
+    auto const egs_pid = it->second;
+    CORSIKA_LOG_DEBUG("position conexObs: {}",
+                      position.getCoordinates(conexObservationCS_));
+
+    auto const coords = position.getCoordinates(conexObservationCS_) / 1_m;
+    double const x = coords[0].magnitude();
+    double const y = coords[1].magnitude();
+
+    double const altitude = ((position - center_).getNorm() - conex::earthRadius) / 1_m;
+    auto const d = position - showerCore_;
+
+    // distance from core to particle projected along shower axis
+    double const slantDistance = -d.dot(showerAxis_.getDirection()) / 1_m;
+
+    // lateral coordinates in CONEX shower frame
+    auto const dShowerPlane = d - d.getParallelProjectionOnto(showerAxis_.getDirection());
+    double const lateralX = dShowerPlane.dot(x_sf_) / 1_m;
+    double const lateralY = dShowerPlane.dot(y_sf_) / 1_m;
+
+    double const slantX = showerAxis_.getProjectedX(position) * (1_cm * 1_cm / 1_g);
+
+    double const time = (t * constants::c - groundDist_) / 1_m;
+
+    // fill u,v,w momentum direction in EGS frame
+    double const u = direction.dot(y_sf_).magnitude();
+    double const v = direction.dot(x_sf_).magnitude();
+    double const w = direction.dot(showerAxis_.getDirection()).magnitude();
+
+    // generation, TO BE CHANGED WHEN WE HAVE THAT INFORMATION AVAILABLE
+    int const latchin = 1;
+
+    double const E = energy / 1_GeV;
+    double const m = mass / 1_GeV;
+
+    CORSIKA_LOG_DEBUG("CONEXhybrid: removing {} {:5e} GeV", egs_pid, energy);
+
+    CORSIKA_LOG_DEBUG("#### parameters to cegs4_() ####");
+    CORSIKA_LOG_DEBUG("egs_pid = {}", egs_pid);
+    CORSIKA_LOG_DEBUG("E = {}", E);
+    CORSIKA_LOG_DEBUG("m = {}", m);
+    CORSIKA_LOG_DEBUG("x = {}", x);
+    CORSIKA_LOG_DEBUG("y = {}", y);
+    CORSIKA_LOG_DEBUG("altitude = {}", altitude);
+    CORSIKA_LOG_DEBUG("slantDistance = {}", slantDistance);
+    CORSIKA_LOG_DEBUG("lateralX = {}", lateralX);
+    CORSIKA_LOG_DEBUG("lateralY = {}", lateralY);
+    CORSIKA_LOG_DEBUG("slantX = {}", slantX);
+    CORSIKA_LOG_DEBUG("time = {}", time);
+    CORSIKA_LOG_DEBUG("u = {}", u);
+    CORSIKA_LOG_DEBUG("v = {}", v);
+    CORSIKA_LOG_DEBUG("w = {}", w);
+
+    ::conex::cxoptl_.dptl[10 - 1] = egs_pid;
+    ::conex::cxoptl_.dptl[4 - 1] = E;
+    ::conex::cxoptl_.dptl[5 - 1] = m;
+    ::conex::cxoptl_.dptl[6 - 1] = x;
+    ::conex::cxoptl_.dptl[7 - 1] = y;
+    ::conex::cxoptl_.dptl[8 - 1] = altitude;
+    ::conex::cxoptl_.dptl[9 - 1] = time;
+    ::conex::cxoptl_.dptl[11 - 1] = weight;
+    ::conex::cxoptl_.dptl[12 - 1] = latchin;
+    ::conex::cxoptl_.dptl[13 - 1] = slantX;
+    ::conex::cxoptl_.dptl[14 - 1] = lateralX;
+    ::conex::cxoptl_.dptl[15 - 1] = lateralY;
+    ::conex::cxoptl_.dptl[16 - 1] = slantDistance;
+    ::conex::cxoptl_.dptl[2 - 1] = u;
+    ::conex::cxoptl_.dptl[1 - 1] = v;
+    ::conex::cxoptl_.dptl[3 - 1] = w;
+
+    int n = 1, i = 1;
+    ::conex::cegs4_(n, i);
+
+    return true;
+  }
+
+  template <typename TOutputE, typename TOutputN>
+  template <typename TStack>
+  inline void CONEXhybrid<TOutputE, TOutputN>::doCascadeEquations(TStack&) {
+
+    ::conex::conexcascade_();
+
+    int nX = ::conex::get_number_of_depth_bins_(); // make sure this works!
+
+    int icut = 1;
+    int icutg = 2;
+    int icute = 3;
+    int icutm = 2;
+    int icuth = 3;
+    int iSec = 0;
+
+    const int maxX = nX;
+
+    auto X = std::make_unique<float[]>(maxX);
+    auto H = std::make_unique<float[]>(maxX);
+    auto D = std::make_unique<float[]>(maxX);
+    auto N = std::make_unique<float[]>(maxX);
+    auto dEdX = std::make_unique<float[]>(maxX);
+    auto Mu = std::make_unique<float[]>(maxX);
+    auto dMu = std::make_unique<float[]>(maxX);
+    auto Photon = std::make_unique<float[]>(maxX);
+    auto Electrons = std::make_unique<float[]>(maxX);
+    auto Hadrons = std::make_unique<float[]>(maxX);
+
+    float EGround[3], fitpars[13];
+
+    ::conex::get_shower_data_(icut, iSec, nX, X[0], N[0], fitpars[0], H[0], D[0]);
+    ::conex::get_shower_edep_(icut, nX, dEdX[0], EGround[0]);
+    ::conex::get_shower_muon_(icutm, nX, Mu[0], dMu[0]);
+    ::conex::get_shower_gamma_(icutg, nX, Photon[0]);
+    ::conex::get_shower_electron_(icute, nX, Electrons[0]);
+    ::conex::get_shower_hadron_(icuth, nX, Hadrons[0]);
+
+    // make sure CONEX binning is same to C8:
+    GrammageType const dX = (X[1] - X[0]) * (1_g / square(1_cm));
+
+    for (int i = 0; i < nX; ++i) {
+      GrammageType const curX = X[i] * (1_g / square(1_cm));
+      SubWriter<TOutputE>::write(curX, curX + dX,
+                                 Code::Unknown, // this is sum of all dEdX
+                                 dEdX[i] * 1_GeV / 1_g * square(1_cm) * dX);
+      SubWriter<TOutputN>::write(curX, curX + dX, Code::Photon, Photon[i]);
+      SubWriter<TOutputN>::write(curX, curX + dX, Code::Proton, Hadrons[i]);
+      SubWriter<TOutputN>::write(curX, curX + dX, Code::Electron, Electrons[i]);
+      SubWriter<TOutputN>::write(curX, curX + dX, Code::MuMinus, Mu[i]);
+    }
+  }
+
+  template <typename TOutputE, typename TOutputN>
+  inline YAML::Node CONEXhybrid<TOutputE, TOutputN>::getConfig() const {
+
+    return YAML::Node();
+  }
+
+} // namespace corsika

corsika/detail/modules/energy_loss/EnergyLoss.inl → corsika/detail/modules/energy_loss/BetheBlochPDG.inl

 /*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
  *
- * See file AUTHORS for a list of contributors.
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
+ * 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/energy_loss/EnergyLoss.hpp>
-
 #include <corsika/framework/core/ParticleProperties.hpp>
 
-#include <corsika/setup/SetupStack.hpp>
-#include <corsika/setup/SetupTrajectory.hpp>
-
 #include <corsika/framework/geometry/Line.hpp>
+#include <corsika/framework/core/Logging.hpp>
 
 #include <cmath>
 #include <fstream>
-#include <iostream>
 #include <limits>
 
-using SetupParticle = corsika::setup::Stack::ParticleType;
-using SetupTrack = corsika::setup::Trajectory;
-
-namespace corsika::energy_loss {
+namespace corsika {
 
-  auto elab2plab = [](corsika::units::si::HEPEnergyType Elab,
-                      corsika::units::si::HEPMassType m) {
-    return sqrt((Elab - m) * (Elab + m));
-  };
+  template <typename TOutput>
+  template <typename... TArgs>
+  inline BetheBlochPDG<TOutput>::BetheBlochPDG(TArgs&&... args)
+      : TOutput(std::forward<TArgs>(args)...) {}
 
-  /**
-   *   PDG2018, passage of particles through matter
-   *
-   * Note, that \f$I_{\mathrm{eff}}\f$ of composite media a determined from \f$ \ln I =
-   * \sum_i a_i \ln(I_i) \f$ where \f$ a_i \f$ is the fraction of the electron population
-   * (\f$\sim Z_i\f$) of the \f$i\f$-th element. This can also be used for shell
-   * corrections or density effects.
-   *
-   * The \f$I_{\mathrm{eff}}\f$ of compounds is not better than a few percent, if not
-   * measured explicitly.
-   *
-   * For shell correction, see Sec 6 of https://www.nap.edu/read/20066/chapter/8#115
-   *
-   */
-  corsika::units::si::HEPEnergyType EnergyLoss::BetheBloch(
-      SetupParticle const& p, corsika::units::si::GrammageType const dX) {
-
-    using namespace corsika::units::si;
+  template <typename TOutput>
+  template <typename TParticle>
+  inline HEPEnergyType BetheBlochPDG<TOutput>::getBetheBloch(TParticle const& p,
+                                                             GrammageType const dX) {
 
     // all these are material constants and have to come through Environment
     // right now: values for nitrogen_D
@@ -69,12 +44,12 @@ namespace corsika::energy_loss {
 
     // this is the Bethe-Bloch coefficiet 4pi N_A r_e^2 m_e c^2
     auto constexpr K = 0.307075_MeV / 1_mol * square(1_cm);
-    HEPEnergyType const E = p.GetEnergy();
-    HEPMassType const m = p.GetMass();
+    HEPEnergyType const E = p.getEnergy();
+    HEPMassType const m = p.getMass();
     double const gamma = E / m;
-    int const Z = p.GetChargeNumber();
+    int const Z = p.getChargeNumber();
     int const Z2 = Z * Z;
-    HEPMassType constexpr me = corsika::Electron::GetMass();
+    HEPMassType constexpr me = Electron::mass;
     auto const m2 = m * m;
     auto constexpr me2 = me * me;
     double const gamma2 = gamma * gamma;
@@ -82,20 +57,20 @@ namespace corsika::energy_loss {
     double const beta2 = (gamma2 - 1) / gamma2; // 1-1/gamma2    (1-1/gamma)*(1+1/gamma);
                                                 // (gamma_2-1)/gamma_2 = (1-1/gamma2);
     double constexpr c2 = 1;                    // HEP convention here c=c2=1
-    std::cout << "BetheBloch beta2=" << beta2 << " gamma2=" << gamma2 << std::endl;
+    CORSIKA_LOG_TRACE("BetheBloch beta2={}, gamma2={}", beta2, gamma2);
     [[maybe_unused]] double const eta2 = beta2 / (1 - beta2);
     HEPMassType const Wmax =
         2 * me * c2 * beta2 * gamma2 / (1 + 2 * gamma * me / m + me2 / m2);
     // approx, but <<1%    HEPMassType const Wmax = 2*me*c2*beta2*gamma2;      for HEAVY
     // PARTICLES Wmax ~ 2me v2 for non-relativistic particles
-    std::cout << "BetheBloch Wmax=" << Wmax << std::endl;
+    CORSIKA_LOG_TRACE("BetheBloch Wmax={}", Wmax);
 
     // Sternheimer parameterization, density corrections towards high energies
     // NOTE/TODO: when Cbar is 0 it needs to be approximated from parameterization ->
     // MISSING
-    std::cout << "BetheBloch p.GetMomentum().GetNorm()/m=" << p.GetMomentum().GetNorm() / m
-         << std::endl;
-    double const x = log10(p.GetMomentum().GetNorm() / m);
+    CORSIKA_LOG_TRACE("BetheBloch p.getMomentum().getNorm()/m={}",
+                      p.getMomentum().getNorm() / m);
+    double const x = log10(p.getMomentum().getNorm() / m);
     double delta = 0;
     if (x >= x1) {
       delta = 2 * (log(10)) * x - Cbar;
@@ -104,7 +79,7 @@ namespace corsika::energy_loss {
     } else if (x < x0) { // and IF conductor (otherwise, this is 0)
       delta = delta0 * pow(100, 2 * (x - x0));
     }
-    std::cout << "BetheBloch delta=" << delta << std::endl;
+    CORSIKA_LOG_TRACE("BetheBloch delta={}", delta);
 
     // with further low energies correction, accurary ~1% down to beta~0.05 (1MeV for p)
 
@@ -124,8 +99,8 @@ namespace corsika::energy_loss {
     // Barkas correction O(Z3) higher-order Born approximation
     // see Appl. Phys. 85 (1999) 1249
     // double A = 1;
-    // if (p.GetPID() == corsika::Code::Nucleus) A = p.GetNuclearA();
-    // double const Erel = (p.GetEnergy()-p.GetMass()) / A / 1_keV;
+    // if (p.getPID() == Code::Nucleus) A = p.getNuclearA();
+    // double const Erel = (p.getEnergy()-p.getMass()) / A / 1_keV;
     // double const Llow = 0.01 * Erel;
     // double const Lhigh = 1.5/pow(Erel, 0.4) + 45000./Zmat * pow(Erel, 1.6);
     // double const barkas = Z * Llow*Lhigh/(Llow+Lhigh); // RU, I think the Z was
@@ -138,146 +113,97 @@ namespace corsika::energy_loss {
     double const y2 = Z * Z * alpha * alpha / beta2;
     double const bloch = -y2 * (1.202 - y2 * (1.042 - 0.855 * y2 + 0.343 * y2 * y2));
 
-    // std::cout << "BetheBloch Erel=" << Erel << " barkas=" << barkas << " bloch=" << bloch <<
-    // std::endl;
-
     double const aux = 2 * me * c2 * beta2 * gamma2 * Wmax / (Ieff * Ieff);
     return -K * Z2 * ZoverA / beta2 *
            (0.5 * log(aux) - beta2 - Cadj / Z - delta / 2 + barkas + bloch) * dX;
   }
 
   // radiation losses according to PDG 2018, ch. 33 ref. [5]
-  corsika::units::si::HEPEnergyType EnergyLoss::RadiationLosses(SetupParticle const& vP,
-                                                                corsika::units::si::GrammageType const vDX) {
-    using namespace corsika::units::si;
+  template <typename TOutput>
+  template <typename TParticle>
+  inline HEPEnergyType BetheBlochPDG<TOutput>::getRadiationLosses(
+      TParticle const& vP, GrammageType const vDX) {
     // simple-minded hard-coded value for b(E) inspired by data from
     // http://pdg.lbl.gov/2018/AtomicNuclearProperties/ for N and O.
     auto constexpr b = 3.0 * 1e-6 * square(1_cm) / 1_g;
-    return -vP.GetEnergy() * b * vDX;
+    return -vP.getEnergy() * b * vDX;
   }
 
-  corsika::units::si::HEPEnergyType EnergyLoss::TotalEnergyLoss(SetupParticle const& vP,
-                                                                corsika::units::si::GrammageType const vDX) {
-    return BetheBloch(vP, vDX) + RadiationLosses(vP, vDX);
+  template <typename TOutput>
+  template <typename TParticle>
+  inline HEPEnergyType BetheBlochPDG<TOutput>::getTotalEnergyLoss(
+      TParticle const& vP, GrammageType const vDX) {
+    return getBetheBloch(vP, vDX) + getRadiationLosses(vP, vDX);
   }
 
-  corsika::EProcessReturn EnergyLoss::DoContinuous(SetupParticle& p,
-                                                   SetupTrack const& t) {
-    using namespace corsika::units::si;
-    if (p.GetChargeNumber() == 0) return corsika::EProcessReturn::eOk;
+  template <typename TOutput>
+  template <typename TParticle>
+  inline ProcessReturn BetheBlochPDG<TOutput>::doContinuous(Step<TParticle>& step,
+                                                            bool const) {
+
+    // if this step was limiting the CORSIKA stepping, the particle is lost
+    /* see Issue https://gitlab.iap.kit.edu/AirShowerPhysics/corsika/-/issues/389
+    if (limitStep) {
+      fillProfile(t, p.getEnergy());
+      p.setEnergy(p.getMass());
+      return ProcessReturn::ParticleAbsorbed;
+    }
+    */
+
+    if (step.getParticlePre().getChargeNumber() == 0) return ProcessReturn::Ok;
 
     GrammageType const dX =
-        p.GetNode()->GetModelProperties().IntegratedGrammage(t, t.GetLength());
-    std::cout << "EnergyLoss " << p.GetPID() << ", z=" << p.GetChargeNumber()
-         << ", dX=" << dX / 1_g * square(1_cm) << "g/cm2" << std::endl;
-    HEPEnergyType dE = TotalEnergyLoss(p, dX);
-    auto E = p.GetEnergy();
-    const auto Ekin = E - p.GetMass();
-    auto Enew = E + dE;
-    std::cout << "EnergyLoss  dE=" << dE / 1_MeV << "MeV, "
-         << " E=" << E / 1_GeV << "GeV,  Ekin=" << Ekin / 1_GeV
-         << ", Enew=" << Enew / 1_GeV << "GeV" << std::endl;
-    auto status = corsika::EProcessReturn::eOk;
-    if (-dE > Ekin) {
-      dE = -Ekin;
-      Enew = p.GetMass();
-      status = corsika::EProcessReturn::eParticleAbsorbed;
-    }
-    p.SetEnergy(Enew);
-    MomentumUpdate(p, Enew);
-    EnergyLossTot_ += dE;
-    FillProfile(p, t, dE);
-    return status;
+        step.getParticlePre().getNode()->getModelProperties().getIntegratedGrammage(
+            step.getStraightTrack());
+    CORSIKA_LOG_TRACE("EnergyLoss pid={}, z={}, dX={} g/cm2",
+                      step.getParticlePre().getPID(),
+                      step.getParticlePre().getChargeNumber(), dX / 1_g * square(1_cm));
+    HEPEnergyType const dE = getTotalEnergyLoss(step.getParticlePre(), dX);
+    //    if (dE > HEPEnergyType::zero())
+    //      dE = -dE;
+    CORSIKA_LOG_TRACE("EnergyLoss  dE={} MeV, Ekin={} GeV, EkinNew={} GeV", dE / 1_MeV,
+                      step.getEkinPre() / 1_GeV, (step.getEkinPre() + dE) / 1_GeV);
+    step.add_dEkin(dE);
+
+    // also send to output
+    TOutput::write(step.getPositionPre(), step.getPositionPost(),
+                   step.getParticlePre().getPID(),
+                   -step.getParticlePre().getWeight() * dE);
+    return ProcessReturn::Ok;
   }
 
-  corsika::units::si::LengthType EnergyLoss::MaxStepLength(
-      SetupParticle const& vParticle, SetupTrack const& vTrack) const {
-    using namespace corsika::units::si;
-    if (vParticle.GetChargeNumber() == 0) {
-      return units::si::meter * std::numeric_limits<double>::infinity();
+  template <typename TOutput>
+  template <typename TParticle, typename TTrajectory>
+  inline LengthType BetheBlochPDG<TOutput>::getMaxStepLength(
+      TParticle const& vParticle, TTrajectory const& vTrack) const {
+    if (vParticle.getChargeNumber() == 0) {
+      return meter * std::numeric_limits<double>::infinity();
     }
 
     auto constexpr dX = 1_g / square(1_cm);
-    auto const dE = -TotalEnergyLoss(vParticle, dX); // dE > 0
-    //~ auto const Ekin = vParticle.GetEnergy() - vParticle.GetMass();
-    auto const maxLoss = 0.01 * vParticle.GetEnergy();
-    auto const maxGrammage = maxLoss / dE * dX;
-
-    return vParticle.GetNode()->GetModelProperties().ArclengthFromGrammage(vTrack,
-                                                                           maxGrammage) *
-           1.0001; // to make sure particle gets absorbed when DoContinuous() is called
+    auto const dEdX = -getTotalEnergyLoss(vParticle, dX) / dX;
+    auto const energy = vParticle.getKineticEnergy();
+    auto const energy_lim =
+        std::max(energy * 0.9, // either 10% relative loss max., or
+                 get_kinetic_energy_propagation_threshold(
+                     vParticle.getPID()) // energy thresholds globally defined
+                                         // for individual particles
+                     * 0.99999 // need to go slightly below global e-cut to assure
+                               // removal in ParticleCut. The 1% does not matter since
+                               // at cut-time the entire energy is removed.
+        );
+    auto const maxGrammage = (energy - energy_lim) / dEdX;
+
+    return vParticle.getNode()->getModelProperties().getArclengthFromGrammage(
+        vTrack, maxGrammage);
   }
 
-  void EnergyLoss::MomentumUpdate(corsika::setup::Stack::ParticleType& vP,
-                                  corsika::units::si::HEPEnergyType Enew) {
-    using namespace corsika::units::si;
-    HEPMomentumType Pnew = elab2plab(Enew, vP.GetMass());
-    auto pnew = vP.GetMomentum();
-    vP.SetMomentum(pnew * Pnew / pnew.GetNorm());
-  }
-
-  void EnergyLoss::FillProfile(SetupParticle const& vP, SetupTrack const& vTrack,
-                               const corsika::units::si::HEPEnergyType dE) {
-
-    using namespace corsika::units::si;
-
-    auto const toStart = vTrack.GetPosition(0) - InjectionPoint_;
-    auto const toEnd = vTrack.GetPosition(1) - InjectionPoint_;
+  template <typename TOutput>
+  inline YAML::Node BetheBlochPDG<TOutput>::getConfig() const {
 
-    auto const v1 = (toStart * 1_Hz).dot(ShowerAxisDirection_);
-    auto const v2 = (toEnd * 1_Hz).dot(ShowerAxisDirection_);
-    corsika::Line const lineToStartBin(InjectionPoint_, ShowerAxisDirection_ * v1);
-    corsika::Line const lineToEndBin(InjectionPoint_, ShowerAxisDirection_ * v2);
-
-    SetupTrack const trajToStartBin(lineToStartBin, 1_s);
-    SetupTrack const trajToEndBin(lineToEndBin, 1_s);
-
-    GrammageType const grammageStart =
-        vP.GetNode()->GetModelProperties().IntegratedGrammage(trajToStartBin,
-                                                              trajToStartBin.GetLength());
-    GrammageType const grammageEnd =
-        vP.GetNode()->GetModelProperties().IntegratedGrammage(trajToEndBin,
-                                                              trajToEndBin.GetLength());
-
-    const int binStart = grammageStart / dX_;
-    const int binEnd = grammageEnd / dX_;
-
-    std::cout << "energy deposit of " << -dE << " between " << grammageStart << " and "
-              << grammageEnd << std::endl;
-
-    auto energyCount = HEPEnergyType::zero();
-
-    auto fill = [&](int bin, GrammageType weight) {
-      const auto dX = grammageEnd - grammageStart;
-      if (dX > dX_threshold_) {
-        auto const increment = -dE * weight / (grammageEnd - grammageStart);
-        Profile_[bin] += increment;
-        energyCount += increment;
-
-        std::cout << "filling bin " << bin << " with weight " << weight << ": "
-                  << increment << std::endl;
-      }
-    };
-
-    // fill longitudinal profile
-    fill(binStart, (1 + binStart) * dX_ - grammageStart);
-    fill(binEnd, grammageEnd - binEnd * dX_);
-
-    if (binStart == binEnd) { fill(binStart, -dX_); }
-
-    for (int bin = binStart + 1; bin < binEnd; ++bin) { fill(bin, dX_); }
-
-    std::cout << "total energy added to histogram: " << energyCount << std::endl;
-  }
-
-  void EnergyLoss::PrintProfile() const {
-    using namespace corsika::units::si;
-    std::ofstream file("EnergyLossProfile.dat");
-    std::cout << "# EnergyLoss PrintProfile  X-bin [g/cm2]  dE/dX [GeV/g/cm2]  " << std::endl;
-    double const deltaX = dX_ / 1_g * square(1_cm);
-    for (auto v : Profile_) {
-      file << v.first * deltaX << " " << v.second / (deltaX * 1_GeV) << std::endl;
-    }
+    YAML::Node node;
+    node["type"] = "BetheBlochPDG";
+    return node;
   }
 
-} // namespace corsika::energy_loss
+} // namespace corsika

corsika/detail/modules/epos/InteractionModel.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/modules/epos/InteractionModel.hpp>
+#include <corsika/modules/epos/EposStack.hpp>
+#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <corsika/framework/utility/COMBoost.hpp>
+#include <corsika/modules/Random.hpp>
+#include <corsika/framework/utility/CorsikaData.hpp>
+
+#include <epos.hpp>
+
+#include <string>
+#include <tuple>
+#include <cmath>
+
+namespace corsika::epos {
+
+  inline InteractionModel::InteractionModel(std::set<Code> vList,
+                                            std::string const& dataPath,
+                                            bool const epos_printout_on)
+      : data_path_(dataPath)
+      , epos_listing_(epos_printout_on) {
+    // initialize Eposlhc
+    corsika::connect_random_stream(RNG_, ::epos::set_rng_function);
+    if (!isInitialized_) {
+      isInitialized_ = true;
+      if (dataPath == "") {
+        data_path_ = (std::string(corsika_data("EPOS").c_str()) + "/").c_str();
+      }
+      initialize();
+      if (vList.empty()) {
+        CORSIKA_LOGGER_DEBUG(logger_,
+                             "set all particles known to CORSIKA stable inside EPOS..");
+        setParticleListStable(get_all_particles());
+      } else {
+        CORSIKA_LOGGER_DEBUG(logger_, "set specific particles stable inside EPOS..");
+        setParticleListStable(vList);
+      }
+    }
+  }
+
+  inline void InteractionModel::setParticleListStable(std::set<Code> vPartList) const {
+    for (auto& p : vPartList) {
+      int const eid = convertToEposRaw(p);
+      if (eid != 0) {
+        // LCOV_EXCL_START
+        // this is only a safeguard against messing up the epos internals by initializing
+        // more than once.
+        unsigned int const n_particles_stable_epos =
+            ::epos::nodcy_.nrnody; // avoid waring -Wsign-compare
+        if (n_particles_stable_epos < ::epos::mxnody) {
+          CORSIKA_LOGGER_TRACE(logger_, "setting {} with EposId={} stable inside EPOS.",
+                               p, eid);
+          ::epos::nodcy_.nrnody = ::epos::nodcy_.nrnody + 1;
+          ::epos::nodcy_.nody[::epos::nodcy_.nrnody - 1] = eid;
+        } else {
+          CORSIKA_LOGGER_ERROR(logger_, "List of stable particles too long for Epos!");
+          throw std::runtime_error("Epos initialization error!");
+        }
+        // LCOV_EXCL_STOP
+      } else {
+        CORSIKA_LOG_TRACE(
+            "particle conversion Corsika-->Epos not known for {}. Using {}. Setting "
+            "unstable in Epos!",
+            p, eid);
+      }
+    }
+    CORSIKA_LOGGER_DEBUG(logger_, "set {} particles stable inside Epos",
+                         ::epos::nodcy_.nrnody);
+  }
+
+  inline bool InteractionModel::isValid(Code const projectileId, Code const targetId,
+                                        HEPEnergyType const sqrtS) const {
+    //! eposlhc only accepts nuclei with X<=A<=Y as targets, or protons aka Hydrogen or
+    //! neutrons (p,n == nucleon)
+    if (!is_nucleus(targetId) && targetId != Code::Neutron && targetId != Code::Proton) {
+      return false;
+    }
+    if (is_nucleus(targetId) && (get_nucleus_A(targetId) >= get_nucleus_A(maxNucleus_))) {
+      return false;
+    }
+    if ((minEnergyCoM_ > sqrtS) || (sqrtS > maxEnergyCoM_)) { return false; }
+    if (!epos::canInteract(projectileId)) { return false; }
+    return true;
+  }
+
+  inline void InteractionModel::initialize() const {
+
+    CORSIKA_LOGGER_DEBUG(logger_, "initializing...");
+
+    // corsika7 ini
+    int iarg = 0;
+    ::epos::aaset_(iarg);
+
+    // debug output settings
+    ::epos::prnt1_.ish = 0;
+    ::epos::prnt3_.iwseed = 0; // 1: printout seeds, 0: off
+    ::epos::files_.ifch = 6;   // output unit, 6: screen
+
+    // dummy set seeds for random number generator in epos. need to fool epos checks...
+    // we will use external generator
+    ::epos::cseed_.seedi = 1;
+    ::epos::cseed_.seedj = 1;
+    ::epos::cseed_.seedc = 1;
+
+    ::epos::enrgy_.egymin = minEnergyCoM_ / 1_GeV; // 6.;
+    ::epos::enrgy_.egymax = maxEnergyCoM_ / 1_GeV; // 2.e6;
+
+    ::epos::lhcparameters_();
+
+    ::epos::hadr6_.isigma = 0; // do not show cross section
+    ::epos::hadr6_.isetcs = 3; /*  !option to obtain pomeron parameters
+      ! 0.....determine parameters but do not use Kfit
+      ! 1.....determine parameters and use Kfit
+      ! else..get from table
+      !         should be sufficiently detailed
+      !          say iclegy1=1,iclegy2=99
+      !         table is always done, more or less detailed!!!
+      !and option to use cross section tables
+      ! 2....tabulation
+      ! 3....simulation
+                               */
+    ::epos::cjinti_.ionudi =
+        1; // !include quasi elastic events but strict calculation of xs
+    ::epos::cjinti_.iorsce = 0; // !color exchange turned on(1) or off(0)
+    ::epos::cjinti_.iorsdf = 3; //  !droplet formation turned on(>0) or off(0)
+    ::epos::cjinti_.iorshh = 0; //    !other hadron-hadron int. turned on(1) or off(0)
+
+    ::epos::othe1_.istore = 0; // do not produce epos output file
+    ::epos::nucl6_.infragm =
+        2; // 0: keep free nucleons in fragmentation,1: one fragment, 2: fragmentation
+
+    ::epos::othe2_.iframe = 11; // cms frame
+
+    // decay settings
+    // activate decays in epos for particles defined by set_stable/set_unstable
+    // ::epos::othe2_.idecay = 0; // no decays in epos
+
+    // set paths to tables in corsika data
+    ::epos::datadir BASE(data_path_);
+    strcpy(::epos::fname_.fnnx, BASE.data);
+    ::epos::nfname_.nfnnx = BASE.length;
+
+    ::epos::datadir TL(data_path_ + "epos.initl");
+    strcpy(::epos::fname_.fnii, TL.data);
+    ::epos::nfname_.nfnii = TL.length;
+
+    ::epos::datadir EV(data_path_ + "epos.iniev");
+    strcpy(::epos::fname_.fnie, EV.data);
+    ::epos::nfname_.nfnie = EV.length;
+
+    ::epos::datadir RJ(data_path_ + "epos.inirj"); // lhcparameters adds ".lhc"
+    strcpy(::epos::fname_.fnrj, RJ.data);
+    ::epos::nfname_.nfnrj = RJ.length;
+
+    ::epos::datadir CS(data_path_ + "epos.inics"); // lhcparameters adds ".lhc"
+    strcpy(::epos::fname_.fncs, CS.data);
+    ::epos::nfname_.nfncs = CS.length;
+
+    // initializes maximum energy and mass
+    initializeEventCoM(
+        maxNucleus_, get_nucleus_A(maxNucleus_), get_nucleus_Z(maxNucleus_), maxNucleus_,
+        get_nucleus_A(maxNucleus_), get_nucleus_Z(maxNucleus_), maxEnergyCoM_);
+  }
+
+  inline void InteractionModel::initializeEventCoM(Code const idBeam, int const iBeamA,
+                                                   int const iBeamZ, Code const idTarget,
+                                                   int const iTargetA, int const iTargetZ,
+                                                   HEPEnergyType const EcmNN) const {
+    CORSIKA_LOGGER_TRACE(logger_,
+                         "initialize event in CoM frame!"
+                         " Ecm={}",
+                         EcmNN);
+    ::epos::lept1_.engy = -1.;
+    ::epos::enrgy_.ecms = -1.;
+    ::epos::enrgy_.elab = -1.;
+    ::epos::enrgy_.ekin = -1.;
+    ::epos::hadr1_.pnll = -1.;
+
+    ::epos::enrgy_.ecms = EcmNN / 1_GeV; // -> c.m.s. frame
+
+    CORSIKA_LOGGER_TRACE(logger_, "inside EPOS: Ecm={}, Elab={}", ::epos::enrgy_.ecms,
+                         ::epos::enrgy_.elab);
+
+    configureParticles(idBeam, iBeamA, iBeamZ, idTarget, iTargetA, iTargetZ);
+    ::epos::ainit_();
+  }
+
+  inline void InteractionModel::configureParticles(Code const idBeam, int const iBeamA,
+                                                   int const iBeamZ, Code const idTarget,
+                                                   int const iTargetA,
+                                                   int const iTargetZ) const {
+    CORSIKA_LOGGER_TRACE(logger_,
+                         "setting "
+                         "Beam={}, "
+                         "BeamA={}, "
+                         "BeamZ={}, "
+                         "Target={}, "
+                         "TargetA={}, "
+                         "TargetZ={} ",
+                         idBeam, iBeamA, iBeamZ, idTarget, iTargetA, iTargetZ);
+
+    if (is_nucleus(idBeam)) {
+      ::epos::hadr25_.idprojin = convertToEposRaw(Code::Proton);
+      ::epos::nucl1_.laproj = iBeamZ;
+      ::epos::nucl1_.maproj = iBeamA;
+    } else {
+      ::epos::hadr25_.idprojin = convertToEposRaw(idBeam);
+      ::epos::nucl1_.laproj = -1;
+      ::epos::nucl1_.maproj = 1;
+    }
+
+    if (is_nucleus(idTarget)) {
+      ::epos::hadr25_.idtargin = convertToEposRaw(Code::Proton);
+      ::epos::nucl1_.matarg = iTargetA;
+      ::epos::nucl1_.latarg = iTargetZ;
+    } else if (idTarget == Code::Proton || idTarget == Code::Hydrogen) {
+      ::epos::hadr25_.idtargin = convertToEposRaw(Code::Proton);
+      ::epos::nucl1_.matarg = 1;
+      ::epos::nucl1_.latarg = -1;
+    } else if (idTarget == Code::Neutron) {
+      ::epos::hadr25_.idtargin = convertToEposRaw(Code::Neutron);
+      ::epos::nucl1_.matarg = 1;
+      ::epos::nucl1_.latarg = -1;
+    }
+
+    CORSIKA_LOGGER_TRACE(logger_,
+                         "inside EPOS: "
+                         "Id beam={}, "
+                         "Z beam={}, "
+                         "A beam={}, "
+                         "XS beam={}, "
+                         "Id target={}, "
+                         "Z target={}, "
+                         "A target={}, "
+                         "XS target={} ",
+                         ::epos::hadr25_.idprojin, ::epos::nucl1_.laproj,
+                         ::epos::nucl1_.maproj, ::epos::had10_.iclpro,
+                         ::epos::hadr25_.idtargin, ::epos::nucl1_.latarg,
+                         ::epos::nucl1_.matarg, ::epos::had10_.icltar);
+  }
+
+  inline InteractionModel::~InteractionModel() {
+    CORSIKA_LOGGER_DEBUG(logger_, "n={} ", count_);
+  }
+
+  inline std::tuple<CrossSectionType, CrossSectionType>
+  InteractionModel::calcCrossSectionCoM(Code const BeamId, int const BeamA,
+                                        int const BeamZ, Code const TargetId,
+                                        int const TargetA, int const TargetZ,
+                                        const HEPEnergyType EnergyCOM) const {
+    CORSIKA_LOGGER_DEBUG(logger_,
+                         "calcCrossSection: input:"
+                         " beamId={}, beamA={}, beamZ={}"
+                         " target={}, targetA={}, targetZ={}"
+                         " Ecm={:4.3f} GeV,",
+                         BeamId, BeamA, BeamZ, TargetId, TargetA, TargetZ,
+                         EnergyCOM / 1_GeV);
+
+    const int iBeam = epos::getEposXSCode(
+        BeamId); // 0 (can not interact, 1: proton-like, 2: pion-like, 3:kaon-like)
+
+    CORSIKA_LOGGER_TRACE(logger_,
+                         "projectile cross section type={} "
+                         "(0: cannot interact, 1:pion, 2:baryon, 3:kaon)",
+                         iBeam);
+    // reset beam particle // (1: pion-like, 2: proton-like, 3:kaon-like)
+    if (iBeam == 1)
+      initializeEventCoM(Code::PiPlus, BeamA, BeamZ, TargetId, TargetA, TargetZ,
+                         EnergyCOM);
+    else if (iBeam == 2)
+      initializeEventCoM(Code::Proton, BeamA, BeamZ, TargetId, TargetA, TargetZ,
+                         EnergyCOM);
+    else if (iBeam == 3)
+      initializeEventCoM(Code::KPlus, BeamA, BeamZ, TargetId, TargetA, TargetZ,
+                         EnergyCOM);
+
+    double sigProd, sigEla = 0;
+    float sigTot1, sigProd1, sigCut1 = 0;
+    if (!is_nucleus(TargetId) && !is_nucleus(BeamId)) {
+      sigProd = ::epos::hadr5_.sigine;
+      sigEla = ::epos::hadr5_.sigela;
+    } else {
+      // calculate from model, SLOW:
+      float sigQEla1 = 0; // target fragmentation/excitation
+      ::epos::crseaaepos_(sigTot1, sigProd1, sigCut1, sigQEla1);
+      sigProd = sigProd1;
+      // sigEla not properly defined here
+    }
+    CORSIKA_LOGGER_DEBUG(logger_,
+                         "calcCrossSectionCoM: output:"
+                         " sigProd={} mb,"
+                         " sigEla={} mb",
+                         sigProd, sigEla);
+
+    return std::make_tuple(sigProd * 1_mb, sigEla * 1_mb);
+  }
+
+  inline std::tuple<CrossSectionType, CrossSectionType>
+  InteractionModel::readCrossSectionTableLab(Code const BeamId, int const BeamA,
+                                             int const BeamZ, Code const TargetId,
+                                             HEPEnergyType const EnergyLab) const {
+    CORSIKA_LOGGER_DEBUG(logger_,
+                         "readCrossSectionTableLab: input: "
+                         "beamId={}, "
+                         "beamA={}, "
+                         "beamZ={} "
+                         "targetId={}, "
+                         "ELab={:12.2f} GeV,",
+                         BeamId, BeamA, BeamZ, TargetId, EnergyLab / 1_GeV);
+
+    // read cross section from epos internal tables
+    int Abeam = 0;
+    float Ekin = -1;
+
+    if (is_nucleus(BeamId)) {
+      Abeam = BeamA;
+      // kinetic energy per nucleon
+      Ekin = (EnergyLab / Abeam - constants::nucleonMass) / 1_GeV;
+    } else {
+      ::epos::hadr2_.idproj = convertToEposRaw(BeamId);
+      int const iBeam = epos::getEposXSCode(
+          BeamId); // 0 (can not interact, 1: pion-like, 2: proton-like, 3:kaon-like)
+      CORSIKA_LOGGER_TRACE(logger_,
+                           "readCrossSectionTableLab: projectile cross section type={} "
+                           "(0: cannot interact, 1:pion, 2:baryon, 3:kaon)",
+                           iBeam);
+
+      ::epos::had10_.iclpro = iBeam;
+      Abeam = 1;
+      Ekin = (EnergyLab - get_mass(BeamId)) / 1_GeV;
+    }
+
+    int Atarget = 1;
+    if (is_nucleus(TargetId)) { Atarget = get_nucleus_A(TargetId); }
+
+    int iMode = 3; // 0: air, >0 not air
+
+    CORSIKA_LOGGER_TRACE(logger_,
+                         "readCrossSectionTableLab: inside Epos "
+                         "beamId={}, beamXS={}",
+                         ::epos::hadr2_.idproj, ::epos::had10_.iclpro);
+
+    CORSIKA_LOGGER_TRACE(logger_,
+                         "readCrossSectionTableLab: calling Epos cross section with"
+                         "Ekin = {}, Abeam = {}, Atarget = {}, iMode = {}",
+                         Ekin, Abeam, Atarget, iMode);
+
+    // cross section from table, FAST
+    float sigProdEpos = ::epos::eposcrse_(Ekin, Abeam, Atarget, iMode);
+    // sig-el from analytic calculation, no fast
+    float sigElaEpos = ::epos::eposelacrse_(Ekin, Abeam, Atarget, iMode);
+
+    CORSIKA_LOGGER_TRACE(logger_,
+                         "readCrossSectionTableLab: result: sigProd = {}, sigEla = {}",
+                         sigProdEpos, sigElaEpos);
+
+    return std::make_tuple(sigProdEpos * 1_mb, sigElaEpos * 1_mb);
+  }
+
+  inline std::tuple<CrossSectionType, CrossSectionType>
+  InteractionModel::getCrossSectionInelEla(Code const projectileId, Code const targetId,
+                                           FourMomentum const& projectileP4,
+                                           FourMomentum const& targetP4) const {
+    auto const sqrtS2 = (projectileP4 + targetP4).getNormSqr();
+    auto const sqrtS = sqrt(sqrtS2);
+
+    if (!isValid(projectileId, targetId, sqrtS)) {
+      return {CrossSectionType::zero(), CrossSectionType::zero()};
+    }
+    HEPEnergyType const Elab = (sqrtS2 - static_pow<2>(get_mass(projectileId)) -
+                                static_pow<2>(get_mass(targetId))) /
+                               (2 * get_mass(targetId));
+    int beamA = 1;
+    int beamZ = 1;
+    if (is_nucleus(projectileId)) {
+      beamA = get_nucleus_A(projectileId);
+      beamZ = get_nucleus_Z(projectileId);
+    }
+
+    CORSIKA_LOGGER_DEBUG(logger_,
+                         "getCrossSectionLab: input:"
+                         " beamId={}, beamA={}, beamZ={}"
+                         " target={}"
+                         " ELab={:4.3f} GeV, sqrtS={}",
+                         projectileId, beamA, beamZ, targetId, Elab / 1_GeV,
+                         sqrtS / 1_GeV);
+    return readCrossSectionTableLab(projectileId, beamA, beamZ, targetId, Elab);
+  }
+
+  template <typename TSecondaryView>
+  inline void InteractionModel::doInteraction(TSecondaryView& view,
+                                              Code const projectileId,
+                                              Code const targetId,
+                                              FourMomentum const& projectileP4,
+                                              FourMomentum const& targetP4) {
+
+    count_ = count_ + 1;
+
+    // define nucleon-nucleon center-of-mass frame
+    auto const projectileP4NN =
+        projectileP4 / (is_nucleus(projectileId) ? get_nucleus_A(projectileId) : 1);
+    auto const targetP4NN =
+        targetP4 / (is_nucleus(targetId) ? get_nucleus_A(targetId) : 1);
+    auto const SNN = (projectileP4NN + targetP4NN).getNormSqr();
+    HEPEnergyType const sqrtSNN = sqrt(SNN);
+    if (!isValid(projectileId, targetId, sqrtSNN)) {
+      throw std::runtime_error("invalid projectile/target/energy combination.");
+    }
+    HEPEnergyType const Elab = (SNN - static_pow<2>(get_mass(projectileId)) -
+                                static_pow<2>(get_mass(targetId))) /
+                               (2 * get_mass(targetId));
+
+    // system of initial-state
+    COMBoost const boost(projectileP4NN, targetP4NN);
+
+    auto const& originalCS = boost.getOriginalCS();
+    auto const& csPrime =
+        boost.getRotatedCS(); // z is along the CM motion (projectile, in Cascade)
+
+    CORSIKA_LOGGER_DEBUG(logger_,
+                         "doInteraction: interaction, projectile id={}, E={}, p3={} ",
+                         projectileId, projectileP4.getTimeLikeComponent(),
+                         projectileP4.getSpaceLikeComponents());
+    CORSIKA_LOGGER_DEBUG(
+        logger_, "doInteraction: projectile per-nucleon ENN={}, p3NN={} ",
+        projectileP4NN.getTimeLikeComponent(), projectileP4NN.getSpaceLikeComponents());
+    CORSIKA_LOGGER_DEBUG(
+        logger_, "doInteraction: interaction, target id={}, E={}, p3={} ", targetId,
+        targetP4.getTimeLikeComponent(), targetP4.getSpaceLikeComponents());
+    CORSIKA_LOGGER_DEBUG(logger_, "doInteraction: target per-nucleon ENN={}, p3NN={} ",
+                         targetP4NN.getTimeLikeComponent(),
+                         targetP4NN.getSpaceLikeComponents());
+    CORSIKA_LOGGER_DEBUG(logger_, "doInteraction: Elab={}, sqrtSNN={} ", Elab, sqrtSNN);
+
+    int beamA = 1;
+    int beamZ = 1;
+    if (is_nucleus(projectileId)) {
+      beamA = get_nucleus_A(projectileId);
+      beamZ = get_nucleus_Z(projectileId);
+      CORSIKA_LOGGER_DEBUG(logger_, "projectile: A={}, Z={} ", beamA, beamZ);
+    }
+
+    // // from corsika7 interface
+    // // NEXLNK-part
+    int targetA = 1;
+    int targetZ = 1;
+    if (is_nucleus(targetId)) {
+      targetA = get_nucleus_A(targetId);
+      targetZ = get_nucleus_Z(targetId);
+      CORSIKA_LOGGER_DEBUG(logger_, "target: A={}, Z={} ", targetA, targetZ);
+    }
+    initializeEventCoM(projectileId, beamA, beamZ, targetId, targetA, targetZ, sqrtSNN);
+
+    // create event
+    int iarg = 1;
+    ::epos::aepos_(iarg);
+    ::epos::afinal_();
+
+    if (epos_listing_) { // LCOV_EXCL_START
+      char nam[9] = "EPOSLHC&";
+      ::epos::alistf_(nam, 9);
+    } // LCOV_EXCL_STOP
+
+    // NSTORE-part
+
+    MomentumVector P_final(originalCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
+    HEPEnergyType E_final = 0_GeV;
+
+    // secondaries
+    EposStack es;
+    CORSIKA_LOGGER_DEBUG(logger_, "number of entries on Epos stack: {}", es.getSize());
+    for (auto& psec : es) {
+      if (!psec.isFinal()) continue;
+
+      auto momentum = psec.getMomentum(csPrime);
+      // transform particle output to frame defined by input 4-momenta
+      auto const P4output = boost.fromCoM(FourVector{psec.getEnergy(), momentum});
+      auto p3output = P4output.getSpaceLikeComponents();
+      p3output.rebase(originalCS); // transform back into standard lab frame
+
+      EposCode const eposId = psec.getPID();
+      Code const pid = epos::convertFromEpos(eposId);
+      HEPEnergyType const mass = get_mass(pid);
+      HEPEnergyType const Ekin = sqrt(p3output.getSquaredNorm() + mass * mass) - mass;
+      CORSIKA_LOGGER_TRACE(logger_,
+                           " id= {}"
+                           " p= {}",
+                           pid, p3output.getComponents() / 1_GeV);
+
+      auto pnew = view.addSecondary(std::make_tuple(pid, Ekin, p3output.normalized()));
+      P_final += pnew.getMomentum();
+      E_final += pnew.getEnergy();
+    }
+    CORSIKA_LOGGER_DEBUG(
+        logger_,
+        "conservation (all GeV): Ecm_final= n/a" /* << Ecm_final / 1_GeV*/
+        ", E_final={} GeV"
+        ", P_final={} GeV"
+        ", no. of particles={}",
+        E_final / 1_GeV, (P_final / 1_GeV).getComponents(), view.getSize());
+  }
+} // namespace corsika::epos

corsika/detail/modules/epos/ParticleConversion.inl

+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <epos.hpp>
+
+namespace corsika::epos {
+
+  inline HEPMassType getEposMass(Code const pCode) {
+    if (is_nucleus(pCode)) throw std::runtime_error("Not suited for Nuclei.");
+    auto sCode = convertToEposRaw(pCode);
+    if (sCode == 0)
+      throw std::runtime_error("getEposMass: unknown particle!");
+    else {
+      float mass2 = 0;
+      ::epos::idmass_(sCode, mass2);
+      return sqrt(mass2) * 1_GeV;
+    }
+  }
+
+  inline PDGCode getEposPDGId(Code const p) {
+    if (!is_nucleus(p)) {
+      int eid = corsika::epos::convertToEposRaw(p);
+      char nxs[4] = "nxs";
+      char pdg[4] = "pdg";
+      return static_cast<PDGCode>(::epos::idtrafo_(nxs, pdg, eid));
+    } else {
+      throw std::runtime_error("Epos id conversion not implemented for nuclei!");
+    }
+  }
+} // namespace corsika::epos

corsika/detail/modules/fluka/InteractionModel.inl

+/*
+ * (c) Copyright 2022 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <algorithm>
+#include <cstdlib>
+#include <vector>
+#include <iterator>
+#include <set>
+#include <utility>
+#include <string_view>
+
+#include <boost/iterator/zip_iterator.hpp>
+#include <Eigen/Dense>
+
+#include <corsika/media/Environment.hpp>
+#include <corsika/media/NuclearComposition.hpp>
+#include <corsika/framework/geometry/FourVector.hpp>
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <corsika/framework/core/EnergyMomentumOperations.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/modules/Random.hpp>
+#include <corsika/modules/fluka/ParticleConversion.hpp>
+
+#include <FLUKA.hpp>
+
+namespace corsika::fluka {
+  inline InteractionModel::InteractionModel(std::set<Code> const& nuccomp)
+      : materials_{genFlukaMaterials(nuccomp)}
+      , cumsgx_{std::make_unique<double[]>(materials_.size() * 3)} {
+    CORSIKA_LOGGER_INFO(logger_, "FLUKA version {}", ::fluka::get_version());
+    for (auto const& [code, matno] : materials_) {
+      CORSIKA_LOGGER_DEBUG(logger_, "FLUKA material initialization: {} -> {}",
+                           get_name(code, full_name{}), matno);
+    }
+
+    if (int const ndmhep = ::fluka::ndmhep_(); ::fluka::nmxhep != ndmhep) {
+      CORSIKA_LOGGER_CRITICAL(logger_, "HEPEVT dimension mismatch. FLUKA reports %d",
+                              ndmhep);
+      throw std::runtime_error{"FLUKA HEPEVT dimension mismatch"};
+    }
+
+    corsika::connect_random_stream(RNG_, ::fluka::set_rng_function);
+  }
+
+  inline bool InteractionModel::isValid(Code projectileID, int material,
+                                        HEPEnergyType /*sqrtS*/) const {
+    if (!fluka::canInteract(projectileID)) {
+      // invalid projectile
+      return false;
+    }
+
+    if (material < 0) {
+      // invalid/uninitialized target
+      return false;
+    }
+
+    // TODO: check validity range of sqrtS
+    return true;
+  }
+
+  inline bool InteractionModel::isValid(Code projectileID, Code targetID,
+                                        HEPEnergyType sqrtS) const {
+    return isValid(projectileID, getMaterialIndex(targetID), sqrtS);
+  }
+
+  inline int InteractionModel::getMaterialIndex(Code targetID) const {
+    auto compare = [=](std::pair<Code, int> const& p) { return p.first == targetID; };
+    if (auto it = std::find_if(materials_.cbegin(), materials_.cend(), compare);
+        it == materials_.cend()) {
+      return -1;
+    } else {
+      return it->second;
+    }
+  }
+
+  inline CrossSectionType InteractionModel::getCrossSection(
+      Code const projectileId, Code const targetId, FourMomentum const& projectileP4,
+      FourMomentum const& targetP4) const {
+    auto const flukaMaterial = getMaterialIndex(targetId);
+
+    HEPEnergyType const sqrtS = (projectileP4 + targetP4).getNorm();
+    if (!isValid(projectileId, flukaMaterial, sqrtS)) { return CrossSectionType::zero(); }
+
+    COMBoost const targetRestBoost{targetP4.getSpaceLikeComponents(), get_mass(targetId)};
+    FourMomentum const projectileLab4mom = targetRestBoost.toCoM(projectileP4);
+    HEPEnergyType const Elab = projectileLab4mom.getTimeLikeComponent();
+    auto constexpr invGeV = 1 / 1_GeV;
+    double const labMomentum =
+        projectileLab4mom.getSpaceLikeComponents().getNorm() * invGeV;
+
+    auto const plab = projectileLab4mom.getSpaceLikeComponents();
+
+    CORSIKA_LOGGER_DEBUG(logger_, fmt::format("Elab = {} GeV", Elab * invGeV));
+    auto const flukaCodeProj =
+        static_cast<FLUKACodeIntType>(convertToFluka(projectileId));
+
+    double const dummyEkin = 0;
+    CrossSectionType const xs = ::fluka::sgmxyz_(&flukaCodeProj, &flukaMaterial,
+                                                 &dummyEkin, &labMomentum, &iflxyz_) *
+                                1_mb;
+    return xs;
+  }
+
+  template <typename TSecondaryView>
+  inline void InteractionModel::doInteraction(TSecondaryView& view,
+                                              Code const projectileId,
+                                              Code const targetId,
+                                              FourMomentum const& projectileP4,
+                                              FourMomentum const& targetP4) {
+
+    auto const flukaCodeProj =
+        static_cast<FLUKACodeIntType>(convertToFluka(projectileId));
+    auto const flukaMaterial = getMaterialIndex(targetId);
+
+    HEPEnergyType const sqrtS = (projectileP4 + targetP4).getNorm();
+    if (!isValid(projectileId, flukaMaterial, sqrtS)) {
+      std::string const errmsg = fmt::format(
+          "Event generation with invalid configuration requested: proj: {}, target: {}",
+          get_name(projectileId, full_name{}), get_name(targetId, full_name{}));
+      CORSIKA_LOGGER_CRITICAL(logger_, errmsg);
+      throw std::runtime_error{errmsg.c_str()};
+    }
+
+    COMBoost const targetRestBoost{targetP4.getSpaceLikeComponents(), get_mass(targetId)};
+    FourMomentum const projectileLab4mom = targetRestBoost.toCoM(projectileP4);
+    auto constexpr invGeV = 1 / 1_GeV;
+
+    auto const plab = projectileLab4mom.getSpaceLikeComponents();
+    auto const& cs = plab.getCoordinateSystem();
+    auto const labMomentum = plab.getNorm();
+    double const labMomentumGeV = labMomentum * invGeV;
+
+    auto const direction = (plab / labMomentum).getComponents().getEigenVector();
+
+    double const dummyEkin = 0;
+    ::fluka::evtxyz_(&flukaCodeProj, &flukaMaterial, &dummyEkin, &labMomentumGeV,
+                     &direction[0], &direction[1], &direction[2], &iflxyz_, cumsgx_.get(),
+                     cumsgx_.get() + materials_.size(),
+                     cumsgx_.get() + materials_.size() * 2);
+
+    for (int i = 0; i < ::fluka::hepevt_.nhep; ++i) {
+      int const status = ::fluka::hepevt_.isthep[i];
+      if (status != 1) // skip non-final-state particles
+        continue;
+
+      auto const pdg = static_cast<corsika::PDGCode>(::fluka::hepevt_.idhep[i]);
+      auto const c8code = corsika::convert_from_PDG(pdg);
+      auto const mom = QuantityVector<hepenergy_d>{
+          Eigen::Map<Eigen::Vector3d>(&::fluka::hepevt_.phep[i][0]) *
+          (1_GeV).magnitude()};
+      auto const pPrime = mom.getNorm();
+      auto const c8mass = corsika::get_mass(c8code);
+
+      auto const fourMomCollisionFrame =
+          FourVector{calculate_total_energy(pPrime, c8mass), MomentumVector{cs, mom}};
+      auto const fourMomOrigFrame = targetRestBoost.fromCoM(fourMomCollisionFrame);
+      auto const momOrigFrame = fourMomOrigFrame.getSpaceLikeComponents();
+      auto const p = momOrigFrame.getNorm();
+
+      view.addSecondary(std::tuple{c8code, corsika::calculate_kinetic_energy(p, c8mass),
+                                   momOrigFrame / p});
+    }
+  }
+
+  inline std::vector<std::pair<Code, int>> InteractionModel::genFlukaMaterials(
+      std::set<Code> const& allElementsInUniverse) {
+
+    /*
+     * We define one material per element/isotope we have in C8. Cross-section averaging
+     * and target isotope selection happen in C8.
+     */
+
+    int const nElements = allElementsInUniverse.size();
+    auto nelmfl = std::make_unique<int[]>(nElements);
+    std::vector<int> izelfl;
+    izelfl.reserve(nElements);
+    auto wfelml = std::make_unique<double[]>(nElements);
+    auto const mxelfl = nElements;
+    double const pptmax = 1e11; // GeV
+    double const ef2dp3 = 0;    // GeV, 0 means default is used
+    double const df2dp3 = -1;   // default
+    bool const lprint = true;
+    auto mtflka = std::make_unique<int[]>(mxelfl);
+    // magic number that FLUKA uses to see if it's the right version
+    char const crvrck[] = "76466879";
+    int const size = 8;
+
+    std::fill(&nelmfl[0], &nelmfl[nElements], 1);
+    std::fill(&wfelml[0], &wfelml[nElements], 1.);
+    std::transform(allElementsInUniverse.cbegin(), allElementsInUniverse.cend(),
+                   std::back_inserter(izelfl), get_nucleus_Z);
+
+    // check if FLUPRO is set
+    if (std::getenv("FLUPRO") == nullptr) {
+      throw std::runtime_error{"FLUPRO not set; required to initialize FLUKA"};
+    }
+
+    // call FLUKA
+    ::fluka::stpxyz_(&nElements, nelmfl.get(), izelfl.data(), wfelml.get(), &mxelfl,
+                     &pptmax, &ef2dp3, &df2dp3, &iflxyz_, &lprint, mtflka.get(), crvrck,
+                     &size);
+
+    // now create & fill vector of (C8 Code, FLUKA mat. no.) pairs
+    std::vector<std::pair<Code, int>> mapping;
+    mapping.reserve(nElements);
+
+    auto it = boost::make_zip_iterator(
+        boost::make_tuple(allElementsInUniverse.begin(), &mtflka[0]));
+    auto const end = boost::make_zip_iterator(
+        boost::make_tuple(allElementsInUniverse.end(), &mtflka[nElements]));
+    for (; it != end; ++it) {
+      boost::tuple<Code const&, int&> tup = *it;
+      mapping.emplace_back(tup.get<0>(), tup.get<1>());
+    }
+
+    return mapping;
+  }
+} // namespace corsika::fluka

corsika/detail/modules/null_model/NullModel.inl

-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * See file AUTHORS for a list of contributors.
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#pragma once
-
-#include <corsika/modules/null_model/NullModel.hpp>
-#include <corsika/setup/SetupStack.hpp>
-#include <corsika/setup/SetupTrajectory.hpp>
-
-namespace corsika::null_model {
-
-  void NullModel::Init() {}
-
-  NullModel::NullModel(corsika::units::si::LengthType maxStepLength)
-      : fMaxStepLength(maxStepLength) {}
-
-  template <>
-  corsika::EProcessReturn NullModel::DoContinuous(corsika::setup::Stack::ParticleType&,
-                                                  corsika::setup::Trajectory&) const {
-    return corsika::EProcessReturn::eOk;
-  }
-
-  template <>
-  units::si::LengthType NullModel::MaxStepLength(corsika::setup::Stack::ParticleType&,
-                                                 corsika::setup::Trajectory&) const {
-    return fMaxStepLength;
-  }
-
-} // namespace corsika::null_model