diff --git a/corsika/detail/framework/geometry/SeparationPlane.inl b/corsika/detail/framework/geometry/SeparationPlane.inl
new file mode 100644
index 0000000000000000000000000000000000000000..38f2c4529bce0797f417a858ace33947fec97aca
--- /dev/null
+++ b/corsika/detail/framework/geometry/SeparationPlane.inl
@@ -0,0 +1,22 @@
+/*
+ * (c) Copyright 2023 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * 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.
+ */
+
+#pragma once
+
+namespace corsika {
+
+ inline SeparationPlane::SeparationPlane(Plane const& plane)
+ : plane_(plane) {}
+
+ inline bool SeparationPlane::contains(Point const& p) const {
+ return !plane_.isAbove(p);
+ }
+
+ inline std::string SeparationPlane::asString() const { return plane_.asString(); }
+
+} // namespace corsika
\ No newline at end of file
diff --git a/corsika/detail/modules/tracking/TrackingLeapFrogCurved.inl b/corsika/detail/modules/tracking/TrackingLeapFrogCurved.inl
index 9f2370a551053d59114365257e09684b408c32f4..579648ac401952421b71ecb109f9d3f80beb974f 100644
--- a/corsika/detail/modules/tracking/TrackingLeapFrogCurved.inl
+++ b/corsika/detail/modules/tracking/TrackingLeapFrogCurved.inl
@@ -345,13 +345,23 @@ namespace corsika {
return tracking_line::Tracking::intersect(particle, plane);
}
+ template <typename TParticle>
+ inline Intersections Tracking::intersect(TParticle const& particle,
+ SeparationPlane const& sepPlane) {
+ return intersect(particle, sepPlane.getPlane());
+ }
+
template <typename TParticle, typename TBaseNodeType>
inline Intersections Tracking::intersect(TParticle const& particle,
TBaseNodeType const& volumeNode) {
Sphere const* sphere = dynamic_cast<Sphere const*>(&volumeNode.getVolume());
if (sphere) { return intersect(particle, *sphere); }
+ SeparationPlane const* sepPlane =
+ dynamic_cast<SeparationPlane const*>(&volumeNode.getVolume());
+ if (sepPlane) { return intersect(particle, *sepPlane); }
throw std::runtime_error(
- "The Volume type provided is not supported in intersect(particle, node)");
+ "The Volume type provided is not supported in "
+ "TrackingLeapFrogCurved::intersect(particle, node)");
}
template <typename TParticle>
diff --git a/corsika/detail/modules/tracking/TrackingStraight.inl b/corsika/detail/modules/tracking/TrackingStraight.inl
index debb95ec9a22bc75dfa593541804dc71a5cd025f..98d58d4408820a8bd1c79fa217e7ac570742c08f 100644
--- a/corsika/detail/modules/tracking/TrackingStraight.inl
+++ b/corsika/detail/modules/tracking/TrackingStraight.inl
@@ -125,10 +125,14 @@ namespace corsika::tracking_line {
return Tracking::intersect<TParticle>(particle, *sphere);
} else if (Box const* box = dynamic_cast<Box const*>(&volumeNode.getVolume()); box) {
return Tracking::intersect<TParticle>(particle, *box);
+ } else if (SeparationPlane const* sepPlane =
+ dynamic_cast<SeparationPlane const*>(&volumeNode.getVolume());
+ sepPlane) {
+ return Tracking::intersect<TParticle>(particle, *sepPlane);
} else {
throw std::runtime_error(
"The Volume type provided is not supported in "
- "Intersect(particle, node)");
+ "TrackingStraight::intersect(particle, node)");
}
}
@@ -149,4 +153,10 @@ namespace corsika::tracking_line {
return Intersections(n.dot(delta) / n_dot_v);
}
+ template <typename TParticle>
+ inline Intersections Tracking::intersect(TParticle const& particle,
+ SeparationPlane const& sepPlane) {
+ return intersect(particle, sepPlane.getPlane());
+ }
+
} // namespace corsika::tracking_line
diff --git a/corsika/framework/geometry/Box.hpp b/corsika/framework/geometry/Box.hpp
index 696ad54555f8418d964e1b7f4f9ab1bea6d68eaf..b4e345465e047ec07c176ba1e54a75fbe09795a4 100644
--- a/corsika/framework/geometry/Box.hpp
+++ b/corsika/framework/geometry/Box.hpp
@@ -16,10 +16,11 @@
namespace corsika {
/**
- * Describes a sphere in space
+ * Describes a box in space
*
* The center point and the orintation of the Box is set by
- * a CoordinateSystemPtr at construction.
+ * a CoordinateSystemPtr at construction and the sides extend
+ * by x, y, z in both directions.
**/
class Box : public IVolume {
diff --git a/corsika/framework/geometry/SeparationPlane.hpp b/corsika/framework/geometry/SeparationPlane.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..82f7b358f1f83476f142ecced381ed680b4b244a
--- /dev/null
+++ b/corsika/framework/geometry/SeparationPlane.hpp
@@ -0,0 +1,43 @@
+/*
+ * (c) Copyright 2023 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * 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.
+ */
+
+#pragma once
+
+#include <corsika/framework/geometry/Plane.hpp>
+#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/framework/geometry/IVolume.hpp>
+
+namespace corsika {
+
+ /**
+ * Describes which exists on one side of a plane
+ *
+ * The separation line is given by a plane where the unit
+ * vector of the plane points outside of the volume
+ **/
+ class SeparationPlane : public IVolume {
+
+ public:
+ SeparationPlane(Plane const& plane);
+
+ ~SeparationPlane() {}
+
+ Plane getPlane() const { return plane_; }
+
+ //! returns true if the Point p is below the plane
+ bool contains(Point const& p) const override;
+
+ std::string asString() const;
+
+ protected:
+ Plane const plane_;
+ };
+
+} // namespace corsika
+
+#include <corsika/detail/framework/geometry/SeparationPlane.inl>
diff --git a/corsika/modules/tracking/TrackingLeapFrogCurved.hpp b/corsika/modules/tracking/TrackingLeapFrogCurved.hpp
index ec31314b489e020a1cbb6a2d5f42d22fe79cd726..2704c5555a7943c257483865d7beb1f8fdf17eaf 100644
--- a/corsika/modules/tracking/TrackingLeapFrogCurved.hpp
+++ b/corsika/modules/tracking/TrackingLeapFrogCurved.hpp
@@ -12,6 +12,7 @@
#include <corsika/framework/geometry/Line.hpp>
#include <corsika/framework/geometry/Plane.hpp>
#include <corsika/framework/geometry/Sphere.hpp>
+#include <corsika/framework/geometry/SeparationPlane.hpp>
#include <corsika/framework/geometry/LeapFrogTrajectory.hpp>
#include <corsika/framework/geometry/Vector.hpp>
#include <corsika/framework/geometry/Intersections.hpp>
@@ -94,6 +95,21 @@ namespace corsika {
template <typename TParticle>
static Intersections intersect(TParticle const& particle, Plane const& plane);
+ /**
+ * find intersection of Separation with Track
+ *
+ * Intersection times of particle are caculated by
+ * and using the intersect-with-plane routine
+ *
+ * @tparam TParticle Type of particle object on stack.
+ * @param particle Particle initial state.
+ * @param sepPlane Separation.
+ * @return Intersections in time units.
+ */
+ template <typename TParticle>
+ static Intersections intersect(TParticle const& particle,
+ SeparationPlane const& sepPlane);
+
static std::string getName() { return "LeapFrog-curved"; }
static std::string getVersion() { return "1.0.0"; }
diff --git a/corsika/modules/tracking/TrackingStraight.hpp b/corsika/modules/tracking/TrackingStraight.hpp
index 68ae17eb1fbbb662ef9bd9c430ff02b4f7d7178e..1bee94a66365179a9c810ddaa3abdedefdac53c0 100644
--- a/corsika/modules/tracking/TrackingStraight.hpp
+++ b/corsika/modules/tracking/TrackingStraight.hpp
@@ -13,6 +13,7 @@
#include <corsika/framework/geometry/Plane.hpp>
#include <corsika/framework/geometry/Sphere.hpp>
#include <corsika/framework/geometry/Box.hpp>
+#include <corsika/framework/geometry/SeparationPlane.hpp>
#include <corsika/framework/geometry/Vector.hpp>
#include <corsika/framework/geometry/StraightTrajectory.hpp>
#include <corsika/framework/geometry/Intersections.hpp>
@@ -56,6 +57,11 @@ namespace corsika::tracking_line {
template <typename TParticle>
static Intersections intersect(TParticle const& particle, Plane const& plane);
+ //! find intersection of SeparationPlane with Track
+ template <typename TParticle>
+ static Intersections intersect(TParticle const& particle,
+ SeparationPlane const& sepPlane);
+
static std::string getName() { return "Tracking-Straight"; }
static std::string getVersion() { return "1.0.0"; }
};
diff --git a/tests/framework/testGeometry.cpp b/tests/framework/testGeometry.cpp
index 8cc52855a21de3b8de2f15b409ff9c33dcb5656c..fcd29c001598f89fb1a90878fc294cdc7c44bd72 100644
--- a/tests/framework/testGeometry.cpp
+++ b/tests/framework/testGeometry.cpp
@@ -11,6 +11,7 @@
#include <cmath>
#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/geometry/Box.hpp>
+#include <corsika/framework/geometry/SeparationPlane.hpp>
#include <corsika/framework/geometry/CoordinateSystem.hpp>
#include <corsika/framework/geometry/Helix.hpp>
#include <corsika/framework/geometry/LeapFrogTrajectory.hpp>
@@ -344,6 +345,41 @@ TEST_CASE("Geometry Box") {
}
}
+TEST_CASE("Geometry SeparationPlane") {
+ CoordinateSystemPtr const& rootCS = get_root_CoordinateSystem();
+ Point const planeCenter = Point(rootCS, {0_m, 0_m, 0_m});
+
+ SECTION("constructor") {
+ DirectionVector const planeNorm{rootCS, {0, 0, 1}};
+ SeparationPlane const sepPlane{Plane(planeCenter, planeNorm)};
+ CHECK(sepPlane.asString() != "");
+ }
+
+ SECTION("isInside") {
+ DirectionVector const planeNorm{rootCS, {0, 0, 1}};
+ SeparationPlane const sepPlane{Plane(planeCenter, planeNorm)};
+
+ CHECK_FALSE(sepPlane.contains(Point(rootCS, {0_m, 0_m, 1_m})));
+ CHECK_FALSE(sepPlane.contains(Point(rootCS, {1_m, 0_m, 1_m})));
+ CHECK_FALSE(sepPlane.contains(Point(rootCS, {-1_m, 0_m, 1_m})));
+ CHECK_FALSE(sepPlane.contains(Point(rootCS, {0_m, 1_m, 1_m})));
+ CHECK_FALSE(sepPlane.contains(Point(rootCS, {0_m, -1_m, 1_m})));
+
+ CHECK(sepPlane.contains(Point(rootCS, {0_m, 0_m, -1_m})));
+ CHECK(sepPlane.contains(Point(rootCS, {1_m, 0_m, -1_m})));
+ CHECK(sepPlane.contains(Point(rootCS, {-1_m, 0_m, -1_m})));
+ CHECK(sepPlane.contains(Point(rootCS, {0_m, 1_m, -1_m})));
+ CHECK(sepPlane.contains(Point(rootCS, {0_m, -1_m, -1_m})));
+ }
+
+ SECTION("getPlane") {
+ DirectionVector const planeNorm{rootCS, {0, 0, 1}};
+ SeparationPlane const sepPlane{Plane(planeCenter, planeNorm)};
+
+ auto planeCheck = sepPlane.getPlane();
+ }
+}
+
TEST_CASE("Geometry Trajectories") {
CoordinateSystemPtr rootCS = get_root_CoordinateSystem();
Point r0(rootCS, {0_m, 0_m, 0_m});
diff --git a/tests/modules/testTracking.cpp b/tests/modules/testTracking.cpp
index 8fabcb6918b5a5b864913655f93db0a3a9193ea2..b66d160dc7b6ed4660320dbdb3e7a538d60e895e 100644
--- a/tests/modules/testTracking.cpp
+++ b/tests/modules/testTracking.cpp
@@ -332,6 +332,75 @@ TEMPLATE_TEST_CASE("TrackingFail", "doesntwork", tracking_leapfrog_curved::Track
CHECK_THROWS(tracking.intersect(particle, dummy));
}
+TEMPLATE_TEST_CASE("TrackingSeparationPlane", "separation_plane",
+ tracking_leapfrog_curved::Tracking,
+ tracking_leapfrog_straight::Tracking, tracking_line::Tracking) {
+
+ logging::set_level(logging::level::info);
+
+ const HEPEnergyType P0 = 10_GeV;
+
+ auto PID = GENERATE(as<Code>{}, Code::MuPlus, Code::MuPlus, Code::Photon);
+ // for algorithms that know magnetic deflections choose: +-50uT, 0uT
+ // otherwise just 0uT
+ auto Bfield = GENERATE(filter(
+ []([[maybe_unused]] MagneticFluxType v) {
+ if constexpr (std::is_same_v<TestType, tracking_line::Tracking>)
+ return v == 0_uT;
+ else
+ return true;
+ },
+ values<MagneticFluxType>({50_uT, 0_uT, -50_uT})));
+
+ SECTION(fmt::format("Tracking PID={}, Bfield={} uT", PID, Bfield / 1_uT)) {
+
+ CORSIKA_LOG_DEBUG(
+ "********************\n TEST algo={} section PID={}, "
+ "Bfield={}uT ",
+ boost::typeindex::type_id<TestType>().pretty_name(), PID, Bfield / 1_uT);
+
+ const int chargeNumber = get_charge_number(PID);
+ LengthType radius = 10_m;
+ int deflect = 0;
+ if (chargeNumber != 0 and Bfield != 0_T) {
+ deflect = 1;
+ LengthType const gyroradius = (convert_HEP_to_SI<MassType::dimension_type>(P0) *
+ constants::c / (abs(get_charge(PID)) * abs(Bfield)));
+ CORSIKA_LOG_DEBUG("Rg={} deflect={}", gyroradius, deflect);
+ radius = gyroradius;
+ }
+
+ auto [env, csPtr, worldPtr] =
+ corsika::setup::testing::setup_environment(Code::Oxygen, Bfield);
+ { [[maybe_unused]] const auto& env_dummy = env; }
+ auto const& cs = *csPtr;
+
+ TestType tracking;
+ Point const center(cs, {0_m, 0_m, 0_m});
+
+ auto [stack, viewPtr] = setup::testing::setup_stack(PID, P0, worldPtr, cs);
+ { [[maybe_unused]] auto& viewPtr_dum = viewPtr; }
+ auto particle = stack->first();
+ // Note: momentum in X-direction
+ // magnetic field in Z-direction
+
+ particle.setPosition(Point(cs, -radius, 0_m, 0_m));
+
+ // put plane where we expect deflection by radius/2
+ Plane const plane(Point(cs, radius * (1 - sqrt(3. / 4.)), 0_m, 0_m),
+ DirectionVector(cs, {-1., 0., 0.}));
+ SeparationPlane const sepPlane{plane};
+
+ Intersections const hit = tracking.intersect(particle, sepPlane);
+
+ CORSIKA_LOG_DEBUG("entry={} exit={}", hit.getEntry(), hit.getExit());
+
+ CHECK(hit.hasIntersections());
+ CHECK(hit.getEntry() / 1_s == Approx(0.00275 * deflect).margin(0.0003));
+ CHECK(hit.getExit() == 1_s * std::numeric_limits<double>::infinity());
+ }
+}
+
TEMPLATE_TEST_CASE("TrackingPlane", "plane", tracking_leapfrog_curved::Tracking,
tracking_leapfrog_straight::Tracking, tracking_line::Tracking) {