From b18d0dc065aa1c8047ba31fcc9dade2d585cc79d Mon Sep 17 00:00:00 2001
From: Maximilian Reininghaus <maximilian.reininghaus@kit.edu>
Date: Tue, 28 May 2019 19:30:08 -0300
Subject: [PATCH] added ObservationPlane

---
 Processes/CMakeLists.txt                      |  1 +
 Processes/ObservationPlane/CMakeLists.txt     | 44 ++++++++++++++++
 .../ObservationPlane/ObservationPlane.cc      | 50 +++++++++++++++++++
 Processes/ObservationPlane/ObservationPlane.h | 49 ++++++++++++++++++
 4 files changed, 144 insertions(+)
 create mode 100644 Processes/ObservationPlane/CMakeLists.txt
 create mode 100644 Processes/ObservationPlane/ObservationPlane.cc
 create mode 100644 Processes/ObservationPlane/ObservationPlane.h

diff --git a/Processes/CMakeLists.txt b/Processes/CMakeLists.txt
index 313dd46f8..758b430e6 100644
--- a/Processes/CMakeLists.txt
+++ b/Processes/CMakeLists.txt
@@ -14,6 +14,7 @@ add_subdirectory (SwitchProcess)
 # continuous physics
 add_subdirectory (EnergyLoss)
 add_subdirectory (TrackWriter)
+add_subdirectory (ObservationPlane)
 # stack processes
 add_subdirectory (StackInspector)
 # secondaries process
diff --git a/Processes/ObservationPlane/CMakeLists.txt b/Processes/ObservationPlane/CMakeLists.txt
new file mode 100644
index 000000000..cb656d312
--- /dev/null
+++ b/Processes/ObservationPlane/CMakeLists.txt
@@ -0,0 +1,44 @@
+set (
+  MODEL_HEADERS
+  ObservationPlane.h
+  )
+
+set (
+  MODEL_SOURCES
+  ObservationPlane.cc
+  )
+
+set (
+  MODEL_NAMESPACE
+  corsika/process/observation_plane
+  )
+
+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_link_libraries (
+  ProcessObservationPlane
+  CORSIKAgeometry
+  CORSIKAprocesssequence
+  )
+
+target_include_directories (
+  ProcessObservationPlane 
+  INTERFACE 
+  $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
+  $<INSTALL_INTERFACE:include/include>
+  )
+
+install (FILES ${MODEL_HEADERS} DESTINATION include/${MODEL_NAMESPACE})
+
+# --------------------
+# code unit testing
+#CORSIKA_ADD_TEST(testObservationPlane)
+#target_link_libraries (
+#  testObservationPlane
+#  ProcessObservationPlane
+#  CORSIKAstackinterface
+#  CORSIKAthirdparty # for catch2
+#)
+
diff --git a/Processes/ObservationPlane/ObservationPlane.cc b/Processes/ObservationPlane/ObservationPlane.cc
new file mode 100644
index 000000000..69c2f1382
--- /dev/null
+++ b/Processes/ObservationPlane/ObservationPlane.cc
@@ -0,0 +1,50 @@
+/*
+ * (c) Copyright 2019 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * 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.
+ */
+
+#include <corsika/process/observation_plane/ObservationPlane.h>
+
+#include <fstream>
+
+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;
+}
+
+corsika::process::EProcessReturn ObservationPlane::DoContinuous(
+    setup::Stack::ParticleType const& vParticle, setup::Trajectory const& vTrajectory) {
+
+  if (fObsPlane.IsAbove(vTrajectory.GetPosition(1.0001))) {
+    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
+                << std::endl;
+
+  return process::EProcessReturn::eParticleAbsorbed;
+}
+
+LengthType ObservationPlane::MaxStepLength(setup::Stack::ParticleType const&,
+                                           setup::Trajectory const& vTrajectory) {
+  if (!fObsPlane.IsAbove(vTrajectory.GetR0())) {
+    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();
+}
diff --git a/Processes/ObservationPlane/ObservationPlane.h b/Processes/ObservationPlane/ObservationPlane.h
new file mode 100644
index 000000000..9ecec9f63
--- /dev/null
+++ b/Processes/ObservationPlane/ObservationPlane.h
@@ -0,0 +1,49 @@
+/*
+ * (c) Copyright 2019 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * 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.
+ */
+
+#ifndef _Processes_ObservationPlane_h_
+#define _Processes_ObservationPlane_h_
+
+#include <corsika/geometry/Plane.h>
+#include <corsika/process/ContinuousProcess.h>
+#include <corsika/setup/SetupStack.h>
+#include <corsika/setup/SetupTrajectory.h>
+#include <corsika/units/PhysicalUnits.h>
+
+#include <fstream>
+
+namespace corsika::process::observation_plane {
+
+  /**
+   * The ObservationPlane writes PDG codes, energies, and distances of particles to the
+   * central point of the plane into its output file. The particles are considered
+   * "absorbed" afterwards.
+   */
+  class ObservationPlane : public corsika::process::ContinuousProcess<ObservationPlane> {
+
+  public:
+    ObservationPlane(geometry::Plane const& vObsPlane, std::string const& vFilename);
+    void Init() {}
+
+    corsika::process::EProcessReturn DoContinuous(
+        corsika::setup::Stack::ParticleType const& vParticle,
+        corsika::setup::Trajectory const& vTrajectory);
+
+    corsika::units::si::LengthType MaxStepLength(
+        corsika::setup::Stack::ParticleType const&,
+        corsika::setup::Trajectory const& vTrajectory);
+
+  private:
+    geometry::Plane const fObsPlane;
+    std::ofstream fOutputStream;
+  };
+} // namespace corsika::process::observation_plane
+
+#endif
-- 
GitLab