Add SeparationPlane volume object

corsika/detail/framework/geometry/SeparationPlane.inl

+/*
+ * (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

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>

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

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 {
 

corsika/framework/geometry/SeparationPlane.hpp

+/*
+ * (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>

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"; }
 

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"; }
   };

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});

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) {