corsika/detail/media/NuclearComposition.inl

 /*
  * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
  *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
  */
 
 #pragma once
@@ -11,7 +10,8 @@
 #include <corsika/framework/core/ParticleProperties.hpp>
 #include <corsika/framework/core/PhysicalUnits.hpp>
 
-#include <corsika/media/WeightProvider.hpp>
+#include <boost/iterator/zip_iterator.hpp>
+#include <boost/iterator/transform_iterator.hpp>
 
 #include <cassert>
 #include <functional>
@@ -22,32 +22,53 @@
 
 namespace corsika {
 
-  NuclearComposition::NuclearComposition(std::vector<Code> const& pComponents,
-                                         std::vector<float> const& pFractions)
+  inline NuclearComposition::NuclearComposition(std::vector<Code> const& pComponents,
+                                                std::vector<double> const& pFractions)
       : numberFractions_(pFractions)
       , components_(pComponents)
-      , avgMassNumber_(std::inner_product(
-            pComponents.cbegin(), pComponents.cend(), pFractions.cbegin(), 0.,
-            std::plus<double>(), [](auto const compID, auto const fraction) -> double {
-              if (is_nucleus(compID)) {
-                return get_nucleus_A(compID) * fraction;
-              } else {
-                return get_mass(compID) / convert_SI_to_HEP(constants::u) * fraction;
-              }
-            })) {
-    assert(pComponents.size() == pFractions.size());
-    auto const sumFractions =
-        std::accumulate(pFractions.cbegin(), pFractions.cend(), 0.f);
-
-    if (!(0.999f < sumFractions && sumFractions < 1.001f)) {
+      , avgMassNumber_(getWeightedSum([](Code const compID) -> double {
+        if (is_nucleus(compID)) {
+          return get_nucleus_A(compID);
+        } else {
+          return get_mass(compID) / convert_SI_to_HEP(constants::u);
+        }
+      })) {
+    if (pComponents.size() != pFractions.size()) {
+      throw std::runtime_error(
+          "Cannot construct NuclearComposition from vectors of different sizes.");
+    }
+    auto const sumFractions = std::accumulate(pFractions.cbegin(), pFractions.cend(), 0.);
+
+    if (!(0.999 < sumFractions && sumFractions < 1.001)) {
       throw std::runtime_error("element fractions do not add up to 1");
     }
     this->updateHash();
   }
 
   template <typename TFunction>
-  inline auto NuclearComposition::getWeightedSum(TFunction const& func) const {
-    using ResultQuantity = decltype(func(*components_.cbegin()));
+  inline auto NuclearComposition::getWeighted(TFunction func) const {
+    using ResultQuantity = decltype(func(std::declval<Code>()));
+    auto const product = [&](auto const compID, auto const fraction) {
+      return func(compID) * fraction;
+    };
+
+    if constexpr (phys::units::is_quantity_v<ResultQuantity>) {
+      std::vector<ResultQuantity> result(components_.size(), ResultQuantity::zero());
+      std::transform(components_.cbegin(), components_.cend(), numberFractions_.cbegin(),
+                     result.begin(), product);
+      return result;
+    } else {
+      std::vector<ResultQuantity> result(components_.size(), ResultQuantity(0));
+      std::transform(components_.cbegin(), components_.cend(), numberFractions_.cbegin(),
+                     result.begin(), product);
+      return result;
+    }
+  } // namespace corsika
+
+  template <typename TFunction>
+  inline auto NuclearComposition::getWeightedSum(TFunction func) const
+      -> decltype(func(std::declval<Code>())) {
+    using ResultQuantity = decltype(func(std::declval<Code>()));
 
     auto const prod = [&](auto const compID, auto const fraction) {
       return func(compID) * fraction;
@@ -68,7 +89,7 @@ namespace corsika {
 
   inline size_t NuclearComposition::getSize() const { return numberFractions_.size(); }
 
-  inline std::vector<float> const& NuclearComposition::getFractions() const {
+  inline std::vector<double> const& NuclearComposition::getFractions() const {
     return numberFractions_;
   }
 
@@ -82,16 +103,28 @@ namespace corsika {
 
   template <class TRNG>
   inline Code NuclearComposition::sampleTarget(std::vector<CrossSectionType> const& sigma,
-                                               TRNG& randomStream) const {
+                                               TRNG&& randomStream) const {
+    if (sigma.size() != numberFractions_.size()) {
+      throw std::runtime_error("incompatible vector sigma as input");
+    }
+
+    auto zip_beg = boost::make_zip_iterator(
+        boost::make_tuple(numberFractions_.cbegin(), sigma.cbegin()));
+    auto zip_end = boost::make_zip_iterator(
+        boost::make_tuple(numberFractions_.cend(), sigma.cend()));
+    using zip_iter_type = decltype(zip_beg);
+
+    auto const mult_func = [](zip_iter_type::value_type const& zipit) -> double {
+      return zipit.get<0>() * zipit.get<1>().magnitude();
+    };
+
+    using transform_iter_type =
+        boost::transform_iterator<decltype(mult_func), zip_iter_type, double, double>;
 
-    assert(sigma.size() == numberFractions_.size());
+    auto trans_beg = transform_iter_type{zip_beg, mult_func};
+    auto trans_end = transform_iter_type{zip_end, mult_func};
 
-    std::discrete_distribution channelDist(
-        WeightProviderIterator<decltype(numberFractions_.begin()),
-                               decltype(sigma.begin())>(numberFractions_.begin(),
-                                                        sigma.begin()),
-        WeightProviderIterator<decltype(numberFractions_.begin()), decltype(sigma.end())>(
-            numberFractions_.end(), sigma.end()));
+    std::discrete_distribution channelDist{trans_beg, trans_end};
 
     auto const iChannel = channelDist(randomStream);
     return components_[iChannel];
@@ -99,11 +132,17 @@ namespace corsika {
 
   // Note: when this class ever modifies its internal data, the hash
   // must be updated, too!
+  // the hash value is important to find tables, etc.
   inline size_t NuclearComposition::getHash() const { return hash_; }
 
+  inline bool NuclearComposition::operator==(NuclearComposition const& v) const {
+    return v.hash_ == hash_;
+  }
+
   inline void NuclearComposition::updateHash() {
     std::vector<std::size_t> hashes;
-    for (float ifrac : this->getFractions()) hashes.push_back(std::hash<float>{}(ifrac));
+    for (double ifrac : this->getFractions())
+      hashes.push_back(std::hash<double>{}(ifrac));
     for (Code icode : this->getComponents())
       hashes.push_back(std::hash<int>{}(static_cast<int>(icode)));
     std::size_t h = std::hash<double>{}(this->getAverageMassNumber());

corsika/detail/media/ShowerAxis.inl

 /*
  * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
  *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
  */
 
 #include <corsika/framework/core/PhysicalUnits.hpp>
@@ -14,9 +13,9 @@
 namespace corsika {
 
   template <typename TEnvModel>
-  ShowerAxis::ShowerAxis(Point const& pStart, Vector<length_d> const& length,
-                         Environment<TEnvModel> const& env, bool const doThrow,
-                         int const steps)
+  inline ShowerAxis::ShowerAxis(Point const& pStart, Vector<length_d> const& length,
+                                Environment<TEnvModel> const& env, bool const doThrow,
+                                int const steps)
       : pointStart_(pStart)
       , length_(length)
       , throw_(doThrow)
@@ -26,9 +25,16 @@ namespace corsika {
       , X_(steps + 1) {
     auto const* const universe = env.getUniverse().get();
 
-    auto rho = [pStart, length, universe](double x) {
+    auto rho = [pStart, length, universe, doThrow](double x) {
       auto const p = pStart + length * x;
       auto const* node = universe->getContainingNode(p);
+      if (!node->hasModelProperties()) {
+        CORSIKA_LOG_CRITICAL(
+            "Unable to construct ShowerAxis. ShowerAxis includes volume "
+            "with no model properties at point {}.",
+            p);
+        if (doThrow) throw std::runtime_error("Unable to construct ShowerAxis.");
+      }
       return node->getModelProperties().getMassDensity(p).magnitude();
     };
 
@@ -58,14 +64,14 @@ namespace corsika {
     assert(std::is_sorted(d_.cbegin(), d_.cend()));
   }
 
-  GrammageType ShowerAxis::getX(LengthType l) const {
+  inline GrammageType ShowerAxis::getX(LengthType l) const {
     double const fractionalBin = l / steplength_;
     int const lower = fractionalBin; // indices of nearest X support points
     double const fraction = fractionalBin - lower;
     unsigned int const upper = lower + 1;
 
     if (fractionalBin < 0) {
-      CORSIKA_LOG_ERROR("cannot extrapolate to points behind point of injection l={} m",
+      CORSIKA_LOG_TRACE("cannot extrapolate to points behind point of injection l={} m",
                         l / 1_m);
       if (throw_) {
         throw std::runtime_error(
@@ -75,7 +81,7 @@ namespace corsika {
     }
 
     if (upper >= X_.size()) {
-      CORSIKA_LOG_ERROR(
+      CORSIKA_LOG_TRACE(
           "shower axis too short, cannot extrapolate (l / max_length_ = {} )",
           l / max_length_);
       if (throw_) {
@@ -84,6 +90,7 @@ namespace corsika {
             l / max_length_);
         throw std::runtime_error(err.c_str());
       }
+      return getMaximumX();
     }
     CORSIKA_LOG_TRACE("showerAxis::X frac={}, fractionalBin={}, lower={}, upper={}",
                       fraction, fractionalBin, lower, upper);
@@ -97,21 +104,21 @@ namespace corsika {
     return X_[upper] * fraction + X_[lower] * (1 - fraction);
   }
 
-  LengthType ShowerAxis::getSteplength() const { return steplength_; }
+  inline LengthType ShowerAxis::getSteplength() const { return steplength_; }
 
-  GrammageType ShowerAxis::getMaximumX() const { return *X_.rbegin(); }
+  inline GrammageType ShowerAxis::getMaximumX() const { return *X_.rbegin(); }
 
-  GrammageType ShowerAxis::getMinimumX() const { return GrammageType::zero(); }
+  inline GrammageType ShowerAxis::getMinimumX() const { return GrammageType::zero(); }
 
-  GrammageType ShowerAxis::getProjectedX(Point const& p) const {
+  inline GrammageType ShowerAxis::getProjectedX(Point const& p) const {
     auto const projectedLength = (p - pointStart_).dot(axis_normalized_);
     return getX(projectedLength);
   }
 
-  Vector<dimensionless_d> const& ShowerAxis::getDirection() const {
+  inline DirectionVector const& ShowerAxis::getDirection() const {
     return axis_normalized_;
   }
 
-  Point const& ShowerAxis::getStart() const { return pointStart_; }
+  inline Point const& ShowerAxis::getStart() const { return pointStart_; }
 
 } // namespace corsika

corsika/detail/media/SlidingPlanarExponential.inl

 /*
  * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
  *
- * This software is distributed under the terms of the 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 <typename T>
-  SlidingPlanarExponential<T>::SlidingPlanarExponential(
-      Point const& p0, MassDensityType rho0, LengthType lambda,
-      NuclearComposition const& nuclComp, LengthType referenceHeight)
-      : BaseExponential<SlidingPlanarExponential<T>>(p0, rho0, lambda)
-      , nuclComp_(nuclComp)
-      , referenceHeight_(referenceHeight) {}
-
-  template <typename T>
-  MassDensityType SlidingPlanarExponential<T>::getMassDensity(Point const& point) const {
-    auto const height =
-        (point - BaseExponential<SlidingPlanarExponential<T>>::getAnchorPoint())
-            .getNorm() -
-        referenceHeight_;
-    return BaseExponential<SlidingPlanarExponential<T>>::getRho0() *
-           exp(BaseExponential<SlidingPlanarExponential<T>>::getInvLambda() * 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 T>
-  NuclearComposition const& SlidingPlanarExponential<T>::getNuclearComposition() const {
+  template <typename TDerived>
+  inline NuclearComposition const&
+  SlidingPlanarExponential<TDerived>::getNuclearComposition() const {
     return nuclComp_;
   }
 
-  template <typename T>
-  GrammageType SlidingPlanarExponential<T>::getIntegratedGrammage(
-      setup::Trajectory const& traj, LengthType l) const {
-    auto const axis = (traj.getPosition(0) -
-                       BaseExponential<SlidingPlanarExponential<T>>::getAnchorPoint())
-                          .normalized();
-    return BaseExponential<SlidingPlanarExponential<T>>::getIntegratedGrammage(traj, 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 <typename T>
-  LengthType SlidingPlanarExponential<T>::getArclengthFromGrammage(
-      setup::Trajectory const& traj, GrammageType const grammage) const {
-    auto const axis = (traj.getPosition(0) -
-                       BaseExponential<SlidingPlanarExponential<T>>::getAnchorPoint())
-                          .normalized();
-    return BaseExponential<SlidingPlanarExponential<T>>::getArclengthFromGrammage(
+  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);
   }
 

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

corsika/detail/media/Universe.inl

 /*
  * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
  *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
  */
 
 #pragma once
@@ -12,10 +11,10 @@
 
 namespace corsika {
 
-  Universe::Universe(CoordinateSystemPtr const& pCS)
+  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

corsika/detail/media/VolumeTreeNode.inl

 /*
  * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
  *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
  */
 
 #pragma once
@@ -14,7 +13,7 @@
 namespace corsika {
 
   template <typename IModelProperties>
-  bool VolumeTreeNode<IModelProperties>::contains(Point const& p) const {
+  inline bool VolumeTreeNode<IModelProperties>::contains(Point const& p) const {
     return geoVolume_->contains(p);
   }
 
@@ -32,7 +31,7 @@ namespace corsika {
    * \class Point \p p, or nullptr iff \p p is not contained in this volume.
    */
   template <typename IModelProperties>
-  VolumeTreeNode<IModelProperties> const*
+  inline VolumeTreeNode<IModelProperties> const*
   VolumeTreeNode<IModelProperties>::getContainingNode(Point const& p) const {
     if (!contains(p)) { return nullptr; }
 
@@ -52,19 +51,37 @@ namespace corsika {
     }
   }
 
+  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(childNodes_.begin(), childNodes_.end(),
-                  [&](auto& v) { v->walk(func); });
+                  [&](auto const& v) { v->walk(func); });
 
     if constexpr (!preorder) { func(*this); };
   }
 
   template <typename IModelProperties>
-  void VolumeTreeNode<IModelProperties>::addChild(
+  inline void VolumeTreeNode<IModelProperties>::addChild(
       typename VolumeTreeNode<IModelProperties>::VTNUPtr pChild) {
     pChild->parentNode_ = this;
     childNodes_.push_back(std::move(pChild));
@@ -74,19 +91,9 @@ namespace corsika {
   }
 
   template <typename IModelProperties>
-  void VolumeTreeNode<IModelProperties>::excludeOverlapWith(
+  inline void VolumeTreeNode<IModelProperties>::excludeOverlapWith(
       typename VolumeTreeNode<IModelProperties>::VTNUPtr const& pNode) {
     excludedNodes_.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)...);
-    }
-  */
 } // namespace corsika

corsika/detail/media/WeightProvider.inl

-/*
- * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#pragma once
-
-#include <vector>
-
-namespace corsika {
-
-  template <class AConstIterator, class BConstIterator>
-  WeightProviderIterator<AConstIterator, BConstIterator>::WeightProviderIterator(
-      AConstIterator a, BConstIterator b)
-      : aIter_(a)
-      , bIter_(b) {}
-
-  template <class AConstIterator, class BConstIterator>
-  typename WeightProviderIterator<AConstIterator, BConstIterator>::value_type
-      WeightProviderIterator<AConstIterator, BConstIterator>::operator*() const {
-    return ((*aIter_) * (*bIter_)).magnitude();
-  }
-
-  template <class AConstIterator, class BConstIterator>
-  WeightProviderIterator<AConstIterator, BConstIterator>&
-  WeightProviderIterator<AConstIterator,
-                         BConstIterator>::operator++() { // prefix ++
-    ++aIter_;
-    ++bIter_;
-    return *this;
-  }
-
-  template <class AConstIterator, class BConstIterator>
-  bool WeightProviderIterator<AConstIterator, BConstIterator>::operator==(
-      WeightProviderIterator other) {
-    return aIter_ == other.aIter_;
-  }
-
-  template <class AConstIterator, class BConstIterator>
-  bool WeightProviderIterator<AConstIterator, BConstIterator>::operator!=(
-      WeightProviderIterator other) {
-    return !(*this == other);
-  }
-
-} // namespace corsika
\ No newline at end of file

corsika/detail/modules/HadronicElasticModel.inl

 /*
  * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
  *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
  */
 
 #pragma once
@@ -16,20 +15,19 @@
 #include <corsika/framework/random/ExponentialDistribution.hpp>
 #include <corsika/framework/utility/COMBoost.hpp>
 
-#include <corsika/setup/SetupStack.hpp>
-
 #include <iomanip>
 #include <iostream>
 
 namespace corsika {
 
-  HadronicElasticInteraction::HadronicElasticInteraction(CrossSectionType x,
-                                                         CrossSectionType y)
+  inline HadronicElasticInteraction::HadronicElasticInteraction(CrossSectionType x,
+                                                                CrossSectionType y)
       : parX_(x)
       , parY_(y) {}
 
-  template <>
-  GrammageType HadronicElasticInteraction::getInteractionLength(SetupParticle const& p) {
+  template <typename TParticle>
+  inline GrammageType HadronicElasticInteraction::getInteractionLength(
+      TParticle const& p) {
     if (p.getPID() == Code::Proton) {
       auto const* currentNode = p.getNode();
       auto const& mediumComposition =
@@ -52,7 +50,7 @@ namespace corsika {
           avgCrossSection += getCrossSection(s) * fractions[i];
         }
 
-        std::cout << "avgCrossSection: " << avgCrossSection / 1_mb << " mb" << std::endl;
+        CORSIKA_LOG_DEBUG("avgCrossSection: {} mb", avgCrossSection / 1_mb);
 
         return avgCrossSection;
       }();
@@ -69,7 +67,7 @@ namespace corsika {
   }
 
   template <typename TParticle>
-  ProcessReturn HadronicElasticInteraction::doInteraction(TParticle& p) {
+  inline ProcessReturn HadronicElasticInteraction::doInteraction(TParticle& p) {
     if (p.getPID() != Code::Proton) { return ProcessReturn::Ok; }
 
     const auto* currentNode = p.getNode();
@@ -116,7 +114,7 @@ namespace corsika {
     auto const s = static_pow<2>(sqrtS);
 
     auto const B = this->B(s);
-    std::cout << B << std::endl;
+    CORSIKA_LOG_DEBUG(B);
 
     ExponentialDistribution tDist(1 / B);
     auto const absT = [&]() {
@@ -133,10 +131,12 @@ namespace corsika {
       return absT;
     }();
 
-    std::cout << "HadronicElasticInteraction: s = " << s * constants::invGeVsq
-              << " GeV²; absT = " << absT * constants::invGeVsq << " GeV² (max./GeV² = "
-              << 4 * constants::invGeVsq * projectileMomentumSquaredNorm << ')'
-              << std::endl;
+    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_);
@@ -158,17 +158,17 @@ namespace corsika {
     return ProcessReturn::Ok;
   }
 
-  HadronicElasticInteraction::inveV2 HadronicElasticInteraction::B(eV2 s) const {
+  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;
-    std::cout << "B(" << s << ") = " << result / constants::invGeVsq << " GeV¯²"
-              << std::endl;
+    CORSIKA_LOG_DEBUG("B({}) = {}  GeV¯²", s, result / constants::invGeVsq);
+
     return result;
   }
 
-  CrossSectionType HadronicElasticInteraction::getCrossSection(
+  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) +
@@ -180,8 +180,8 @@ namespace corsika {
         static_pow<2>(sigmaTotal) /
         (16 * constants::pi * convert_HEP_to_SI<CrossSectionType::dimension_type>(B(s)));
 
-    std::cout << "HEM sigmaTot = " << sigmaTotal / 1_mb << " mb" << std::endl;
-    std::cout << "HEM sigmaElastic = " << sigmaElastic / 1_mb << " mb" << std::endl;
+    CORSIKA_LOG_DEBUG("HEM sigmaTot = {} mb", sigmaTotal / 1_mb);
+    CORSIKA_LOG_DEBUG("HEM sigmaElastic = {} mb", sigmaElastic / 1_mb);
     return sigmaElastic;
   }
 

corsika/detail/modules/LongitudinalProfile.inl

 /*
  * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
  *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
  */
 
 #include <corsika/framework/core/ParticleProperties.hpp>
 #include <corsika/framework/core/Logging.hpp>
-
+#include <corsika/framework/core/Step.hpp>
 #include <corsika/modules/LongitudinalProfile.hpp>
 
-#include <corsika/setup/SetupStack.hpp>
-#include <corsika/setup/SetupTrajectory.hpp>
-
 #include <cmath>
 #include <iomanip>
 #include <limits>
+#include <utility>
 
 namespace corsika {
 
-  LongitudinalProfile::LongitudinalProfile(ShowerAxis const& shower_axis, GrammageType dX)
-      : dX_(dX)
-      , shower_axis_{shower_axis}
-      , profiles_{static_cast<unsigned int>(shower_axis.getMaximumX() / dX_) + 1} {}
-
-  template <typename TParticle, typename TTrack>
-  ProcessReturn LongitudinalProfile::doContinuous(TParticle const& vP,
-                                                  TTrack const& vTrack) {
-    auto const pid = vP.getPID();
-
-    GrammageType const grammageStart = shower_axis_.getProjectedX(vTrack.getPosition(0));
-    GrammageType const grammageEnd = shower_axis_.getProjectedX(vTrack.getPosition(1));
+  template <typename TOutput>
+  template <typename... TArgs>
+  inline LongitudinalProfile<TOutput>::LongitudinalProfile(TArgs&&... args)
+      : TOutput(std::forward<TArgs>(args)...) {}
 
-    CORSIKA_LOG_INFO(
-        "pos1={} m, pos2={}, X={} g/cm2", vTrack.getPosition(0).getCoordinates() / 1_m,
-        vTrack.getPosition(1).getCoordinates() / 1_m, grammageStart / 1_g * square(1_cm));
-
-    const int binStart = std::ceil(grammageStart / dX_);
-    const int binEnd = std::floor(grammageEnd / dX_);
-
-    for (int b = binStart; b <= binEnd; ++b) {
-      if (pid == Code::Gamma) {
-        profiles_.at(b)[ProfileIndex::Gamma]++;
-      } else if (pid == Code::Positron) {
-        profiles_.at(b)[ProfileIndex::Positron]++;
-      } else if (pid == Code::Electron) {
-        profiles_.at(b)[ProfileIndex::Electron]++;
-      } else if (pid == Code::MuPlus) {
-        profiles_.at(b)[ProfileIndex::MuPlus]++;
-      } else if (pid == Code::MuMinus) {
-        profiles_.at(b)[ProfileIndex::MuMinus]++;
-      } else if (is_hadron(pid)) {
-        profiles_.at(b)[ProfileIndex::Hadron]++;
-      }
-    }
+  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;
   }
 
-  void LongitudinalProfile::save(std::string const& filename, const int width,
-                                 const int precision) {
-    std::ofstream f{filename};
-    f << "# X / g·cm¯², gamma, e+, e-, mu+, mu-, all hadrons" << std::endl;
-    for (size_t b = 0; b < profiles_.size(); ++b) {
-      f << std::setprecision(5) << std::setw(11) << b * (dX_ / (1_g / 1_cm / 1_cm));
-      for (auto const& N : profiles_.at(b)) {
-        f << std::setw(width) << std::setprecision(precision) << std::scientific << N;
-      }
-      f << std::endl;
-    }
+  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 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.
  */
 
-#include <corsika/modules/ObservationPlane.hpp>
-
-#include <fstream>
-
 namespace corsika {
 
-  ObservationPlane::ObservationPlane(Plane const& obsPlane, DirectionVector const& x_axis,
-                                     std::string const& filename, bool deleteOnHit)
-      : plane_(obsPlane)
-      , outputStream_(filename)
-      , deleteOnHit_(deleteOnHit)
-      , energy_ground_(0_GeV)
-      , count_ground_(0)
+  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_)) {
-    outputStream_ << "#PDG code, energy / eV, x distance / m, y distance / m"
-                  << std::endl;
-  }
-
-  ProcessReturn ObservationPlane::doContinuous(
-      corsika::setup::Stack::particle_type& particle,
-      corsika::setup::Trajectory& trajectory) {
-    TimeType const timeOfIntersection =
-        (plane_.getCenter() - trajectory.getPosition(0)).dot(plane_.getNormal()) /
-        trajectory.getVelocity(0).dot(plane_.getNormal());
-
-    if (timeOfIntersection < TimeType::zero()) { return ProcessReturn::Ok; }
-
-    if (plane_.isAbove(trajectory.getPosition(0)) ==
-        plane_.isAbove(trajectory.getPosition(1))) {
-      return ProcessReturn::Ok;
+      , 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;
     }
 
-    const auto energy = particle.getEnergy();
-    auto const displacement = trajectory.getPosition(1) - plane_.getCenter();
+    HEPEnergyType const kineticEnergy = step.getEkinPost();
+    Point const pointOfIntersection = step.getPositionPost();
+    Vector const displacement = pointOfIntersection - plane_.getCenter();
+    DirectionVector const direction = step.getDirectionPost();
 
-    outputStream_ << static_cast<int>(get_PDG(particle.getPID())) << ' ' << energy / 1_eV
-                  << ' ' << displacement.dot(xAxis_) / 1_m << ' '
-                  << displacement.dot(yAxis_) / 1_m
-                  << (trajectory.getPosition(1) - plane_.getCenter()).getNorm() / 1_m
-                  << std::endl;
+    // 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_) {
-      count_ground_++;
-      energy_ground_ += energy;
-      particle.erase();
       return ProcessReturn::ParticleAbsorbed;
     } else {
       return ProcessReturn::Ok;
     }
-  }
+  } // namespace corsika
 
-  LengthType ObservationPlane::getMaxStepLength(
-      corsika::setup::Stack::particle_type const& vParticle,
-      corsika::setup::Trajectory const& trajectory) {
+  template <typename TTracking, typename TOutput>
+  template <typename TParticle, typename TTrajectory>
+  inline LengthType ObservationPlane<TTracking, TOutput>::getMaxStepLength(
+      TParticle const& particle, TTrajectory const& trajectory) {
 
-    int chargeNumber;
-    if (is_nucleus(vParticle.getPID())) {
-      chargeNumber = vParticle.getNuclearZ();
-    } else {
-      chargeNumber = get_charge_number(vParticle.getPID());
-    }
-    auto const* currentLogicalVolumeNode = vParticle.getNode();
-    auto magneticfield = currentLogicalVolumeNode->getModelProperties().getMagneticField(
-        vParticle.getPosition());
-    auto direction = trajectory.getVelocity(0).normalized();
-
-    if (chargeNumber != 0 &&
-        abs(plane_.getNormal().dot(
-            trajectory.getLine().getVelocity().cross(magneticfield))) *
-                1_s / 1_m / 1_T >
-            1e-6) {
-      auto const* currentLogicalVolumeNode = vParticle.getNode();
-      auto magneticfield =
-          currentLogicalVolumeNode->getModelProperties().getMagneticField(
-              vParticle.getPosition());
-      auto k =
-          chargeNumber * constants::c * 1_eV / (vParticle.getMomentum().getNorm() * 1_V);
-
-      if (direction.dot(plane_.getNormal()) * direction.dot(plane_.getNormal()) -
-              (plane_.getNormal().dot(trajectory.getPosition(0) - plane_.getCenter()) *
-               plane_.getNormal().dot(direction.cross(magneticfield)) * 2 * k) <
-          0) {
-        return std::numeric_limits<double>::infinity() * 1_m;
-      }
-
-      LengthType MaxStepLength1 =
-          (sqrt(direction.dot(plane_.getNormal()) * direction.dot(plane_.getNormal()) -
-                (plane_.getNormal().dot(trajectory.getPosition(0) - plane_.getCenter()) *
-                 plane_.getNormal().dot(direction.cross(magneticfield)) * 2 * k)) -
-           direction.dot(plane_.getNormal()) / direction.getNorm()) /
-          (plane_.getNormal().dot(direction.cross(magneticfield)) * k);
-
-      LengthType MaxStepLength2 =
-          (-sqrt(direction.dot(plane_.getNormal()) * direction.dot(plane_.getNormal()) -
-                 (plane_.getNormal().dot(trajectory.getPosition(0) - plane_.getCenter()) *
-                  plane_.getNormal().dot(direction.cross(magneticfield)) * 2 * k)) -
-           direction.dot(plane_.getNormal()) / direction.getNorm()) /
-          (plane_.getNormal().dot(direction.cross(magneticfield)) * k);
-
-      if (MaxStepLength1 <= 0_m && MaxStepLength2 <= 0_m) {
-        return std::numeric_limits<double>::infinity() * 1_m;
-      } else if (MaxStepLength1 <= 0_m || MaxStepLength2 < MaxStepLength1) {
-        std::cout << " steplength to obs plane 2: " << MaxStepLength2 << std::endl;
-        return MaxStepLength2 *
-               (direction + direction.cross(magneticfield) * MaxStepLength2 * k / 2)
-                   .getNorm() *
-               1.001;
-      } else if (MaxStepLength2 <= 0_m || MaxStepLength1 < MaxStepLength2) {
-        std::cout << " steplength to obs plane 1: " << MaxStepLength1 << std::endl;
-        return MaxStepLength1 *
-               (direction + direction.cross(magneticfield) * MaxStepLength2 * k / 2)
-                   .getNorm() *
-               1.001;
-      }
-    }
+    CORSIKA_LOG_TRACE("getMaxStepLength, particle={}, pos={}, dir={}, plane={}",
+                      particle.asString(), particle.getPosition(),
+                      particle.getDirection(), plane_.asString());
 
-    TimeType const timeOfIntersection =
-        (plane_.getCenter() - trajectory.getPosition(0)).dot(plane_.getNormal()) /
-        trajectory.getVelocity(0).dot(plane_.getNormal());
+    auto const intersection = TTracking::intersect(particle, plane_);
 
-    if (timeOfIntersection < TimeType::zero()) {
+    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();
-
-    auto const pointOfIntersection = trajectory.getPosition(fractionOfIntersection);
-    auto dist = (trajectory.getPosition(0) - pointOfIntersection).getNorm() * 1.0001;
-    CORSIKA_LOG_TRACE("ObservationPlane w/o magnetic field: getMaxStepLength l={} m",
-                      dist / 1_m);
-    return dist;
-  }
-
-  void ObservationPlane::showResults() const {
-    CORSIKA_LOG_INFO(
-        " ******************************\n"
-        " ObservationPlane: \n"
-        " energy in ground (GeV)     :  {}\n"
-        " no. of particles in ground :  {}\n"
-        " ******************************",
-        energy_ground_ / 1_GeV, count_ground_);
+    CORSIKA_LOG_TRACE("ObservationPlane: getMaxStepLength dist={} m, pos={}",
+                      trajectory.getLength(fractionOfIntersection) / 1_m,
+                      trajectory.getPosition(fractionOfIntersection));
+    return trajectory.getLength(fractionOfIntersection);
   }
 
-  void ObservationPlane::reset() {
-    energy_ground_ = 0_GeV;
-    count_ground_ = 0;
+  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 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/FourVector.hpp>
 #include <corsika/modules/OnShellCheck.hpp>
+#include <corsika/framework/core/Logging.hpp>
 
 namespace corsika {
 
-  OnShellCheck::OnShellCheck(const double vMassTolerance, const double vEnergyTolerance,
-                             const bool vError)
+  inline OnShellCheck::OnShellCheck(double const vMassTolerance,
+                                    double const vEnergyTolerance, bool const vError)
       : mass_tolerance_(vMassTolerance)
       , energy_tolerance_(vEnergyTolerance)
       , throw_error_(vError) {
-    std::cout << "OnShellCheck: mass tolerance is set to " << mass_tolerance_ * 100 << "%"
-              << std::endl
-              << "              energy tolerance is set to " << energy_tolerance_ * 100
-              << "%" << std::endl;
+    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);
   }
 
-  OnShellCheck::~OnShellCheck() {
-    std::cout << "OnShellCheck: summary" << std::endl
-              << " particles shifted: " << int(count_) << std::endl;
-    if (count_)
-      std::cout << " average energy shift (%): " << average_shift_ / count_ * 100.
-                << std::endl
-                << " max. energy shift (%): " << max_shift_ * 100. << std::endl;
+  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>
-  void OnShellCheck::doSecondaries(TView& vS) {
+  inline void OnShellCheck::doSecondaries(TView& vS) {
     for (auto& p : vS) {
       auto const pid = p.getPID();
-      if (!is_hadron(pid) || is_nucleus(pid)) continue;
+      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);
@@ -51,33 +51,36 @@ namespace corsika {
         count_ = count_ + 1;
         average_shift_ += abs(e_shift_relative);
         if (abs(e_shift_relative) > max_shift_) max_shift_ = abs(e_shift_relative);
-        std::cout << "OnShellCheck: shift particle mass for " << pid << std::endl
-                  << std::setw(40) << std::setfill(' ')
-                  << "corsika mass (GeV): " << m_corsika / 1_GeV << std::endl
-                  << std::setw(40) << std::setfill(' ')
-                  << "kinetic mass (GeV): " << m_kinetic / 1_GeV << std::endl
-                  << std::setw(40) << std::setfill(' ')
-                  << "m_kin-m_cor (GeV): " << m_err_abs / 1_GeV << std::endl
-                  << std::setw(40) << std::setfill(' ')
-                  << "mass tolerance (GeV): " << (m_corsika * mass_tolerance_) / 1_GeV
-                  << std::endl;
+        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_) {
-          std::cout << "OnShellCheck: warning! shifted particle energy by "
-                    << e_shift_relative * 100 << " %" << std::endl;
-          if (throw_error_)
+          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
-        std::cout << "OnShellCheck: particle mass for " << pid << " OK" << std::endl;
+      } 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 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/ParticleCut.hpp>
+#include <corsika/framework/core/Logging.hpp>
 
 namespace corsika {
 
-  ParticleCut::ParticleCut(const HEPEnergyType eCut, bool em, bool inv)
-      : energy_cut_(eCut)
-      , doCutEm_(em)
-      , doCutInv_(inv)
-      , energy_(0_GeV)
-      , em_energy_(0_GeV)
-      , em_count_(0)
-      , inv_energy_(0_GeV)
-      , inv_count_(0) {}
-
-  template <typename TParticle>
-  bool ParticleCut::isBelowEnergyCut(TParticle const& vP) const {
-    auto const energyLab = vP.getEnergy();
-    // nuclei
-    if (vP.getPID() == Code::Nucleus) {
-      // calculate energy per nucleon
-      auto const ElabNuc = energyLab / vP.getNuclearA();
-      return (ElabNuc < energy_cut_);
-    } else {
-      return (energyLab < energy_cut_);
+  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);
   }
 
-  bool ParticleCut::isEmParticle(Code vCode) const {
-    // FOR NOW: switch
-    switch (vCode) {
-      case Code::Gamma:
-      case Code::Electron:
-      case Code::Positron:
-        return true;
-      default:
-        return false;
+  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);
   }
 
-  bool ParticleCut::isInvisible(Code vCode) const {
-    switch (vCode) {
-      case Code::NuE:
-      case Code::NuEBar:
-      case Code::NuMu:
-      case Code::NuMuBar:
-        return true;
+  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");
+  }
 
-      default:
-        return false;
+  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 TParticle>
-  bool ParticleCut::checkCutParticle(const TParticle& particle) {
-
-    const Code pid = particle.getPID();
-    HEPEnergyType energy = particle.getEnergy();
-    CORSIKA_LOG_DEBUG(fmt::format("ParticleCut: checking {}, E= {} GeV, EcutTot={} GeV",
-                                  pid, energy / 1_GeV,
-                                  (em_energy_ + inv_energy_ + energy_) / 1_GeV));
-    if (doCutEm_ && isEmParticle(pid)) {
-      CORSIKA_LOG_DEBUG("removing em. particle...");
-      em_energy_ += energy;
-      em_count_ += 1;
-      return true;
-    } else if (doCutInv_ && isInvisible(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...");
-      inv_energy_ += energy;
-      inv_count_ += 1;
       return true;
-    } else if (isBelowEnergyCut(particle)) {
+    } else if (isBelowEnergyCut(pid, kine_energy)) {
       CORSIKA_LOG_DEBUG("removing low en. particle...");
-      energy_ += energy;
       return true;
-    } else if (particle.getTime() > 10_ms) {
+    } else if (timePost > 10_ms) {
       CORSIKA_LOG_DEBUG("removing OLD particle...");
-      energy_ += energy;
       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
   }
 
-  void ParticleCut::doSecondaries(corsika::setup::StackView& vS) {
+  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()) {
-      if (checkCutParticle(particle)) { particle.erase(); }
+      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);
   }
 
-  ProcessReturn ParticleCut::doContinuous(corsika::setup::Stack::particle_type& particle,
-                                          corsika::setup::Trajectory const&) {
-    CORSIKA_LOG_TRACE("ParticleCut::DoContinuous");
-    if (checkCutParticle(particle)) {
+  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");
-      particle.erase();
       // signal to upstream code that this particle was deleted
       return ProcessReturn::ParticleAbsorbed;
     }
     return ProcessReturn::Ok;
   }
 
-  void ParticleCut::showResults() {
-    CORSIKA_LOG_INFO(
-        " ******************************\n"
-        " energy in em.  component (GeV): {} \n "
-        " no. of em.  particles injected: {} \n "
-        " energy in inv. component (GeV): {} \n "
-        " no. of inv. particles injected: {} \n "
-        " energy below particle cut (GeV): {} \n"
-        " ******************************",
-        em_energy_ / 1_GeV, em_count_, inv_energy_ / 1_GeV, inv_count_, energy_ / 1_GeV);
+  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);
+    }
   }
 
-  void ParticleCut::reset() {
-    em_energy_ = 0_GeV;
-    em_count_ = 0;
-    inv_energy_ = 0_GeV;
-    inv_count_ = 0;
-    energy_ = 0_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 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,8 +13,6 @@
 #include <corsika/framework/core/PhysicalUnits.hpp>
 #include <corsika/framework/geometry/RootCoordinateSystem.hpp>
 
-#include <corsika/setup/SetupTrajectory.hpp>
-
 #include <chrono>
 #include <iomanip>
 #include <iostream>
@@ -24,56 +21,97 @@
 namespace corsika {
 
   template <typename TStack>
-  StackInspector<TStack>::StackInspector(const int vNStep, const bool vReportStack,
-                                         const HEPEnergyType vE0)
+  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()) {}
+      , StartTime_(std::chrono::system_clock::now())
+      , energyPostInit_(HEPEnergyType::zero())
+      , timePostInit_(std::chrono::system_clock::now()) {}
 
   template <typename TStack>
-  StackInspector<TStack>::~StackInspector() {}
+  inline StackInspector<TStack>::~StackInspector() {}
 
   template <typename TStack>
-  void StackInspector<TStack>::doStack(const TStack& vS) {
+  inline void StackInspector<TStack>::doStack(TStack const& vS) {
 
     [[maybe_unused]] int i = 0;
     HEPEnergyType Etot = 0_GeV;
 
-    for (const auto& iterP : vS) {
+    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);
-        std::cout << "StackInspector: i=" << std::setw(5) << std::fixed << (i++)
-                  << ", id=" << std::setw(30) << iterP.getPID() << " E=" << std::setw(15)
-                  << std::scientific << (E / 1_GeV) << " GeV, "
-                  << " pos=" << pos << " node = " << iterP.getNode();
-        if (iterP.getPID() == Code::Nucleus)
-          std::cout << " nuc_ref=" << iterP.getNucleusRef();
-        std::cout << std::endl;
+        CORSIKA_LOG_INFO(
+            "StackInspector: i= {:5d}"
+            ", id= {:^15}"
+            " E={:15e} GeV, "
+            " pos= {}"
+            " node = {}",
+            (i++), iterP.getPID(), (E / 1_GeV), pos, fmt::ptr(iterP.getNode()));
       }
     }
 
-    auto const now = std::chrono::system_clock::now();
-    const std::chrono::duration<double> elapsed_seconds = now - StartTime_;
-    std::time_t const now_time = std::chrono::system_clock::to_time_t(now);
-    auto const dE = E0_ - Etot;
-    if (dE < dE_threshold_) return;
-    double const progress = dE / E0_;
-
-    double const eta_seconds = elapsed_seconds.count() / progress;
-    std::time_t const eta_time = std::chrono::system_clock::to_time_t(
-        StartTime_ + std::chrono::seconds((int)eta_seconds));
-
-    std::cout << "StackInspector: "
-              << " time=" << std::put_time(std::localtime(&now_time), "%T")
-              << ", running=" << elapsed_seconds.count() << " seconds"
-              << " (" << std::setw(3) << int(progress * 100) << "%)"
-              << ", nStep=" << getStep() << ", stackSize=" << vS.getSize()
-              << ", Estack=" << Etot / 1_GeV << " GeV"
-              << ", ETA=" << std::put_time(std::localtime(&eta_time), "%T") << std::endl;
+    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
\ No newline at end of file
+} // 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 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/TrackWriter.hpp>
-
 #include <corsika/framework/core/ParticleProperties.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
 
-#include <corsika/setup/SetupStack.hpp>
-#include <corsika/setup/SetupTrajectory.hpp>
-
-#include <iomanip>
 #include <limits>
 
 namespace corsika {
 
-  TrackWriter::TrackWriter(std::string const& filename)
-      : filename_(filename) {
-    using namespace std::string_literals;
+  template <typename TOutput>
+  inline TrackWriter<TOutput>::TrackWriter() {}
 
-    file_.open(filename_);
-    file_
-        << "# PID, E / eV, start coordinates / m, displacement vector to end / m, steplength / m "s
-        << '\n';
-  }
+  template <typename TOutput>
+  template <typename TParticle>
+  inline ProcessReturn TrackWriter<TOutput>::doContinuous(Step<TParticle> const& step,
+                                                          bool const) {
 
-  template <typename TParticle, typename TTrack>
-  ProcessReturn TrackWriter::doContinuous(const TParticle& vP, const TTrack& vT) {
-    auto const start = vT.getPosition(0).getCoordinates();
-    auto const delta = vT.getPosition(1).getCoordinates() - start;
-    auto const pdg = static_cast<int>(get_PDG(vP.getPID()));
-
-    // clang-format off
-    file_ << std::setw(7) << pdg
-          << std::setw(width_) << std::scientific << std::setprecision(precision_) << vP.getEnergy() / 1_eV
-          << std::setw(width_) << std::scientific << std::setprecision(precision_) << start[0] / 1_m 
-          << std::setw(width_) << std::scientific << std::setprecision(precision_) << start[1] / 1_m
-          << std::setw(width_) << std::scientific << std::setprecision(precision_) << start[2] / 1_m
-          << std::setw(width_) << std::scientific << std::setprecision(precision_) << delta[0] / 1_m
-          << std::setw(width_) << std::scientific << std::setprecision(precision_) << delta[1] / 1_m
-          << std::setw(width_) << std::scientific << std::setprecision(precision_) << delta[2] / 1_m 
-          << std::setw(width_) << std::scientific << std::setprecision(precision_) << delta.getNorm() / 1_m
-          << '\n';
-    // clang-format on
+    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>
-  LengthType TrackWriter::getMaxStepLength(const TParticle&, const 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/TrackingLine.inl

-/*
- * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#pragma once
-
-#include <corsika/framework/geometry/Point.hpp>
-#include <corsika/framework/geometry/QuantityVector.hpp>
-#include <corsika/framework/geometry/Sphere.hpp>
-#include <corsika/framework/geometry/Vector.hpp>
-#include <corsika/media/Environment.hpp>
-
-#include <limits>
-#include <stdexcept>
-#include <utility>
-
-namespace corsika::tracking_line {
-
-  std::optional<std::pair<TimeType, TimeType>> TimeOfIntersection(
-      corsika::Line const& line, corsika::Sphere const& sphere) {
-    auto const delta = line.getStartPoint() - sphere.getCenter();
-    auto const v = line.getVelocity();
-    auto const vSqNorm =
-        v.getSquaredNorm(); // todo: get rid of this by having V0 normalized always
-    auto const R = sphere.getRadius();
-
-    auto const vDotDelta = v.dot(delta);
-    auto const discriminant =
-        vDotDelta * vDotDelta - vSqNorm * (delta.getSquaredNorm() - R * R);
-
-    if (discriminant.magnitude() > 0) {
-      auto const sqDisc = sqrt(discriminant);
-      auto const invDenom = 1 / vSqNorm;
-      return std::make_pair((-vDotDelta - sqDisc) * invDenom,
-                            (-vDotDelta + sqDisc) * invDenom);
-    } else {
-      return {};
-    }
-  }
-
-  TimeType getTimeOfIntersection(Line const& vLine, Plane const& vPlane) {
-
-    auto const delta = vPlane.getCenter() - vLine.getStartPoint();
-    auto const v = vLine.getVelocity();
-    auto const n = vPlane.getNormal();
-    auto const c = n.dot(v);
-
-    if (c.magnitude() == 0) {
-      return std::numeric_limits<TimeType::value_type>::infinity() * 1_s;
-    } else {
-      return n.dot(delta) / c;
-    }
-  }
-
-} // namespace corsika::tracking_line

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 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.
  */
 
 #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 <iomanip>
+#include <numeric>
 #include <utility>
 
 namespace corsika {
 
-  CONEXhybrid::CONEXhybrid(Point center, ShowerAxis const& showerAxis,
-                           LengthType groundDist, LengthType injectionHeight,
-                           HEPEnergyType primaryEnergy, PDGCode primaryPDG)
-      : center_{center}
+  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();
@@ -35,19 +45,11 @@ namespace corsika {
             make_translation(c8cs, translation.getComponents(c8cs));
         auto const transformCS = make_rotationToZ(intermediateCS, translation);
 
-        std::cout << "translation C8/CONEX obs: " << translation.getComponents()
-                  << std::endl;
+        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!
-        std::cout << transform.matrix() << std::endl << std::endl;
-        std::cout << CoordinateSystem::getTransformation(intermediateCS, c8cs).matrix()
-                  << std::endl
-                  << std::endl;
-        std::cout << CoordinateSystem::getTransformation(intermediateCS2, intermediateCS)
-                         .matrix()
-                  << std::endl;
         */
         return transformCS;
       })}
@@ -63,29 +65,51 @@ namespace corsika {
       })}
       , y_sf_{showerAxis_.getDirection().cross(x_sf_)} {
 
-    std::cout << "x_sf (conexObservationCS): " << x_sf_.getComponents(conexObservationCS_)
-              << std::endl;
-    std::cout << "x_sf (C8): " << x_sf_.getComponents(center.getCoordinateSystem())
-              << std::endl;
-
-    std::cout << "y_sf (conexObservationCS): " << y_sf_.getComponents(conexObservationCS_)
-              << std::endl;
-    std::cout << "y_sf (C8): " << y_sf_.getComponents(center.getCoordinateSystem())
-              << std::endl;
-
-    std::cout << "showerAxisDirection (conexObservationCS): "
-              << showerAxis_.getDirection().getComponents(conexObservationCS_)
-              << std::endl;
-    std::cout << "showerAxisDirection (C8): "
-              << showerAxis_.getDirection().getComponents(center.getCoordinateSystem())
-              << std::endl;
-
-    std::cout << "showerCore (conexObservationCS): "
-              << showerCore_.getCoordinates(conexObservationCS_) << std::endl;
-    std::cout << "showerCore (C8): "
-              << showerCore_.getCoordinates(center.getCoordinateSystem()) << std::endl;
-
-    int randomSeeds[3] = {1234, 0, 0}; // will be overwritten later??
+    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.
@@ -100,8 +124,10 @@ namespace corsika {
                         particleListMode,
 #endif
                         configPath.c_str(), configPath.size());
+  }
 
-    double eprima = primaryEnergy / 1_GeV;
+  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}};
@@ -114,23 +140,29 @@ namespace corsika {
                             showerAxisConex.getX().magnitude()) *
                  180 / M_PI;
 
-    std::cout << "theta (deg) = " << theta << "; phi (deg) = " << phi << std::endl;
+    CORSIKA_LOG_DEBUG(
+        "theta (deg) = {}"
+        "; phi (deg) = {}",
+        theta, phi);
 
-    int ipart = static_cast<int>(primaryPDG);
-    auto rng = RNGManager::getInstance().getRandomStream("cascade");
+    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
 
-    std::array<int, 3> ioseed{static_cast<int>(rng()), static_cast<int>(rng()),
-                              static_cast<int>(rng())};
+    // 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 xminp = injectionHeight / 1_m;
+    double eprima = primaryEnergy_ / 1_GeV;
+    double xminp = injectionHeight_ / 1_m;
 
     ::conex::conexrun_(ipart, eprima, theta, phi, xminp, dimpact, ioseed.data());
   }
 
-  void CONEXhybrid::doSecondaries(setup::StackView& vS) {
+  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();
@@ -142,9 +174,10 @@ namespace corsika {
     }
   }
 
-  bool CONEXhybrid::addParticle(Code pid, HEPEnergyType energy, HEPEnergyType mass,
-                                Point const& position, DirectionVector const& direction,
-                                TimeType t) {
+  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; });
@@ -152,8 +185,8 @@ namespace corsika {
 
     // EM particle
     auto const egs_pid = it->second;
-    std::cout << "position conexObs: " << position.getCoordinates(conexObservationCS_)
-              << std::endl;
+    CORSIKA_LOG_DEBUG("position conexObs: {}",
+                      position.getCoordinates(conexObservationCS_));
 
     auto const coords = position.getCoordinates(conexObservationCS_) / 1_m;
     double const x = coords[0].magnitude();
@@ -179,32 +212,29 @@ namespace corsika {
     double const v = direction.dot(x_sf_).magnitude();
     double const w = direction.dot(showerAxis_.getDirection()).magnitude();
 
-    double const weight = 1; // NEEDS TO BE CHANGED WHEN WE HAVE WEIGHTS!
-
     // 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;
 
-    std::cout << "CONEXhybrid: removing " << egs_pid << " " << std::scientific << energy
-              << " GeV" << std::endl;
-
-    std::cout << "#### parameters to cegs4_() ####" << std::endl;
-    std::cout << "egs_pid = " << egs_pid << std::endl;
-    std::cout << "E = " << E << std::endl;
-    std::cout << "m = " << m << std::endl;
-    std::cout << "x = " << x << std::endl;
-    std::cout << "y = " << y << std::endl;
-    std::cout << "altitude = " << altitude << std::endl;
-    std::cout << "slantDistance = " << slantDistance << std::endl;
-    std::cout << "lateralX = " << lateralX << std::endl;
-    std::cout << "lateralY = " << lateralY << std::endl;
-    std::cout << "slantX = " << slantX << std::endl;
-    std::cout << "time = " << time << std::endl;
-    std::cout << "u = " << u << std::endl;
-    std::cout << "v = " << v << std::endl;
-    std::cout << "w = " << w << std::endl;
+    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;
@@ -229,7 +259,9 @@ namespace corsika {
     return true;
   }
 
-  void CONEXhybrid::solveCE() {
+  template <typename TOutputE, typename TOutputN>
+  template <typename TStack>
+  inline void CONEXhybrid<TOutputE, TOutputN>::doCascadeEquations(TStack&) {
 
     ::conex::conexcascade_();
 
@@ -251,7 +283,7 @@ namespace corsika {
     auto dEdX = std::make_unique<float[]>(maxX);
     auto Mu = std::make_unique<float[]>(maxX);
     auto dMu = std::make_unique<float[]>(maxX);
-    auto Gamma = 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);
 
@@ -260,35 +292,29 @@ namespace corsika {
     ::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, Gamma[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]);
 
-    std::ofstream file{"conex_output.txt"};
-    file << fmt::format("#{:>8} {:>13} {:>13} {:>13} {:>13} {:>13} {:>13} {:>13}\n", "X",
-                        "N", "dEdX", "Mu", "dMu", "Gamma", "El", "Had");
+    // 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) {
-      file << fmt::format(
-          " {:>8.2f} {:>13.3} {:>13.3} {:>13.3} {:>13.3} {:>13.3} {:>13.3} {:>13.3}\n",
-          X[i], N[i], dEdX[i], Mu[i], dMu[i], Gamma[i], Electrons[i], Hadrons[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]);
     }
-    file.close();
-
-    std::ofstream fitout{"conex_fit.txt"};
-    fitout << fitpars[1 - 1] << " # log10(eprima/eV)" << std::endl;
-    fitout << fitpars[2 - 1] << " # theta" << std::endl;
-    fitout << fitpars[3 - 1] << " # X1 (first interaction)" << std::endl;
-    fitout << fitpars[4 - 1] << " # Nmax" << std::endl;
-    fitout << fitpars[5 - 1] << " # X0" << std::endl;
-    fitout << fitpars[6 - 1] << " # P1" << std::endl;
-    fitout << fitpars[7 - 1] << " # P2" << std::endl;
-    fitout << fitpars[8 - 1] << " # P3" << std::endl;
-    fitout << fitpars[9 - 1] << " # chi^2 / sqrt(Nmax)" << std::endl;
-    fitout << fitpars[10 - 1] << " # Xmax" << std::endl;
-    fitout << fitpars[11 - 1] << " # phi" << std::endl;
-    fitout << fitpars[12 - 1] << " # inelasticity 1st int." << std::endl;
-    fitout << fitpars[13 - 1] << " # ???" << std::endl;
-    fitout.close();
+  }
+
+  template <typename TOutputE, typename TOutputN>
+  inline YAML::Node CONEXhybrid<TOutputE, TOutputN>::getConfig() const {
+
+    return YAML::Node();
   }
 
 } // namespace corsika

corsika/detail/modules/energy_loss/BetheBlochPDG.inl

 /*
  * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
  *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
  */
 
 #pragma once
 
 #include <corsika/framework/core/ParticleProperties.hpp>
 
-#include <corsika/setup/SetupStack.hpp>
-#include <corsika/setup/SetupTrajectory.hpp>
-
 #include <corsika/framework/geometry/Line.hpp>
 #include <corsika/framework/core/Logging.hpp>
 
@@ -22,19 +18,15 @@
 
 namespace corsika {
 
-  auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
-    return sqrt((Elab - m) * (Elab + m));
-  };
-
-  BetheBlochPDG::BetheBlochPDG(ShowerAxis const& shower_axis, HEPEnergyType emCut)
-      : dX_(10_g / square(1_cm)) // profile binning
-      , dX_threshold_(0.0001_g / square(1_cm))
-      , shower_axis_(shower_axis)
-      , emCut_(emCut)
-      , profile_(int(shower_axis.getMaximumX() / dX_) + 1) {}
+  template <typename TOutput>
+  template <typename... TArgs>
+  inline BetheBlochPDG<TOutput>::BetheBlochPDG(TArgs&&... args)
+      : TOutput(std::forward<TArgs>(args)...) {}
 
-  HEPEnergyType BetheBlochPDG::getBetheBloch(setup::Stack::particle_type const& p,
-                                             GrammageType const dX) {
+  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
@@ -65,18 +57,18 @@ namespace corsika {
     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
-    CORSIKA_LOG_DEBUG("BetheBloch beta2={}, gamma2={}", beta2, gamma2);
+    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
-    CORSIKA_LOG_DEBUG("BetheBloch Wmax={}", Wmax);
+    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
-    CORSIKA_LOG_DEBUG("BetheBloch 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;
@@ -87,7 +79,7 @@ namespace corsika {
     } else if (x < x0) { // and IF conductor (otherwise, this is 0)
       delta = delta0 * pow(100, 2 * (x - x0));
     }
-    CORSIKA_LOG_DEBUG("BetheBloch delta={}", delta);
+    CORSIKA_LOG_TRACE("BetheBloch delta={}", delta);
 
     // with further low energies correction, accurary ~1% down to beta~0.05 (1MeV for p)
 
@@ -127,135 +119,91 @@ namespace corsika {
   }
 
   // radiation losses according to PDG 2018, ch. 33 ref. [5]
-  HEPEnergyType BetheBlochPDG::getRadiationLosses(setup::Stack::particle_type const& vP,
-                                                  GrammageType const vDX) {
+  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;
   }
 
-  HEPEnergyType BetheBlochPDG::getTotalEnergyLoss(setup::Stack::particle_type const& vP,
-                                                  GrammageType const 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);
   }
 
-  ProcessReturn BetheBlochPDG::doContinuous(setup::Stack::particle_type& p,
-                                            setup::Trajectory const& t) {
-    if (p.getChargeNumber() == 0) return ProcessReturn::Ok;
+  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().getIntegratedGrammage(t, t.getLength());
-    CORSIKA_LOG_DEBUG("EnergyLoss pid={}, z={}, dX={} g/cm2", p.getPID(),
-                      p.getChargeNumber(), dX / 1_g * square(1_cm));
-    HEPEnergyType dE = getTotalEnergyLoss(p, dX);
-    auto E = p.getEnergy();
-    const auto Ekin = E - p.getMass();
-    auto Enew = E + dE;
-    CORSIKA_LOG_DEBUG("EnergyLoss  dE={} MeV, E={} GeV, Ekin={} GeV, Enew={} GeV",
-                      dE / 1_MeV, E / 1_GeV, Ekin / 1_GeV, Enew / 1_GeV);
-    p.setEnergy(Enew);
-    updateMomentum(p, Enew);
-    fillProfile(t, dE);
+        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;
   }
 
-  LengthType BetheBlochPDG::getMaxStepLength(setup::Stack::particle_type const& vParticle,
-                                             setup::Trajectory const& vTrack) const {
+  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 dEdX = -getTotalEnergyLoss(vParticle, dX) / dX; // dE > 0
-    //~ auto const Ekin = vParticle.getEnergy() - vParticle.getMass();
-
-    // in any case: never go below 0.99*emCut_ This needs to be
-    // slightly smaller than emCut_ since, either this Step is limited
-    // by energy_lim, then the particle is stopped in a very short
-    // range (before doing anythin else) and is then removed
-    // instantly. The exact position where it reaches emCut is not
-    // important, the important fact is that its E_kin is zero
-    // afterwards.
-    //
-    const auto energy = vParticle.getEnergy();
-    auto energy_lim = std::max(0.9 * energy, 0.99 * emCut_);
-
+    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 BetheBlochPDG::updateMomentum(corsika::setup::Stack::particle_type& vP,
-                                     HEPEnergyType Enew) {
-    HEPMomentumType Pnew = elab2plab(Enew, vP.getMass());
-    auto pnew = vP.getMomentum();
-    vP.setMomentum(pnew * Pnew / pnew.getNorm());
-  }
-
-  void BetheBlochPDG::fillProfile(setup::Trajectory const& vTrack,
-                                  const HEPEnergyType dE) {
-
-    GrammageType const grammageStart = shower_axis_.getProjectedX(vTrack.getPosition(0));
-    GrammageType const grammageEnd = shower_axis_.getProjectedX(vTrack.getPosition(1));
-    const auto deltaX = grammageEnd - grammageStart;
-
-    int binStart = grammageStart / dX_;
-    if (binStart < 0) return;
-    int binEnd = grammageEnd / dX_;
-    if (binEnd > int(profile_.size() - 1)) return;
-    if (deltaX < dX_threshold_) return;
-
-    CORSIKA_LOG_DEBUG("energy deposit of -dE={} between {} and {}", -dE, grammageStart,
-                      grammageEnd);
-
-    auto energyCount = HEPEnergyType::zero();
+  template <typename TOutput>
+  inline YAML::Node BetheBlochPDG<TOutput>::getConfig() const {
 
-    auto fill = [&](const int bin, const double weight) {
-      auto const increment = -dE * weight;
-      profile_[bin] += increment;
-      energyCount += increment;
-
-      CORSIKA_LOG_DEBUG("filling bin {} with weight {} : {} ", bin, weight, increment);
-    };
-
-    // fill longitudinal profile
-    if (binStart == binEnd) {
-      fill(binStart, 1);
-    } else {
-      fill(binStart, ((1 + binStart) * dX_ - grammageStart) / deltaX);
-      fill(binEnd, (grammageEnd - binEnd * dX_) / deltaX);
-      for (int bin = binStart + 1; bin < binEnd; ++bin) { fill(bin, 1); }
-    }
-
-    CORSIKA_LOG_DEBUG("total energy added to histogram: {} ", energyCount);
-  }
-
-  void BetheBlochPDG::printProfile() const {
-    std::ofstream file("EnergyLossProfile.dat");
-    file << "# EnergyLoss profile" << std::endl
-         << "# lower X bin edge [g/cm2]  dE/dX [GeV/g/cm2]\n";
-    double const deltaX = dX_ / 1_g * square(1_cm);
-    for (size_t i = 0; i < profile_.size(); ++i) {
-      file << std::scientific << std::setw(15) << i * deltaX << std::setw(15)
-           << profile_.at(i) / (deltaX * 1_GeV) << '\n';
-    }
-    file.close();
-  }
-
-  HEPEnergyType BetheBlochPDG::getTotal() const {
-    return std::accumulate(profile_.cbegin(), profile_.cend(), HEPEnergyType::zero());
+    YAML::Node node;
+    node["type"] = "BetheBlochPDG";
+    return node;
   }
 
-  void BetheBlochPDG::showResults() const {
-    CORSIKA_LOG_INFO(
-        " ******************************\n"
-        " PROCESS::ContinuousProcess: \n"
-        " energy lost dE (GeV)      : {}\n  ",
-        energy_lost_ / 1_GeV);
-  }
-
-  void BetheBlochPDG::reset() { energy_lost_ = 0_GeV; }
-
 } // namespace corsika