corsika/detail/media/HomogeneousMedium.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

corsika/detail/media/InhomogeneousMedium.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

corsika/detail/media/LayeredSphericalAtmosphereBuilder.hpp

 /*
  * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
  *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
  */
 
 #pragma once

corsika/detail/media/LayeredSphericalAtmosphereBuilder.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

corsika/detail/media/LinearApproximationIntegrator.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

corsika/detail/media/MediumPropertyModel.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

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>
@@ -46,7 +46,7 @@ namespace corsika {
   }
 
   template <typename TFunction>
-  inline auto NuclearComposition::getWeighted(TFunction const& func) const {
+  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;
@@ -66,7 +66,7 @@ namespace corsika {
   } // namespace corsika
 
   template <typename TFunction>
-  inline auto NuclearComposition::getWeightedSum(TFunction const& func) const
+  inline auto NuclearComposition::getWeightedSum(TFunction func) const
       -> decltype(func(std::declval<Code>())) {
     using ResultQuantity = decltype(func(std::declval<Code>()));
 
@@ -108,9 +108,23 @@ namespace corsika {
       throw std::runtime_error("incompatible vector sigma as input");
     }
 
-    std::discrete_distribution channelDist(
-        WeightProviderIterator(numberFractions_.begin(), sigma.begin()),
-        WeightProviderIterator(numberFractions_.end(), sigma.end()));
+    auto zip_beg = boost::make_zip_iterator(
+        boost::make_tuple(numberFractions_.cbegin(), sigma.cbegin()));
+    auto zip_end = boost::make_zip_iterator(
+        boost::make_tuple(numberFractions_.cend(), sigma.cend()));
+    using zip_iter_type = decltype(zip_beg);
+
+    auto const mult_func = [](zip_iter_type::value_type const& zipit) -> double {
+      return zipit.get<0>() * zipit.get<1>().magnitude();
+    };
+
+    using transform_iter_type =
+        boost::transform_iterator<decltype(mult_func), zip_iter_type, double, double>;
+
+    auto trans_beg = transform_iter_type{zip_beg, mult_func};
+    auto trans_end = transform_iter_type{zip_end, mult_func};
+
+    std::discrete_distribution channelDist{trans_beg, trans_end};
 
     auto const iChannel = channelDist(randomStream);
     return components_[iChannel];

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>
@@ -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();
     };
 
@@ -65,7 +71,7 @@ namespace corsika {
     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_) {

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

corsika/detail/media/SlidingPlanarTabular.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

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

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

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
@@ -52,13 +51,31 @@ 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>
-  inline 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); };
   }

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>
-  inline WeightProviderIterator<AConstIterator, BConstIterator>::WeightProviderIterator(
-      AConstIterator a, BConstIterator b)
-      : aIter_(a)
-      , bIter_(b) {}
-
-  template <class AConstIterator, class BConstIterator>
-  inline typename WeightProviderIterator<AConstIterator, BConstIterator>::value_type
-      WeightProviderIterator<AConstIterator, BConstIterator>::operator*() const {
-    return ((*aIter_) * (*bIter_)).magnitude();
-  }
-
-  template <class AConstIterator, class BConstIterator>
-  inline WeightProviderIterator<AConstIterator, BConstIterator>&
-  WeightProviderIterator<AConstIterator,
-                         BConstIterator>::operator++() { // prefix ++
-    ++aIter_;
-    ++bIter_;
-    return *this;
-  }
-
-  template <class AConstIterator, class BConstIterator>
-  inline bool WeightProviderIterator<AConstIterator, BConstIterator>::operator==(
-      WeightProviderIterator other) {
-    return aIter_ == other.aIter_;
-  }
-
-  template <class AConstIterator, class BConstIterator>
-  inline bool WeightProviderIterator<AConstIterator, BConstIterator>::operator!=(
-      WeightProviderIterator other) {
-    return !(*this == other);
-  }
-
-} // namespace corsika

corsika/detail/modules/HadronicElasticModel.inl

 /*
  * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
  *
- * This software is distributed under the terms of the 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

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 <cmath>
 #include <iomanip>
 #include <limits>
+#include <utility>
 
 namespace corsika {
 
-  inline 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>
-  inline ProcessReturn LongitudinalProfile::doContinuous(TParticle const& vP,
-                                                         TTrack const& vTrack,
-                                                         bool const) {
-    auto const pid = vP.getPID();
-
-    GrammageType const grammageStart = shower_axis_.getProjectedX(vTrack.getPosition(0));
-    GrammageType const grammageEnd = shower_axis_.getProjectedX(vTrack.getPosition(1));
-
-    CORSIKA_LOG_DEBUG("longprof: pos1={} m, pos2={}, X1={} g/cm2, X2={} g/cm2",
-                      vTrack.getPosition(0).getCoordinates() / 1_m,
-                      vTrack.getPosition(1).getCoordinates() / 1_m,
-                      grammageStart / 1_g * square(1_cm),
-                      grammageEnd / 1_g * square(1_cm));
+  template <typename TOutput>
+  template <typename... TArgs>
+  inline LongitudinalProfile<TOutput>::LongitudinalProfile(TArgs&&... args)
+      : TOutput(std::forward<TArgs>(args)...) {}
 
-    // Note: particle may go also "upward", thus, grammageEnd<grammageStart
-    int const binStart = std::ceil(grammageStart / dX_);
-    int const binEnd = std::floor(grammageEnd / dX_);
-
-    for (int b = binStart; b <= binEnd; ++b) {
-      if (pid == Code::Photon) {
-        profiles_.at(b)[ProfileIndex::Photon]++;
-      } 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]++;
-      } else if (is_neutrino(pid)) {
-        profiles_.at(b)[ProfileIndex::Invisible]++;
-      }
-    }
+  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;
   }
 
-  inline void LongitudinalProfile::save(std::string const& filename, const int width,
-                                        const int precision) {
-    CORSIKA_LOG_DEBUG("Write longprof to {}", filename);
-    std::ofstream f{filename};
-    f << "# X / g·cm¯², photon, e+, e-, mu+, mu-, all hadrons, neutrinos" << 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.
  */
 
 namespace corsika {
 
   template <typename TTracking, typename TOutput>
+  template <typename... TArgs>
   ObservationPlane<TTracking, TOutput>::ObservationPlane(Plane const& obsPlane,
                                                          DirectionVector const& x_axis,
-                                                         bool deleteOnHit)
-      : plane_(obsPlane)
-      , deleteOnHit_(deleteOnHit)
-      , energy_ground_(0_GeV)
-      , count_ground_(0)
+                                                         bool const deleteOnHit,
+                                                         TArgs&&... args)
+      : TOutput(std::forward<TArgs>(args)...)
+      , plane_(obsPlane)
       , xAxis_(x_axis.normalized())
-      , yAxis_(obsPlane.getNormal().cross(xAxis_)) {}
+      , yAxis_(obsPlane.getNormal().cross(xAxis_))
+      , deleteOnHit_(deleteOnHit) {}
 
   template <typename TTracking, typename TOutput>
-  template <typename TParticle, typename TTrajectory>
+  template <typename TParticle>
   inline ProcessReturn ObservationPlane<TTracking, TOutput>::doContinuous(
-      TParticle& particle, TTrajectory& step, bool const stepLimit) {
+      Step<TParticle>& step, bool const stepLimit) {
     /*
        The current step did not yet reach the ObservationPlane, do nothing now and wait:
      */
@@ -33,7 +33,7 @@ namespace corsika {
       if (deleteOnHit_) {
         // since this is basically a bug, it cannot be tested LCOV_EXCL_START
         LengthType const check =
-            (particle.getPosition() - plane_.getCenter()).dot(plane_.getNormal());
+            (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.");
@@ -45,24 +45,26 @@ namespace corsika {
         return ProcessReturn::Ok;
     }
 
-    HEPEnergyType const energy = particle.getEnergy();
-    Point const pointOfIntersection = step.getPosition(1);
+    HEPEnergyType const kineticEnergy = step.getEkinPost();
+    Point const pointOfIntersection = step.getPositionPost();
     Vector const displacement = pointOfIntersection - plane_.getCenter();
+    DirectionVector const direction = step.getDirectionPost();
 
     // add our particles to the output file stream
-    this->write(particle.getPID(), energy, displacement.dot(xAxis_),
-                displacement.dot(yAxis_), particle.getTime());
+    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;
       return ProcessReturn::ParticleAbsorbed;
     } else {
       return ProcessReturn::Ok;
     }
-  }
+  } // namespace corsika
 
   template <typename TTracking, typename TOutput>
   template <typename TParticle, typename TTrajectory>
@@ -77,7 +79,7 @@ namespace corsika {
 
     TimeType const timeOfIntersection = intersection.getEntry();
     CORSIKA_LOG_TRACE("timeOfIntersection={}", timeOfIntersection);
-    if (timeOfIntersection < TimeType::zero()) {
+    if (timeOfIntersection <= TimeType::zero()) {
       return std::numeric_limits<double>::infinity() * 1_m;
     }
     if (timeOfIntersection > trajectory.getDuration()) {
@@ -90,17 +92,6 @@ namespace corsika {
     return trajectory.getLength(fractionOfIntersection);
   }
 
-  template <typename TTracking, typename TOutput>
-  inline void ObservationPlane<TTracking, TOutput>::showResults() const {
-    CORSIKA_LOG_INFO(
-        " ******************************\n"
-        " ObservationPlane: \n"
-        " energy an ground (GeV)     :  {}\n"
-        " no. of particles at ground :  {}\n"
-        " ******************************",
-        energy_ground_ / 1_GeV, count_ground_);
-  }
-
   template <typename TTracking, typename TOutput>
   inline YAML::Node ObservationPlane<TTracking, TOutput>::getConfig() const {
     using namespace units::si;
@@ -110,7 +101,9 @@ namespace corsika {
 
     // basic info
     node["type"] = "ObservationPlane";
-    node["units"] = "m"; // add default units for values
+    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()};
@@ -144,10 +137,4 @@ namespace corsika {
     return node;
   }
 
-  template <typename TTracking, typename TOutput>
-  inline void ObservationPlane<TTracking, TOutput>::reset() {
-    energy_ground_ = 0_GeV;
-    count_ground_ = 0;
-  }
-
 } // 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

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
@@ -12,168 +11,136 @@
 
 namespace corsika {
 
-  inline ParticleCut::ParticleCut(HEPEnergyType const eEleCut,
-                                  HEPEnergyType const ePhoCut,
-                                  HEPEnergyType const eHadCut, HEPEnergyType const eMuCut,
-                                  bool const inv)
-      : doCutEm_(false)
-      , doCutInv_(inv)
-      , energy_cut_(0_GeV)
-      , energy_timecut_(0_GeV)
-      , energy_emcut_(0_GeV)
-      , energy_invcut_(0_GeV)
-      , em_count_(0)
-      , inv_count_(0)
-      , energy_count_() {
-    for (auto p : get_all_particles()) {
+  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_threshold(p, eHadCut);
+        set_kinetic_energy_propagation_threshold(p, eHadCut);
       else if (is_muon(p))
-        set_kinetic_energy_threshold(p, eMuCut);
+        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_threshold(p, eEleCut);
+        set_kinetic_energy_propagation_threshold(p, eEleCut);
       else if (p == Code::Photon)
-        set_kinetic_energy_threshold(p, ePhoCut);
+        set_kinetic_energy_propagation_threshold(p, ePhoCut);
     }
-    set_kinetic_energy_threshold(Code::Nucleus, eHadCut);
+    set_kinetic_energy_propagation_threshold(Code::Nucleus, eHadCut);
     CORSIKA_LOG_DEBUG(
         "setting kinetic energy thresholds: electrons = {} GeV, photons = {} GeV, "
         "hadrons = {} GeV, "
         "muons = {} GeV",
-        eEleCut / 1_GeV, ePhoCut / 1_GeV, eHadCut / 1_GeV, eMuCut / 1_GeV);
+        "tau = {} GeV", eEleCut / 1_GeV, ePhoCut / 1_GeV, eHadCut / 1_GeV, eMuCut / 1_GeV,
+        eTauCut / 1_GeV);
   }
 
-  inline ParticleCut::ParticleCut(HEPEnergyType const eHadCut, HEPEnergyType const eMuCut,
-                                  bool const inv)
-      : doCutEm_(true)
-      , doCutInv_(inv)
-      , energy_cut_(0_GeV)
-      , energy_timecut_(0_GeV)
-      , energy_emcut_(0_GeV)
-      , energy_invcut_(0_GeV)
-      , em_count_(0)
-      , inv_count_(0)
-      , energy_count_(0) {
-
+  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()) {
-      if (is_hadron(p))
-        set_kinetic_energy_threshold(p, eHadCut);
-      else if (is_muon(p))
-        set_kinetic_energy_threshold(p, eMuCut);
+      set_kinetic_energy_propagation_threshold(p, eCut);
     }
-    set_kinetic_energy_threshold(Code::Nucleus, eHadCut);
-    CORSIKA_LOG_DEBUG(
-        "setting thresholds: hadrons = {} GeV, "
-        "muons = {} GeV",
-        eHadCut / 1_GeV, eMuCut / 1_GeV);
+    set_kinetic_energy_propagation_threshold(Code::Nucleus, eCut);
+    CORSIKA_LOG_DEBUG("setting kinetic energy threshold {} GeV", eCut / 1_GeV);
   }
 
-  inline ParticleCut::ParticleCut(HEPEnergyType const eCut, bool const em, bool const inv)
-      : doCutEm_(em)
-      , doCutInv_(inv)
-      , energy_cut_(0_GeV)
-      , energy_timecut_(0_GeV)
-      , energy_emcut_(0_GeV)
-      , energy_invcut_(0_GeV)
-      , em_count_(0)
-      , inv_count_(0)
-      , energy_count_(0) {
-    for (auto p : get_all_particles()) set_kinetic_energy_threshold(p, eCut);
-    set_kinetic_energy_threshold(Code::Nucleus, eCut);
-    CORSIKA_LOG_DEBUG("setting kinetic energy threshold for all particles to {} GeV",
-                      eCut / 1_GeV);
-  }
-
-  inline ParticleCut::ParticleCut(
-      std::unordered_map<Code const, HEPEnergyType const> const& eCuts, bool const em,
-      bool const inv)
-      : doCutEm_(em)
-      , doCutInv_(inv)
-      , energy_cut_(0_GeV)
-      , energy_timecut_(0_GeV)
-      , energy_emcut_(0_GeV)
-      , energy_invcut_(0_GeV)
-      , em_count_(0)
-      , inv_count_(0)
-      , energy_count_(0) {
-    set_kinetic_energy_thresholds(eCuts);
+  template <typename TOutput>
+  template <typename... TArgs>
+  inline ParticleCut<TOutput>::ParticleCut(
+      std::unordered_map<Code const, HEPEnergyType const> const& eCuts, bool const inv,
+      TArgs&&... args)
+      : TOutput(std::forward<TArgs>(args)...)
+      , doCutInv_(inv) {
+    for (auto const& cut : eCuts) {
+      set_kinetic_energy_propagation_threshold(cut.first, cut.second);
+    }
     CORSIKA_LOG_DEBUG("setting threshold particles individually");
   }
 
-  template <typename TParticle>
-  inline bool ParticleCut::isBelowEnergyCut(TParticle const& vP) const {
-    auto const energyLab = vP.getKineticEnergy();
-    auto const pid = vP.getPID();
+  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 < calculate_kinetic_energy_threshold(pid));
+      return (ElabNuc < get_kinetic_energy_propagation_threshold(pid));
     } else {
-      return (energyLab < calculate_kinetic_energy_threshold(pid));
+      return (energyLab < get_kinetic_energy_propagation_threshold(pid));
     }
   }
 
-  inline bool ParticleCut::isInvisible(Code const& vCode) const {
-    return is_neutrino(vCode);
-  }
+  template <typename TOutput>
+  inline bool ParticleCut<TOutput>::checkCutParticle(Code const pid,
+                                                     HEPEnergyType const kine_energy,
+                                                     TimeType const timePost) const {
 
-  template <typename TParticle>
-  inline bool ParticleCut::checkCutParticle(TParticle const& particle) {
-
-    Code const pid = particle.getPID();
-    HEPEnergyType const kine_energy = particle.getKineticEnergy();
-    HEPEnergyType const energy = particle.getEnergy();
+    HEPEnergyType const energy = kine_energy + get_mass(pid);
     CORSIKA_LOG_DEBUG(
-        "ParticleCut: checking {} ({}), E_kin= {} GeV, EcutTot={} GeV, E={} GeV, m={} "
+        "ParticleCut: checking {} ({}), E_kin= {} GeV, E={} GeV, m={} "
         "GeV",
-        pid, particle.getPDG(), kine_energy / 1_GeV,
-        (energy_emcut_ + energy_invcut_ + energy_cut_ + energy_timecut_) / 1_GeV,
-        energy / 1_GeV, particle.getMass() / 1_GeV);
-    CORSIKA_LOG_DEBUG("p={}", particle.asString());
-    if (doCutEm_ && is_em(pid)) {
-      CORSIKA_LOG_DEBUG("removing em. particle...");
-      energy_emcut_ += kine_energy;
-      em_count_ += 1;
-      energy_event_ += kine_energy;
-      return true;
-    } else if (doCutInv_ && is_neutrino(pid)) {
+        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...");
-      energy_invcut_ += kine_energy;
-      inv_count_ += 1;
-      energy_event_ += kine_energy;
       return true;
-    } else if (isBelowEnergyCut(particle)) {
+    } else if (isBelowEnergyCut(pid, kine_energy)) {
       CORSIKA_LOG_DEBUG("removing low en. particle...");
-      energy_cut_ += kine_energy;
-      energy_count_ += 1;
-      energy_event_ += kine_energy;
       return true;
-    } else if (particle.getTime() > 10_ms) {
+    } else if (timePost > 10_ms) {
       CORSIKA_LOG_DEBUG("removing OLD particle...");
-      energy_timecut_ += kine_energy;
-      energy_event_ += kine_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
   }
 
+  template <typename TOutput>
   template <typename TStackView>
-  inline void ParticleCut::doSecondaries(TStackView& vS) {
-    energy_event_ = 0_GeV; // per event counting for printout
+  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);
+    CORSIKA_LOG_DEBUG("Event cut: {} GeV", energy_event / 1_GeV);
   }
 
-  template <typename TParticle, typename TTrajectory>
-  inline ProcessReturn ParticleCut::doContinuous(TParticle& particle, TTrajectory const&,
-                                                 bool 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");
       // signal to upstream code that this particle was deleted
       return ProcessReturn::ParticleAbsorbed;
@@ -181,34 +148,40 @@ namespace corsika {
     return ProcessReturn::Ok;
   }
 
-  inline void ParticleCut::printThresholds() {
-    for (auto p : get_all_particles()) {
-      auto const Eth = calculate_kinetic_energy_threshold(p);
-      CORSIKA_LOG_INFO("kinetic energy threshold for particle {} is {} GeV", p,
-                       Eth / 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);
     }
   }
 
-  inline void ParticleCut::showResults() {
-    CORSIKA_LOG_INFO(
-        "\n ******************************\n "
-        " kinetic energy removed by cut of electromagnetic (GeV): {} (number: {})\n "
-        " kinetic energy removed by cut of invisible (GeV): {} (number: {})\n "
-        " kinetic energy removed by kinetic energy cut (GeV): {} (number: {}) \n "
-        " kinetic energy removed by time cut (GeV): {} \n"
-        " ******************************",
-        energy_emcut_ / 1_GeV, em_count_, energy_invcut_ / 1_GeV, inv_count_,
-        energy_cut_ / 1_GeV, energy_count_, energy_timecut_ / 1_GeV);
-  }
-
-  inline void ParticleCut::reset() {
-    energy_emcut_ = 0_GeV;
-    em_count_ = 0;
-    energy_invcut_ = 0_GeV;
-    inv_count_ = 0;
-    energy_cut_ = 0_GeV;
-    energy_count_ = 0;
-    energy_timecut_ = 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