From 4b1691d935c7a334e1ac55d665e78b6fdf335e41 Mon Sep 17 00:00:00 2001
From: Maximilian Reininghaus <maximilian.reininghaus@kit.edu>
Date: Tue, 4 Dec 2018 22:00:49 +0100
Subject: [PATCH] implemented propagation to next volume boundary for line
trajectory
---
Documentation/Examples/cascade_example.cc | 7 +++-
Environment/Environment.h | 2 +
Environment/VolumeTreeNode.h | 10 +++--
Processes/TrackingLine/TrackingLine.h | 45 +++++++++++++++++++---
Processes/TrackingLine/testTrackingLine.cc | 9 +++++
5 files changed, 62 insertions(+), 11 deletions(-)
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index 9b664046b..dd776a1f3 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/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);
diff --git a/Environment/Environment.h b/Environment/Environment.h
index 0bc2fa112..65de554e7 100644
--- a/Environment/Environment.h
+++ b/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();
}
diff --git a/Environment/VolumeTreeNode.h b/Environment/VolumeTreeNode.h
index 470bb8aa0..f5779f9ad 100644
--- a/Environment/VolumeTreeNode.h
+++ b/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; }
diff --git a/Processes/TrackingLine/TrackingLine.h b/Processes/TrackingLine/TrackingLine.h
index 386215801..99fe8272a 100644
--- a/Processes/TrackingLine/TrackingLine.h
+++ b/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);
}
};
diff --git a/Processes/TrackingLine/testTrackingLine.cc b/Processes/TrackingLine/testTrackingLine.cc
index a07101096..050b1b932 100644
--- a/Processes/TrackingLine/testTrackingLine.cc
+++ b/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));
+ }
}
--
GitLab