"transparent" observation plane

Processes/ObservationPlane/CMakeLists.txt

@@ -16,7 +16,7 @@ set (
 add_library (ProcessObservationPlane STATIC ${MODEL_SOURCES})
 CORSIKA_COPY_HEADERS_TO_NAMESPACE (ProcessObservationPlane ${MODEL_NAMESPACE} ${MODEL_HEADERS})
 
-# target dependencies on other libraries (also the header onlys)
+#target dependencies on other libraries(also the header onlys)
 target_link_libraries (
   ProcessObservationPlane
   CORSIKAgeometry
@@ -32,8 +32,8 @@ target_include_directories (
 
 install (FILES ${MODEL_HEADERS} DESTINATION include/${MODEL_NAMESPACE})
 
-# --------------------
-# code unit testing
+#-- -- -- -- -- -- -- -- -- --
+#code unit testing
 CORSIKA_ADD_TEST(testObservationPlane)
 target_link_libraries (
   testObservationPlane
@@ -41,4 +41,3 @@ target_link_libraries (
   CORSIKAstackinterface
   CORSIKAthirdparty # for catch2
 )
-

Processes/ObservationPlane/ObservationPlane.cc

@@ -15,36 +15,48 @@
 using namespace corsika::process::observation_plane;
 using namespace corsika::units::si;
 
-ObservationPlane::ObservationPlane(geometry::Plane const& vObsPlane,
-                                   std::string const& vFilename)
-    : fObsPlane(vObsPlane)
-    , fOutputStream(vFilename) {
-  fOutputStream << "#PDG code, energy / eV, distance to center / m" << std::endl;
+ObservationPlane::ObservationPlane(geometry::Plane const& obsPlane,
+                                   std::string const& filename, bool deleteOnHit)
+    : plane_(obsPlane)
+    , outputStream_(filename)
+    , deleteOnHit_(deleteOnHit) {
+  outputStream_ << "#PDG code, energy / eV, distance to center / m" << std::endl;
 }
 
 corsika::process::EProcessReturn ObservationPlane::DoContinuous(
-    setup::Stack::ParticleType const& vParticle, setup::Trajectory const& vTrajectory) {
+    setup::Stack::ParticleType const& particle, setup::Trajectory const& trajectory) {
+  TimeType const timeOfIntersection =
+      (plane_.GetCenter() - trajectory.GetR0()).dot(plane_.GetNormal()) /
+      trajectory.GetV0().dot(plane_.GetNormal());
 
-  if (fObsPlane.IsAbove(vTrajectory.GetPosition(1.0001))) {
+  if (timeOfIntersection < TimeType::zero()) { return process::EProcessReturn::eOk; }
+
+  if (plane_.IsAbove(trajectory.GetR0()) == plane_.IsAbove(trajectory.GetPosition(1))) {
     return process::EProcessReturn::eOk;
   }
 
-  fOutputStream << static_cast<int>(particles::GetPDG(vParticle.GetPID())) << ' '
-                << vParticle.GetEnergy() * (1 / 1_eV) << ' '
-                << (vTrajectory.GetPosition(1) - fObsPlane.GetCenter()).norm() / 1_m
+  outputStream_ << static_cast<int>(particles::GetPDG(particle.GetPID())) << ' '
+                << particle.GetEnergy() * (1 / 1_eV) << ' '
+                << (trajectory.GetPosition(1) - plane_.GetCenter()).norm() / 1_m
                 << std::endl;
 
-  return process::EProcessReturn::eParticleAbsorbed;
+  if (deleteOnHit_) {
+    return process::EProcessReturn::eParticleAbsorbed;
+  } else {
+    return process::EProcessReturn::eOk;
+  }
 }
 
 LengthType ObservationPlane::MaxStepLength(setup::Stack::ParticleType const&,
-                                           setup::Trajectory const& vTrajectory) {
-  if (!fObsPlane.IsAbove(vTrajectory.GetR0())) {
+                                           setup::Trajectory const& trajectory) {
+  TimeType const timeOfIntersection =
+      (plane_.GetCenter() - trajectory.GetR0()).dot(plane_.GetNormal()) /
+      trajectory.GetV0().dot(plane_.GetNormal());
+
+  if (timeOfIntersection < TimeType::zero()) {
     return std::numeric_limits<double>::infinity() * 1_m;
   }
 
-  auto const pointOfIntersection = vTrajectory.GetPosition(
-      (fObsPlane.GetCenter() - vTrajectory.GetR0()).dot(fObsPlane.GetNormal()) /
-      vTrajectory.GetV0().dot(fObsPlane.GetNormal()));
-  return (vTrajectory.GetR0() - pointOfIntersection).norm();
+  auto const pointOfIntersection = trajectory.GetPosition(timeOfIntersection);
+  return (trajectory.GetR0() - pointOfIntersection).norm() * 1.0001;
 }

Processes/ObservationPlane/ObservationPlane.h

@@ -29,7 +29,7 @@ namespace corsika::process::observation_plane {
   class ObservationPlane : public corsika::process::ContinuousProcess<ObservationPlane> {
 
   public:
-    ObservationPlane(geometry::Plane const& vObsPlane, std::string const& vFilename);
+    ObservationPlane(geometry::Plane const&, std::string const&, bool = true);
     void Init() {}
 
     corsika::process::EProcessReturn DoContinuous(
@@ -41,8 +41,9 @@ namespace corsika::process::observation_plane {
         corsika::setup::Trajectory const& vTrajectory);
 
   private:
-    geometry::Plane const fObsPlane;
-    std::ofstream fOutputStream;
+    geometry::Plane const plane_;
+    std::ofstream outputStream_;
+    bool const deleteOnHit_;
   };
 } // namespace corsika::process::observation_plane
 

Processes/ObservationPlane/testObservationPlane.cc

@@ -41,7 +41,7 @@ TEST_CASE("ContinuousProcess interface", "[proccesses][observation_plane]") {
   Vector<units::si::SpeedType::dimension_type> vec(rootCS, 0_m / second, 0_m / second,
                                                    -units::constants::c);
   Line line(start, vec);
-  Trajectory<Line> track(line, 10_m / units::constants::c);
+  Trajectory<Line> track(line, 12_m / units::constants::c);
 
   // setup particle stack, and add primary particle
   setup::Stack stack;
@@ -64,14 +64,14 @@ TEST_CASE("ContinuousProcess interface", "[proccesses][observation_plane]") {
 
     Plane const obsPlane(Point(rootCS, {0_m, 0_m, 0_m}),
                          Vector<dimensionless_d>(rootCS, {0., 0., 1.}));
-    ObservationPlane obs(obsPlane, "particles.dat");
+    ObservationPlane obs(obsPlane, "particles.dat", true);
 
     obs.Init();
-    [[maybe_unused]] const LengthType length = obs.MaxStepLength(particle, track);
-    [[maybe_unused]] const process::EProcessReturn ret =
-        obs.DoContinuous(particle, track);
+    const LengthType length = obs.MaxStepLength(particle, track);
+    const process::EProcessReturn ret = obs.DoContinuous(particle, track);
 
-    SECTION("steplength") { REQUIRE(length == 10_m); }
+    REQUIRE(length / 10_m == Approx(1).margin(1e-4));
+    REQUIRE(ret == process::EProcessReturn::eParticleAbsorbed);
 
     /*
     SECTION("horizontal plane") {
@@ -81,16 +81,18 @@ TEST_CASE("ContinuousProcess interface", "[proccesses][observation_plane]") {
     */
   }
 
-  SECTION("inclined plane") {
+  SECTION("inclined plane") {}
+
+  SECTION("transparent plane") {
     Plane const obsPlane(Point(rootCS, {0_m, 0_m, 0_m}),
-                         Vector<dimensionless_d>(rootCS, {1., 1., 0.5}));
-    ObservationPlane obs(obsPlane, "particles.dat");
+                         Vector<dimensionless_d>(rootCS, {0., 0., 1.}));
+    ObservationPlane obs(obsPlane, "particles.dat", false);
 
     obs.Init();
-    [[maybe_unused]] const LengthType length = obs.MaxStepLength(particle, track);
-    [[maybe_unused]] const process::EProcessReturn ret =
-        obs.DoContinuous(particle, track);
+    const LengthType length = obs.MaxStepLength(particle, track);
+    const process::EProcessReturn ret = obs.DoContinuous(particle, track);
 
-    SECTION("steplength") { REQUIRE(length == 12_m); }
+    REQUIRE(length / 10_m == Approx(1).margin(1e-4));
+    REQUIRE(ret == process::EProcessReturn::eOk);
   }
 }