corsika/detail/setup/SetupStack.inl

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

corsika/detail/stack/NuclearStackExtension.inl

-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#pragma once
-
-#include <corsika/framework/core/ParticleProperties.hpp>
-#include <corsika/framework/core/PhysicalUnits.hpp>
-#include <corsika/framework/stack/Stack.hpp>
-#include <corsika/framework/geometry/Point.hpp>
-#include <corsika/framework/geometry/Vector.hpp>
-#include <corsika/stack/VectorStack.hpp>
-
-#include <algorithm>
-#include <tuple>
-#include <vector>
-
-namespace corsika::nuclear_stack {
-
-  template <template <typename> class InnerParticleInterface,
-            typename StackIteratorInterface>
-  inline void NuclearParticleInterface<InnerParticleInterface, StackIteratorInterface>::
-      setParticleData(particle_data_type const& v) {
-
-    if (std::get<0>(v) == Code::Nucleus) {
-      std::ostringstream err;
-      err << "NuclearStackExtension: no A and Z specified for new Nucleus!";
-      throw std::runtime_error(err.str());
-    }
-
-    super_type::setParticleData(v);
-    setNucleusRef(-1); // this is not a nucleus
-  }
-
-  template <template <typename> class InnerParticleInterface,
-            typename StackIteratorInterface>
-  inline void NuclearParticleInterface<InnerParticleInterface, StackIteratorInterface>::
-      setParticleData(altenative_particle_data_type const& v) {
-    const unsigned short A = std::get<5>(v);
-    const unsigned short Z = std::get<6>(v);
-    if (std::get<0>(v) != Code::Nucleus || A == 0 || Z == 0) {
-      std::ostringstream err;
-      err << "NuclearStackExtension: no A and Z specified for new Nucleus!";
-      throw std::runtime_error(err.str());
-    }
-    setNucleusRef(
-        super_type::getStackData().getNucleusNextRef()); // store this nucleus data ref
-    setNuclearA(A);
-    setNuclearZ(Z);
-    super_type::setParticleData(particle_data_type{
-        std::get<0>(v), std::get<1>(v), std::get<2>(v), std::get<3>(v), std::get<4>(v)});
-  }
-
-  template <template <typename> class InnerParticleInterface,
-            typename StackIteratorInterface>
-  inline void NuclearParticleInterface<InnerParticleInterface, StackIteratorInterface>::
-      setParticleData(super_type& p, particle_data_type const& v) {
-    if (std::get<0>(v) == Code::Nucleus) {
-      std::ostringstream err;
-      err << "NuclearStackExtension: no A and Z specified for new Nucleus!";
-      throw std::runtime_error(err.str());
-    }
-
-    super_type::setParticleData(
-        p, particle_data_type{std::get<0>(v), std::get<1>(v), std::get<2>(v),
-                              std::get<3>(v), std::get<4>(v)});
-
-    setNucleusRef(-1); // this is not a nucleus
-  }
-
-  template <template <typename> class InnerParticleInterface,
-            typename StackIteratorInterface>
-  inline void NuclearParticleInterface<InnerParticleInterface, StackIteratorInterface>::
-      setParticleData(super_type& p, altenative_particle_data_type const& v) {
-
-    const unsigned short A = std::get<5>(v);
-    const unsigned short Z = std::get<6>(v);
-
-    if (std::get<0>(v) != Code::Nucleus || A == 0 || Z == 0) {
-      std::ostringstream err;
-      err << "NuclearStackExtension: no A and Z specified for new Nucleus!";
-      throw std::runtime_error(err.str());
-    }
-
-    setNucleusRef(
-        super_type::getStackData().getNucleusNextRef()); // store this nucleus data ref
-    setNuclearA(A);
-    setNuclearZ(Z);
-    super_type::setParticleData(
-        p, particle_data_type{std::get<0>(v), std::get<1>(v), std::get<2>(v),
-                              std::get<3>(v), std::get<4>(v)});
-  }
-
-  template <template <typename> class InnerParticleInterface,
-            typename StackIteratorInterface>
-  inline std::string NuclearParticleInterface<InnerParticleInterface,
-                                              StackIteratorInterface>::asString() const {
-    return fmt::format(
-        "{}, nuc({})", super_type::asString(),
-        (isNucleus() ? fmt::format("A={}, Z={}", getNuclearA(), getNuclearZ()) : "n/a"));
-  }
-
-  template <template <typename> class InnerParticleInterface,
-            typename StackIteratorInterface>
-  inline HEPMassType NuclearParticleInterface<InnerParticleInterface,
-                                              StackIteratorInterface>::getMass() const {
-    if (super_type::getPID() == Code::Nucleus)
-      return get_nucleus_mass(getNuclearA(), getNuclearZ());
-    return super_type::getMass();
-  }
-
-  template <template <typename> class InnerParticleInterface,
-            typename StackIteratorInterface>
-  inline int16_t NuclearParticleInterface<
-      InnerParticleInterface, StackIteratorInterface>::getChargeNumber() const {
-    if (super_type::getPID() == Code::Nucleus) return getNuclearZ();
-    return super_type::getChargeNumber();
-  }
-
-  template <typename InnerStackImpl>
-  inline int NuclearStackExtensionImpl<InnerStackImpl>::getNucleusNextRef() {
-    nuclearA_.push_back(0);
-    nuclearZ_.push_back(0);
-    return nuclearA_.size() - 1;
-  }
-
-  template <typename InnerStackImpl>
-  inline int NuclearStackExtensionImpl<InnerStackImpl>::getNucleusRef(
-      const unsigned int i) const {
-    if (nucleusRef_[i] >= 0) return nucleusRef_[i];
-    std::ostringstream err;
-    err << "NuclearStackExtension: no nucleus at ref=" << i;
-    throw std::runtime_error(err.str());
-  }
-
-  template <typename InnerStackImpl>
-  inline void NuclearStackExtensionImpl<InnerStackImpl>::copy(const unsigned int i1,
-                                                              const unsigned int i2) {
-    // index range check
-    if (i1 >= getSize() || i2 >= getSize()) {
-      std::ostringstream err;
-      err << "NuclearStackExtension: trying to access data beyond size of stack!";
-      throw std::runtime_error(err.str());
-    }
-    // copy internal particle data p[i2] = p[i1]
-    super_type::copy(i1, i2);
-    // check if any of p[i1] or p[i2] was a Code::Nucleus
-    const int ref1 = nucleusRef_[i1];
-    const int ref2 = nucleusRef_[i2];
-    if (ref2 < 0) {
-      if (ref1 >= 0) {
-        // i1 is nucleus, i2 is not
-        nucleusRef_[i2] = getNucleusNextRef();
-        nuclearA_[nucleusRef_[i2]] = nuclearA_[ref1];
-        nuclearZ_[nucleusRef_[i2]] = nuclearZ_[ref1];
-      } else {
-        // neither i1 nor i2 are nuclei
-      }
-    } else {
-      if (ref1 >= 0) {
-        // both are nuclei, i2 is overwritten with nucleus i1
-        // fNucleusRef stays the same, but A and Z data is overwritten
-        nuclearA_[ref2] = nuclearA_[ref1];
-        nuclearZ_[ref2] = nuclearZ_[ref1];
-      } else {
-        // i2 is overwritten with non-nucleus i1
-        nucleusRef_[i2] = -1;                       // flag as non-nucleus
-        nuclearA_.erase(nuclearA_.cbegin() + ref2); // remove data for i2
-        nuclearZ_.erase(nuclearZ_.cbegin() + ref2); // remove data for i2
-        const int n = nucleusRef_.size(); // update fNucleusRef: indices above ref2
-        // must be decremented by 1
-        for (int i = 0; i < n; ++i) {
-          if (nucleusRef_[i] > ref2) { nucleusRef_[i] -= 1; }
-        }
-      }
-    }
-  }
-
-  template <typename InnerStackImpl>
-  inline void NuclearStackExtensionImpl<InnerStackImpl>::clear() {
-    super_type::clear();
-    nucleusRef_.clear();
-    nuclearA_.clear();
-    nuclearZ_.clear();
-  }
-
-  template <typename InnerStackImpl>
-  inline void NuclearStackExtensionImpl<InnerStackImpl>::swap(const unsigned int i1,
-                                                              const unsigned int i2) {
-    // index range check
-    if (i1 >= getSize() || i2 >= getSize()) {
-      std::ostringstream err;
-      err << "NuclearStackExtension: trying to access data beyond size of stack!";
-      throw std::runtime_error(err.str());
-    }
-    // swap original particle data
-    super_type::swap(i1, i2);
-    // swap corresponding nuclear reference data
-    std::swap(nucleusRef_[i2], nucleusRef_[i1]);
-  }
-
-  template <typename InnerStackImpl>
-  inline void NuclearStackExtensionImpl<InnerStackImpl>::incrementSize() {
-    super_type::incrementSize();
-    nucleusRef_.push_back(-1);
-  }
-
-  template <typename InnerStackImpl>
-  inline void NuclearStackExtensionImpl<InnerStackImpl>::decrementSize() {
-    super_type::decrementSize();
-    if (nucleusRef_.size() > 0) {
-      const int ref = nucleusRef_.back();
-      nucleusRef_.pop_back();
-      if (ref >= 0) {
-        nuclearA_.erase(nuclearA_.begin() + ref);
-        nuclearZ_.erase(nuclearZ_.begin() + ref);
-        const int n = nucleusRef_.size();
-        for (int i = 0; i < n; ++i) {
-          if (nucleusRef_[i] >= ref) { nucleusRef_[i] -= 1; }
-        }
-      }
-    }
-  }
-
-} // namespace corsika::nuclear_stack

corsika/detail/stack/VectorStack.inl

 /*
  * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
  *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
  */
 
 #pragma once
@@ -24,67 +23,95 @@ namespace corsika {
 
   template <typename StackIteratorInterface>
   inline void ParticleInterface<StackIteratorInterface>::setParticleData(
-      std::tuple<Code, HEPEnergyType, MomentumVector, Point, TimeType> const& v) {
+      particle_data_type const& v) {
     this->setPID(std::get<0>(v));
-    this->setEnergy(std::get<1>(v));
-    this->setMomentum(std::get<2>(v));
+    this->setKineticEnergy(std::get<1>(v));
+    this->setDirection(std::get<2>(v));
     this->setPosition(std::get<3>(v));
     this->setTime(std::get<4>(v));
   }
 
   template <typename StackIteratorInterface>
   inline void ParticleInterface<StackIteratorInterface>::setParticleData(
-      ParticleInterface<StackIteratorInterface> const&,
-      std::tuple<Code, HEPEnergyType, MomentumVector, Point, TimeType> const& v) {
+      ParticleInterface<StackIteratorInterface> const& parent,
+      secondary_data_type const& v) {
     this->setPID(std::get<0>(v));
-    this->setEnergy(std::get<1>(v));
-    this->setMomentum(std::get<2>(v));
-    this->setPosition(std::get<3>(v));
-    this->setTime(std::get<4>(v));
+    this->setKineticEnergy(std::get<1>(v));
+    this->setDirection(std::get<2>(v));
+    this->setPosition(parent.getPosition()); // position
+    this->setTime(parent.getTime());         // parent time
+  }
+
+  template <typename StackIteratorInterface>
+  inline void ParticleInterface<StackIteratorInterface>::setParticleData(
+      ParticleInterface<StackIteratorInterface> const& parent,
+      secondary_extended_data_type const& v) {
+    this->setPID(std::get<0>(v));
+    this->setKineticEnergy(std::get<1>(v));
+    this->setDirection(std::get<2>(v));
+    this->setPosition(parent.getPosition() + std::get<3>(v)); // + position
+    this->setTime(parent.getTime() + std::get<4>(v));         // + parent time
+  }
+
+  template <typename StackIteratorInterface>
+  inline std::string ParticleInterface<StackIteratorInterface>::asString() const {
+    return fmt::format("particle: i={}, PID={}, Ekin={}GeV", super_type::getIndex(),
+                       get_name(this->getPID()), this->getKineticEnergy() / 1_GeV);
   }
 
   inline void VectorStackImpl::clear() {
     dataPID_.clear();
-    dataE_.clear();
-    momentum_.clear();
+    dataEkin_.clear();
+    direction_.clear();
     position_.clear();
     time_.clear();
   }
 
-  inline void VectorStackImpl::copy(size_t i1, size_t i2) {
+  inline void VectorStackImpl::copy(size_t const i1, size_t const i2) {
+    // index range check
+    if (i1 >= getSize() || i2 >= getSize()) {
+      std::ostringstream err;
+      err << "VectorStackImpl: trying to access data beyond size of stack !";
+      throw std::runtime_error(err.str());
+    }
     dataPID_[i2] = dataPID_[i1];
-    dataE_[i2] = dataE_[i1];
-    momentum_[i2] = momentum_[i1];
+    dataEkin_[i2] = dataEkin_[i1];
+    direction_[i2] = direction_[i1];
     position_[i2] = position_[i1];
     time_[i2] = time_[i1];
   }
 
-  inline void VectorStackImpl::swap(size_t i1, size_t i2) {
+  inline void VectorStackImpl::swap(size_t const i1, size_t const i2) {
+    // index range check
+    if (i1 >= getSize() || i2 >= getSize()) {
+      std::ostringstream err;
+      err << "VectorStackImpl: trying to access data beyond size of stack !";
+      throw std::runtime_error(err.str());
+    }
     std::swap(dataPID_[i2], dataPID_[i1]);
-    std::swap(dataE_[i2], dataE_[i1]);
-    std::swap(momentum_[i2], momentum_[i1]);
+    std::swap(dataEkin_[i2], dataEkin_[i1]);
+    std::swap(direction_[i2], direction_[i1]);
     std::swap(position_[i2], position_[i1]);
     std::swap(time_[i2], time_[i1]);
   }
 
   inline void VectorStackImpl::incrementSize() {
     dataPID_.push_back(Code::Unknown);
-    dataE_.push_back(0 * electronvolt);
+    dataEkin_.push_back(0 * electronvolt);
 
     CoordinateSystemPtr const& dummyCS = get_root_CoordinateSystem();
 
-    momentum_.push_back(
-        MomentumVector(dummyCS, {0 * electronvolt, 0 * electronvolt, 0 * electronvolt}));
+    direction_.push_back(DirectionVector(dummyCS, {0, 0, 0}));
 
     position_.push_back(Point(dummyCS, {0 * meter, 0 * meter, 0 * meter}));
     time_.push_back(0 * second);
   }
 
   inline void VectorStackImpl::decrementSize() {
-    if (dataE_.size() > 0) {
+    if (dataEkin_.size() > 0) {
       dataPID_.pop_back();
-      dataE_.pop_back();
-      momentum_.pop_back();
+      dataEkin_.pop_back();
+      direction_.pop_back();
       position_.pop_back();
       time_.pop_back();
     }

corsika/detail/stack/WeightStackExtension.inl

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

corsika/detail/stack/history/HistoryObservationPlane.hpp

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

corsika/detail/stack/history/HistoryObservationPlane.inl

 /*
  * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
  *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
  */
 
 #include <corsika/logging/Logging.h>
@@ -13,15 +12,18 @@
 
 namespace corsika::history {
 
-  inline HistoryObservationPlane::HistoryObservationPlane(setup::Stack const& stack,
-                                                          Plane const& obsPlane,
-                                                          bool deleteOnHit)
+  template <typename TStack>
+  inline HistoryObservationPlane<TStack> HistoryObservationPlane(TStack const& stack,
+                                                                 Plane const& obsPlane,
+                                                                 bool deleteOnHit)
       : stack_{stack}
       , plane_{obsPlane}
       , deleteOnHit_{deleteOnHit} {}
 
-  inline ProcessReturn HistoryObservationPlane::DoContinuous(
-      setup::Stack::ParticleType const& particle, setup::Trajectory const& trajectory) {
+  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());
@@ -45,8 +47,10 @@ namespace corsika::history {
     }
   }
 
-  inline LengthType HistoryObservationPlane::MaxStepLength(
-      setup::Stack::ParticleType const&, setup::Trajectory const& trajectory) {
+  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());
@@ -59,8 +63,10 @@ namespace corsika::history {
     return (trajectory.getR0() - pointOfIntersection).norm() * 1.0001;
   }
 
-  inline void HistoryObservationPlane::fillHistoryHistogram(
-      setup::Stack::ParticleType const& muon) {
+  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};

corsika/framework/analytics/ClassTimer.hpp

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

corsika/framework/analytics/FunctionTimer.hpp

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

corsika/framework/analytics/Timer.hpp

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

corsika/framework/core/Cascade.hpp

 /*
  * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
  *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
  */
 
 #pragma once
 
-#include <corsika/framework/process/ProcessReturn.hpp>
+#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/framework/core/Logging.hpp>
 
-/*  see Issue 161, we need to include SetupStack only because we need
-    to globally define StackView. This is clearly not nice and should
-    be changed, when possible. It might be that StackView needs to be
-    templated in Cascade, but this would be even worse... so we don't
-    do that until it is really needed.
- */
-#include <corsika/setup/SetupStack.hpp>
+#include <corsika/stack/history/HistoryStackExtension.hpp>
 
 #include <cassert>
 #include <cmath>
 #include <limits>
 #include <type_traits>
 
-/**
- * The cascade namespace assembles all objects needed to simulate full particles cascades.
- */
-
 namespace corsika {
 
   /**
@@ -46,35 +38,31 @@ namespace corsika {
    * <b>TTracking</b> must be a class according to the
    * TrackingInterface providing the functions:
    *
-   * <code>
-   * auto getTrack(Particle const& p)</auto>,
+   * @code
+   * auto getTrack(particle_type const& p)</auto>,
    * with the return type <code>geometry::Trajectory<Line>
-   * </code>
+   * @endcode
    *
    * <b>TProcessList</b> must be a ProcessSequence.   *
    * <b>Stack</b> is the storage object for particle data, i.e. with
-   * Particle class type <code>Stack::ParticleType</code>
-   *
-   *
+   * particle class type `Stack::particle_type`.
    */
-  template <typename TTracking, typename TProcessList, typename TOutput, typename TStack,
-            /*
-             TStackView is needed as explicit template parameter because
-             of issue 161 and the
-             inability of clang to understand "stack::MakeView" so far.
-            */
-            typename TStackView = corsika::setup::StackView>
+  template <typename TTracking, typename TProcessList, typename TOutput, typename TStack>
   class Cascade {
 
-    typedef typename TStack::particle_type Particle;
-    typedef std::remove_pointer_t<decltype(((Particle*)nullptr)->getNode())>
-        VolumeTreeNode;
-    typedef typename VolumeTreeNode::IModelProperties MediumInterface;
+    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:
     /**
-     * \group constructors
-     * \{
+     * @name constructors
+     * @{
      * Cascade class cannot be default constructed, but needs a valid
      * list of physics processes for configuration at construct time.
      */
@@ -83,25 +71,13 @@ namespace corsika {
     Cascade(Cascade&&) = default;
     ~Cascade() = default;
     Cascade& operator=(Cascade const&) = default;
-    Cascade(Environment<MediumInterface> const& env, TTracking& tr, TProcessList& pl,
-            TOutput& out, TStack& stack)
-        : environment_(env)
-        , tracking_(tr)
-        , sequence_(pl)
-        , output_(out)
-        , stack_(stack) {
-      CORSIKA_LOG_INFO(c8_ascii_);
-      CORSIKA_LOG_INFO("Tracking algorithm: {} (version {})", TTracking::getName(),
-                       TTracking::getVersion());
-      if constexpr (TStackView::has_event) {
-        CORSIKA_LOG_INFO("Stack - with full cascade HISTORY.");
-      }
-    }
-    //! \}
+    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
+     * position.
      */
     void setNodes();
 
@@ -115,9 +91,18 @@ namespace corsika {
      * 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
@@ -129,19 +114,23 @@ namespace corsika {
      * New particles produced in one step are subject to further
      * processing, e.g. thinning, etc.
      */
-    void step(Particle& vParticle);
+    void step(particle_type& vParticle);
 
-    ProcessReturn decay(TStackView& view);
-    ProcessReturn interaction(TStackView& view);
-    void setEventType(TStackView& view, history::EventType);
+    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<MediumInterface> const& environment_;
+    Environment<medium_interface_type> const& environment_;
     TTracking& tracking_;
     TProcessList& sequence_;
     TOutput& output_;
     TStack& stack_;
-    default_prng_type& rng_ = RNGManager::getInstance().getRandomStream("cascade");
+    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:

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 GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
  */
 
 /**
@@ -22,23 +21,31 @@
 
 // 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 DEBUG
+#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>
@@ -49,7 +56,8 @@ namespace corsika {
   /*
    * The default pattern for CORSIKA8 loggers.
    */
-  const std::string default_pattern{"[%n:%^%-8l%$] %v"};
+  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"};
 
   /**
@@ -91,7 +99,6 @@ namespace corsika {
    */
   static inline std::shared_ptr<spdlog::logger> corsika_logger =
       get_logger("corsika", true);
-  // corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
 
   // many of these free functions are special to the logging
   // infrastructure so we hide them in the corsika::logging namespace.
@@ -103,7 +110,7 @@ namespace corsika {
     /**
      * Set the default log level for all *newly* created loggers.
      *
-     *  @param name    The minimum log level required to print.
+     *  @param minlevel The minimum log level required to print.
      *
      */
     auto set_default_level(level::level_enum const minlevel) -> void;
@@ -149,3 +156,4 @@ namespace corsika {
 } // 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 GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
  */
 
 /**
-   @file ParticleProperties.hpp
-
-   Interface to particle properties
+ *   @file ParticleProperties.hpp
+ *
+ * Interface to particle properties.
  */
 
 #pragma once
@@ -27,99 +26,214 @@
 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
-     particles, and each individual particles 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
+   * @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
+   * @{
+   */
 
-     The names, relations and properties of all particles known to CORSIKA 8 are listed
-     below.
+  /**
+   * @enum Code
+   *
+   * The Code enum is the actual place to define CORSIKA 8 particle codes.
+   */
+  enum class Code : std::int32_t;
 
-     @addtogroup Particles
-     @{
+  /**
+   * @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;
 
-  /** The Code enum is the actual place to define CORSIKA 8 particle codes. */
-  enum class Code : int16_t;
+  /**
+   * Internal integer type for enum PDGCode.
+   */
+  typedef std::underlying_type<PDGCode>::type PDGCodeIntType;
+} // namespace corsika
 
-  /** Specifically for PDG ids */
-  enum class PDGCode : int32_t;
+// data arrays, etc., as generated automatically
+#include <corsika/framework/core/GeneratedParticleProperties.inc>
 
-  using CodeIntType = std::underlying_type<Code>::type;
-  using PDGCodeType = std::underlying_type<PDGCode>::type;
+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
-  HEPEnergyType constexpr get_energy_threshold(
-      Code const); //!< get energy threshold below which the particle is discarded, by
-                   //!< default set to particle mass
-  void constexpr set_energy_threshold(
-      Code const, HEPEnergyType const); //!< set energy threshold below which the particle
-                                        //!< is discarded
-
-  inline void set_energy_threshold(std::pair<Code const, HEPEnergyType const> p) {
-    set_energy_threshold(p.first, p.second);
-  }
-  inline void set_energy_thresholds(
-      std::unordered_map<Code const, HEPEnergyType const> const& eCuts) {
-    for (auto v : eCuts) set_energy_threshold(v);
-  }
+
+  /**
+   * 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
-  TimeType constexpr get_lifetime(Code const);     //!< lifetime
+  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);
 
-  //! true iff the particle is a hard-coded nucleus or Code::Nucleus
+  /**
+   * Checks if Code corresponds to a nucleus.
+   *
+   * @return true if nucleus.
+   * @return false  if not nucleus.
+   */
   bool constexpr is_nucleus(Code const);
-  bool constexpr is_hadron(Code const); //!< true iff particle is hadron
-  bool constexpr is_em(Code const); //!< true iff particle is electron, positron or gamma
-  bool constexpr is_muon(Code const);     //!< true iff particle is mu+ or mu-
-  bool constexpr is_neutrino(Code const); //!< true iff particle is (anti-) neutrino
-  int constexpr get_nucleus_A(
+
+  /**
+   * 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
-  int constexpr get_nucleus_Z(
+
+  /**
+   * 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
 
-  //! returns mass of (A,Z) nucleus, disregarding binding energy
-  HEPMassType get_nucleus_mass(unsigned int const, unsigned int const);
+  /**
+   * @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);
 
-  //! convert PDG code to CORSIKA 8 internal 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();
 
-  //! the output stream operator for human-readable particle codes
+  /**
+   * @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
 
-// data arrays, etc., as generated automatically
-#include <corsika/framework/core/GeneratedParticleProperties.inc>
-
 #include <corsika/detail/framework/core/ParticleProperties.inl>
 
 // constants in namespaces-like static classes, generated automatically

corsika/framework/core/PhysicalConstants.hpp

@@ -2,9 +2,8 @@
  * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
  *
  *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
  */
 
 #pragma once
@@ -35,6 +34,10 @@ namespace corsika::constants {
   // 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};
 
@@ -65,7 +68,8 @@ namespace corsika::constants {
    */
   namespace EarthRadius {
     static constexpr auto Mean{6'371'000 * meter};
-    static constexpr auto Eqautorial{6'378'137 * 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

corsika/framework/core/PhysicalGeometry.hpp

 /*
  * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
  *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
  */
 
 #pragma once
@@ -23,18 +22,34 @@ namespace corsika {
   /**
    * A 3D vector defined in a specific coordinate system with units HEPMomentumType
    **/
-  typedef Vector<hepmomentum_d> MomentumVector;
+  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".
    **/
-  typedef Vector<dimensionless_d> DirectionVector;
+  using DirectionVector = Vector<dimensionless_d>;
 
   /**
    * A 3D vector defined in a specific coordinate system with units "velocity_t".
    *
    **/
-  typedef Vector<SpeedType::dimension_type> VelocityVector;
+  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 GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
  */
 
 #pragma once
@@ -92,10 +91,16 @@ namespace corsika::units::si {
       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() {

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 GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
  */
 
 #pragma once

corsika/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>