diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc index 9b664046b1c4fc51dc1e7ceae7795247cf881120..dd776a1f326186d39b8a95ef1fad315b0b26be91 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 0bc2fa11273b615edaec5089ee5527ecd2808b47..65de554e79b4aa105939a26625de6d4878c1d7a5 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 470bb8aa0c301fc9494417282d9bb32aad181927..f5779f9ad339528e56ba4a5f14078af5afcca8e2 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 386215801676a144dd28a8ab123e12e15981324e..99fe8272ac09ad9a95c590d72939392daea5e92b 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 a07101096e13eda22daec046fb3390df03dc9882..050b1b932a52bacf83d77b216f7cdd440e3f53eb 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)); + } }