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>
@@ -24,12 +23,13 @@ namespace corsika {
       : TOutput(std::forward<TArgs>(args)...) {}
 
   template <typename TOutput>
-  template <typename TParticle, typename TTrack>
+  template <typename TParticle>
   inline ProcessReturn LongitudinalProfile<TOutput>::doContinuous(
-      TParticle const& particle, TTrack const& track, bool const) {
+      Step<TParticle> const& step, bool const) {
 
-    auto const pid = particle.getPID();
-    this->write(track, pid, 1.0); // weight hardcoded so far
+    auto const pid = step.getParticlePre().getPID();
+    this->write(step.getPositionPre(), step.getPositionPost(), pid,
+                step.getParticlePre().getWeight()); // weight hardcoded so far
     return ProcessReturn::Ok;
   }
 

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 {
@@ -21,9 +20,9 @@ namespace corsika {
       , 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:
      */
@@ -34,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.");
@@ -46,14 +45,17 @@ 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
-    double const weight = 1.; // particle.getWeight()
-    this->write(particle.getPID(), energy, displacement.dot(xAxis_),
-                displacement.dot(yAxis_), 0_m, weight);
+    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_);
 
@@ -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()) {
@@ -100,6 +102,8 @@ namespace corsika {
     // 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()};

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
@@ -17,19 +16,23 @@ namespace corsika {
   inline ParticleCut<TOutput>::ParticleCut(HEPEnergyType const eEleCut,
                                            HEPEnergyType const ePhoCut,
                                            HEPEnergyType const eHadCut,
-                                           HEPEnergyType const eMuCut, bool const inv,
+                                           HEPEnergyType const eMuCut,
+                                           HEPEnergyType const eTauCut, bool const inv,
                                            TArgs&&... outputArgs)
       : TOutput(std::forward<TArgs>(outputArgs)...)
       , cut_electrons_(eEleCut)
       , cut_photons_(ePhoCut)
-      , cut_muons_(eMuCut)
       , cut_hadrons_(eHadCut)
+      , cut_muons_(eMuCut)
+      , cut_tau_(eTauCut)
       , doCutInv_(inv) {
-    for (auto p : get_all_particles()) {
+    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)
@@ -40,7 +43,8 @@ namespace corsika {
         "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);
   }
 
   template <typename TOutput>
@@ -70,10 +74,8 @@ namespace corsika {
   }
 
   template <typename TOutput>
-  template <typename TParticle>
-  inline bool ParticleCut<TOutput>::isBelowEnergyCut(TParticle const& vP) const {
-    auto const energyLab = vP.getKineticEnergy();
-    auto const pid = vP.getPID();
+  inline bool ParticleCut<TOutput>::isBelowEnergyCut(
+      Code const pid, HEPEnergyType const energyLab) const {
     // nuclei
     if (is_nucleus(pid)) {
       // calculate energy per nucleon
@@ -85,25 +87,22 @@ namespace corsika {
   }
 
   template <typename TOutput>
-  template <typename TParticle>
-  inline bool ParticleCut<TOutput>::checkCutParticle(TParticle const& particle) {
+  inline bool ParticleCut<TOutput>::checkCutParticle(Code const pid,
+                                                     HEPEnergyType const kine_energy,
+                                                     TimeType const timePost) const {
 
-    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, E={} GeV, m={} "
         "GeV",
-        pid, particle.getPDG(), kine_energy / 1_GeV, energy / 1_GeV,
-        particle.getMass() / 1_GeV);
-    CORSIKA_LOG_DEBUG("p={}", particle.asString());
+        pid, get_PDG(pid), kine_energy / 1_GeV, energy / 1_GeV, get_mass(pid) / 1_GeV);
     if (doCutInv_ && is_neutrino(pid)) {
       CORSIKA_LOG_DEBUG("removing inv. particle...");
       return true;
-    } else if (isBelowEnergyCut(particle)) {
+    } else if (isBelowEnergyCut(pid, kine_energy)) {
       CORSIKA_LOG_DEBUG("removing low en. particle...");
       return true;
-    } else if (particle.getTime() > 10_ms) {
+    } else if (timePost > 10_ms) {
       CORSIKA_LOG_DEBUG("removing OLD particle...");
       return true;
     } else {
@@ -120,9 +119,10 @@ namespace corsika {
     HEPEnergyType energy_event = 0_GeV; // per event counting for printout
     auto particle = vS.begin();
     while (particle != vS.end()) {
-      if (checkCutParticle(particle)) {
-        HEPEnergyType Ekin = particle.getKineticEnergy();
-        this->write(particle.getPosition(), particle.getPID(), Ekin);
+      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
@@ -131,12 +131,16 @@ namespace corsika {
   }
 
   template <typename TOutput>
-  template <typename TParticle, typename TTrajectory>
-  inline ProcessReturn ParticleCut<TOutput>::doContinuous(TParticle& particle,
-                                                          TTrajectory const&,
+  template <typename TParticle>
+  inline ProcessReturn ParticleCut<TOutput>::doContinuous(Step<TParticle>& step,
                                                           bool const) {
-    if (checkCutParticle(particle)) {
-      this->write(particle.getPosition(), particle.getPID(), particle.getKineticEnergy());
+    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;
@@ -152,7 +156,8 @@ namespace corsika {
     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 hadros is {} 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_) {
@@ -171,6 +176,7 @@ namespace corsika {
     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;

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
@@ -27,7 +26,9 @@ namespace corsika {
       : 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>
   inline StackInspector<TStack>::~StackInspector() {}
@@ -54,54 +55,62 @@ namespace corsika {
       }
     }
 
-    std::chrono::system_clock::time_point const now = std::chrono::system_clock::now();
-    std::chrono::duration<double> const elapsed_seconds = now - StartTime_; // seconds
-    auto const dE = E0_ - Etot;
-    if (dE < dE_threshold_) return;
-    double const progress = dE / E0_;
+    if ((E0_ - Etot) < dE_threshold_) return;
 
-    // for printout
-    std::time_t const now_time = std::chrono::system_clock::to_time_t(now);
-    std::time_t const start_time = std::chrono::system_clock::to_time_t(StartTime_);
+    // 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
 
-    if (progress > 0) {
-
-      double const eta_seconds = elapsed_seconds.count() / progress;
-      std::chrono::system_clock::time_point const eta =
-          StartTime_ + std::chrono::seconds((int)eta_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 eta_time = std::chrono::system_clock::to_time_t(eta);
+      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);
 
-      int const yday0 = std::localtime(&start_time)->tm_yday;
-      int const yday1 = std::localtime(&eta_time)->tm_yday;
+      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;
 
-      CORSIKA_LOG_INFO(
+      std::stringstream ETA_string;
+      ETA_string << std::put_time(std::localtime(&ETA_time), "%T %Z");
+
+      CORSIKA_LOG_DEBUG(
           "StackInspector: "
           " time={}"
-          ", running={} seconds"
+          ", running={:.1f} seconds"
           " ( {:.1f}%)"
           ", nStep={}"
           ", stackSize={}"
-          ", Estack={} GeV"
+          ", Estack={:.1f} GeV"
           ", ETA={}{}",
-          std::put_time(std::localtime(&now_time), "%T"), elapsed_seconds.count(),
+          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)),
-          std::put_time(std::localtime(&eta_time), "%T"));
-    } else {
-      CORSIKA_LOG_INFO(
-          "StackInspector: "
-          " time={}"
-          ", running={} seconds"
-          " ( {:.1f}%)"
-          ", nStep={}"
-          ", stackSize={}"
-          ", Estack={} GeV"
-          ", ETA={}{}",
-          std::put_time(std::localtime(&now_time), "%T"), 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();
+      }
     }
   }
 

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