implemented propagation to next volume boundary for line trajectory

Documentation/Examples/cascade_example.cc

@@ -246,10 +246,13 @@ double s_rndm_(int&) {
 
 int main() {
   Environment env;
-  auto& universe = env.GetUniverse();
+  //~ auto& universe = env.GetUniverse();
+  auto& universe = *(env.GetUniverse());
 
-  auto theMedium = Environment::CreateNode<Sphere>(
+  auto const theMedium = Environment::CreateNode<Sphere>(
       Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m}, 100_km);
+      
+  universe.AddChild(std::move(theMedium));
 
   tracking_line::TrackingLine<setup::Stack> tracking(env);
   stack_inspector::StackInspector<setup::Stack> p0(true);

Environment/Environment.h

@@ -23,6 +23,8 @@ namespace corsika::environment {
     using IEnvironmentModel = corsika::setup::IEnvironmentModel;
 
     auto& GetUniverse() { return universe; }
+    auto const& GetUniverse() const { return universe; }
+
     auto const& GetCoordinateSystem() const {
       return corsika::geometry::RootCoordinateSystem::GetInstance().GetRootCS();
     }

Environment/VolumeTreeNode.h

@@ -65,7 +65,7 @@ namespace corsika::environment {
       }
     }
 
-    void addChild(VTNUPtr pChild) {
+    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
@@ -73,11 +73,15 @@ namespace corsika::environment {
       // later and the caller won't notice.
     }
 
-    void excludeOverlapWith(VTNUPtr const& pNode) {
+    void ExcludeOverlapWith(VTNUPtr const& pNode) {
       fExcludedNodes.push_back(pNode.get());
     }
 
-    auto GetParent() const { return fParentNode; };
+    auto* GetParent() const { return fParentNode; };
+
+    auto const& GetChildNodes() const { return fChildNodes; }
+
+    auto const& GetExcludedNodes() const { return fExcludedNodes; }
 
     auto const& GetVolume() const { return *fGeoVolume; }
 

Processes/TrackingLine/TrackingLine.h

@@ -11,9 +11,9 @@
 #include <corsika/setup/SetupStack.h>
 #include <corsika/setup/SetupTrajectory.h>
 
-#include <optional>
-
+#include <algorithm>
 #include <iostream>
+#include <optional>
 
 using namespace corsika;
 
@@ -29,8 +29,7 @@ namespace corsika::process {
 
     public:
       std::optional<corsika::units::si::TimeType> TimeOfIntersection(
-          geometry::Trajectory<corsika::geometry::Line> const& traj,
-          geometry::Sphere const& sphere) {
+          corsika::geometry::Line const& traj, geometry::Sphere const& sphere) {
         using namespace corsika::units::si;
         auto const& cs = fEnvironment.GetCoordinateSystem();
         geometry::Point const origin(cs, 0_m, 0_m, 0_m);
@@ -60,12 +59,46 @@ namespace corsika::process {
       TrackingLine(corsika::environment::Environment const& pEnv)
           : fEnvironment(pEnv) {}
       void Init() {}
+
       auto GetTrack(Particle& p) {
         using namespace corsika::units::si;
         geometry::Vector<SpeedType::dimension_type> const velocity =
             p.GetMomentum() / p.GetEnergy() * corsika::units::si::constants::cSquared;
-        geometry::Line traj(p.GetPosition(), velocity);
-        return geometry::Trajectory<corsika::geometry::Line>(traj, 100_ns);
+
+        auto const currentPosition = p.GetPosition();
+        geometry::Line line(currentPosition, velocity);
+
+        auto const* currentVolumeNode =
+            fEnvironment.GetUniverse()->GetContainingNode(currentPosition);
+        auto const& children = currentVolumeNode->GetChildNodes();
+        auto const& excluded = currentVolumeNode->GetExcludedNodes();
+
+        std::vector<TimeType> intersectionTimes;
+
+        auto addIfIntersects = [&](auto& vtn) {
+          auto const& volume = vtn.GetVolume();
+          auto const& sphere = dynamic_cast<geometry::Sphere const&>(
+              volume); // for the moment we are a bit bold here and assume
+                       // everything is a sphere, crashes with exception if not
+
+          if (auto opt = TimeOfIntersection(line, sphere); opt.has_value())
+            intersectionTimes.push_back(*opt);
+        };
+
+        for (auto const& child : children) { addIfIntersects(*child); }
+
+        for (auto const* child : excluded) { addIfIntersects(*child); }
+
+        addIfIntersects(*currentVolumeNode);
+
+        auto const minIter =
+            std::min_element(intersectionTimes.cbegin(), intersectionTimes.cend());
+
+        if (minIter == intersectionTimes.cend()) {
+          throw std::string("no intersection with anything!");
+        }
+
+        return geometry::Trajectory<corsika::geometry::Line>(line, *minIter);
       }
     };
 

Processes/TrackingLine/testTrackingLine.cc

@@ -55,4 +55,13 @@ TEST_CASE("TrackingLine") {
     auto const optNoIntersection = tracking.TimeOfIntersection(traj, Sphere(Point(cs, {5_m, 0_m, 10_m}), 1_m));
     REQUIRE_FALSE(optNoIntersection.has_value());
   }
+  
+  SECTION("maximally possible propagation") {
+      auto& universe = *(env.GetUniverse());
+      
+      auto theMedium = corsika::environment::Environment::CreateNode<Sphere>(
+      Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m}, 100_km);
+      
+      universe.AddChild(std::move(theMedium));
+  }
 }