CMakeLists.txt

@@ -62,6 +62,7 @@ find_package (Eigen3 REQUIRED)
 add_subdirectory (ThirdParty)
 #add_subdirectory (Utilities)
 add_subdirectory (Framework)
+add_subdirectory (Environment)
 add_subdirectory (Stack)
 add_subdirectory (Setup)
 add_subdirectory (Processes)

Documentation/Examples/cascade_example.cc

-
 /**
  * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
  *
@@ -17,22 +16,33 @@
 #include <corsika/setup/SetupStack.h>
 #include <corsika/setup/SetupTrajectory.h>
 
+#include <corsika/environment/Environment.h>
+#include <corsika/environment/HomogeneousMedium.h>
+#include <corsika/environment/NuclearComposition.h>
+
 #include <corsika/random/RNGManager.h>
 
 #include <corsika/cascade/SibStack.h>
 #include <corsika/cascade/sibyll2.3c.h>
+#include <corsika/geometry/Sphere.h>
 #include <corsika/process/sibyll/ParticleConversion.h>
 
 #include <corsika/units/PhysicalUnits.h>
+
+#include <iostream>
+#include <limits>
+#include <typeinfo>
+
 using namespace corsika;
 using namespace corsika::process;
 using namespace corsika::units;
 using namespace corsika::particles;
 using namespace corsika::random;
+using namespace corsika::geometry;
+using namespace corsika::environment;
 
-#include <iostream>
-#include <typeinfo>
 using namespace std;
+using namespace corsika::units::si;
 
 static int fCount = 0;
 
@@ -42,7 +52,6 @@ public:
 
   template <typename Particle, typename Track>
   double MinStepLength(Particle& p, Track&) const {
-
     // beam particles for sibyll : 1, 2, 3 for p, pi, k
     // read from cross section code table
     int kBeam = 1;
@@ -94,7 +103,6 @@ public:
 
   template <typename Particle, typename Track, typename Stack>
   EProcessReturn DoContinuous(Particle&, Track&, Stack&) const {
-    // corsika::utls::ignore(p);
     return EProcessReturn::eOk;
   }
 
@@ -186,9 +194,8 @@ public:
     //    rnd_ini_();
 
     corsika::random::RNGManager& rmng = corsika::random::RNGManager::GetInstance();
-    ;
-    const std::string str_name = "s_rndm";
-    rmng.RegisterRandomStream(str_name);
+
+    rmng.RegisterRandomStream("s_rndm");
 
     // //    corsika::random::RNG srng;
     // auto srng = rmng.GetRandomStream("s_rndm");
@@ -236,9 +243,28 @@ double s_rndm_(int&) {
   return rmng() / (double)rmng.max();
 }
 
+class A {};
+class B : public A {};
+
 int main() {
+  corsika::environment::Environment env; // dummy environment
+  auto& universe = *(env.GetUniverse());
+
+  auto theMedium = corsika::environment::Environment::CreateNode<Sphere>(
+      Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m},
+      1_km * std::numeric_limits<double>::infinity());
+
+  using MyHomogeneousModel =
+      corsika::environment::HomogeneousMedium<corsika::environment::IMediumModel>;
+  theMedium->SetModelProperties<MyHomogeneousModel>(
+      1_g / (1_m * 1_m * 1_m),
+      corsika::environment::NuclearComposition(
+          std::vector<corsika::particles::Code>{corsika::particles::Code::Proton},
+          std::vector<float>{1.}));
+
+  universe.AddChild(std::move(theMedium));
 
-  tracking_line::TrackingLine<setup::Stack> tracking;
+  tracking_line::TrackingLine<setup::Stack> tracking(env);
   stack_inspector::StackInspector<setup::Stack> p0(true);
   ProcessSplit p1;
   const auto sequence = p0 + p1;
@@ -248,7 +274,7 @@ int main() {
 
   stack.Clear();
   auto particle = stack.NewParticle();
-  EnergyType E0 = 100_GeV;
+  EnergyType E0 = 100_TeV;
   particle.SetEnergy(E0);
   particle.SetPID(Code::Proton);
   EAS.Init();

Environment/CMakeLists.txt

+set (
+  ENVIRONMENT_HEADERS
+  VolumeTreeNode.h
+  IMediumModel.h
+  NuclearComposition.h
+  HomogeneousMedium.h
+  Environment.h
+  )
+
+set (
+  ENVIRONMENT_NAMESPACE
+  corsika/environment
+  )
+
+add_library (CORSIKAenvironment INTERFACE)
+CORSIKA_COPY_HEADERS_TO_NAMESPACE (CORSIKAenvironment ${ENVIRONMENT_NAMESPACE} ${ENVIRONMENT_HEADERS})
+
+# target dependencies on other libraries (also the header onlys)
+target_link_libraries (
+  CORSIKAenvironment
+  INTERFACE
+  CORSIKAgeometry
+  CORSIKAparticles
+  CORSIKAunits
+  )
+
+target_include_directories (
+  CORSIKAenvironment
+  INTERFACE
+  $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
+  $<INSTALL_INTERFACE:include>
+  )
+
+install (
+  TARGETS CORSIKAenvironment
+  LIBRARY DESTINATION lib
+  ARCHIVE DESTINATION lib
+  PUBLIC_HEADER DESTINATION include/${ENVIRONMENT_NAMESPACE}
+  )
+
+# --------------------
+# code unit testing
+add_executable (testEnvironment testEnvironment.cc)
+
+target_link_libraries (
+  testEnvironment
+  CORSIKAenvironment
+  CORSIKAthirdparty # for catch2
+  )
+
+add_test (NAME testGeometry COMMAND testGeometry)
+

Environment/Environment.h

+/**
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * 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.
+ */
+
+#ifndef _include_Environment_h
+#define _include_Environment_h
+
+#include <corsika/environment/IMediumModel.h>
+#include <corsika/environment/VolumeTreeNode.h>
+#include <corsika/geometry/RootCoordinateSystem.h>
+#include <corsika/setup/SetupEnvironment.h>
+
+namespace corsika::environment {
+  struct Universe : public corsika::geometry::Volume {
+    bool Contains(corsika::geometry::Point const&) const override { return true; }
+  };
+
+  // template <typename IEnvironmentModel>
+  class Environment {
+  public:
+    Environment()
+        : fUniverse(std::make_unique<VolumeTreeNode<IEnvironmentModel>>(
+              std::make_unique<Universe>())) {}
+
+    using IEnvironmentModel = corsika::setup::IEnvironmentModel;
+
+    auto& GetUniverse() { return fUniverse; }
+    auto const& GetUniverse() const { return fUniverse; }
+
+    auto const& GetCoordinateSystem() const {
+      return corsika::geometry::RootCoordinateSystem::GetInstance().GetRootCS();
+    }
+
+    // factory method for creation of VolumeTreeNodes
+    template <class VolumeType, typename... Args>
+    static auto CreateNode(Args&&... args) {
+      static_assert(std::is_base_of_v<corsika::geometry::Volume, VolumeType>,
+                    "unusable type provided, needs to be derived from "
+                    "\"corsika::geometry::Volume\"");
+
+      return std::make_unique<VolumeTreeNode<IEnvironmentModel>>(
+          std::make_unique<VolumeType>(std::forward<Args>(args)...));
+    }
+
+  private:
+    VolumeTreeNode<IEnvironmentModel>::VTNUPtr fUniverse;
+  };
+
+} // namespace corsika::environment
+
+#endif

Environment/HomogeneousMedium.h

+/**
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * 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.
+ */
+
+#ifndef _include_HomogeneousMedium_h_
+#define _include_HomogeneousMedium_h_
+
+#include <corsika/environment/NuclearComposition.h>
+#include <corsika/geometry/Line.h>
+#include <corsika/geometry/Point.h>
+#include <corsika/geometry/Trajectory.h>
+#include <corsika/particles/ParticleProperties.h>
+#include <corsika/units/PhysicalUnits.h>
+
+/**
+ * a homogeneous medium
+ */
+
+namespace corsika::environment {
+
+  template <class T>
+  class HomogeneousMedium : public T {
+    corsika::units::si::MassDensityType const fDensity;
+    NuclearComposition const fNuclComp;
+
+  public:
+    HomogeneousMedium(corsika::units::si::MassDensityType pDensity,
+                      NuclearComposition pNuclComp)
+        : fDensity(pDensity)
+        , fNuclComp(pNuclComp){};
+
+    corsika::units::si::MassDensityType GetMassDensity(
+        corsika::geometry::Point const&) const override {
+      return fDensity;
+    }
+    NuclearComposition const& GetNuclearComposition() const override { return fNuclComp; }
+
+    corsika::units::si::GrammageType IntegratedGrammage(
+        corsika::geometry::Trajectory<corsika::geometry::Line> const& pTraj,
+        corsika::units::si::TimeType pTo) const override {
+      using namespace corsika::units::si;
+      return pTraj.ArcLength(0_s, pTo) * fDensity;
+    }
+
+    corsika::units::si::TimeType FromGrammage(
+        corsika::geometry::Trajectory<corsika::geometry::Line> const& pTraj,
+        corsika::units::si::GrammageType pGrammage) const override {
+      return pTraj.TimeFromArclength(pGrammage / fDensity);
+    }
+  };
+
+} // namespace corsika::environment
+#endif

Environment/IMediumModel.h

+#ifndef _include_IMediumModel_h
+#define _include_IMediumModel_h
+
+#include <corsika/environment/NuclearComposition.h>
+#include <corsika/geometry/Line.h>
+#include <corsika/geometry/Point.h>
+#include <corsika/geometry/Trajectory.h>
+#include <corsika/units/PhysicalUnits.h>
+
+namespace corsika::environment {
+
+  class IMediumModel {
+  public:
+    virtual ~IMediumModel() = default;
+
+    virtual corsika::units::si::MassDensityType GetMassDensity(
+        corsika::geometry::Point const&) const = 0;
+
+    // todo: think about the mixin inheritance of the trajectory vs the BaseTrajectory
+    // approach for now, only lines are supported
+    virtual corsika::units::si::GrammageType IntegratedGrammage(
+        corsika::geometry::Trajectory<corsika::geometry::Line> const&,
+        corsika::units::si::TimeType) const = 0;
+
+    virtual corsika::units::si::TimeType FromGrammage(
+        corsika::geometry::Trajectory<corsika::geometry::Line> const&,
+        corsika::units::si::GrammageType) const = 0;
+
+    virtual NuclearComposition const& GetNuclearComposition() const = 0;
+  };
+
+} // namespace corsika::environment
+
+#endif

Environment/NuclearComposition.h

+#ifndef _include_NuclearComposition_h
+#define _include_NuclearComposition_h
+
+#include <corsika/particles/ParticleProperties.h>
+#include <numeric>
+#include <stdexcept>
+#include <vector>
+
+namespace corsika::environment {
+  class NuclearComposition {
+    std::vector<float> const fNumberFractions; //!< relative fractions of number density
+    std::vector<corsika::particles::Code> const
+        fComponents; //!< particle codes of consitutents
+
+  public:
+    NuclearComposition(std::vector<corsika::particles::Code> pComponents,
+                       std::vector<float> pFractions)
+        : fNumberFractions(pFractions)
+        , fComponents(pComponents) {
+      auto const sumFractions =
+          std::accumulate(pFractions.cbegin(), pFractions.cend(), 0.f);
+
+      if (!(0.999f < sumFractions && sumFractions < 1.001f)) {
+        throw std::runtime_error("element fractions do not add up to 1");
+      }
+    }
+
+    auto size() const { return fNumberFractions.size(); }
+
+    auto const& GetFractions() const { return fNumberFractions; }
+    auto const& GetComponents() const { return fComponents; }
+  };
+
+} // namespace corsika::environment
+
+#endif

Environment/VolumeTreeNode.h

+/**
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * 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.
+ */
+
+#ifndef _include_VolumeTreeNode_H
+#define _include_VolumeTreeNode_H
+
+#include <corsika/geometry/Volume.h>
+#include <memory>
+#include <vector>
+
+namespace corsika::environment {
+
+  class Empty {}; //<! intended for usage as default template argument
+
+  template <typename IModelProperties = Empty>
+  class VolumeTreeNode {
+  public:
+    using VTNUPtr = std::unique_ptr<VolumeTreeNode<IModelProperties>>;
+    using IMPSharedPtr = std::shared_ptr<IModelProperties>;
+    using VolUPtr = std::unique_ptr<corsika::geometry::Volume>;
+
+    VolumeTreeNode(VolUPtr pVolume = nullptr)
+        : fGeoVolume(std::move(pVolume)) {}
+
+    bool Contains(corsika::geometry::Point const& p) const {
+      return fGeoVolume->Contains(p);
+    }
+
+    VolumeTreeNode<IModelProperties> const* Excludes(
+        corsika::geometry::Point const& p) const {
+      auto exclContainsIter =
+          std::find_if(fExcludedNodes.cbegin(), fExcludedNodes.cend(),
+                       [&](auto const& s) { return bool(s->Contains(p)); });
+
+      return exclContainsIter != fExcludedNodes.cend() ? *exclContainsIter : nullptr;
+    }
+
+    /** returns a pointer to the sub-VolumeTreeNode which is "responsible" for the given
+     * \class Point \arg p, or nullptr iff \arg p is not contained in this volume.
+     */
+    VolumeTreeNode<IModelProperties> const* GetContainingNode(
+        corsika::geometry::Point const& p) const {
+      if (!Contains(p)) { return nullptr; }
+
+      if (auto const childContainsIter =
+              std::find_if(fChildNodes.cbegin(), fChildNodes.cend(),
+                           [&](auto const& s) { return bool(s->Contains(p)); });
+          childContainsIter == fChildNodes.cend()) // not contained in any of the children
+      {
+        if (auto const exclContainsIter = Excludes(p)) // contained in any excluded nodes
+        {
+          return exclContainsIter->GetContainingNode(p);
+        } else {
+          return this;
+        }
+      } else {
+        return (*childContainsIter)->GetContainingNode(p);
+      }
+    }
+
+    void AddChild(VTNUPtr pChild) {
+      pChild->fParentNode = this;
+      fChildNodes.push_back(std::move(pChild));
+      // It is a bad idea to return an iterator to the inserted element
+      // because it might get invalidated when the vector needs to grow
+      // later and the caller won't notice.
+    }
+
+    void ExcludeOverlapWith(VTNUPtr const& pNode) {
+      fExcludedNodes.push_back(pNode.get());
+    }
+
+    auto* GetParent() const { return fParentNode; };
+
+    auto const& GetChildNodes() const { return fChildNodes; }
+
+    auto const& GetExcludedNodes() const { return fExcludedNodes; }
+
+    auto const& GetVolume() const { return *fGeoVolume; }
+
+    auto const& GetModelProperties() const { return *fModelProperties; }
+
+    template <typename TModelProperties, typename... Args>
+    auto SetModelProperties(Args&&... args) {
+      static_assert(std::is_base_of_v<IModelProperties, TModelProperties>,
+                    "unusable type provided");
+
+      fModelProperties = std::make_shared<TModelProperties>(std::forward<Args>(args)...);
+      return fModelProperties;
+    }
+
+    void SetModelProperties(IMPSharedPtr ptr) { fModelProperties = ptr; }
+
+    template <class MediumType, typename... Args>
+    static auto CreateMedium(Args&&... args) {
+      static_assert(std::is_base_of_v<IMediumModel, MediumType>,
+                    "unusable type provided, needs to be derived from \"IMediumModel\"");
+
+      return std::make_shared<MediumType>(std::forward<Args>(args)...);
+    }
+
+  private:
+    std::vector<VTNUPtr> fChildNodes;
+    std::vector<VolumeTreeNode<IModelProperties> const*> fExcludedNodes;
+    VolumeTreeNode<IModelProperties> const* fParentNode = nullptr;
+    VolUPtr fGeoVolume;
+    IMPSharedPtr fModelProperties;
+  };
+
+} // namespace corsika::environment
+
+#endif

Environment/testEnvironment.cc

+/**
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * 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.
+ */
+
+#define CATCH_CONFIG_MAIN // This tells Catch to provide a main() - only do this in one
+                          // cpp file
+
+#include <corsika/environment/HomogeneousMedium.h>
+#include <corsika/environment/IMediumModel.h>
+#include <corsika/environment/NuclearComposition.h>
+#include <corsika/environment/VolumeTreeNode.h>
+#include <corsika/particles/ParticleProperties.h>
+#include <catch2/catch.hpp>
+
+using namespace corsika::geometry;
+using namespace corsika::environment;
+using namespace corsika::units::si;
+
+TEST_CASE("HomogeneousMedium") {
+  NuclearComposition const protonComposition(
+      std::vector<corsika::particles::Code>{corsika::particles::Code::Proton},
+      std::vector<float>{1.f});
+  HomogeneousMedium<IMediumModel> const medium(19.2_g / cube(1_cm), protonComposition);
+}

Framework/Cascade/Cascade.h

@@ -19,8 +19,6 @@
 
 #include <corsika/setup/SetupTrajectory.h>
 
-using namespace corsika::units::si;
-
 namespace corsika::cascade {
 
   template <typename Tracking, typename ProcessList, typename Stack>

Framework/Cascade/testCascade.cc

-
 /**
  * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
  *
@@ -9,6 +8,10 @@
  * the license.
  */
 
+#include <limits>
+
+#include <corsika/environment/Environment.h>
+
 #include <corsika/cascade/Cascade.h>
 
 #include <corsika/process/ProcessSequence.h>
@@ -36,6 +39,7 @@ using namespace corsika::geometry;
 
 #include <iostream>
 using namespace std;
+using namespace corsika::units::si;
 
 static int fCount = 0;
 
@@ -75,8 +79,15 @@ private:
 };
 
 TEST_CASE("Cascade", "[Cascade]") {
-
-  tracking_line::TrackingLine<setup::Stack> tracking;
+  corsika::environment::Environment env; // dummy environment
+  auto& universe = *(env.GetUniverse());
+  auto const radius = 1_m * std::numeric_limits<double>::infinity();
+  ;
+  auto theMedium = corsika::environment::Environment::CreateNode<Sphere>(
+      Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m}, radius);
+  universe.AddChild(std::move(theMedium));
+
+  tracking_line::TrackingLine<setup::Stack> tracking(env);
 
   stack_inspector::StackInspector<setup::Stack> p0(true);
   ProcessSplit p1;

Framework/Geometry/BaseTrajectory.h

@@ -42,7 +42,10 @@ namespace corsika::geometry {
      * parameterized by \arg t1 and \arg t2. Requires \arg t2 > \arg t1.
      */
 
-    virtual LengthType GetDistance(corsika::units::si::TimeType t1,
+    virtual corsika::units::si::TimeType TimeFromArclength(
+        corsika::units::si::LengthType) const = 0;
+
+    virtual LengthType ArcLength(corsika::units::si::TimeType t1,
                                    corsika::units::si::TimeType t2) const = 0;
 
     virtual corsika::units::si::TimeType GetDuration(

Framework/Geometry/CMakeLists.txt

@@ -10,6 +10,7 @@ set (
   Point.h
   Line.h
   Sphere.h
+  Volume.h
   CoordinateSystem.h
   RootCoordinateSystem.h
   Helix.h

Framework/Geometry/CoordinateSystem.cc

@@ -10,6 +10,7 @@
  */
 
 #include <corsika/geometry/CoordinateSystem.h>
+#include <stdexcept>
 
 using namespace corsika::geometry;
 
@@ -40,7 +41,7 @@ EigenTransform CoordinateSystem::GetTransformation(CoordinateSystem const& pFrom
     commonBase = a;
 
   } else {
-    throw std::string("no connection between coordinate systems found!");
+    throw std::runtime_error("no connection between coordinate systems found!");
   }
 
   EigenTransform t = EigenTransform::Identity();

Framework/Geometry/CoordinateSystem.h

@@ -15,6 +15,7 @@
 #include <corsika/geometry/QuantityVector.h>
 #include <corsika/units/PhysicalUnits.h>
 #include <Eigen/Dense>
+#include <stdexcept>
 
 typedef Eigen::Transform<double, 3, Eigen::Affine> EigenTransform;
 typedef Eigen::Translation<double, 3> EigenTranslation;
@@ -60,7 +61,7 @@ namespace corsika::geometry {
 
     auto rotate(QuantityVector<phys::units::length_d> axis, double angle) const {
       if (axis.eVector.isZero()) {
-        throw std::string("null-vector given as axis parameter");
+        throw std::runtime_error("null-vector given as axis parameter");
       }
 
       EigenTransform const rotation{Eigen::AngleAxisd(angle, axis.eVector.normalized())};
@@ -71,7 +72,7 @@ namespace corsika::geometry {
     auto translateAndRotate(QuantityVector<phys::units::length_d> translation,
                             QuantityVector<phys::units::length_d> axis, double angle) {
       if (axis.eVector.isZero()) {
-        throw std::string("null-vector given as axis parameter");
+        throw std::runtime_error("null-vector given as axis parameter");
       }
 
       EigenTransform const transf{Eigen::AngleAxisd(angle, axis.eVector.normalized()) *

Framework/Geometry/Helix.h

@@ -58,10 +58,15 @@ namespace corsika::geometry {
 
     auto GetRadius() const { return radius; }
 
-    LengthType GetDistanceBetween(corsika::units::si::TimeType t1,
-                                  corsika::units::si::TimeType t2) const {
+    corsika::units::si::LengthType ArcLength(
+      corsika::units::si::TimeType t1, corsika::units::si::TimeType t2) const {
       return (vPar + vPerp).norm() * (t2 - t1);
     }
+
+    corsika::units::si::TimeType TimeFromArclength(
+        corsika::units::si::LengthType t) const {
+      return t / (vPar + vPerp).norm();
+    }
   };
 
 } // namespace corsika::geometry

Framework/Geometry/Line.h

@@ -32,11 +32,18 @@ namespace corsika::geometry {
 
     Point GetPosition(corsika::units::si::TimeType t) const { return r0 + v0 * t; }
 
-    LengthType GetDistanceBetween(corsika::units::si::TimeType t1,
-                                  corsika::units::si::TimeType t2) const {
-      // assert(t2 >= t1);
+    LengthType ArcLength(corsika::units::si::TimeType t1,
+                         corsika::units::si::TimeType t2) const {
       return v0.norm() * (t2 - t1);
     }
+
+    corsika::units::si::TimeType TimeFromArclength(
+        corsika::units::si::LengthType t) const {
+      return t / v0.norm();
+    }
+
+    auto GetR0() const { return r0; }
+    auto GetV0() const { return v0; }
   };
 
 } // namespace corsika::geometry

Framework/Geometry/Sphere.h

@@ -13,23 +13,27 @@
 #define _include_SPHERE_H_
 
 #include <corsika/geometry/Point.h>
+#include <corsika/geometry/Volume.h>
 #include <corsika/units/PhysicalUnits.h>
 
 namespace corsika::geometry {
 
-  class Sphere {
-    Point center;
-    LengthType const radius;
+  class Sphere : public Volume {
+    Point const fCenter;
+    LengthType const fRadius;
 
   public:
     Sphere(Point const& pCenter, LengthType const pRadius)
-        : center(pCenter)
-        , radius(pRadius) {}
+        : fCenter(pCenter)
+        , fRadius(pRadius) {}
 
-    //! returns true if the Point \a p is within the sphere
-    auto Contains(Point const& p) const {
-      return radius * radius > (center - p).squaredNorm();
+    //! returns true if the Point p is within the sphere
+    bool Contains(Point const& p) const override {
+      return fRadius * fRadius > (fCenter - p).squaredNorm();
     }
+
+    auto& GetCenter() const { return fCenter; }
+    auto GetRadius() const { return fRadius; }
   };
 
 } // namespace corsika::geometry

Framework/Geometry/Trajectory.h

@@ -13,6 +13,7 @@
 #define _include_TRAJECTORY_H
 
 #include <corsika/units/PhysicalUnits.h>
+#include <corsika/geometry/Point.h>
 
 using corsika::units::si::LengthType;
 using corsika::units::si::TimeType;
@@ -25,7 +26,7 @@ namespace corsika::geometry {
     corsika::units::si::TimeType fTimeLength;
 
   public:
-    using T::GetDistanceBetween;
+    using T::ArcLength;
     using T::GetPosition;
 
     Trajectory(T const& theT, corsika::units::si::TimeType timeLength)
@@ -36,14 +37,14 @@ namespace corsika::geometry {
       return fTraj.GetPosition(t + fTStart);
       }*/
 
-    Point GetPosition(const double u) const { return T::GetPosition(fTimeLength * u); }
+    Point GetPosition(double u) const { return T::GetPosition(fTimeLength * u); }
 
     TimeType GetDuration() const { return fTimeLength; }
 
-    LengthType GetDistance(const corsika::units::si::TimeType t) const {
+    LengthType GetDistance(corsika::units::si::TimeType t) const {
       assert(t > fTimeLength);
       assert(t >= 0 * corsika::units::si::second);
-      return T::DistanceBetween(0, t);
+      return T::ArcLength(0, t);
     }
   };
 

Framework/Geometry/Vector.h

@@ -181,6 +181,17 @@ namespace corsika::geometry {
                                                              bareResult);
       }
     }
+
+    template <typename dim2>
+    auto dot(Vector<dim2> pV) const {
+      auto const c1 = GetComponents().eVector;
+      auto const c2 = pV.GetComponents(*BaseVector<dim>::cs).eVector;
+      auto const bareResult = c1.dot(c2);
+
+      using ProdQuantity = phys::units::detail::Product<dim, dim2, double, double>;
+
+      return ProdQuantity(phys::units::detail::magnitude_tag, bareResult);
+    }
   };
 
 } // namespace corsika::geometry