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 + inline Intersections Tracking::intersect(TParticle const& particle, + SeparationPlane const& sepPlane) { + return intersect(particle, sepPlane.getPlane()); + } + template inline Intersections Tracking::intersect(TParticle const& particle, TBaseNodeType const& volumeNode) { Sphere const* sphere = dynamic_cast(&volumeNode.getVolume()); if (sphere) { return intersect(particle, *sphere); } + SeparationPlane const* sepPlane = + dynamic_cast(&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 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(particle, *sphere); } else if (Box const* box = dynamic_cast(&volumeNode.getVolume()); box) { return Tracking::intersect(particle, *box); + } else if (SeparationPlane const* sepPlane = + dynamic_cast(&volumeNode.getVolume()); + sepPlane) { + return Tracking::intersect(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 + 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 +#include +#include + +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 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 #include #include +#include #include #include #include @@ -94,6 +95,21 @@ namespace corsika { template 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 + 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 #include #include +#include #include #include #include @@ -56,6 +57,11 @@ namespace corsika::tracking_line { template static Intersections intersect(TParticle const& particle, Plane const& plane); + //! find intersection of SeparationPlane with Track + template + 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 #include #include +#include #include #include #include @@ -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::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) + return v == 0_uT; + else + return true; + }, + values({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().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(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::infinity()); + } +} + TEMPLATE_TEST_CASE("TrackingPlane", "plane", tracking_leapfrog_curved::Tracking, tracking_leapfrog_straight::Tracking, tracking_line::Tracking) {