corsika.dox

-/**
-@mainpage CORSIKA air shower simulations 
-
-Source code documentation and reference guide for the CORSIKA8
-software framework for air shower simulations.
-
-
-*/

corsika/detail/framework/analytics/ClassTimer.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.
+ */
+
+// 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>
+
+namespace corsika {
+
+  // Common
+
+  template <typename TClass, typename TRet, typename... TArgs,
+            TRet (TClass::*TFuncPtr)(TArgs...), typename TTimer>
+  inline ClassTimer<TRet (TClass::*)(TArgs...), TFuncPtr, TTimer>::ClassTimer(TClass& obj)
+      : ClassTimerImpl<TClass, TTimer>(obj) {}
+
+  template <typename TClass, typename TRet, typename... TArgs,
+            TRet (TClass::*TFuncPtr)(TArgs...), typename TTimer>
+  inline TRet ClassTimer<TRet (TClass::*)(TArgs...), TFuncPtr, TTimer>::call(
+      TArgs... args) {
+    this->start_ = ClassTimerImpl<TClass, TTimer>::clock_type::now();
+    auto tmp = (this->obj_.*TFuncPtr)(std::forward<TArgs>(args)...);
+    this->timeDiff_ = std::chrono::duration_cast<
+        typename ClassTimerImpl<TClass, TTimer>::duration_type>(
+        ClassTimerImpl<TClass, TTimer>::clock_type::now() - this->start_);
+    return tmp;
+  }
+
+  // Specialisation 1
+
+  template <typename TClass, typename... TArgs, void (TClass::*TFuncPtr)(TArgs...),
+            typename TTimer>
+  inline ClassTimer<void (TClass::*)(TArgs...), TFuncPtr, TTimer>::ClassTimer(TClass& obj)
+      : ClassTimerImpl<TClass, TTimer>(obj) {}
+
+  template <typename TClass, typename... TArgs, void (TClass::*TFuncPtr)(TArgs...),
+            typename TTimer>
+  inline void ClassTimer<void (TClass::*)(TArgs...), TFuncPtr, TTimer>::call(
+      TArgs... args) {
+    this->start_ = ClassTimerImpl<TClass, TTimer>::clock_type::now();
+    (this->obj_.*TFuncPtr)(std::forward<TArgs>(args)...);
+    this->timeDiff_ = std::chrono::duration_cast<
+        typename ClassTimerImpl<TClass, TTimer>::duration_type>(
+        ClassTimerImpl<TClass, TTimer>::clock_type::now() - this->start_);
+    return;
+  }
+
+  /// Specialisation 2
+
+  template <typename TClass, typename TRet, typename... TArgs,
+            TRet (TClass::*TFuncPtr)(TArgs...) const, typename TTimer>
+  inline ClassTimer<TRet (TClass::*)(TArgs...) const, TFuncPtr, TTimer>::ClassTimer(
+      TClass& obj)
+      : ClassTimerImpl<TClass, TTimer>(obj) {}
+
+  template <typename TClass, typename TRet, typename... TArgs,
+            TRet (TClass::*TFuncPtr)(TArgs...) const, typename TTimer>
+  inline TRet ClassTimer<TRet (TClass::*)(TArgs...) const, TFuncPtr, TTimer>::call(
+      TArgs... args) {
+    this->start_ = ClassTimerImpl<TClass, TTimer>::clock_type::now();
+    auto tmp = (this->obj_.*TFuncPtr)(std::forward<TArgs>(args)...);
+    this->timeDiff_ = std::chrono::duration_cast<
+        typename ClassTimerImpl<TClass, TTimer>::duration_type>(
+        ClassTimerImpl<TClass, TTimer>::clock_type::now() - this->start_);
+    return tmp;
+  }
+
+  /// Specialisation 3
+  template <typename TClass, typename... TArgs, void (TClass::*TFuncPtr)(TArgs...) const,
+            typename TTimer>
+  inline ClassTimer<void (TClass::*)(TArgs...) const, TFuncPtr, TTimer>::ClassTimer(
+      TClass& obj)
+      : ClassTimerImpl<TClass, TTimer>(obj) {}
+
+  template <typename TClass, typename... TArgs, void (TClass::*TFuncPtr)(TArgs...) const,
+            typename TTimer>
+  inline void ClassTimer<void (TClass::*)(TArgs...) const, TFuncPtr, TTimer>::call(
+      TArgs... args) {
+    this->start_ = ClassTimerImpl<TClass, TTimer>::clock_type::now();
+    (this->obj_.*TFuncPtr)(std::forward<TArgs>(args)...);
+    this->timeDiff_ = std::chrono::duration_cast<
+        typename ClassTimerImpl<TClass, TTimer>::duration_type>(
+        ClassTimerImpl<TClass, TTimer>::clock_type::now() - this->start_);
+    return;
+  }
+
+} // namespace corsika

corsika/detail/framework/analytics/FunctionTimer.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 <chrono>
+#include <utility>
+
+namespace corsika {
+
+  template <typename TFunc, typename TTime>
+  inline FunctionTimer<TFunc, TTime>::FunctionTimer(TFunc f)
+      : function_(f) {}
+
+  template <typename TFunc, typename TTime>
+  template <typename... TArgs>
+  inline auto FunctionTimer<TFunc, TTime>::operator()(TArgs&&... args)
+      -> std::invoke_result_t<TFunc, TArgs...> {
+    this->startTimer();
+    auto tmp = function_(std::forward<TArgs>(args)...);
+    this->stopTimer();
+    return tmp;
+  }
+
+} // namespace corsika

corsika/detail/framework/core/Cascade.inl

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/core/Step.hpp>
+
+#include <corsika/framework/process/ProcessReturn.hpp>
+#include <corsika/framework/process/ContinuousProcessStepLength.hpp>
+#include <corsika/framework/process/ContinuousProcessIndex.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/utility/COMBoost.hpp>
+
+#include <corsika/media/Environment.hpp>
+#include <corsika/media/NuclearComposition.hpp>
+
+#include <cassert>
+#include <cmath>
+#include <iostream>
+#include <limits>
+#include <type_traits>
+
+namespace corsika {
+
+  template <typename TTracking, typename TProcessList, typename TOutput, typename TStack>
+  inline Cascade<TTracking, TProcessList, TOutput, TStack>::Cascade(
+      Environment<medium_interface_type> const& env, TTracking& tr, TProcessList& pl,
+      TOutput& out, TStack& stack)
+      : environment_(env)
+      , tracking_(tr)
+      , sequence_(pl)
+      , output_(out)
+      , stack_(stack)
+      , forceInteraction_(false)
+      , forceDecay_(false) {
+    CORSIKA_LOG_INFO(c8_ascii_);
+    CORSIKA_LOG_INFO("This is CORSIKA {}.{}.{}.{}", CORSIKA_RELEASE_NUMBER,
+                     CORSIKA_MAJOR_NUMBER, CORSIKA_MINOR_NUMBER, CORSIKA_PATCH_NUMBER);
+    CORSIKA_LOG_INFO(
+        "The C8 author list can be found at: "
+        "https://gitlab.iap.kit.edu/AirShowerPhysics/corsika/-/wikis/"
+        "Current-CORSIKA-8-author-list");
+    CORSIKA_LOG_INFO("Tracking algorithm: {} (version {})", TTracking::getName(),
+                     TTracking::getVersion());
+    if constexpr (stack_view_type::has_event) {
+      CORSIKA_LOG_INFO("Stack - with full cascade HISTORY.");
+    }
+  }
+
+  template <typename TTracking, typename TProcessList, typename TOutput, typename TStack>
+  inline void Cascade<TTracking, TProcessList, TOutput, TStack>::run() {
+
+    // trigger the start of the outputs for this shower
+    output_.startOfShower();
+
+    setNodes(); // put each particle on stack in correct environment volume
+
+    while (!stack_.isEmpty()) {
+
+      sequence_.initCascadeEquations();
+
+      while (!stack_.isEmpty()) {
+        CORSIKA_LOG_TRACE("Stack: {}", stack_.asString());
+        count_++;
+        auto pNext = stack_.getNextParticle();
+
+        CORSIKA_LOG_TRACE(
+            "============== next particle : count={}, pid={}"
+            ", stack entries={}"
+            ", stack deleted={}",
+            count_, pNext.getPID(), stack_.getEntries(), stack_.getErased());
+
+        step(pNext);
+        sequence_.doStack(stack_);
+      }
+
+      // do cascade equations, which can put new particles on Stack,
+      // thus, the double loop
+      sequence_.doCascadeEquations(stack_);
+    }
+
+    // indicate end of shower
+    output_.endOfShower();
+  }
+
+  template <typename TTracking, typename TProcessList, typename TOutput, typename TStack>
+  inline void Cascade<TTracking, TProcessList, TOutput, TStack>::forceInteraction() {
+    forceInteraction_ = true;
+    if (forceDecay_) {
+      CORSIKA_LOG_ERROR("Cannot set forceInteraction when forceDecay is already set");
+      throw std::runtime_error(
+          "Cannot set forceInteraction when forceDecay is already set");
+    }
+  }
+
+  template <typename TTracking, typename TProcessList, typename TOutput, typename TStack>
+  inline void Cascade<TTracking, TProcessList, TOutput, TStack>::forceDecay() {
+    forceDecay_ = true;
+    if (forceInteraction_) {
+      CORSIKA_LOG_ERROR("Cannot set forceDecay when forceInteraction is already set");
+      throw std::runtime_error(
+          "Cannot set forceDecay when forceInteraction is already set");
+    }
+  }
+
+  template <typename TTracking, typename TProcessList, typename TOutput, typename TStack>
+  inline void Cascade<TTracking, TProcessList, TOutput, TStack>::step(
+      particle_type& particle) {
+
+    // determine the volume where the particle is (last) known to be
+    auto const* currentLogicalNode = particle.getNode();
+
+    // assert that particle stays outside void Universe if it has no
+    // model properties set
+    assert((currentLogicalNode != &*environment_.getUniverse() ||
+            environment_.getUniverse()->hasModelProperties()) &&
+           "FATAL: The environment model has no valid properties set!");
+
+    NuclearComposition const& composition =
+        currentLogicalNode->getModelProperties().getNuclearComposition();
+
+    // determine projectile
+    HEPEnergyType const Elab = particle.getEnergy();
+    FourMomentum const projectileP4{Elab, particle.getMomentum()};
+
+    // determine combined full inelastic cross section of the particles in the material
+    auto const targetMomentum = MomentumVector{
+        particle.getMomentum().getCoordinateSystem(), {0_GeV, 0_GeV, 0_GeV}};
+
+    auto const xs_function = [&](Code const targetId) -> CrossSectionType {
+      FourMomentum const targetP4{get_mass(targetId), targetMomentum};
+      return sequence_.getCrossSection(particle, targetId, targetP4);
+    };
+
+    CrossSectionType const total_cx_pre = composition.getWeightedSum(xs_function);
+
+    if (forceInteraction_) {
+      CORSIKA_LOG_TRACE("forced interaction!");
+      forceInteraction_ = false; // just one (first) interaction
+      stack_view_type secondaries(particle);
+      interaction(secondaries, projectileP4, composition, total_cx_pre);
+      sequence_.doSecondaries(secondaries);
+      particle.erase(); // primary particle is done
+      return;
+    }
+
+    if (forceDecay_) {
+      CORSIKA_LOG_TRACE("forced decay!");
+      forceDecay_ = false; // just one decay
+      stack_view_type secondaries(particle);
+      decay(secondaries, sequence_.getInverseLifetime(particle));
+      if (secondaries.getSize() == 1 && secondaries.getProjectile().getPID() ==
+                                            secondaries.getNextParticle().getPID()) {
+        throw std::runtime_error(
+            fmt::format("Particle {} decays into itself!",
+                        get_name(secondaries.getProjectile().getPID())));
+      }
+      sequence_.doSecondaries(secondaries);
+      particle.erase(); // primary particle is done
+      return;
+    }
+
+    // calculate interaction length in medium
+    GrammageType const total_lambda =
+        (composition.getAverageMassNumber() * constants::u) / total_cx_pre;
+
+    // sample random exponential step length in grammage
+    ExponentialDistribution expDist{total_lambda};
+    GrammageType const next_interact = expDist(rng_);
+
+    CORSIKA_LOG_DEBUG("total_lambda={} g/cm2, next_interact={} g/cm2",
+                      double(total_lambda / 1_g * 1_cm * 1_cm),
+                      double(next_interact / 1_g * 1_cm * 1_cm));
+
+    // determine combined total inverse decay time
+    InverseTimeType const total_inv_lifetime_pre = sequence_.getInverseLifetime(particle);
+
+    // sample random exponential decay time
+    ExponentialDistribution expDistDecay(1 / total_inv_lifetime_pre);
+    TimeType const next_decay = expDistDecay(rng_);
+
+    CORSIKA_LOG_DEBUG("total_lifetime={} ns, next_decay={} ns",
+                      (1 / total_inv_lifetime_pre) / 1_ns, next_decay / 1_ns);
+
+    // convert next_decay from time to length [m]
+    LengthType const distance_decay = next_decay * particle.getMomentum().getNorm() /
+                                      particle.getEnergy() * constants::c;
+
+    // determine geometric tracking
+    auto [track, nextVol] = tracking_.getTrack(particle);
+    auto geomMaxLength = track.getLength(1);
+
+    // convert next_step from grammage to length
+    LengthType const distance_interact =
+        currentLogicalNode->getModelProperties().getArclengthFromGrammage(track,
+                                                                          next_interact);
+
+    // determine the maximum geometric step length
+    ContinuousProcessStepLength const continuousMaxStep =
+        sequence_.getMaxStepLength(particle, track);
+    LengthType const continuous_max_dist = continuousMaxStep;
+
+    // take minimum of geometry, interaction, decay for next step
+    LengthType const min_discrete = std::min(distance_interact, distance_decay);
+    LengthType const min_non_continuous = std::min(min_discrete, geomMaxLength);
+    LengthType const min_distance = std::min(min_non_continuous, continuous_max_dist);
+
+    bool const isContinuous = continuous_max_dist < min_non_continuous;
+
+    // inform ContinuousProcesses (if applicable) that it is responsible for step-limit
+    // this would become simpler if we follow the idea of Max to enumerate ALL types of
+    // processes. Then non-continuous are included and no further logic is needed to
+    // distinguish between continuous and non-continuous limit.
+    auto const limitingId = isContinuous ? continuousMaxStep : ContinuousProcessIndex{};
+    // // the current step IS limited by a known continuous process
+
+    CORSIKA_LOG_DEBUG(
+        "transport particle by : {} m "
+        "Medium transition after: {} m "
+        "Decay after: {} m "
+        "Interaction after: {} m "
+        "Continuous limit: {} m ",
+        min_distance / 1_m, geomMaxLength / 1_m, distance_decay / 1_m,
+        distance_interact / 1_m, continuous_max_dist / 1_m);
+
+    // move particle along the trajectory to new position
+    // also update momentum/direction/time
+    track.setLength(min_distance);
+
+    Step step{particle, track};
+
+    // apply all continuous processes on particle + track
+    if (sequence_.doContinuous(step, limitingId) == ProcessReturn::ParticleAbsorbed) {
+      CORSIKA_LOG_DEBUG("Cascade: delete absorbed particle PID={} E={} GeV",
+                        particle.getPID(), particle.getEnergy() / 1_GeV);
+      if (particle.isErased()) {
+        CORSIKA_LOG_WARN(
+            "Particle marked as Absorbed in doContinuous, but prematurely erased. This "
+            "may be bug. Check.");
+      } else {
+        particle.erase();
+      }
+      return; // particle is gone -> return
+    }
+    particle.setTime(step.getTimePost());
+    particle.setPosition(step.getPositionPost());
+    particle.setDirection(step.getDirectionPost());
+    particle.setKineticEnergy(step.getEkinPost());
+
+    if (isContinuous) {
+      return; // there is nothing further, step is finished
+    }
+
+    CORSIKA_LOG_DEBUG("discrete process before geometric limit ? {}",
+                      ((min_distance < geomMaxLength) ? "yes" : "no"));
+
+    if (geomMaxLength < min_discrete) {
+      // geometric / tracking limit
+
+      if (nextVol != currentLogicalNode) {
+        // boundary crossing, step is limited by volume boundary
+
+        CORSIKA_LOG_DEBUG("volume boundary crossing to {}", fmt::ptr(nextVol));
+
+        if (nextVol == environment_.getUniverse().get()) {
+          CORSIKA_LOG_DEBUG(
+              "particle left physics world, is now in unknown space -> delete");
+          particle.erase();
+        }
+        particle.setNode(nextVol);
+        /*
+          doBoundary may delete the particle (or not)
+
+          caveat: any changes to particle, or even the production
+          of new secondaries is currently not passed to ParticleCut,
+          thus, particles outside the desired phase space may be produced.
+
+          \todo: this must be fixed.
+        */
+
+        sequence_.doBoundaryCrossing(particle, *currentLogicalNode, *nextVol);
+        return; // step finished
+      }
+
+      CORSIKA_LOG_DEBUG("step limit reached (e.g. deflection). nothing further happens.");
+
+      // final sanity check, no actions
+      {
+        auto const* numericalNodeAfterStep =
+            environment_.getUniverse()->getContainingNode(particle.getPosition());
+        CORSIKA_LOG_TRACE(
+            "Geometry check: numericalNodeAfterStep={} currentLogicalNode={}",
+            fmt::ptr(numericalNodeAfterStep), fmt::ptr(currentLogicalNode));
+        if (numericalNodeAfterStep != currentLogicalNode) {
+          CORSIKA_LOG_DEBUG(
+              "expect to be in node currentLogicalNode={} but are in "
+              "numericalNodeAfterStep={}. Continue, but without guarantee.",
+              fmt::ptr(currentLogicalNode), fmt::ptr(numericalNodeAfterStep));
+        }
+      }
+      // we did not cross any volume boundary
+
+      // step length limit
+      return;
+    }
+
+    // interaction or decay to happen in this step
+    // the outcome of decay or interaction MAY be a) new particles in
+    // secondaries, b) the projectile particle deleted (or
+    // changed)
+
+    stack_view_type secondaries{particle};
+
+    /*
+      Create SecondaryView object on Stack. The data container
+      remains untouched and identical, and 'projectile' is identical
+      to 'particle' above this line. However,
+      projectile.addSecondaries populate the SecondaryView, which can
+      then be used afterwards for further processing. Thus: it is
+      important to use projectile/view (and not particle) for Interaction,
+      and Decay!
+    */
+
+    FourMomentum const projectileP4Post{particle.getEnergy(), particle.getMomentum()};
+
+    bool eraseParticle =
+        false; // only erase original particle if it decayed or interacted
+
+    if (distance_interact < distance_decay) {
+      eraseParticle = isInteracted(
+          interaction(secondaries, projectileP4Post, composition, total_cx_pre));
+    } else {
+      [[maybe_unused]] auto projectile = secondaries.getProjectile();
+      if (decay(secondaries, total_inv_lifetime_pre) == ProcessReturn::Decayed) {
+        eraseParticle = true;
+        if (secondaries.getSize() == 1 &&
+            projectile.getPID() == secondaries.getNextParticle().getPID()) {
+          throw std::runtime_error(fmt::format("Particle {} decays into itself!",
+                                               get_name(projectile.getPID())));
+        }
+      }
+    }
+
+    if (eraseParticle) {
+      // doSecondaries() makes sense only if there was an actual event
+      sequence_.doSecondaries(secondaries);
+      particle.erase();
+    }
+  }
+
+  template <typename TTracking, typename TProcessList, typename TOutput, typename TStack>
+  inline ProcessReturn Cascade<TTracking, TProcessList, TOutput, TStack>::decay(
+      stack_view_type& view, InverseTimeType initial_inv_decay_time) {
+    CORSIKA_LOG_DEBUG("decay");
+
+    // one option is that decay_time is now larger (less
+    // probability for decay) than it was before the step, thus,
+    // no decay might actually occur and is allowed
+
+    UniformRealDistribution<InverseTimeType> uniDist(initial_inv_decay_time);
+    const auto sample_process = uniDist(rng_);
+
+    auto const returnCode = sequence_.selectDecay(view, sample_process);
+    if (returnCode != ProcessReturn::Decayed) {
+      CORSIKA_LOG_ERROR("Particle {} did not decay!",
+                        get_name(view.getProjectile().getPID()));
+    }
+    setEventType(view, history::EventType::Decay);
+    return returnCode;
+  }
+
+  template <typename TTracking, typename TProcessList, typename TOutput, typename TStack>
+  inline ProcessReturn Cascade<TTracking, TProcessList, TOutput, TStack>::interaction(
+      stack_view_type& view, FourMomentum const& projectileP4,
+      NuclearComposition const& composition,
+      CrossSectionType const initial_cross_section) {
+
+    CORSIKA_LOG_DEBUG("collide");
+
+    // one option is that cross section is now smaller (less
+    // probability for collision) than it was before the step, thus,
+    // no interaction might actually occur and is allowed
+
+    UniformRealDistribution<CrossSectionType> uniDist(initial_cross_section);
+    CrossSectionType const sample_process_by_cx = uniDist(rng_);
+    auto const returnCode = sequence_.selectInteraction(view, projectileP4, composition,
+                                                        rng_, sample_process_by_cx);
+    if (returnCode != ProcessReturn::Interacted) {
+      CORSIKA_LOG_DEBUG("Particle did not interact!");
+    }
+    setEventType(view, history::EventType::Interaction);
+    return returnCode;
+  }
+
+  template <typename TTracking, typename TProcessList, typename TOutput, typename TStack>
+  inline void Cascade<TTracking, TProcessList, TOutput, TStack>::setNodes() {
+    std::for_each(stack_.begin(), stack_.end(), [&](auto& p) {
+      auto const* numericalNode =
+          environment_.getUniverse()->getContainingNode(p.getPosition());
+      p.setNode(numericalNode);
+    });
+  }
+
+  template <typename TTracking, typename TProcessList, typename TOutput, typename TStack>
+  inline void Cascade<TTracking, TProcessList, TOutput, TStack>::setEventType(
+      stack_view_type& view, [[maybe_unused]] history::EventType eventType) {
+    if constexpr (stack_view_type::has_event) {
+      for (auto&& sec : view) { sec.getEvent()->setEventType(eventType); }
+    }
+  }
+
+} // namespace corsika

corsika/detail/framework/core/Logging.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 <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 {
+
+  // many of these free functions are special to the logging
+  // infrastructure so we hide them in the corsika::logging namespace.
+  namespace logging {
+
+    /*
+     * The default pattern for CORSIKA8 loggers.
+     */
+    inline auto set_default_level(level::level_enum const minlevel) -> void {
+      spdlog::set_level(minlevel);
+    }
+
+    template <typename TLogger>
+    inline auto add_source_info(TLogger& logger) -> void {
+      logger->set_pattern(source_pattern);
+    }
+
+    template <typename TLogger>
+    inline auto reset_pattern(TLogger& logger) -> void {
+      logger->set_pattern(default_pattern);
+    }
+
+  } // namespace logging
+
+  inline std::shared_ptr<spdlog::logger> create_logger(std::string const& name,
+                                                       bool const defaultlog) {
+
+    // create the logger
+    // this is currently a colorized multi-threading safe logger
+    auto logger = spdlog::stdout_color_mt(name);
+
+    // set the default C8 format
+#if (!defined(_GLIBCXX_USE_CXX11_ABI) || _GLIBCXX_USE_CXX11_ABI == 1)
+    logger->set_pattern(default_pattern);
+#else
+    // special case: gcc from the software collections devtoolset
+    std::string dp(default_pattern);
+    logger->set_pattern(dp);
+#endif
+
+    // if defaultlog is True, we set this as the default spdlog logger.
+    if (defaultlog) { spdlog::set_default_logger(logger); }
+
+    return logger;
+  }
+
+  inline std::shared_ptr<spdlog::logger> get_logger(std::string const& name,
+                                                    bool const defaultlog) {
+
+    // attempt to get the logger from the registry
+    auto logger = spdlog::get(name);
+
+    // weg found the logger, so just return it
+    if (logger) {
+      return logger;
+    } else { // logger was not found so create it
+      return create_logger(name, defaultlog);
+    }
+  }
+
+} // namespace corsika

corsika/detail/framework/core/ParticleProperties.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 <cmath>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/core/Logging.hpp>
+
+namespace corsika {
+
+  inline HEPEnergyType get_kinetic_energy_propagation_threshold(Code const code) {
+    if (is_nucleus(code)) return particle::detail::threshold_nuclei;
+    return particle::detail::propagation_thresholds[static_cast<CodeIntType>(code)];
+  }
+
+  inline void set_kinetic_energy_propagation_threshold(Code const code,
+                                                       HEPEnergyType const val) {
+    if (is_nucleus(code))
+      particle::detail::threshold_nuclei = val;
+    else
+      particle::detail::propagation_thresholds[static_cast<CodeIntType>(code)] = val;
+  }
+
+  inline HEPMassType constexpr get_mass(Code const code) {
+    if (is_nucleus(code)) { return get_nucleus_mass(code); }
+    return particle::detail::masses[static_cast<CodeIntType>(code)];
+  }
+
+  inline HEPEnergyType get_energy_production_threshold(Code const p) {
+    return particle::detail::production_thresholds[static_cast<CodeIntType>(p)];
+  }
+
+  inline void set_energy_production_threshold(Code const p, HEPEnergyType const val) {
+    particle::detail::production_thresholds[static_cast<CodeIntType>(p)] = val;
+  }
+
+  inline bool constexpr is_charged(Code const c) { return get_charge_number(c) != 0; }
+
+  inline bool constexpr is_nucleus(Code const code) { return code >= Code::Nucleus; }
+
+  inline Code constexpr get_nucleus_code(size_t const A,
+                                         size_t const Z) { // 10LZZZAAAI
+    if (Z > A) { throw std::runtime_error("Z cannot be larger than A in nucleus."); }
+    return static_cast<Code>(static_cast<CodeIntType>(Code::Nucleus) + Z * 10000 +
+                             A * 10);
+  }
+
+  inline size_t constexpr get_nucleus_Z(Code const code) {
+    return (static_cast<CodeIntType>(code) % static_cast<CodeIntType>(Code::Nucleus)) /
+           10000;
+  }
+
+  inline size_t constexpr get_nucleus_A(Code const code) {
+    return (static_cast<CodeIntType>(code) % 10000) / 10;
+  }
+
+  inline PDGCode constexpr get_PDG(Code const code) {
+    if (code < Code::Nucleus) {
+      return particle::detail::pdg_codes[static_cast<CodeIntType>(code)];
+    }
+    size_t const Z = get_nucleus_Z(code);
+    size_t const A = get_nucleus_A(code);
+    return static_cast<PDGCode>(static_cast<CodeIntType>(Code::Nucleus) + Z * 10000 +
+                                A * 10); // 10LZZZAAAI
+  }
+
+  inline PDGCode constexpr get_PDG(unsigned int const A, unsigned int const Z) {
+    return PDGCode(1000000000 + Z * 10000 + A * 10);
+  }
+
+  inline int16_t constexpr get_charge_number(Code const code) {
+    if (is_nucleus(code)) return get_nucleus_Z(code);
+    return particle::detail::electric_charges[static_cast<CodeIntType>(code)];
+  }
+
+  inline ElectricChargeType constexpr get_charge(Code const code) {
+    if (code == Code::Nucleus)
+      throw std::runtime_error("charge of particle::Nucleus undefined");
+    return get_charge_number(code) * constants::e;
+  }
+
+  inline std::string_view constexpr get_name(Code const code) {
+    if (is_nucleus(code)) { return "nucleus"; }
+    return particle::detail::names[static_cast<CodeIntType>(code)];
+  }
+
+  inline std::string get_name(Code code, full_name) {
+    if (is_nucleus(code)) {
+      return fmt::format("nucleus ({},{})", get_nucleus_A(code), get_nucleus_Z(code));
+    }
+
+    return std::string{get_name(code)};
+  }
+
+  inline TimeType constexpr get_lifetime(Code const p) {
+    return particle::detail::lifetime[static_cast<CodeIntType>(p)] * second;
+  }
+
+  inline bool constexpr is_hadron(Code const code) {
+    if (is_nucleus(code)) return true;
+    return particle::detail::isHadron[static_cast<CodeIntType>(code)];
+  }
+
+  inline bool constexpr is_em(Code const c) {
+    return c == Code::Electron || c == Code::Positron || c == Code::Photon;
+  }
+
+  inline bool constexpr is_muon(Code const c) {
+    return c == Code::MuPlus || c == Code::MuMinus;
+  }
+
+  inline bool constexpr is_neutrino(Code const c) {
+    return c == Code::NuE || c == Code::NuMu || c == Code::NuTau || c == Code::NuEBar ||
+           c == Code::NuMuBar || c == Code::NuTauBar;
+  }
+
+  inline std::ostream& operator<<(std::ostream& stream, corsika::Code const code) {
+    return stream << get_name(code);
+  }
+
+  inline Code convert_from_PDG(PDGCode const p) {
+    static_assert(particle::detail::conversionArray.size() % 2 == 1);
+    // this will fail, for the strange case where the maxPDG is negative...
+    int constexpr maxPDG{(particle::detail::conversionArray.size() - 1) >> 1};
+    auto const k = static_cast<PDGCodeIntType>(p);
+    if (std::abs(k) <= maxPDG) {
+      return particle::detail::conversionArray[k + maxPDG];
+    } else {
+      if (1000000000 <= k && k <= 1009999990) { // nucleus (no L or I)
+        int const Z = (k - 1000000000) / 10000;
+        int const A = (k - 1000000000 - 10000 * Z) / 10;
+        return get_nucleus_code(A, Z);
+      }
+      return particle::detail::conversionMap.at(p);
+    }
+  }
+
+  inline HEPMassType constexpr get_nucleus_mass(Code const code) {
+    unsigned int const A = get_nucleus_A(code);
+    unsigned int const Z = get_nucleus_Z(code);
+    return get_nucleus_mass(A, Z);
+  }
+
+  inline HEPMassType constexpr get_nucleus_mass(unsigned int const A,
+                                                unsigned int const Z) {
+    return get_mass(Code::Proton) * Z + (A - Z) * get_mass(Code::Neutron);
+  }
+
+  inline std::string get_nucleus_name(Code const code) {
+    size_t const A = get_nucleus_A(code);
+    size_t const Z = get_nucleus_Z(code);
+    return fmt::format("Nucleus_A{}_Z{}", A, Z);
+  }
+
+  inline std::initializer_list<Code> constexpr get_all_particles() {
+    return particle::detail::all_particles;
+  }
+
+} // namespace corsika

corsika/detail/framework/core/SpdlogSpecializations.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 <fmt/format.h>
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <boost/filesystem/path.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <Eigen/Dense>
+
+//-----------------------------------
+// STD
+//-----------------------------------
+namespace std {
+
+  auto inline format_as(std::_Put_time<char> const& arg) {
+    std::ostringstream os;
+    os << arg;
+    return os.str();
+  }
+
+} // namespace std
+
+//-----------------------------------
+// CORSIKA
+//-----------------------------------
+namespace corsika {
+
+  // formatters for particle codes declared on ParticleProperties.hpp
+  auto inline format_as(Code code) { return get_name(code); }
+  auto inline format_as(PDGCode code) { return fmt::underlying(code); }
+
+  template <typename Type>
+  auto inline format_as(Type const& arg) {
+    std::ostringstream os;
+    os << arg;
+    return os.str();
+  }
+
+} // namespace corsika
+
+//-----------------------------------
+// boost::filesystem
+//-----------------------------------
+namespace boost::filesystem {
+
+  auto inline format_as(path const& fname) { return fname.string().c_str(); }
+
+} // namespace boost::filesystem
+
+//-----------------------------------
+// phys::units
+//-----------------------------------
+namespace phys::units {
+
+  template <typename Dimensions>
+  auto inline format_as(phys::units::quantity<Dimensions> const& arg) {
+
+    return io::to_string(arg);
+  }
+
+} // namespace phys::units
+
+//----------------------------------
+// Eigen
+//----------------------------------
+namespace Eigen {
+
+  template <typename Scalar, int M, int N>
+  auto inline format_as(Matrix<Scalar, M, N> const& arg) {
+
+    std::ostringstream os;
+    os << arg;
+    return os.str();
+  }
+
+  template <typename Matrix1, typename Matrix2>
+  auto inline format_as(Product<Matrix1, Matrix2> const& arg) {
+
+    std::ostringstream os;
+    os << arg;
+    return os.str();
+  }
+
+} // namespace Eigen

corsika/detail/framework/geometry/BaseVector.inl

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/geometry/CoordinateSystem.hpp>
+#include <corsika/framework/geometry/QuantityVector.hpp>
+
+namespace corsika {
+
+  template <typename TDimension>
+  inline CoordinateSystemPtr BaseVector<TDimension>::getCoordinateSystem() const {
+    return cs_;
+  }
+
+  template <typename TDimension>
+  inline QuantityVector<TDimension> const& BaseVector<TDimension>::getQuantityVector()
+      const {
+    return quantityVector_;
+  }
+
+  template <typename TDimension>
+  inline QuantityVector<TDimension>& BaseVector<TDimension>::getQuantityVector() {
+    return quantityVector_;
+  }
+
+} // namespace corsika

corsika/detail/framework/geometry/Box.inl

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+namespace corsika {
+  inline Box::Box(CoordinateSystemPtr cs, LengthType const x, LengthType const y,
+                  LengthType const z)
+      : center_(Point(cs, {0_m, 0_m, 0_m}))
+      , cs_(cs)
+      , x_(x)
+      , y_(y)
+      , z_(z) {}
+
+  inline Box::Box(CoordinateSystemPtr cs, LengthType const side)
+      : center_(Point(cs, {0_m, 0_m, 0_m}))
+      , cs_(cs)
+      , x_(side / 2)
+      , y_(side / 2)
+      , z_(side / 2) {}
+
+  inline bool Box::contains(Point const& p) const {
+    if ((abs(p.getX(cs_)) < x_) && (abs(p.getY(cs_)) < y_) && (abs(p.getZ(cs_)) < z_))
+      return true;
+    else
+      return false;
+  }
+
+  inline std::string Box::asString() const {
+    std::ostringstream txt;
+    txt << "center=" << center_ << ", x-axis=" << DirectionVector{cs_, {1, 0, 0}}
+        << ", y-axis: " << DirectionVector{cs_, {0, 1, 0}}
+        << ", z-axis: " << DirectionVector{cs_, {0, 0, 1}};
+    return txt.str();
+  }
+
+  template <typename TDim>
+  inline void Box::rotate(QuantityVector<TDim> const& axis, double const angle) {
+    cs_ = make_rotation(cs_, axis, angle);
+  }
+
+} // namespace corsika
\ No newline at end of file

corsika/detail/framework/geometry/CoordinateSystem.inl

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/geometry/RootCoordinateSystem.hpp>
+#include <corsika/framework/geometry/QuantityVector.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+
+#include <Eigen/Dense>
+
+#include <memory>
+#include <stdexcept>
+
+namespace corsika {
+
+  inline CoordinateSystemPtr CoordinateSystem::getReferenceCS() const {
+    return referenceCS_;
+  }
+
+  inline EigenTransform const& CoordinateSystem::getTransform() const { return transf_; }
+
+  inline bool CoordinateSystem::operator==(CoordinateSystem const& cs) const {
+    return referenceCS_ == cs.referenceCS_ && transf_.matrix() == cs.transf_.matrix();
+  }
+
+  inline bool CoordinateSystem::operator!=(CoordinateSystem const& cs) const {
+    return !(cs == *this);
+  }
+
+  /// find transformation between two CS, using most optimal common base
+  inline EigenTransform get_transformation(CoordinateSystem const& pFrom,
+                                           CoordinateSystem const& pTo) {
+    CoordinateSystem const* a{&pFrom};
+    CoordinateSystem const* b{&pTo};
+
+    while (a != b && b) {
+
+      // traverse pFrom
+      a = &pFrom;
+      while (a != b && a) { a = a->getReferenceCS().get(); }
+
+      if (a == b) break;
+
+      b = b->getReferenceCS().get();
+    }
+
+    if (a != b || a == nullptr) {
+      throw std::runtime_error("no connection between coordinate systems found!");
+    }
+
+    CoordinateSystem const* commonBase = a;
+    CoordinateSystem const* p = &pFrom;
+    EigenTransform t = EigenTransform::Identity();
+    while ((*p) != (*commonBase)) {
+      t = p->getTransform() * t;
+      p = p->getReferenceCS().get();
+    }
+
+    p = &pTo;
+
+    while (*p != *commonBase) {
+      t = t * p->getTransform().inverse(Eigen::TransformTraits::Isometry);
+      p = p->getReferenceCS().get();
+    }
+
+    return t;
+  }
+
+  inline CoordinateSystemPtr make_translation(CoordinateSystemPtr const& cs,
+                                              QuantityVector<length_d> const& vector) {
+    EigenTransform const translation{EigenTranslation(vector.getEigenVector())};
+    return CoordinateSystemPtr{new CoordinateSystem(cs, translation)};
+  }
+
+  template <typename TDim>
+  inline CoordinateSystemPtr make_rotationToZ(CoordinateSystemPtr const& cs,
+                                              Vector<TDim> const& vVec) {
+    auto const vVecComp = vVec.getComponents(cs);
+    if (vVecComp.getX().magnitude() == 0 && vVecComp.getY().magnitude() == 0 &&
+        vVecComp.getZ().magnitude() == 0) {
+      return cs;
+    }
+
+    auto const a = vVecComp.normalized().getEigenVector();
+    auto const a1 = a(0), a2 = a(1), a3 = a(2);
+
+    Eigen::Matrix3d A, B;
+
+    if (a3 > 0) {
+      auto const c = 1 / (1 + a3);
+      A << 1, 0, a1,                      // comment to prevent clang-format
+          0, 1, a2,                       // .
+          -a1, -a2, 1;                    // .
+      B << -a1 * a1 * c, -a1 * a2 * c, 0, // .
+          -a1 * a2 * c, -a2 * a2 * c, 0,  // .
+          0, 0, -(a1 * a1 + a2 * a2) * c; // .
+
+    } else {
+      auto const c = 1 / (1 - a3);
+      A << 1, 0, a1,                      // .
+          0, -1, a2,                      // .
+          a1, -a2, -1;                    // .
+      B << -a1 * a1 * c, +a1 * a2 * c, 0, // .
+          -a1 * a2 * c, +a2 * a2 * c, 0,  // .
+          0, 0, (a1 * a1 + a2 * a2) * c;  // .
+    }
+
+    return CoordinateSystemPtr{new CoordinateSystem{cs, EigenTransform{A + B}}};
+  }
+
+  template <typename TDim>
+  inline CoordinateSystemPtr make_rotation(CoordinateSystemPtr const& cs,
+                                           QuantityVector<TDim> const& axis,
+                                           double const angle) {
+    if (axis.getEigenVector().isZero()) {
+      throw std::runtime_error("null-vector given as axis parameter");
+    }
+
+    EigenTransform const rotation{
+        Eigen::AngleAxisd(angle, axis.getEigenVector().normalized())};
+
+    return CoordinateSystemPtr{new CoordinateSystem{cs, rotation}};
+  }
+
+  template <typename TDim>
+  inline CoordinateSystemPtr make_translationAndRotation(
+      CoordinateSystemPtr const& cs, QuantityVector<length_d> const& translation,
+      QuantityVector<TDim> const& axis, double const angle) {
+    if (axis.getEigenVector().isZero()) {
+      throw std::runtime_error("null-vector given as axis parameter");
+    }
+
+    EigenTransform const transf{
+        Eigen::AngleAxisd(angle, axis.getEigenVector().normalized()) *
+        EigenTranslation(translation.getEigenVector())};
+
+    return CoordinateSystemPtr{new CoordinateSystem{cs, transf}};
+  }
+
+} // namespace corsika

corsika/detail/framework/geometry/FourVector.inl

+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <type_traits>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/geometry/Vector.hpp>
+
+namespace corsika {
+
+  template <typename TTimeType, typename TSpaceVecType>
+  inline TTimeType FourVector<TTimeType, TSpaceVecType>::getTimeLikeComponent() const {
+    return timeLike_;
+  }
+
+  template <typename TTimeType, typename TSpaceVecType>
+  inline TSpaceVecType& FourVector<TTimeType, TSpaceVecType>::getSpaceLikeComponents() {
+    return spaceLike_;
+  }
+
+  template <typename TTimeType, typename TSpaceVecType>
+  inline TSpaceVecType const&
+  FourVector<TTimeType, TSpaceVecType>::getSpaceLikeComponents() const {
+    return spaceLike_;
+  }
+
+  template <typename TTimeType, typename TSpaceVecType>
+  inline typename FourVector<TTimeType, TSpaceVecType>::norm_square_type
+  FourVector<TTimeType, TSpaceVecType>::getNormSqr() const {
+    return getTimeSquared() - spaceLike_.getSquaredNorm();
+  }
+
+  template <typename TTimeType, typename TSpaceVecType>
+  inline typename FourVector<TTimeType, TSpaceVecType>::norm_type
+  FourVector<TTimeType, TSpaceVecType>::getNorm() const {
+
+    return sqrt(abs(getNormSqr()));
+  }
+
+  template <typename TTimeType, typename TSpaceVecType>
+  inline bool FourVector<TTimeType, TSpaceVecType>::isTimelike() const {
+    return getTimeSquared() < spaceLike_.getSquaredNorm();
+  }
+
+  template <typename TTimeType, typename TSpaceVecType>
+  inline bool FourVector<TTimeType, TSpaceVecType>::isSpacelike() const {
+    return getTimeSquared() > spaceLike_.getSquaredNorm();
+  }
+
+  template <typename TTimeType, typename TSpaceVecType>
+  inline FourVector<TTimeType, TSpaceVecType>& FourVector<TTimeType, TSpaceVecType>::
+  operator+=(FourVector const& b) {
+    timeLike_ += b.timeLike_;
+    spaceLike_ += b.spaceLike_;
+
+    return *this;
+  }
+
+  template <typename TTimeType, typename TSpaceVecType>
+  inline FourVector<TTimeType, TSpaceVecType>& FourVector<TTimeType, TSpaceVecType>::
+  operator-=(FourVector const& b) {
+    timeLike_ -= b.timeLike_;
+    spaceLike_ -= b.spaceLike_;
+    return *this;
+  }
+
+  template <typename TTimeType, typename TSpaceVecType>
+  inline FourVector<TTimeType, TSpaceVecType>& FourVector<TTimeType, TSpaceVecType>::
+  operator*=(double const b) {
+    timeLike_ *= b;
+    spaceLike_ *= b;
+    return *this;
+  }
+
+  template <typename TTimeType, typename TSpaceVecType>
+  inline FourVector<TTimeType, TSpaceVecType>& FourVector<TTimeType, TSpaceVecType>::
+  operator/=(double const b) {
+    timeLike_ /= b;
+    spaceLike_.getComponents() /= b;
+    return *this;
+  }
+
+  template <typename TTimeType, typename TSpaceVecType>
+  inline FourVector<TTimeType, TSpaceVecType>& FourVector<TTimeType, TSpaceVecType>::
+  operator/(double const b) {
+    *this /= b;
+    return *this;
+  }
+
+  template <typename TTimeType, typename TSpaceVecType>
+  inline typename FourVector<TTimeType, TSpaceVecType>::norm_type
+      FourVector<TTimeType, TSpaceVecType>::operator*(FourVector const& b) {
+    if constexpr (std::is_same<time_type, decltype(std::declval<space_type>() / meter *
+                                                   second)>::value)
+      return timeLike_ * b.timeLike_ * constants::cSquared - spaceLike_.norm();
+    else
+      return timeLike_ * timeLike_ - spaceLike_.norm();
+  }
+
+  template <typename TTimeType, typename TSpaceVecType>
+  inline typename FourVector<TTimeType, TSpaceVecType>::norm_square_type
+  FourVector<TTimeType, TSpaceVecType>::getTimeSquared() const {
+    if constexpr (std::is_same<time_type, decltype(std::declval<space_type>() / meter *
+                                                   second)>::value)
+      return timeLike_ * timeLike_ * constants::cSquared;
+    else
+      return timeLike_ * timeLike_;
+  }
+
+  template <typename TTimeType, typename TSpaceVecType>
+  inline std::ostream& operator<<(
+      std::ostream& os, corsika::FourVector<TTimeType, TSpaceVecType> const qv) {
+
+    os << '(' << qv.timeLike_ << ", " << qv.spaceLike_ << ") ";
+    return os;
+  }
+
+} // namespace corsika

corsika/detail/framework/geometry/Helix.inl

+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/geometry/Vector.hpp>
+#include <cmath>
+
+namespace corsika {
+
+  inline LengthType Helix::getRadius() const { return radius_; }
+
+  inline Point Helix::getPosition(TimeType const t) const {
+    return r0_ + vPar_ * t +
+           (vPerp_ * (std::cos(omegaC_ * t) - 1) + uPerp_ * std::sin(omegaC_ * t)) /
+               omegaC_;
+  }
+
+  inline Point Helix::getPositionFromArclength(LengthType const l) const {
+    return getPosition(getTimeFromArclength(l));
+  }
+
+  inline LengthType Helix::getArcLength(TimeType const t1, TimeType const t2) const {
+    return (vPar_ + vPerp_).getNorm() * (t2 - t1);
+  }
+
+  inline TimeType Helix::getTimeFromArclength(LengthType const l) const {
+    return l / (vPar_ + vPerp_).getNorm();
+  }
+
+} // namespace corsika

corsika/detail/framework/geometry/LeapFrogTrajectory.inl

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/geometry/Line.hpp>
+#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/framework/geometry/PhysicalGeometry.hpp>
+#include <corsika/framework/utility/QuarticSolver.hpp>
+
+namespace corsika {
+
+  inline Line LeapFrogTrajectory::getLine() const {
+    auto D = getPosition(1) - getPosition(0);
+    auto d = D.getNorm();
+    auto v = initialVelocity_;
+    if (d > 1_um) { // if trajectory is ultra-short, we do not
+                    // re-calculate velocity, just use initial
+                    // value. Otherwise, this is numerically unstable
+      v = D / d * getVelocity(0).getNorm();
+    }
+    return Line(getPosition(0), v);
+  }
+
+  inline Point LeapFrogTrajectory::getPosition(double const u) const {
+    if (u == 0) return initialPosition_;
+    Point const position = initialPosition_ + initialVelocity_ * timeStep_ * u / 2;
+    VelocityVector const velocity =
+        initialVelocity_ + initialVelocity_.cross(magneticfield_) * timeStep_ * u * k_;
+    return position + velocity * timeStep_ * u / 2;
+  }
+
+  inline VelocityVector LeapFrogTrajectory::getVelocity(double const u) const {
+    return getDirection(u) * initialVelocity_.getNorm();
+  }
+
+  inline DirectionVector LeapFrogTrajectory::getDirection(double const u) const {
+    if (u == 0) return initialDirection_;
+
+    return (initialDirection_ +
+            initialDirection_.cross(magneticfield_) * timeStep_ * u * k_)
+        .normalized();
+  }
+
+  inline TimeType LeapFrogTrajectory::getDuration(double const u) const {
+    TimeType const step = timeStep_ * u;
+    double const correction = 1;
+    // the eventual (delta-L to L) correction factor is:
+    //    (initialDirection_ + initialDirection_.cross(magneticfield_) * step *
+    //    k_).getNorm();
+    return step / 2 * (correction + 1);
+  }
+
+  template <typename Particle>
+  inline TimeType LeapFrogTrajectory::getTime(Particle const& particle,
+                                              double const u) const {
+    return particle.getTime() + getDuration(u);
+  }
+
+  inline LengthType LeapFrogTrajectory::getLength(double const u) const {
+    return getDuration(u) * initialVelocity_.getNorm();
+  }
+
+  inline void LeapFrogTrajectory::setLength(LengthType const limit) {
+    if (initialVelocity_.getNorm() == SpeedType::zero()) setDuration(0_s);
+    setDuration(limit / initialVelocity_.getNorm());
+  }
+
+  inline void LeapFrogTrajectory::setDuration(TimeType const limit) {
+    /*
+    initial attempt to calculate delta-L from assumed full-leap-frog-length L:
+
+    Note: often return 0. Not good enough yet.
+
+    LengthType const L = initialVelocity_.getNorm() * limit; // distance
+    double const a = (initialVelocity_.cross(magneticfield_) * k_).getSquaredNorm() / 4 /
+                     square(1_m) * static_pow<4>(1_s);
+    double const d = L * initialVelocity_.getNorm() / square(1_m) * 1_s;
+    double const e = -square(L) / square(1_m);
+    std::vector<double> solutions = solve_quartic_real(a, 0, 0, d, e);
+    CORSIKA_LOG_DEBUG("setDuration limit={} L={} solution={}", limit, L,
+                      fmt::join(solutions, ", "));
+    */
+    timeStep_ = limit;
+  }
+
+} // namespace corsika

corsika/detail/framework/geometry/Line.inl

+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/geometry/Vector.hpp>
+
+namespace corsika {
+
+  inline Point Line::getPosition(TimeType const t) const {
+    return start_point_ + velocity_ * t;
+  }
+
+  inline VelocityVector const& Line::getVelocity(TimeType const) const {
+    return velocity_;
+  }
+
+  inline Point Line::getPositionFromArclength(LengthType const l) const {
+    return start_point_ + velocity_.normalized() * l;
+  }
+
+  inline LengthType Line::getArcLength(TimeType const t1, TimeType const t2) const {
+    return velocity_.getNorm() * (t2 - t1);
+  }
+
+  inline TimeType Line::getTimeFromArclength(LengthType const t) const {
+    return t / velocity_.getNorm();
+  }
+
+  inline Point const& Line::getStartPoint() const { return start_point_; }
+
+  inline DirectionVector Line::getDirection() const { return velocity_.normalized(); }
+
+  inline VelocityVector const& Line::getVelocity() const { return velocity_; }
+
+} // namespace corsika

corsika/detail/framework/geometry/Path.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 <deque>
+
+namespace corsika {
+
+  inline Path::Path(Point const& point) { points_.push_front(point); }
+
+  inline Path::Path(std::deque<Point> const& points)
+      : points_(points) {
+    int dequesize_ = points.size();
+    if (dequesize_ == 0 || dequesize_ == 1) {
+      length_ = LengthType::zero();
+    } else if (dequesize_ == 2) {
+      length_ = (points.back() - points.front()).getNorm();
+    } else {
+      for (auto point = points.begin(); point != points.end() - 1; ++point) {
+        auto point_next = *(point + 1);
+        auto point_now = *(point);
+        length_ += (point_next - point_now).getNorm();
+      }
+    }
+  }
+
+  inline void Path::addToEnd(Point const& point) {
+    length_ += (point - points_.back()).getNorm();
+    points_.push_back(point);
+  }
+
+  inline void Path::removeFromEnd() {
+    auto lastpoint_ = points_.back();
+    points_.pop_back();
+    int dequesize_ = points_.size();
+    if (dequesize_ == 0 || dequesize_ == 1) {
+      length_ = LengthType::zero();
+    } else if (dequesize_ == 2) {
+      length_ = (points_.back() - points_.front()).getNorm();
+    } else {
+      length_ -= (lastpoint_ - points_.back()).getNorm();
+    }
+  }
+
+  inline LengthType Path::getLength() const { return length_; }
+
+  inline Point const& Path::getStart() const { return points_.front(); }
+
+  inline Point const& Path::getEnd() const { return points_.back(); }
+
+  inline Point const& Path::getPoint(std::size_t const index) const {
+    return points_.at(index);
+  }
+
+  inline Path::const_iterator Path::begin() const { return points_.cbegin(); }
+
+  inline Path::const_iterator Path::end() const { return points_.cend(); }
+
+  inline Path::iterator Path::begin() { return points_.begin(); }
+
+  inline Path::iterator Path::end() { return points_.end(); }
+
+  inline int Path::getNSegments() const { return points_.size() - 1; }
+
+} // namespace corsika
\ No newline at end of file

corsika/detail/framework/geometry/Plane.inl

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/framework/geometry/Vector.hpp>
+
+#include <sstream>
+
+namespace corsika {
+
+  inline bool Plane::isAbove(Point const& vP) const {
+    return normal_.dot(vP - center_) > LengthType::zero();
+  }
+
+  inline LengthType Plane::getDistanceTo(Point const& vP) const {
+    return (normal_ * (vP - center_).dot(normal_)).getNorm();
+  }
+
+  inline Point const& Plane::getCenter() const { return center_; }
+
+  inline DirectionVector const& Plane::getNormal() const { return normal_; }
+
+  inline std::string Plane::asString() const {
+    std::ostringstream txt;
+    txt << "center=" << center_ << ", normal=" << normal_;
+    return txt.str();
+  }
+
+} // namespace corsika

corsika/detail/framework/geometry/Point.inl

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/geometry/BaseVector.hpp>
+#include <corsika/framework/geometry/QuantityVector.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/geometry/Vector.hpp>
+
+namespace corsika {
+
+  inline QuantityVector<length_d> const& Point::getCoordinates() const {
+    return BaseVector<length_d>::getQuantityVector();
+  }
+
+  inline QuantityVector<length_d>& Point::getCoordinates() {
+    return BaseVector<length_d>::getQuantityVector();
+  }
+
+  inline LengthType Point::getX(CoordinateSystemPtr const& pCS) const {
+    return getCoordinates(pCS).getX();
+  }
+
+  inline LengthType Point::getY(CoordinateSystemPtr const& pCS) const {
+    return getCoordinates(pCS).getY();
+  }
+
+  inline LengthType Point::getZ(CoordinateSystemPtr const& pCS) const {
+    return getCoordinates(pCS).getZ();
+  }
+
+  /// this always returns a QuantityVector as triple
+  inline QuantityVector<length_d> Point::getCoordinates(
+      CoordinateSystemPtr const& pCS) const {
+    CoordinateSystemPtr const& cs = BaseVector<length_d>::getCoordinateSystem();
+    if (pCS == cs) {
+      return BaseVector<length_d>::getQuantityVector();
+    } else {
+      return QuantityVector<length_d>(
+          get_transformation(*cs.get(), *pCS.get()) *
+          BaseVector<length_d>::getQuantityVector().eigenVector_);
+    }
+  }
+
+  /// this always returns a QuantityVector as triple
+  inline QuantityVector<length_d>& Point::getCoordinates(CoordinateSystemPtr const& pCS) {
+    if (pCS != BaseVector<length_d>::getCoordinateSystem()) { rebase(pCS); }
+    return BaseVector<length_d>::getQuantityVector();
+  }
+
+  inline void Point::rebase(CoordinateSystemPtr const& pCS) {
+    BaseVector<length_d>::setQuantityVector(QuantityVector<length_d>(
+        get_transformation(*BaseVector<length_d>::getCoordinateSystem().get(),
+                           *pCS.get()) *
+        BaseVector<length_d>::getQuantityVector().eigenVector_));
+    BaseVector<length_d>::setCoordinateSystem(pCS);
+  }
+
+  inline Point Point::operator+(Vector<length_d> const& pVec) const {
+    CoordinateSystemPtr const& cs = BaseVector<length_d>::getCoordinateSystem();
+    return Point(cs, getCoordinates() + pVec.getComponents(cs));
+  }
+
+  inline Point Point::operator-(Vector<length_d> const& pVec) const {
+    CoordinateSystemPtr const& cs = BaseVector<length_d>::getCoordinateSystem();
+    return Point(cs, getCoordinates() - pVec.getComponents(cs));
+  }
+
+  inline Vector<length_d> Point::operator-(Point const& pB) const {
+    CoordinateSystemPtr const& cs = BaseVector<length_d>::getCoordinateSystem();
+    return Vector<length_d>(cs, getCoordinates() - pB.getCoordinates(cs));
+  }
+  inline std::ostream& operator<<(std::ostream& os, corsika::Point const& p) {
+    auto const& qv = p.getCoordinates();
+    os << qv << " (ref:" << fmt::ptr(p.getCoordinateSystem()) << ")";
+    return os;
+  }
+
+  inline LengthType distance(Point const& p1, Point const& p2) {
+    return (p1 - p2).getNorm();
+  }
+
+} // namespace corsika

corsika/detail/framework/geometry/QuantityVector.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 <Eigen/Dense>
+
+#include <iostream>
+#include <utility>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+
+namespace corsika {
+
+  template <typename TDimension>
+  inline typename QuantityVector<TDimension>::quantity_type QuantityVector<TDimension>::
+  operator[](size_t const index) const {
+    return quantity_type(phys::units::detail::magnitude_tag, eigenVector_[index]);
+  }
+
+  template <typename TDimension>
+  inline typename QuantityVector<TDimension>::quantity_type
+  QuantityVector<TDimension>::getX() const {
+    return (*this)[0];
+  }
+
+  template <typename TDimension>
+  inline typename QuantityVector<TDimension>::quantity_type
+  QuantityVector<TDimension>::getY() const {
+    return (*this)[1];
+  }
+
+  template <typename TDimension>
+  inline typename QuantityVector<TDimension>::quantity_type
+  QuantityVector<TDimension>::getZ() const {
+    return (*this)[2];
+  }
+
+  template <typename TDimension>
+  inline typename QuantityVector<TDimension>::quantity_type
+  QuantityVector<TDimension>::getNorm() const {
+    return quantity_type(phys::units::detail::magnitude_tag, eigenVector_.norm());
+  }
+
+  template <typename TDimension>
+  inline typename QuantityVector<TDimension>::quantity_square_type
+  QuantityVector<TDimension>::getSquaredNorm() const {
+    using QuantitySquared =
+        decltype(std::declval<quantity_type>() * std::declval<quantity_type>());
+    return QuantitySquared(phys::units::detail::magnitude_tag,
+                           eigenVector_.squaredNorm());
+  }
+
+  template <typename TDimension>
+  inline QuantityVector<TDimension> QuantityVector<TDimension>::operator+(
+      QuantityVector<TDimension> const& pQVec) const {
+    return QuantityVector<TDimension>(eigenVector_ + pQVec.eigenVector_);
+  }
+
+  template <typename TDimension>
+  inline QuantityVector<TDimension> QuantityVector<TDimension>::operator-(
+      QuantityVector<TDimension> const& pQVec) const {
+    return QuantityVector<TDimension>(eigenVector_ - pQVec.eigenVector_);
+  }
+
+  template <typename TDimension>
+  template <typename TScalarDim>
+  inline auto QuantityVector<TDimension>::operator*(
+      phys::units::quantity<TScalarDim, double> const p) const {
+    using ResQuantity =
+        phys::units::detail::Product<TScalarDim, TDimension, double, double>;
+
+    if constexpr (std::is_same<ResQuantity, double>::value) // result dimensionless, not
+                                                            // a "quantity_type" anymore
+    {
+      return QuantityVector<phys::units::dimensionless_d>(eigenVector_ * p.magnitude());
+    } else {
+      return QuantityVector<typename ResQuantity::dimension_type>(eigenVector_ *
+                                                                  p.magnitude());
+    }
+  }
+
+  template <typename TDimension>
+  template <typename TScalarDim>
+  inline auto QuantityVector<TDimension>::operator/(
+      phys::units::quantity<TScalarDim, double> const p) const {
+    return (*this) * (1 / p);
+  }
+
+  template <typename TDimension>
+  inline auto QuantityVector<TDimension>::operator*(double const p) const {
+    return QuantityVector<TDimension>(eigenVector_ * p);
+  }
+
+  template <typename TDimension>
+  inline auto QuantityVector<TDimension>::operator/(double const p) const {
+    return QuantityVector<TDimension>(eigenVector_ / p);
+  }
+
+  template <typename TDimension>
+  inline auto& QuantityVector<TDimension>::operator/=(double const p) {
+    eigenVector_ /= p;
+    return *this;
+  }
+
+  template <typename TDimension>
+  inline auto& QuantityVector<TDimension>::operator*=(double const p) {
+    eigenVector_ *= p;
+    return *this;
+  }
+
+  template <typename TDimension>
+  inline auto& QuantityVector<TDimension>::operator+=(
+      QuantityVector<TDimension> const& pQVec) {
+    eigenVector_ += pQVec.eigenVector_;
+    return *this;
+  }
+
+  template <typename TDimension>
+  inline auto& QuantityVector<TDimension>::operator-=(
+      QuantityVector<TDimension> const& pQVec) {
+    eigenVector_ -= pQVec.eigenVector_;
+    return *this;
+  }
+
+  template <typename TDimension>
+  inline auto QuantityVector<TDimension>::operator-() const {
+    return QuantityVector<TDimension>(-eigenVector_);
+  }
+
+  template <typename TDimension>
+  inline auto QuantityVector<TDimension>::normalized() const {
+    return QuantityVector<TDimension>(eigenVector_.normalized());
+  }
+
+  template <typename TDimension>
+  inline auto QuantityVector<TDimension>::operator==(
+      QuantityVector<TDimension> const& p) const {
+    return eigenVector_ == p.eigenVector_;
+  }
+
+  template <typename TDimension>
+  inline std::ostream& operator<<(std::ostream& os,
+                                  corsika::QuantityVector<TDimension> const& qv) {
+    using quantity_type = phys::units::quantity<TDimension, double>;
+
+    os << '(' << qv.eigenVector_(0) << ' ' << qv.eigenVector_(1) << ' '
+       << qv.eigenVector_(2) << ") "
+       << phys::units::to_unit_symbol<TDimension, double>(
+              quantity_type(phys::units::detail::magnitude_tag, 1));
+    return os;
+  }
+
+} // namespace corsika

corsika/detail/framework/geometry/SeparationPlane.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.
+ */
+
+#pragma once
+
+namespace corsika {
+
+  inline SeparationPlane::SeparationPlane(Plane const& plane)
+      : plane_(plane) {}
+
+  inline bool SeparationPlane::contains(Point const& p) const {
+    return !plane_.isAbove(p);
+  }
+
+  inline std::string SeparationPlane::asString() const { return plane_.asString(); }
+
+} // namespace corsika
\ No newline at end of file

corsika/detail/framework/geometry/Sphere.inl

+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+
+namespace corsika {
+
+  inline bool Sphere::contains(Point const& p) const {
+    return radius_ * radius_ > (center_ - p).getSquaredNorm();
+  }
+
+  inline Point const& Sphere::getCenter() const { return center_; }
+
+  inline void Sphere::setCenter(Point const& p) { center_ = p; }
+
+  inline LengthType Sphere::getRadius() const { return radius_; }
+
+  inline void Sphere::setRadius(LengthType const r) { radius_ = r; }
+
+  inline CoordinateSystemPtr const Sphere::getCoordinateSystem() const {
+    return center_.getCoordinateSystem();
+  }
+
+  inline std::string Sphere::asString() const {
+    std::ostringstream txt;
+    txt << "center=" << center_ << ", radius=" << radius_;
+    return txt.str();
+  }
+
+} // namespace corsika