From af464a9df28c51cb9ae13d2b29e9a93710bb2ebe Mon Sep 17 00:00:00 2001
From: Maximilian Reininghaus <maximilian.reininghaus@tu-dortmund.de>
Date: Wed, 5 Dec 2018 01:01:50 +0100
Subject: [PATCH] propagate only in forward direction, note that testCascade is
failing currently
---
Documentation/Examples/cascade_example.cc | 5 ++-
Environment/Environment.h | 12 +++++--
Framework/Cascade/testCascade.cc | 7 +++-
Framework/Geometry/Vector.h | 2 +-
Processes/StackInspector/StackInspector.cc | 1 +
Processes/TrackingLine/TrackingLine.h | 15 ++++++---
Processes/TrackingLine/testTrackingLine.cc | 39 +++++++++++++++++++---
7 files changed, 65 insertions(+), 16 deletions(-)
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index dd776a1f..2203a2cb 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -245,11 +245,10 @@ double s_rndm_(int&) {
}
int main() {
- Environment env;
- //~ auto& universe = env.GetUniverse();
+ corsika::environment::Environment env; // dummy environment
auto& universe = *(env.GetUniverse());
- auto const theMedium = Environment::CreateNode<Sphere>(
+ auto theMedium = corsika::environment::Environment::CreateNode<Sphere>(
Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m}, 100_km);
universe.AddChild(std::move(theMedium));
diff --git a/Environment/Environment.h b/Environment/Environment.h
index 65de554e..285659a6 100644
--- a/Environment/Environment.h
+++ b/Environment/Environment.h
@@ -17,13 +17,19 @@
#include <corsika/setup/SetupEnvironment.h>
namespace corsika::environment {
+ struct Universe : public corsika::geometry::Volume {
+ bool Contains(corsika::geometry::Point const&) const override {return true;}
+ };
class Environment {
public:
+ Environment() : fUniverse(std::make_unique<VolumeTreeNode<IEnvironmentModel>>(
+ std::make_unique<Universe>())) {}
+
using IEnvironmentModel = corsika::setup::IEnvironmentModel;
- auto& GetUniverse() { return universe; }
- auto const& GetUniverse() const { return universe; }
+ auto& GetUniverse() { return fUniverse; }
+ auto const& GetUniverse() const { return fUniverse; }
auto const& GetCoordinateSystem() const {
return corsika::geometry::RootCoordinateSystem::GetInstance().GetRootCS();
@@ -41,7 +47,7 @@ namespace corsika::environment {
}
private:
- VolumeTreeNode<IEnvironmentModel>::VTNUPtr universe;
+ VolumeTreeNode<IEnvironmentModel>::VTNUPtr fUniverse;
};
} // namespace corsika::environment
diff --git a/Framework/Cascade/testCascade.cc b/Framework/Cascade/testCascade.cc
index 5017d5ac..5419838f 100644
--- a/Framework/Cascade/testCascade.cc
+++ b/Framework/Cascade/testCascade.cc
@@ -78,8 +78,13 @@ private:
TEST_CASE("Cascade", "[Cascade]") {
corsika::environment::Environment env; // dummy environment
+ auto& universe = *(env.GetUniverse());
+ auto const radius = 20_km;
+ 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);
+ tracking_line::TrackingLine<setup::Stack> tracking(env);
stack_inspector::StackInspector<setup::Stack> p0(true);
ProcessSplit p1;
diff --git a/Framework/Geometry/Vector.h b/Framework/Geometry/Vector.h
index 904b7039..94394b90 100644
--- a/Framework/Geometry/Vector.h
+++ b/Framework/Geometry/Vector.h
@@ -190,7 +190,7 @@ namespace corsika::geometry {
using ProdQuantity = phys::units::detail::Product<dim, dim2, double, double>;
- return ProdQuantity(bareResult);
+ return ProdQuantity(phys::units::detail::magnitude_tag, bareResult);
}
};
diff --git a/Processes/StackInspector/StackInspector.cc b/Processes/StackInspector/StackInspector.cc
index b0acefc4..1986c83b 100644
--- a/Processes/StackInspector/StackInspector.cc
+++ b/Processes/StackInspector/StackInspector.cc
@@ -12,6 +12,7 @@
#include <corsika/geometry/RootCoordinateSystem.h>
#include <corsika/process/stack_inspector/StackInspector.h>
#include <corsika/units/PhysicalUnits.h>
+#include <corsika/particles/ParticleProperties.h>
#include <corsika/logging/Logger.h>
diff --git a/Processes/TrackingLine/TrackingLine.h b/Processes/TrackingLine/TrackingLine.h
index 99fe8272..a53fbc6c 100644
--- a/Processes/TrackingLine/TrackingLine.h
+++ b/Processes/TrackingLine/TrackingLine.h
@@ -14,6 +14,7 @@
#include <algorithm>
#include <iostream>
#include <optional>
+#include <utility>
using namespace corsika;
@@ -28,7 +29,7 @@ namespace corsika::process {
corsika::environment::Environment const& fEnvironment;
public:
- std::optional<corsika::units::si::TimeType> TimeOfIntersection(
+ std::optional<std::pair<corsika::units::si::TimeType, corsika::units::si::TimeType>> TimeOfIntersection(
corsika::geometry::Line const& traj, geometry::Sphere const& sphere) {
using namespace corsika::units::si;
auto const& cs = fEnvironment.GetCoordinateSystem();
@@ -49,8 +50,8 @@ namespace corsika::process {
//~ std::cout << "beta: " << beta << std::endl;
if (discriminant.magnitude() > 0) {
- auto const t = (-alpha - sqrt(discriminant)) / (2 * v0.squaredNorm());
- return t;
+ (-alpha - sqrt(discriminant)) / (2 * v0.squaredNorm());
+ return std::make_pair((-alpha - sqrt(discriminant)) / (2 * v0.squaredNorm()), (-alpha + sqrt(discriminant)) / (2 * v0.squaredNorm()));
} else {
return {};
}
@@ -82,7 +83,13 @@ namespace corsika::process {
// everything is a sphere, crashes with exception if not
if (auto opt = TimeOfIntersection(line, sphere); opt.has_value())
- intersectionTimes.push_back(*opt);
+ {
+ auto const [t1, t2] = *opt;
+ if (t1.magnitude() >= 0)
+ intersectionTimes.push_back(t1);
+ else if (t2.magnitude() >= 0)
+ intersectionTimes.push_back(t2);
+ }
};
for (auto const& child : children) { addIfIntersects(*child); }
diff --git a/Processes/TrackingLine/testTrackingLine.cc b/Processes/TrackingLine/testTrackingLine.cc
index 050b1b93..ef44f239 100644
--- a/Processes/TrackingLine/testTrackingLine.cc
+++ b/Processes/TrackingLine/testTrackingLine.cc
@@ -32,12 +32,27 @@ using namespace corsika::geometry;
using namespace std;
using namespace corsika::units::si;
+struct DummyParticle {
+ EnergyType fEnergy;
+ Vector<momentum_d> fMomentum;
+ Point fPosition;
+
+ DummyParticle(EnergyType pEnergy, Vector<momentum_d> pMomentum, Point pPosition) : fEnergy(pEnergy), fMomentum(pMomentum), fPosition(pPosition) {}
+
+ auto GetEnergy() const { return fEnergy; }
+ auto GetMomentum() const { return fMomentum; }
+ auto GetPosition() const { return fPosition; }
+};
+
+struct DummyStack {
+ using ParticleType = DummyParticle;
+};
TEST_CASE("TrackingLine") {
corsika::environment::Environment env; // dummy environment
auto const& cs = env.GetCoordinateSystem();
- tracking_line::TrackingLine<setup::Stack> tracking(env);
+ tracking_line::TrackingLine<DummyStack> tracking(env);
SECTION("intersection with sphere") {
Point const origin(cs, {0_m, 0_m, 0_m});
@@ -50,7 +65,10 @@ TEST_CASE("TrackingLine") {
auto const opt = tracking.TimeOfIntersection(traj, Sphere(Point(cs, {0_m, 0_m, 10_m}), 1_m));
REQUIRE(opt.has_value());
- REQUIRE(opt.value() / 9_s == Approx(1));
+
+ auto [t1, t2] = opt.value();
+ REQUIRE(t1 / 9_s == Approx(1));
+ REQUIRE(t2 / 11_s == Approx(1));
auto const optNoIntersection = tracking.TimeOfIntersection(traj, Sphere(Point(cs, {5_m, 0_m, 10_m}), 1_m));
REQUIRE_FALSE(optNoIntersection.has_value());
@@ -59,9 +77,22 @@ TEST_CASE("TrackingLine") {
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);
+ //~ std::cout << env.GetUniverse().get() << std::endl;
+
+ DummyParticle p(1_J, Vector<momentum_d>(cs, 0*kilogram*meter/second, 0*kilogram*meter/second, 1*kilogram*meter/second), Point(cs, 0_m, 0_m,0_m));
+ auto const radius = 20_m;
+
+ auto theMedium = corsika::environment::Environment::CreateNode<Sphere>(
+ Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m}, radius);
universe.AddChild(std::move(theMedium));
+
+ Point const origin(cs, {0_m, 0_m, 0_m});
+ Vector<corsika::units::si::SpeedType::dimension_type> v(cs, 0_m/second, 0_m/second, 1_m/second);
+ Line line(origin, v);
+
+ auto const traj = tracking.GetTrack(p);
+
+ REQUIRE((traj.GetPosition(1.) - Point(cs, 0_m, 0_m, radius)).GetComponents(cs).norm().magnitude() == Approx(0).margin(1e-4));
}
}
--
GitLab