corsika/detail/stack/history/HistoryObservationPlane.inl

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#include <corsika/logging/Logging.h>
+#include <corsika/stack/history/HistoryObservationPlane.hpp>
+
+#include <boost/histogram/ostream.hpp>
+
+namespace corsika::history {
+
+  template <typename TStack>
+  inline HistoryObservationPlane<TStack> HistoryObservationPlane(TStack const& stack,
+                                                                 Plane const& obsPlane,
+                                                                 bool deleteOnHit)
+      : stack_{stack}
+      , plane_{obsPlane}
+      , deleteOnHit_{deleteOnHit} {}
+
+  template <typename TStack>
+  template <typename TParticle, typename TTrajectory>
+  inline ProcessReturn HistoryObservationPlane<TStack> DoContinuous(
+      TParticle const& particle, TTrajectory const& trajectory) {
+    TimeType const timeOfIntersection =
+        (plane_.getCenter() - trajectory.getR0()).dot(plane_.getNormal()) /
+        trajectory.getV0().dot(plane_.getNormal());
+
+    if (timeOfIntersection < TimeType::zero()) { return ProcessReturn::Ok; }
+
+    if (plane_.isAbove(trajectory.getR0()) == plane_.isAbove(trajectory.getPosition(1))) {
+      return ProcessReturn::Ok;
+    }
+
+    CORSIKA_LOG_DEBUG(fmt::format("HistoryObservationPlane: Particle detected: pid={}",
+                                  particle.getPID()));
+
+    auto const pid = particle.getPID();
+    if (particles::isMuon(pid)) { fillHistoryHistogram(particle); }
+
+    if (deleteOnHit_) {
+      return ProcessReturn::ParticleAbsorbed;
+    } else {
+      return ProcessReturn::Ok;
+    }
+  }
+
+  template <typename TStack>
+  template <typename TParticle, typename TTrajectory>
+  inline LengthType HistoryObservationPlane<TStack> MaxStepLength(
+      TParticle const&, TTrajectory const& trajectory) {
+    TimeType const timeOfIntersection =
+        (plane_.getCenter() - trajectory.getR0()).dot(plane_.getNormal()) /
+        trajectory.getV0().dot(plane_.getNormal());
+
+    if (timeOfIntersection < TimeType::zero()) {
+      return std::numeric_limits<double>::infinity() * 1_m;
+    }
+
+    auto const pointOfIntersection = trajectory.getPosition(timeOfIntersection);
+    return (trajectory.getR0() - pointOfIntersection).norm() * 1.0001;
+  }
+
+  template <typename TStack>
+  template <typename TParticle>
+  inline void HistoryObservationPlane<TStack> fillHistoryHistogram(
+      TParticle const& muon) {
+    double const muon_energy = muon.getEnergy() / 1_GeV;
+
+    int genctr{0};
+    Event const* event = muon.getEvent().get();
+    while (event) {
+      auto const projectile = stack_.cfirst() + event->projectileIndex();
+      if (event->eventType() == EventType::Interaction) {
+        genctr++;
+        double const projEnergy = projectile.getEnergy() / 1_GeV;
+        int const pdg = static_cast<int>(particles::getPDG(projectile.getPID()));
+
+        histogram_(muon_energy, projEnergy, pdg);
+      }
+      event = event->parentEvent().get(); // projectile.getEvent().get();
+    }
+  }
+} // namespace corsika::history

corsika/framework/analytics/ClassTimer.hpp

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+// Another possibility:
+// https://en.wikibooks.org/wiki/More_C%2B%2B_Idioms/Execute-Around_Pointer
+// for a more global approach
+//
+// In this case here only a single function is measured via member function pointer.
+
+#pragma once
+
+#include <chrono>
+#include <utility>
+
+#include <corsika/framework/analytics/Timer.hpp>
+
+namespace corsika {
+
+  template <class TClass, typename TTimer>
+  class ClassTimerImpl : public TTimer {
+    static_assert(
+        is_timer_v<TTimer>,
+        "TTimer is not a timer!"); // Better
+                                   // https://en.cppreference.com/w/cpp/language/constraints
+                                   // but not available in C++17
+
+  protected:
+    /// Reference to the class object on which the function should be called
+    TClass& obj_;
+
+  public:
+    ClassTimerImpl(TClass& obj)
+        : obj_(obj){};
+  };
+
+  /// Measure the runtime of a single class function
+  /**
+   * @tparam TClassFunc Type of the member function pointer that should be wrapped
+   * @tparam TFunc Actual function of the type defined in TClass
+   */
+  template <typename TClassFunc, TClassFunc TFunc,
+            typename TTimer =
+                Timer<std::chrono::high_resolution_clock, std::chrono::microseconds>>
+  class ClassTimer;
+
+  /// Measure the runtime of a single class function
+  /** Specialisation to capture exact information about the composition of the member
+   * function pointer used.
+   *
+   *  This class wrapes a single function and allowes the measureing of its runtime if it
+   * called via the "call(...)" function
+   *
+   * @tparam TClass Class of the function that should be wrapped
+   * @tparam TRet   Return value of the wrapped function
+   * @tparam TArgs  Arguments passed to the wrapped function
+   * @tparam TFuncPtr Actual function of the type defined by TRet
+   * TClass::TFuncPtr(TArgs...)
+   */
+  template <typename TClass, typename TRet, typename... TArgs,
+            TRet (TClass::*TFuncPtr)(TArgs...), typename TTimer>
+  class ClassTimer<TRet (TClass::*)(TArgs...), TFuncPtr, TTimer>
+      : public ClassTimerImpl<TClass, TTimer> {
+  private:
+  public:
+    ClassTimer(TClass& obj);
+
+    /// Executes the wrapped function
+    /** This function executes and measure the runtime of the wrapped function with the
+     * highest precision available (high_resolution_clock).
+     *
+     * @param args Arguments are perfect forwarded to the wrapped function.
+     * @return Returns the return value of the wrapped function. This value get copied
+     * during the process and therefore must be copy constructible!
+     */
+    TRet call(TArgs... args);
+  };
+
+  /// Specialisation for member functions without return value
+  template <typename TClass, typename... TArgs, void (TClass::*TFuncPtr)(TArgs...),
+            typename TTimer>
+  class ClassTimer<void (TClass::*)(TArgs...), TFuncPtr, TTimer>
+      : public ClassTimerImpl<TClass, TTimer> {
+  public:
+    ClassTimer(TClass& obj);
+
+    void call(TArgs... args);
+  };
+
+  /// Specialisation for const member functions
+  template <typename TClass, typename TRet, typename... TArgs,
+            TRet (TClass::*TFuncPtr)(TArgs...) const, typename TTimer>
+  class ClassTimer<TRet (TClass::*)(TArgs...) const, TFuncPtr, TTimer>
+      : public ClassTimerImpl<TClass, TTimer> {
+  public:
+    ClassTimer(TClass& obj);
+
+    TRet call(TArgs... args);
+  };
+
+  /// Specialisation for const member functions without return value
+  template <typename TClass, typename... TArgs, void (TClass::*TFuncPtr)(TArgs...) const,
+            typename TTimer>
+  class ClassTimer<void (TClass::*)(TArgs...) const, TFuncPtr, TTimer>
+      : public ClassTimerImpl<TClass, TTimer> {
+  public:
+    ClassTimer(TClass& obj);
+
+    void call(TArgs... args);
+  };
+
+} // namespace corsika
+
+#include <corsika/detail/framework/analytics/ClassTimer.inl>
\ No newline at end of file

corsika/framework/analytics/FunctionTimer.hpp

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+#pragma once
+
+#include <chrono>
+#include <utility>
+#include <type_traits>
+
+#include <corsika/framework/analytics/Timer.hpp>
+
+namespace corsika {
+
+  /// Wraps and measures the runtime of a single function type object
+  /**
+   *
+   * @tparam TFunc funtion pointer that should be wrapped
+   * @tparam TClock type of the clock that should be used for measurements, default is
+   * high_resolution_clock
+   * @tparam TDuration type of std::duration to measure the elapsed time, default is
+   * microseconds
+   */
+  template <typename TFunc, class TTimer = Timer<std::chrono::high_resolution_clock,
+                                                 std::chrono::microseconds>>
+  class FunctionTimer : public TTimer {
+    static_assert(
+        is_timer_v<TTimer>,
+        "TTimer is not a timer!"); // Better
+                                   // https://en.cppreference.com/w/cpp/language/constraints
+                                   // but not available in C++17
+  public:
+    /** Constructs the wrapper with the given functionpointer
+     *  @param f  Function or functor whose runtime should be measured
+     **/
+    FunctionTimer(TFunc f);
+
+    /** Functor for calling the wrapped function
+     *  This functor calls the wrapped function and measures the elapsed time between call
+     *and return. The return value needs to be copy constructible.
+     *  @tparam TArgs Parameter types that are forwarded to the function. The use of
+     *correct types is the responsibility of the user, no checks are done.
+     *  @param args Arguments are forwarded to the wrapped function. This method does not
+     *support overloaded function resolution.
+     *  @return The return value of the wrapped function is temporarily copied and then
+     *returned by value
+     **/
+    template <typename... TArgs>
+    auto operator()(TArgs&&... args) -> std::invoke_result_t<TFunc, TArgs...>;
+
+  private:
+    TFunc function_;
+  };
+
+} // namespace corsika
+
+#include <corsika/detail/framework/analytics/FunctionTimer.inl>
\ No newline at end of file

corsika/framework/analytics/Timer.hpp

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <chrono>
+#include <utility>
+
+#include <type_traits>
+
+namespace corsika {
+
+  template <typename TClock = std::chrono::high_resolution_clock,
+            typename TDuration = std::chrono::microseconds>
+  class Timer {
+  protected:
+    /// Default clock used for time measurement
+    using clock_type = TClock;
+
+    /// Internal resolution of the time measurement
+    using duration_type = TDuration;
+
+    /// Startpoint of time measurement
+    typename clock_type::time_point start_;
+
+    /// Measured runtime of the function
+    duration_type timeDiff_;
+
+  public:
+    Timer()
+        : timeDiff_(0){};
+
+    /// Start the timer
+    void startTimer() { start_ = clock_type::now(); }
+
+    /// Stop the timer
+    void stopTimer() {
+      timeDiff_ = std::chrono::duration_cast<duration_type>(clock_type::now() - start_);
+    }
+
+    /**
+     * Returns the runtime of the last call to the wrapped function
+     * @return Returns the measured runtime of the wrapped function/functor in the unit
+     * given by TDuration
+     **/
+    duration_type getTime() const { return timeDiff_; }
+  };
+
+  std::false_type is_timer_impl(...);
+  template <typename T, typename U>
+  std::true_type is_timer_impl(Timer<T, U> const volatile&);
+
+  template <typename T>
+  constexpr bool is_timer_v =
+      std::is_same_v<decltype(is_timer_impl(std::declval<T&>())), std::true_type>;
+
+} // namespace corsika

corsika/framework/core/Cascade.hpp

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/corsika.hpp>
+
+#include <corsika/framework/process/ProcessReturn.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/core/Logging.hpp>
+#include <corsika/framework/random/ExponentialDistribution.hpp>
+#include <corsika/framework/random/RNGManager.hpp>
+#include <corsika/framework/random/UniformRealDistribution.hpp>
+#include <corsika/framework/stack/SecondaryView.hpp>
+#include <corsika/framework/geometry/FourVector.hpp>
+
+#include <corsika/media/Environment.hpp>
+
+#include <corsika/stack/history/HistoryStackExtension.hpp>
+
+#include <cassert>
+#include <cmath>
+#include <limits>
+#include <type_traits>
+
+namespace corsika {
+
+  /**
+   *
+   * The Cascade class is constructed from template arguments making
+   * it very versatile. Via the template arguments physics models are
+   * plugged into the cascade simulation.
+   *
+   * <b>TTracking</b> must be a class according to the
+   * TrackingInterface providing the functions:
+   *
+   * @code
+   * auto getTrack(particle_type const& p)</auto>,
+   * with the return type <code>geometry::Trajectory<Line>
+   * @endcode
+   *
+   * <b>TProcessList</b> must be a ProcessSequence.   *
+   * <b>Stack</b> is the storage object for particle data, i.e. with
+   * particle class type `Stack::particle_type`.
+   */
+  template <typename TTracking, typename TProcessList, typename TOutput, typename TStack>
+  class Cascade {
+
+    typedef typename TStack::stack_view_type stack_view_type;
+
+    typedef typename TStack::particle_type particle_type;
+
+    typedef std::remove_pointer_t<decltype(std::declval<particle_type>().getNode())>
+        volume_tree_node_type;
+
+    typedef typename volume_tree_node_type::IModelProperties medium_interface_type;
+
+  public:
+    /**
+     * @name constructors
+     * @{
+     * Cascade class cannot be default constructed, but needs a valid
+     * list of physics processes for configuration at construct time.
+     */
+    Cascade() = delete;
+    Cascade(Cascade const&) = default;
+    Cascade(Cascade&&) = default;
+    ~Cascade() = default;
+    Cascade& operator=(Cascade const&) = default;
+    Cascade(Environment<medium_interface_type> const& env, TTracking& tr,
+            TProcessList& pl, TOutput& out, TStack& stack);
+    //! @}
+
+    /**
+     * set the nodes for all particles on the stack according to their numerical
+     * position.
+     */
+    void setNodes();
+
+    /**
+     * The Run function is the main simulation loop, which processes
+     * particles from the Stack until the Stack is empty.
+     */
+    void run();
+
+    /**
+     * Force an interaction of the top particle of the stack at its current position.
+     * Note that setNodes() or an equivalent procedure needs to be called first if you
+     * want to call forceInteraction() for the primary interaction.
+     * Incompatible with forceDecay()
+     */
+    void forceInteraction();
+
+    /**
+     * Force an decay of the top particle of the stack at its current position.
+     * Note that setNodes() or an equivalent procedure needs to be called first if you
+     * want to call forceDecay() for the primary interaction.
+     * Incompatible with forceInteraction()
+     */
+    void forceDecay();
+
+  private:
+    /**
+     * The Step function is executed for each particle from the
+     * stack. It will calcualte geometric transport of the particles,
+     * and apply continuous and stochastic processes to it, which may
+     * lead to energy losses, scattering, absorption, decays and the
+     * production of secondary particles.
+     *
+     * New particles produced in one step are subject to further
+     * processing, e.g. thinning, etc.
+     */
+    void step(particle_type& vParticle);
+
+    ProcessReturn decay(stack_view_type& view, InverseTimeType initial_inv_decay_time);
+    ProcessReturn interaction(stack_view_type& view, FourMomentum const& projectileP4,
+                              NuclearComposition const& composition,
+                              CrossSectionType const initial_cross_section);
+    void setEventType(stack_view_type& view, history::EventType);
+
+    // data members
+    Environment<medium_interface_type> const& environment_;
+    TTracking& tracking_;
+    TProcessList& sequence_;
+    TOutput& output_;
+    TStack& stack_;
+    default_prng_type& rng_ = RNGManager<>::getInstance().getRandomStream("cascade");
+    bool forceInteraction_;
+    bool forceDecay_;
+    unsigned int count_ = 0;
+
+    // but this here temporarily. Should go into dedicated file later:
+    const char* c8_ascii_ =
+        R"V0G0N(
+  ,ad8888ba,     ,ad8888ba,    88888888ba    ad88888ba   88  88      a8P          db              ad88888ba   
+ d8"'    `"8b   d8"'    `"8b   88      "8b  d8"     "8b  88  88    ,88'          d88b            d8"     "8b  
+d8'            d8'        `8b  88      ,8P  Y8,          88  88  ,88"           d8'`8b           Y8a     a8P  
+88             88          88  88aaaaaa8P'  `Y8aaaaa,    88  88,d88'           d8'  `8b           "Y8aaa8P"   
+88             88          88  88""""88'      `"""""8b,  88  8888"88,         d8YaaaaY8b          ,d8"""8b,   
+Y8,            Y8,        ,8P  88    `8b            `8b  88  88P   Y8b       d8""""""""8b        d8"     "8b  
+ Y8a.    .a8P   Y8a.    .a8P   88     `8b   Y8a     a8P  88  88     "88,    d8'        `8b       Y8a     a8P  
+  `"Y8888Y"'     `"Y8888Y"'    88      `8b   "Y88888P"   88  88       Y8b  d8'          `8b       "Y88888P"
+       )V0G0N";
+  };
+
+} // namespace corsika
+
+#include <corsika/detail/framework/core/Cascade.inl>

corsika/framework/core/EnergyMomentumOperations.hpp

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/core/PhysicalGeometry.hpp>
+#include <corsika/framework/core/PhysicalConstants.hpp>
+
+/**
+ * \file EnergyMomentumOperations.hpp
+ *
+ * Relativistic energy momentum calculations.
+ */
+
+namespace corsika {
+
+  namespace detail {
+    using HEPEnergyTypeSqr = decltype(1_GeV * 1_GeV);
+  }
+
+  /**
+   * \f[m^2=E_{tot}^2-p^2\f]
+   *
+   * @param E total energy.
+   * @param p particle momentum.
+   * @return HEPEnergyType-squared
+   */
+  auto constexpr calculate_mass_sqr(HEPEnergyType const E, HEPMomentumType const p) {
+    return (E - p) * (E + p);
+  }
+
+  /**
+   * \f[m=sqrt(E_{tot}^2-p^2) \f]
+   *
+   * @param E total energy.
+   * @param p particle momentum.
+   * @return HEPEnergyType
+   */
+  HEPEnergyType constexpr calculate_mass(HEPEnergyType const E, HEPMomentumType const p) {
+    return sqrt(calculate_mass_sqr(E, p));
+  }
+
+  /**
+   * \f[p^2=E_{tot}^2-m^2\f]
+   *
+   * @param E total energy.
+   * @param m particle mass.
+   * @return HEPEnergyType-squared
+   */
+  auto constexpr calculate_momentum_sqr(HEPEnergyType const E, HEPMassType const m) {
+    return (E - m) * (E + m);
+  }
+
+  /**
+   * \f[p=sqrt(E_{tot}^2-m^2) \f]
+   *
+   * @param E total energy.
+   * @param m particle mass.
+   * @return HEPEnergyType
+   */
+  HEPEnergyType constexpr calculate_momentum(HEPEnergyType const E, HEPMassType const m) {
+    return sqrt(calculate_momentum_sqr(E, m));
+  }
+
+  /**
+   * \f[E_{tot}^2=p^2+m^2\f]
+   *
+   * @param p momentum.
+   * @param m particle mass.
+   * @return HEPEnergyType-squared
+   */
+  auto constexpr calculate_total_energy_sqr(HEPMomentumType const p,
+                                            HEPMassType const m) {
+    return p * p + m * m;
+  }
+
+  /**
+   * \f[E_{tot}=sqrt(p^2+m^2)\f]
+   *
+   * @param p momentum.
+   * @param m particle mass.
+   * @return HEPEnergyType
+   */
+  HEPEnergyType constexpr calculate_total_energy(HEPMomentumType const p,
+                                                 HEPMassType const m) {
+    return sqrt(calculate_total_energy_sqr(p, m));
+  }
+
+  /**
+   * \f[E_{kin}=sqrt(p^2+m^2) - m \f]
+   *
+   * @param p momentum.
+   * @param m particle mass.
+   * @return HEPEnergyType
+   */
+  HEPEnergyType constexpr calculate_kinetic_energy(HEPMomentumType const p,
+                                                   HEPMassType const m) {
+    return calculate_total_energy(p, m) - m;
+  }
+
+  /**
+   * \f[E_{lab}=(\sqrt{s}^2 - m_{proj}^2 - m_{targ}^2) / (2 m_{targ}) \f]
+   *
+   * @param p momentum.
+   * @param m particle mass.
+   * @return HEPEnergyType
+   */
+  HEPEnergyType constexpr calculate_lab_energy(detail::HEPEnergyTypeSqr sqrtS_sqr,
+                                               HEPMassType const m_proj,
+                                               HEPMassType const m_targ) {
+    return (sqrtS_sqr - static_pow<2>(m_proj) - static_pow<2>(m_targ)) / (2 * m_targ);
+  }
+
+  /**
+   * \f[E_{com}=sqrt{2 * m_{proj} * m_{targ} * E_{lab} + m_{proj}^2 + m_{targ}^2} \f]
+   *
+   * @param E lab. energy.
+   * @param m particle mass.
+   * @return HEPEnergyType
+   */
+  HEPEnergyType constexpr calculate_com_energy(HEPEnergyType Elab,
+                                               HEPMassType const m_proj,
+                                               HEPMassType const m_targ) {
+    return sqrt(2 * Elab * m_targ + static_pow<2>(m_proj) + static_pow<2>(m_targ));
+  }
+
+} // namespace corsika

corsika/framework/core/Logging.hpp

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+/**
+ *  @file Logging.hpp
+ *
+ * CORSIKA8 logging utilities.
+ *
+ * See testLogging.cpp for a complete set of examples for
+ * how the logging functions should be used.
+ */
+
+#pragma once
+
+// Configure some behaviour of sdlog.
+// This must be done before spdlog is included.
+
+// use the coarse system clock. This is *much* faster
+// but introduces a timestamp error of O(10 ms) which is fine for us.
+#ifndef SPDLOG_CLOCK_COARSE
+#define SPDLOG_CLOCK_COARSE
+#endif
+// do not create a default logger (we provide our own "corsika" logger)
+#ifndef SPDLOG_DISABLE_DEFAULT_LOGGER
+#define SPDLOG_DISABLE_DEFAULT_LOGGER
+#endif
+// use __PRETTY_FUNCTION__ instead of __FUNCTION__ where
+// printing function names in trace statements. This is much
+// nicer than __FUNCTION__ under GCC/clang.
+#ifndef SPDLOG_FUNCTION
+#define SPDLOG_FUNCTION __PRETTY_FUNCTION__
+#endif
+
+// if this is a Debug build, include debug messages in objects
+#ifdef _C8_DEBUG_
+// trace is the highest level of logging (ALL messages will be printed)
+#ifndef SPDLOG_ACTIVE_LEVEL
+#define SPDLOG_ACTIVE_LEVEL SPDLOG_LEVEL_TRACE
+#endif
+#else // otherwise, remove everything but "error" and worse messages
+#ifndef SPDLOG_ACTIVE_LEVEL
+#define SPDLOG_ACTIVE_LEVEL SPDLOG_LEVEL_DEBUG
+#endif
+#endif
+
+#include <spdlog/fmt/ostr.h> // will output whenerver a streaming operator is found
+#include <spdlog/sinks/stdout_color_sinks.h>
+#include <spdlog/spdlog.h>
+
+namespace corsika {
+
+  /*
+   * The default pattern for CORSIKA8 loggers.
+   */
+  const std::string minimal_pattern{"[%n:%^%-8l%$] %v"};
+  const std::string default_pattern{"[%n:%^%-8l%$(%s:%#)] %v"};
+  const std::string source_pattern{"[%n:%^%-8l%$(%s:%!:%#)] %v"};
+
+  /**
+   * Create a new C8-style logger.
+   *
+   * Use this if you are explicitly (and can guarantee) that you
+   * are creating a logger for the first time. It is recommended
+   * that for regular usage, the `get_logger` function is used instead
+   * as that will also create the logger if it has not yet been created.
+   *
+   * Calling `create_logger` twice to create the same logger will
+   * result in an spdlog duplicate exception.
+   *
+   * @param name           The unique name of the logger.
+   * @param defaultlog     If True, set this as the default logger.
+   * @returns              The constructed and formatted logger.
+   */
+  std::shared_ptr<spdlog::logger> create_logger(std::string const& name,
+                                                bool const defaultlog = false);
+
+  /**
+   * Get a smart pointer to an existing logger.
+   *
+   * This should be the default method for code to obtain a
+   * logger. If the logger *does not* exist, it is *created* and
+   * returned to the caller.
+   *
+   * This should be preferred over `create_logger`.
+   *
+   * @param name    The name of the logger to get.
+   * @param defaultlog   If True, make this the default logger.
+   * @returns              The constructed and formatted logger.
+   */
+  std::shared_ptr<spdlog::logger> get_logger(std::string const& name,
+                                             bool const defaultlog = false);
+
+  /**
+   * The default "corsika" logger.
+   */
+  static inline std::shared_ptr<spdlog::logger> corsika_logger =
+      get_logger("corsika", true);
+
+  // many of these free functions are special to the logging
+  // infrastructure so we hide them in the corsika::logging namespace.
+  namespace logging {
+
+    // bring spdlog into the corsika::logging namespace
+    using namespace spdlog;
+
+    /**
+     * Set the default log level for all *newly* created loggers.
+     *
+     *  @param minlevel The minimum log level required to print.
+     *
+     */
+    auto set_default_level(level::level_enum const minlevel) -> void;
+
+    /**
+     * Add the source (filename, line no) info to the logger.
+     *
+     * @param logger  The logger to set the level of.
+     *
+     */
+    template <typename TLogger>
+    auto add_source_info(TLogger& logger) -> void;
+
+    /**
+     * Reset the logging pattern to the default.
+     *
+     * @param logger  The logger to set the level of.
+     *
+     */
+    template <typename TLogger>
+    auto reset_pattern(TLogger& logger) -> void;
+
+  } // namespace logging
+
+// define our macro-style loggers
+// these use the default "corsika" logger
+#define CORSIKA_LOG_TRACE SPDLOG_TRACE
+#define CORSIKA_LOG_DEBUG SPDLOG_DEBUG
+#define CORSIKA_LOG_INFO SPDLOG_INFO
+#define CORSIKA_LOG_WARN SPDLOG_WARN
+#define CORSIKA_LOG_ERROR SPDLOG_ERROR
+#define CORSIKA_LOG_CRITICAL SPDLOG_CRITICAL
+
+// and the specific logger versions
+// these take a logger instance as their first argument
+#define CORSIKA_LOGGER_TRACE SPDLOG_LOGGER_TRACE
+#define CORSIKA_LOGGER_DEBUG SPDLOG_LOGGER_DEBUG
+#define CORSIKA_LOGGER_INFO SPDLOG_LOGGER_INFO
+#define CORSIKA_LOGGER_WARN SPDLOG_LOGGER_WARN
+#define CORSIKA_LOGGER_ERROR SPDLOG_LOGGER_ERROR
+#define CORSIKA_LOGGER_CRITICAL SPDLOG_LOGGER_CRITICAL
+
+} // namespace corsika
+
+#include <corsika/detail/framework/core/Logging.inl>
+#include <corsika/detail/framework/core/SpdlogSpecializations.inl>

corsika/framework/core/ParticleProperties.hpp

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+/**
+ *   @file ParticleProperties.hpp
+ *
+ * Interface to particle properties.
+ */
+
+#pragma once
+
+#include <array>
+#include <cstdint>
+#include <cmath>
+#include <iosfwd>
+#include <string_view>
+#include <type_traits>
+#include <unordered_map>
+
+#include <corsika/framework/core/PhysicalUnits.hpp>
+
+namespace corsika {
+
+  /**
+   * @defgroup Particles Particle Properties
+   *
+   * The properties of all particles are saved in static and flat
+   * arrays. There is a enum corsika::Code to identify each
+   * particle, and each individual particle has its own static class,
+   * which can be used to retrieve its physical properties.
+   *
+   * The properties of all elementary particles are accessible here. The data
+   * are taken from the Pythia ParticleData.xml file.
+   *
+   * Particle data can be accessed via global function in namespace corsika, or via
+   * static classes for each particle type. These classes all have the interface (example
+   * for the class corsika::Electron):
+   *
+   * @code{.cpp}
+   *   static constexpr Code code{Code::Electron};
+   *   static constexpr Code anti_code{Code::Positron};
+   *   static constexpr HEPMassType mass{corsika::get_mass(code)};
+   *   static constexpr ElectricChargeType charge{corsika::get_charge(code)};
+   *   static constexpr int charge_number{corsika::get_charge_number(code)};
+   *   static constexpr std::string_view name{corsika::get_name(code)};
+   *   static constexpr bool is_nucleus{corsika::is_nucleus(code)};
+   * @endcode
+   *
+   * The names, relations and properties of all particles known to CORSIKA 8 are listed
+   * below.
+   *
+   * **Note** on energy threshold on particle production as well as particle propagation.
+   * The functions:
+   * @code {.cpp}
+   * HEPEnergyType constexpr get_energy_production_threshold(Code const);
+   * void constexpr set_energy_production_threshold(Code const, HEPEnergyType const);
+   * @endcode
+   * can be used to tune the transition where explicit production of new particles, e.g.
+   * in Bremsstrahlung, is simulated versus a continuous handling of low-energy particles
+   * as generic energy losses. The default value for all particle types is 1 MeV.
+   *
+   * Furthermore, the functions:
+   * @code {.cpp}
+   * HEPEnergyType constexpr get_kinetic_energy_propagation_threshold(Code const);
+   * void constexpr set_kinetic_energy_propagation_threshold(Code const, HEPEnergyType
+   *                                                         const);
+   * @endcode
+   * are used to discard low energy particle during tracking. The default value for all
+   * particle types is 1 GeV.
+   *
+   * @addtogroup Particles
+   * @{
+   */
+
+  /**
+   * @enum Code
+   *
+   * The Code enum is the actual place to define CORSIKA 8 particle codes.
+   */
+  enum class Code : std::int32_t;
+
+  /**
+   * @enum PDGCode
+   *
+   * Specifically for PDG ids.
+   */
+  enum class PDGCode : std::int32_t;
+  /**
+   * Internal integer type for enum Code.
+   */
+  typedef std::underlying_type<Code>::type CodeIntType;
+
+  /**
+   * Internal integer type for enum PDGCode.
+   */
+  typedef std::underlying_type<PDGCode>::type PDGCodeIntType;
+} // namespace corsika
+
+// data arrays, etc., as generated automatically
+#include <corsika/framework/core/GeneratedParticleProperties.inc>
+
+namespace corsika {
+
+  // forward declarations to be used in GeneratedParticleProperties
+  struct full_name {}; //!< tag class for get_name()
+
+  int16_t constexpr get_charge_number(Code const);     //!< electric charge in units of e
+  ElectricChargeType constexpr get_charge(Code const); //!< electric charge
+  HEPMassType constexpr get_mass(Code const);          //!< mass
+
+  /**
+   * Get the kinetic energy propagation threshold.
+   *
+   * Particles are tracked only above the kinetic energy propagation threshold. Below
+   * this, they are discarded and removed. Sensible default values must be configured for
+   * a simulation.
+   */
+  HEPEnergyType get_kinetic_energy_propagation_threshold(Code const);
+
+  /**
+   * Set the kinetic energy propagation threshold object.
+   */
+  void set_kinetic_energy_propagation_threshold(Code const, HEPEnergyType const);
+
+  /**
+   * Get the particle production energy threshold.
+   *
+   * The (total) energy below which a particle is only handled stoachastically (no
+   * production below this energy). This is for example important for stochastic discrete
+   * Bremsstrahlung versus low-energy Bremsstrahlung as part of continuous energy losses.
+   */
+  HEPEnergyType get_energy_production_threshold(Code const); //!<
+
+  /**
+   * Set the particle production energy threshold in total energies.
+   */
+  void set_energy_production_threshold(Code const, HEPEnergyType const);
+
+  //! Particle code according to PDG, "Monte Carlo Particle Numbering Scheme"
+  PDGCode constexpr get_PDG(Code const);
+  PDGCode constexpr get_PDG(unsigned int const A, unsigned int const Z);
+  std::string_view constexpr get_name(Code const); //!< name of the particle as string
+  std::string get_name(Code,
+                       full_name); //!< get name of particle, including (A,Z) for nuclei
+  TimeType constexpr get_lifetime(Code const); //!< lifetime
+
+  bool constexpr is_hadron(Code const); //!< true if particle is hadron
+  bool constexpr is_em(Code const); //!< true if particle is electron, positron or photon
+  bool constexpr is_muon(Code const);     //!< true if particle is mu+ or mu-
+  bool constexpr is_neutrino(Code const); //!< true if particle is (anti-) neutrino
+  bool constexpr is_charged(Code const);  //!< true if particle is charged
+
+  /**
+   * @brief Creates the Code for a nucleus of type 10LZZZAAAI.
+   *
+   * @return internal nucleus Code
+   */
+  Code constexpr get_nucleus_code(size_t const A, size_t const Z);
+
+  /**
+   * Checks if Code corresponds to a nucleus.
+   *
+   * @return true if nucleus.
+   * @return false  if not nucleus.
+   */
+  bool constexpr is_nucleus(Code const);
+
+  /**
+   * Get the mass number A for nucleus.
+   *
+   * @return int size of nucleus.
+   */
+  size_t constexpr get_nucleus_A(
+      Code const); //!< returns A for hard-coded nucleus, otherwise 0
+
+  /**
+   * Get the charge number Z for nucleus.
+   *
+   * @return int charge of nucleus.
+   */
+  size_t constexpr get_nucleus_Z(
+      Code const); //!< returns Z for hard-coded nucleus, otherwise 0
+
+  /**
+   * @brief Calculates the mass of nucleus.
+   *
+   * @return HEPMassType the mass of (A,Z) nucleus, disregarding binding energy.
+   */
+  HEPMassType constexpr get_nucleus_mass(Code const code);
+
+  /**
+   * @brief Calculates the mass of nucleus.
+   *
+   * @return HEPMassType the mass of (A,Z) nucleus, disregarding binding energy.
+   */
+  HEPMassType constexpr get_nucleus_mass(unsigned int const A, unsigned int const Z);
+
+  /**
+   * @brief Get the nucleus name.
+   *
+   * @param code
+   * @return std::string_view
+   */
+  inline std::string get_nucleus_name(Code const code);
+
+  /**
+   * @brief convert PDG code to CORSIKA 8 internal code.
+   *
+   * @return Code internal code.
+   */
+  Code convert_from_PDG(PDGCode const);
+
+  /**
+   * @brief Returns list of all non-nuclei particles.
+   *
+   * @return std::initializer_list<Code> constexpr
+   */
+  std::initializer_list<Code> constexpr get_all_particles();
+
+  /**
+   * @brief Code output operator.
+   *
+   * The output stream operator for human-readable particle codes.
+   *
+   * @return std::ostream&
+   */
+  std::ostream& operator<<(std::ostream&, corsika::Code);
+
+  /** @}*/
+
+} // namespace corsika
+
+#include <corsika/detail/framework/core/ParticleProperties.inl>
+
+// constants in namespaces-like static classes, generated automatically
+#include <corsika/framework/core/GeneratedParticleClasses.inc>

corsika/framework/core/PhysicalConstants.hpp

+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/core/PhysicalUnits.hpp>
+
+/**
+ * \file PhysicalConstants.hpp
+ *
+ * Constants are defined with static units, based on the package
+ * (namespace) phys::units, imported in PhysicsUnits.hpp
+ *
+ * \namespace corsika::constants
+ *
+ * Physical and mathematical constants with units.
+ */
+
+namespace corsika::constants {
+
+  using namespace phys::units;
+
+  // acceleration of free-fall, standard
+  constexpr quantity<acceleration_d> g_sub_n{Rep(9.80665L) * meter / square(second)};
+
+  // Avogadro constant
+  constexpr quantity<dimensions<0, 0, 0, 0, 0, -1>> N_sub_A{Rep(6.02214199e+23L) / mole};
+
+  // elementary charge
+  constexpr quantity<electric_charge_d> e{Rep(1.6021766208e-19L) * coulomb};
+
+  // vacuum permittivity
+  constexpr quantity<dimensions<-3, -1, 4, 2>> epsilonZero{Rep(8.8541878128e-12L) *
+                                                           farad / meter};
+
+  // electronvolt
+  // constexpr quantity<hepenergy_d> eV{e / coulomb * joule};
+
+  // Planck constant
+  constexpr quantity<dimensions<2, 1, -1>> h{Rep(6.62606876e-34L) * joule * second};
+  constexpr quantity<dimensions<2, 1, -1>> hBar{h / (2 * M_PI)};
+
+  // speed of light in a vacuum
+  constexpr quantity<speed_d> c{Rep(299792458L) * meter / second};
+  constexpr auto cSquared = c * c;
+
+  // hbar * c
+  constexpr quantity<dimensions<1, 0, 0, 0, 0, 0, 0, 1>> hBarC{
+      Rep(1.973'269'78e-7L) * electronvolt * meter}; // from RPP 2018
+
+  auto constexpr invGeVsq = 1e-18 / (electronvolt * electronvolt);
+
+  // unified atomic mass unit
+  constexpr quantity<mass_d> u{Rep(1.6605402e-27L) * kilogram};
+
+  auto constexpr nucleonMass = 0.5 * (0.93827 + 0.93957) * 1e9 * electronvolt;
+
+  // molar gas constant
+  auto constexpr R = Rep(8.314'459'8) * joule / (mole * kelvin);
+
+  /**
+   * A namespace containing various Earth radii.
+   */
+  namespace EarthRadius {
+    static constexpr auto Mean{6'371'000 * meter};
+    static constexpr auto Geomagnetic_reference{6'371'200 * meter};
+    static constexpr auto Equatorial{6'378'137 * meter};
+    static constexpr auto Polar{6'356'752 * meter};
+    static constexpr auto PolarCurvature{6'399'593 * meter};
+  } // namespace EarthRadius
+
+  // etc.
+
+} // namespace corsika::constants

corsika/framework/core/PhysicalGeometry.hpp

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/geometry/Vector.hpp>
+
+/**
+ * \file PhysicalUnits.hpp
+ *
+ * Import and extend the phys::units package. The SI units are also imported into the
+ * `\namespace corsika`, since they are used everywhere as integral part of the framework.
+ */
+
+namespace corsika {
+
+  /**
+   * A 3D vector defined in a specific coordinate system with units HEPMomentumType
+   **/
+  using MomentumVector = Vector<hepmomentum_d>;
+
+  /**
+   * A 3D vector defined in a specific coordinate system with no units. But, note, this is
+   * not automatically normaliyed! It is not a "NormalVector".
+   **/
+  using DirectionVector = Vector<dimensionless_d>;
+
+  /**
+   * A 3D vector defined in a specific coordinate system with units "velocity_t".
+   *
+   **/
+  using VelocityVector = Vector<SpeedType::dimension_type>;
+
+  /**
+   * A 3D vector defined in a specific coordinate system with units "length_t".
+   *
+   **/
+  using LengthVector = Vector<length_d>;
+
+  /**
+   * A 3D vector defined in a specific coordinate system with units ElectricFieldType
+   **/
+  typedef Vector<ElectricFieldType::dimension_type> ElectricFieldVector;
+
+  /**
+   * A 3D vector defined in a specific coordinate system with units VectorPotentialType
+   **/
+  typedef Vector<VectorPotentialType::dimension_type> VectorPotential;
+
+} // namespace corsika

corsika/framework/core/PhysicalUnits.hpp

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+// the templated static-unit package we use:
+#include <corsika/framework/units/io.hpp>
+#include <corsika/framework/units/quantity.hpp>
+
+#include <corsika/framework/core/PhysicalConstants.hpp>
+
+/**
+ * \file PhysicalUnits.hpp
+ *
+ * Import and extend the phys::units package. The SI units are also imported into the
+ * `\namespace corsika`, since they are used everywhere as integral part of the framework.
+ */
+
+/*
+  It is essentially a bug of the phys_units package to define the
+  operator<< not in the same namespace as the types it is working
+  on. This breaks ADL (argument-dependent lookup). Here we "fix" this:
+ */
+namespace phys::units {
+  using phys::units::io::operator<<;
+} // namespace phys::units
+
+/**
+ * \namespace corsika::units
+ *
+ * Extension of the phys::units package.
+ *
+ */
+
+namespace corsika::units {
+  template <int N, typename T>
+  auto constexpr static_pow([[maybe_unused]] T x) {
+    if constexpr (N == 0) {
+      return 1;
+    } else if constexpr (N > 0) {
+      return x * static_pow<N - 1, T>(x);
+    } else {
+      return 1 / static_pow<-N, T>(x);
+    }
+  }
+} // namespace corsika::units
+
+/**
+ * \namespace corsika::units::si
+ *
+ * SI units as used mainly in CORSIKA8 as basedline.
+ *
+ */
+
+namespace corsika::units::si {
+
+  using namespace phys::units;
+  using namespace phys::units::literals;
+  using namespace phys::units::io;
+
+  using phys::units::io::operator<<;
+
+  /// defining momentum you suckers
+  /// dimensions, i.e. composition in base SI dimensions
+  using hepmomentum_d = phys::units::hepenergy_d;
+  using hepmass_d = phys::units::hepenergy_d;
+
+  /// defining cross section as area
+  using sigma_d = phys::units::area_d;
+
+  /// add the unit-types
+  using DimensionlessType = phys::units::quantity<phys::units::dimensionless_d, double>;
+  using LengthType = phys::units::quantity<phys::units::length_d, double>;
+  using TimeType = phys::units::quantity<phys::units::time_interval_d, double>;
+  using SpeedType = phys::units::quantity<phys::units::speed_d, double>;
+  using FrequencyType = phys::units::quantity<phys::units::frequency_d, double>;
+  using ElectricChargeType =
+      phys::units::quantity<phys::units::electric_charge_d, double>;
+  using HEPEnergyType = phys::units::quantity<phys::units::hepenergy_d, double>;
+  using MassType = phys::units::quantity<phys::units::mass_d, double>;
+  using HEPMassType = phys::units::quantity<hepmass_d, double>;
+  using MassDensityType = phys::units::quantity<phys::units::mass_density_d, double>;
+  using GrammageType = phys::units::quantity<phys::units::dimensions<-2, 1, 0>, double>;
+  using HEPMomentumType = phys::units::quantity<hepmomentum_d, double>;
+  using CrossSectionType = phys::units::quantity<area_d, double>;
+  using InverseLengthType =
+      phys::units::quantity<phys::units::dimensions<-1, 0, 0>, double>;
+  using InverseTimeType =
+      phys::units::quantity<phys::units::dimensions<0, 0, -1>, double>;
+  using InverseMassDensityType =
+      phys::units::quantity<phys::units::dimensions<3, -1, 0>, double>;
+  using InverseGrammageType =
+      phys::units::quantity<phys::units::dimensions<2, -1, 0>, double>;
+  using MagneticFluxType =
+      phys::units::quantity<phys::units::magnetic_flux_density_d, double>;
+  using ElectricFieldType =
+      phys::units::quantity<phys::units::dimensions<1, 1, -3, -1>, double>;
+  using VectorPotentialType =
+      phys::units::quantity<phys::units::dimensions<1, 1, -2, -1>, double>;
+
+  template <typename DimFrom, typename DimTo>
+  auto constexpr conversion_factor_HEP_to_SI() {
+    static_assert(DimFrom::dim1 == 0 && DimFrom::dim2 == 0 && DimFrom::dim3 == 0 &&
+                      DimFrom::dim4 == 0 && DimFrom::dim5 == 0 && DimFrom::dim6 == 0 &&
+                      DimFrom::dim7 == 0,
+                  "must be a pure HEP type");
+
+    static_assert(
+        DimTo::dim4 == 0 && DimTo::dim5 == 0 && DimTo::dim6 == 0 && DimTo::dim7 == 0,
+        "conversion possible only into L, M, T dimensions");
+
+    int constexpr e = DimFrom::dim8; // HEP dim.
+
+    int constexpr l = DimTo::dim1; // SI length dim.
+    int constexpr m = DimTo::dim2; // SI mass dim.
+    int constexpr t = DimTo::dim3; // SI time dim.
+
+    int constexpr p = m;
+    int constexpr q = -m - t;
+    static_assert(q == l + e - 2 * m, "HEP/SI dimension mismatch!");
+
+    return static_pow<-e>(constants::hBarC) * static_pow<p>(constants::hBar) *
+           static_pow<q>(constants::c);
+  }
+
+  template <typename DimFrom>
+  auto constexpr conversion_factor_SI_to_HEP() {
+    static_assert(DimFrom::dim4 == 0 && DimFrom::dim5 == 0 && DimFrom::dim6 == 0 &&
+                      DimFrom::dim7 == 0 && DimFrom::dim8 == 0,
+                  "must be pure L, M, T type");
+
+    int constexpr l = DimFrom::dim1; // SI length dim.
+    int constexpr m = DimFrom::dim2; // SI mass dim.
+    int constexpr t = DimFrom::dim3; // SI time dim.
+
+    int constexpr p = -m;
+    int constexpr q = m + t;
+    int constexpr e = m - t - l;
+
+    return static_pow<e>(constants::hBarC) * static_pow<p>(constants::hBar) *
+           static_pow<q>(constants::c);
+  }
+
+  template <typename DimTo, typename DimFrom>
+  auto constexpr convert_HEP_to_SI(quantity<DimFrom> q) {
+    return conversion_factor_HEP_to_SI<DimFrom, DimTo>() * q;
+  }
+
+  template <typename DimFrom>
+  auto constexpr convert_SI_to_HEP(quantity<DimFrom> q) {
+    return conversion_factor_SI_to_HEP<DimFrom>() * q;
+  }
+
+} // end namespace corsika::units::si
+
+/**
+ * @file PhysicalUnits
+ *
+ */
+
+namespace phys {
+  namespace units {
+    namespace literals {
+
+      /**
+       * Define new _XeV literals, alowing 10_GeV in the code.
+       * Define new _barn literal
+       */
+
+      QUANTITY_DEFINE_SCALING_LITERALS(eV, hepenergy_d, 1)
+
+      QUANTITY_DEFINE_SCALING_LITERALS(b, corsika::units::si::sigma_d,
+                                       magnitude(corsika::constants::barn))
+
+    } // namespace literals
+  }   // namespace units
+} // namespace phys
+
+// import into main \namespace corsika here:
+namespace corsika {
+  using namespace units;
+  using namespace units::si;
+  using namespace phys::units;
+  using namespace phys::units::literals;
+  using namespace phys::units::io;
+  using phys::units::io::operator<<;
+} // namespace corsika

corsika/framework/core/Step.hpp

+/*
+ * (c) Copyright 2022 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/core/PhysicalGeometry.hpp>
+#include <corsika/framework/geometry/Vector.hpp>
+#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/framework/geometry/StraightTrajectory.hpp>
+
+namespace corsika {
+
+  struct DeltaParticleState {
+  public:
+    DeltaParticleState(HEPEnergyType dEkin, TimeType dT, LengthVector const& ds,
+                       DirectionVector const& du)
+        : delta_Ekin_{dEkin}
+        , delta_time_{dT}
+        , displacement_{ds}
+        , delta_direction_{du} {}
+
+    HEPEnergyType delta_Ekin_ = HEPEnergyType::zero();
+    TimeType delta_time_ = TimeType::zero();
+    LengthVector displacement_;
+    DirectionVector delta_direction_;
+  };
+
+  template <typename TParticle>
+  class Step {
+  public:
+    template <typename TTrajectory>
+    Step(TParticle const& particle, TTrajectory const& track)
+        : particlePreStep_{particle} //~ , track_{track}
+        , diff_{HEPEnergyType::zero(), track.getDuration(1),
+                track.getPosition(1) - particle.getPosition(),
+                track.getDirection(1) - particle.getDirection()}
+
+    {}
+
+    HEPEnergyType const& add_dEkin(HEPEnergyType dEkin) {
+      diff_.delta_Ekin_ += dEkin;
+      return diff_.delta_Ekin_;
+    }
+
+    TimeType const& add_dt(TimeType dt) {
+      diff_.delta_time_ += dt;
+      return diff_.delta_time_;
+    }
+
+    LengthVector const& add_displacement(LengthVector const& dis) {
+      diff_.displacement_ += dis;
+      return diff_.displacement_;
+    }
+
+    DirectionVector const& add_dU(DirectionVector const& du) {
+      diff_.delta_direction_ += du;
+      return diff_.delta_direction_;
+    }
+
+    // getters for difference
+
+    HEPEnergyType getDiffEkin() const { return diff_.delta_Ekin_; }
+
+    TimeType getDiffT() const { return diff_.delta_time_; };
+
+    DirectionVector const& getDiffDirection() const { return diff_.delta_direction_; }
+
+    LengthVector const& getDisplacement() const { return diff_.displacement_; }
+
+    //! alias for getDisplacement()
+    LengthVector const& getDiffPosition() const { return getDisplacement(); }
+
+    // getters for absolute
+    TParticle const& getParticlePre() const { return particlePreStep_; }
+
+    HEPEnergyType getEkinPre() const { return getParticlePre().getKineticEnergy(); }
+
+    HEPEnergyType getEkinPost() const { return getEkinPre() + getDiffEkin(); }
+
+    TimeType getTimePre() const { return getParticlePre().getTime(); }
+
+    TimeType getTimePost() const { return getTimePre() + getDiffT(); }
+
+    DirectionVector const getDirectionPre() const {
+      return getParticlePre().getDirection();
+    }
+
+    VelocityVector getVelocityVector() const { return getDisplacement() / getDiffT(); }
+
+    StraightTrajectory getStraightTrack() const {
+      Line const line(getPositionPre(), getVelocityVector());
+      StraightTrajectory track(line, getDiffT());
+      return track;
+    }
+
+    DirectionVector getDirectionPost() const {
+      return (getDirectionPre() + getDiffDirection()).normalized();
+    }
+
+    Point const& getPositionPre() const { return getParticlePre().getPosition(); }
+
+    Point getPositionPost() const { return getPositionPre() + getDisplacement(); }
+
+  private:
+    TParticle const& particlePreStep_;
+    //~ TTrajectory const& track_; // TODO: perhaps remove
+    DeltaParticleState diff_;
+  };
+
+} // namespace corsika

corsika/framework/geometry/BaseTrajectory.hpp

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/geometry/Line.hpp>
+#include <corsika/framework/geometry/Point.hpp>
+
+namespace corsika {
+
+  /**
+   *
+   * A Trajectory is a description of a momvement of an object in
+   * three-dimensional space that describes the trajectory (connection
+   * between two Points in space), as well as the direction of motion
+   * at any given point.
+   *
+   * A Trajectory has a start `0` and an end `1`, where
+   * e.g. getPosition(0) returns the start point and getDirection(1)
+   * the direction of motion at the end. Values outside 0...1 are not
+   * defined.
+   *
+   * A Trajectory has a length in [m], getLength, a duration in [s], getDuration.
+   *
+   * Note: so far it is assumed that the speed (d|vec{r}|/dt) between
+   * start and end does not change and is constant for the entire
+   * Trajectory.
+   *
+   **/
+
+  class BaseTrajectory {
+
+  public:
+    virtual Point getPosition(double const u) const = 0;
+
+    virtual VelocityVector getVelocity(double const u) const = 0;
+
+    virtual DirectionVector getDirection(double const u) const = 0;
+
+    virtual TimeType getDuration(double const u = 1) const = 0;
+
+    virtual LengthType getLength(double const u = 1) const = 0;
+
+    virtual void setLength(LengthType const limit) = 0;
+
+    virtual void setDuration(TimeType const limit) = 0;
+  };
+
+} // namespace corsika

corsika/framework/geometry/BaseVector.hpp

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/geometry/CoordinateSystem.hpp>
+#include <corsika/framework/geometry/QuantityVector.hpp>
+
+#include <memory>
+
+namespace corsika {
+
+  /*!
+   * Common base class for Vector and Point.
+   *
+   * This holds a QuantityVector and a CoordinateSystem. The
+   * BaseVector manages resources for many geometry objects.
+   *
+   */
+  template <typename TDimension>
+  class BaseVector {
+
+  public:
+    BaseVector(CoordinateSystemPtr const& pCS, QuantityVector<TDimension> const& pQVector)
+        : quantityVector_(pQVector)
+        , cs_(pCS) {}
+
+    BaseVector() = delete; // we only want to creat initialized
+                           // objects
+    BaseVector(BaseVector const&) = default;
+    BaseVector(BaseVector&& a) = default;
+    BaseVector& operator=(BaseVector const&) = default;
+    ~BaseVector() = default;
+
+    CoordinateSystemPtr getCoordinateSystem() const;
+    void setCoordinateSystem(CoordinateSystemPtr const& cs) { cs_ = cs; }
+
+  protected:
+    QuantityVector<TDimension> const& getQuantityVector() const;
+    QuantityVector<TDimension>& getQuantityVector();
+    void setQuantityVector(QuantityVector<TDimension> const& v) { quantityVector_ = v; }
+
+  private:
+    QuantityVector<TDimension> quantityVector_;
+    CoordinateSystemPtr cs_;
+  };
+
+} // namespace corsika
+
+#include <corsika/detail/framework/geometry/BaseVector.inl>

corsika/framework/geometry/Box.hpp

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/framework/geometry/Line.hpp>
+#include <corsika/framework/geometry/IVolume.hpp>
+
+namespace corsika {
+
+  /**
+   * Describes a box in space
+   *
+   *  The center point and the orintation of the Box is set by
+   *  a CoordinateSystemPtr at construction and the sides extend
+   *  by x, y, z in both directions.
+   **/
+  class Box : public IVolume {
+
+  public:
+    Box(CoordinateSystemPtr cs, LengthType const x, LengthType const y,
+        LengthType const z);
+
+    Box(CoordinateSystemPtr cs, LengthType const side);
+
+    //! returns true if the Point p is within the sphere
+    bool contains(Point const& p) const override;
+
+    Point const& getCenter() const { return center_; };
+    CoordinateSystemPtr const getCoordinateSystem() const { return cs_; }
+
+    LengthType const getX() const { return x_; }
+    LengthType const getY() const { return y_; }
+    LengthType const getZ() const { return z_; }
+
+    std::string asString() const;
+
+    template <typename TDim>
+    void rotate(QuantityVector<TDim> const& axis, double const angle);
+
+  protected:
+    Point center_;
+    CoordinateSystemPtr cs_; // local coordinate system with center_ in coordinate (0, 0,
+                             // 0) and user defined orientation
+    LengthType x_;
+    LengthType y_;
+    LengthType z_;
+  };
+
+} // namespace corsika
+
+#include <corsika/detail/framework/geometry/Box.inl>

corsika/framework/geometry/CoordinateSystem.hpp

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+/**
+ * \file CoordinateSystem.hpp
+ **/
+
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/geometry/QuantityVector.hpp>
+#include <corsika/framework/core/Logging.hpp>
+
+#include <Eigen/Dense>
+#include <stdexcept>
+#include <memory>
+
+namespace corsika {
+
+  typedef Eigen::Transform<double, 3, Eigen::Affine> EigenTransform;
+  typedef Eigen::Translation<double, 3> EigenTranslation;
+
+  template <typename T>
+  class Vector; // fwd decl
+
+  class CoordinateSystem; // fwd decl
+  /**
+   * To refer to CoordinateSystems, only the CoordinateSystemPtr must be used.
+   */
+  using CoordinateSystemPtr = std::shared_ptr<CoordinateSystem const>;
+
+  /// this is the only way to create ONE unique root CS
+  CoordinateSystemPtr const& get_root_CoordinateSystem();
+
+  /**
+   * Creates new CoordinateSystemPtr by translation along \a vector
+   */
+  CoordinateSystemPtr make_translation(CoordinateSystemPtr const& cs,
+                                       QuantityVector<length_d> const& vector);
+
+  /**
+   * creates a new CoordinateSystem in which vVec points in direction of the new z-axis,
+   * \a vVec
+   */
+  template <typename TDim>
+  CoordinateSystemPtr make_rotationToZ(CoordinateSystemPtr const& cs,
+                                       Vector<TDim> const& vVec);
+
+  /**
+   * creates a new CoordinateSystem, rotated around axis by angle.
+   */
+  template <typename TDim>
+  CoordinateSystemPtr make_rotation(CoordinateSystemPtr const& cs,
+                                    QuantityVector<TDim> const& axis, double const angle);
+
+  /**
+   * creates a new CoordinateSystem, translated by \a translation and rotated around \a
+   * axis by \a angle.
+   */
+  template <typename TDim>
+  CoordinateSystemPtr make_translationAndRotation(
+      CoordinateSystemPtr const& cs, QuantityVector<length_d> const& translation,
+      QuantityVector<TDim> const& axis, double const angle);
+
+  /**
+   * A class to store the reference coordinate system for a geometric object
+   *
+   * A CoordinateSystem can only be created in reference and relative
+   * to other CoordinateSystems. Thus, the geometric
+   * transformation between all CoordinateSystems is always known and stored.
+   *
+   * The static (singleton) function \ref make_root_CoordinateSystem is
+   * the only way to create and access the global top-level
+   * CoordinateSystem obect. CoordinateSystem objects should be
+   * *abosulte* *only* handled in their form of CoordinateSystemPtr,
+   * which are shared_ptr that handle the lifetime of the entire
+   * CoordinateSystem hirarchy.
+   *
+   * Thus, new CoordinateSystem are only be created (via
+   * CoordinateSystemPtr) by transforing existing CoordinateSystem
+   * using: \ref make_rotationToZ, \ref make_rotation, or \ref
+   * make_translationAndRotation, see below.
+   *
+   * Warning: As a consequence, never try to access, modify, copy, the raw
+   * CoordinateSystem objects directly, this will almost certainly result in undefined
+   * behaviour. Only access, copy, handle them via CoordinateSystemPtr.
+   */
+
+  class CoordinateSystem {
+
+    /**
+     * Constructor only from referenceCS, given the transformation matrix transf
+     */
+    CoordinateSystem(CoordinateSystemPtr const& referenceCS, EigenTransform const& transf)
+        : referenceCS_(referenceCS)
+        , transf_(transf) {}
+
+    /**
+     * for creating the root CS
+     */
+    CoordinateSystem()
+        : referenceCS_(nullptr)
+        , transf_(EigenTransform::Identity()) {}
+
+  public:
+    // default resource allocation
+    CoordinateSystem(CoordinateSystem const&) = delete;
+    CoordinateSystem(CoordinateSystem&&) = delete;
+    CoordinateSystem& operator=(CoordinateSystem const& pCS) =
+        delete; // avoid making copies
+    ~CoordinateSystem() = default;
+
+    /**
+     * Checks, if this is the unique ROOT CS
+     */
+    bool isRoot() const { return !referenceCS_; }
+
+    CoordinateSystemPtr getReferenceCS() const;
+
+    EigenTransform const& getTransform() const;
+
+    bool operator==(CoordinateSystem const&) const;
+    bool operator!=(CoordinateSystem const&) const;
+
+  protected:
+    /**
+     * \name Friends
+     * Manipulation and creation functions.
+     * \{
+     **/
+
+    friend CoordinateSystemPtr const& get_root_CoordinateSystem();
+
+    friend CoordinateSystemPtr make_translation(CoordinateSystemPtr const& cs,
+                                                QuantityVector<length_d> const& vector);
+    template <typename TDim>
+    friend CoordinateSystemPtr make_rotationToZ(CoordinateSystemPtr const& cs,
+                                                Vector<TDim> const& vVec);
+    template <typename TDim>
+    friend CoordinateSystemPtr make_rotation(CoordinateSystemPtr const& cs,
+                                             QuantityVector<TDim> const& axis,
+                                             double const angle);
+    template <typename TDim>
+    friend CoordinateSystemPtr make_translationAndRotation(
+        CoordinateSystemPtr const& cs, QuantityVector<length_d> const& translation,
+        QuantityVector<TDim> const& axis, double const angle);
+
+    /** \} **/
+
+  private:
+    CoordinateSystemPtr referenceCS_;
+    EigenTransform transf_;
+  };
+
+  /**
+   * Transformation matrix from one reference system to another.
+   *
+   * returns the transformation matrix necessary to transform primitives with coordinates
+   * in \a pFrom to \a pTo, e.g.
+   * \f$ \vec{v}^{\text{(to)}} = \mathcal{M} \vec{v}^{\text{(from)}} \f$
+   * (\f$ \vec{v}^{(.)} \f$ denotes the coordinates/components of the component in
+   * the indicated CoordinateSystem).
+   *
+   * \todo make this a protected member of CoordinateSystem
+   */
+  EigenTransform get_transformation(CoordinateSystem const& c1,
+                                    CoordinateSystem const& c2);
+
+} // namespace corsika
+
+#include <corsika/detail/framework/geometry/CoordinateSystem.inl>

corsika/framework/geometry/FourVector.hpp

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/core/PhysicalGeometry.hpp>
+#include <type_traits>
+
+/**
+ * @file FourVector.hpp
+ * @author Ralf Ulrich
+ * @brief General FourVector object.
+ * @date 2021-10-16
+ */
+
+namespace corsika {
+
+  /**
+   * Description of physical four-vectors
+   *
+   * FourVector fully supports units, e.g. E in [GeV/c] and p in [GeV],
+   * or also t in [s] and r in [m], etc.
+   *
+   * However, for HEP applications it is also possible to use E and p
+   * both in [GeV].
+   *
+   * Thus, the input units of time-like and space-like coordinates
+   * must either be idential (e.g. GeV) or scaled by "c" as in
+   * [E/c]=[p].
+   *
+   * The FourVector can return its squared-norm \ref getNormSqr and its
+   * norm \ref getNorm, whereas norm is sqrt(abs(norm-squared)). The
+   * physical units are always calculated and returned properly.
+   *
+   * FourVector can also return if it is TimeLike, SpaceLike or PhotonLike.
+   *
+   * When a FourVector is initialized with a lvalue references,
+   * e.g. as `FourVector<TimeType&, Vector<length_d>&>`, references
+   * are also used as internal data types, which should lead to
+   * complete disappearance of the FourVector class during
+   * optimization.
+   */
+
+  template <typename TTimeType, typename TSpaceVecType>
+  class FourVector {
+
+  public:
+    using space_vec_type = typename std::decay<TSpaceVecType>::type;
+    using space_type = typename space_vec_type::quantity_type;
+    using time_type = typename std::decay<TTimeType>::type;
+
+    // check the types and the physical units here:
+    static_assert(std::is_same<time_type, space_type>::value ||
+                      std::is_same<time_type, decltype(std::declval<space_type>() /
+                                                       meter * second)>::value,
+                  "Units of time-like and space-like coordinates must either be idential "
+                  "(e.g. GeV) or [E/c]=[p]");
+
+    using norm_type = space_type;
+    using norm_square_type =
+        decltype(std::declval<norm_type>() * std::declval<norm_type>());
+
+  public:
+    // resource management
+    FourVector() = default;
+    FourVector(FourVector&&) = default;
+    FourVector(FourVector const&) = default;
+    FourVector& operator=(const FourVector&) = default;
+    ~FourVector() = default;
+
+    FourVector(TTimeType const& eT, TSpaceVecType const& eS)
+        : timeLike_(eT)
+        , spaceLike_(eS) {}
+
+    /**
+     * @return timeLike_
+     */
+    TTimeType getTimeLikeComponent() const;
+
+    /**
+     * @return spaceLike_
+     */
+    TSpaceVecType& getSpaceLikeComponents();
+
+    /**
+     *
+     * @return spaceLike_;
+     */
+    TSpaceVecType const& getSpaceLikeComponents() const;
+
+    /**
+     *
+     * @return \f$p_0^2 - \vec{p}^2\f$
+     */
+    norm_square_type getNormSqr() const;
+
+    /**
+     *
+     * @return \f$\sqrt(p_0^2 - \vec{p}^2)\f$
+     */
+    norm_type getNorm() const;
+
+    /**
+     * \todo FIXME: a better alternative would be to define an enumeration
+     * enum { SpaceLike =-1, TimeLike, LightLike } V4R_Category;
+     * and a method called  V4R_Category GetCategory() const;
+     * RU: then you have to decide in the constructor which avoids "lazyness"
+     **/
+    ///\return if \f$|p_0|>|\vec{p}|\f$
+    bool isTimelike() const;
+    ///\return if \f$|p_0|<|\vec{p}|\f$
+    bool isSpacelike() const;
+
+    /**
+     * \name Math operators (class members)
+     * @{
+     */
+    FourVector& operator+=(FourVector const&);
+    FourVector& operator-=(FourVector const&);
+    FourVector& operator*=(double const);
+    FourVector& operator/=(double const);
+    FourVector& operator/(double const);
+
+    /**
+     * Scalar product of two FourVectors.
+     *
+     *  Note that the product between two 4-vectors assumes that you use
+     *  the same "c" convention for both. Only the LHS vector is checked
+     *  for this. You cannot mix different conventions due to
+     *  unit-checking.
+     */
+    norm_type operator*(FourVector const& b);
+
+    /** @} */
+
+  protected:
+    // the data members
+    TTimeType timeLike_;
+    TSpaceVecType spaceLike_;
+
+    /**
+     * \name Free math operators
+     * We need to define them inline here since we do not want to
+     * implement them as template functions. They are valid only for
+     * the specific types as defined right here.
+     *
+     * Note, these are "free function" (even if they don't look as
+     *  such). Thus, even if the input object uses internal references
+     *  for storage, the free math operators, of course, must provide
+     *  value-copies.
+     * @{
+     *
+     */
+    friend FourVector<time_type, space_vec_type> operator+(FourVector const& a,
+                                                           FourVector const& b) {
+      return FourVector<time_type, space_vec_type>(a.timeLike_ + b.timeLike_,
+                                                   a.spaceLike_ + b.spaceLike_);
+    }
+
+    friend FourVector<time_type, space_vec_type> operator-(FourVector const& a,
+                                                           FourVector const& b) {
+      return FourVector<time_type, space_vec_type>(a.timeLike_ - b.timeLike_,
+                                                   a.spaceLike_ - b.spaceLike_);
+    }
+
+    friend FourVector<time_type, space_vec_type> operator*(FourVector const& a,
+                                                           double const b) {
+      return FourVector<time_type, space_vec_type>(a.timeLike_ * b, a.spaceLike_ * b);
+    }
+
+    friend FourVector<time_type, space_vec_type> operator/(FourVector const& a,
+                                                           double const b) {
+      return FourVector<time_type, space_vec_type>(a.timeLike_ / b, a.spaceLike_ / b);
+    }
+
+    /** @} */
+
+  private:
+    /**
+     * This function is there to automatically remove the eventual
+     * extra factor of "c" for the time-like quantity.
+     */
+    norm_square_type getTimeSquared() const;
+  };
+
+  /**
+   * streaming operator
+   */
+
+  template <typename TTimeType, typename TSpaceVecType>
+  std::ostream& operator<<(std::ostream& os,
+                           corsika::FourVector<TTimeType, TSpaceVecType> const& qv);
+
+  /**
+   * @typedef FourMomentum A FourVector with HEPEnergyType and MomentumVector.
+   */
+  typedef FourVector<HEPEnergyType, MomentumVector> FourMomentum;
+
+} // namespace corsika
+
+#include <corsika/detail/framework/geometry/FourVector.inl>

corsika/framework/geometry/Helix.hpp

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/core/PhysicalGeometry.hpp>
+#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/framework/geometry/Vector.hpp>
+
+#include <cmath>
+
+namespace corsika {
+
+  /*!
+   * Defines a helical path
+   *
+   * A Helix is defined by the cyclotron frequency \f$ \omega_c \f$, the initial
+   * Point r0 and
+   * the velocity vectors \f$ \vec{v}_{\parallel} \f$ and \f$ \vec{v}_{\perp} \f$
+   * denoting the projections of the initial velocity \f$ \vec{v}_0 \f$ parallel
+   * and perpendicular to the axis \f$ \vec{B} \f$, respectively, i.e.
+   * \f{align*}{
+        \vec{v}_{\parallel} &= \frac{\vec{v}_0 \cdot \vec{B}}{\vec{B}^2} \vec{B} \\
+        \vec{v}_{\perp} &= \vec{v}_0 - \vec{v}_{\parallel}
+     \f}
+   */
+
+  class Helix {
+
+  public:
+    Helix(Point const& pR0, FrequencyType pOmegaC, VelocityVector const& pvPar,
+          VelocityVector const& pvPerp)
+        : r0_(pR0)
+        , omegaC_(pOmegaC)
+        , vPar_(pvPar)
+        , vPerp_(pvPerp)
+        , uPerp_(vPerp_.cross(vPar_.normalized()))
+        , radius_(pvPar.getNorm() / abs(pOmegaC)) {}
+
+    LengthType getRadius() const;
+
+    Point getPosition(TimeType const t) const;
+
+    VelocityVector getVelocity(TimeType const t) const;
+
+    Point getPositionFromArclength(LengthType const l) const;
+
+    LengthType getArcLength(TimeType const t1, TimeType const t2) const;
+
+    TimeType getTimeFromArclength(LengthType const l) const;
+
+  private:
+    Point r0_;             ///! origin of helix, but this is in the center of the
+                           ///! "cylinder" on which the helix rotates
+    FrequencyType omegaC_; ///! speed of angular rotation
+    VelocityVector vPar_;  ///! speed along direction of "cylinder"
+    VelocityVector vPerp_, uPerp_;
+    LengthType radius_;
+  };
+
+} // namespace corsika
+
+#include <corsika/detail/framework/geometry/Helix.inl>

corsika/framework/geometry/IVolume.hpp

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/geometry/Point.hpp>
+
+namespace corsika {
+
+  class IVolume {
+
+  public:
+    //! returns true if the Point p is within the volume
+    virtual bool contains(Point const& p) const = 0;
+
+    virtual ~IVolume() = default;
+  };
+
+} // namespace corsika

corsika/framework/geometry/Intersections.hpp

+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/core/PhysicalUnits.hpp>
+
+#include <map> // for pair
+
+namespace corsika {
+
+  /**
+   *
+   * Container to store and return a list of intersections of a
+   * trajectory with a geometric volume objects in space.
+   *
+   **/
+
+  class Intersections {
+
+    Intersections(const Intersections&) = delete;
+    Intersections(Intersections&&) = delete;
+    Intersections& operator=(const Intersections&) = delete;
+
+  public:
+    Intersections()
+        : has_intersections_(false) {}
+
+    Intersections(TimeType&& t1, TimeType&& t2)
+        : has_intersections_(true)
+        , intersections_(std::make_pair(t1, t2)) {}
+
+    Intersections(TimeType&& t)
+        : has_intersections_(true)
+        , intersections_(std::make_pair(
+              t, std::numeric_limits<TimeType::value_type>::infinity() * second)) {}
+
+    bool hasIntersections() const { return has_intersections_; }
+
+    ///! where did the trajectory currently enter the volume
+    TimeType getEntry() const {
+      if (has_intersections_)
+        return intersections_.first;
+      else
+        return std::numeric_limits<TimeType::value_type>::infinity() * second;
+    }
+
+    ///! where did the trajectory currently exit the volume
+    TimeType getExit() const {
+      if (has_intersections_)
+        return intersections_.second;
+      else
+        return std::numeric_limits<TimeType::value_type>::infinity() * second;
+    }
+
+  private:
+    bool has_intersections_;
+    std::pair<TimeType, TimeType> intersections_;
+  };
+
+} // namespace corsika