From 26e1489355d19548c481d23639a4db5678c2f859 Mon Sep 17 00:00:00 2001
From: Maximilian Reininghaus <maximilian.reininghaus@tu-dortmund.de>
Date: Fri, 12 Jun 2020 23:59:39 +0200
Subject: [PATCH] added LongitudinalProfile process

---
 Processes/CMakeLists.txt                      |  1 +
 Processes/LongitudinalProfile/CMakeLists.txt  | 62 +++++++++++++++++
 .../LongitudinalProfile.cc                    | 69 +++++++++++++++++++
 .../LongitudinalProfile/LongitudinalProfile.h | 60 ++++++++++++++++
 4 files changed, 192 insertions(+)
 create mode 100644 Processes/LongitudinalProfile/CMakeLists.txt
 create mode 100644 Processes/LongitudinalProfile/LongitudinalProfile.cc
 create mode 100644 Processes/LongitudinalProfile/LongitudinalProfile.h

diff --git a/Processes/CMakeLists.txt b/Processes/CMakeLists.txt
index 6e93793ec..d9083f299 100644
--- a/Processes/CMakeLists.txt
+++ b/Processes/CMakeLists.txt
@@ -13,6 +13,7 @@ add_subdirectory (UrQMD)
 
 # continuous physics
 add_subdirectory (EnergyLoss)
+add_subdirectory (LongitudinalProfile)
 add_subdirectory (TrackWriter)
 add_subdirectory (ObservationPlane)
 # stack processes
diff --git a/Processes/LongitudinalProfile/CMakeLists.txt b/Processes/LongitudinalProfile/CMakeLists.txt
new file mode 100644
index 000000000..33e0df66d
--- /dev/null
+++ b/Processes/LongitudinalProfile/CMakeLists.txt
@@ -0,0 +1,62 @@
+set (
+  MODEL_SOURCES
+  LongitudinalProfile.cc
+  )
+
+set (
+  MODEL_HEADERS
+  LongitudinalProfile.h
+  )
+
+set (
+  MODEL_NAMESPACE
+  corsika/process/longitudinal_profile
+  )
+
+add_library (ProcessLongitudinalProfile STATIC ${MODEL_SOURCES})
+CORSIKA_COPY_HEADERS_TO_NAMESPACE (ProcessLongitudinalProfile ${MODEL_NAMESPACE} ${MODEL_HEADERS})
+
+set_target_properties (
+  ProcessLongitudinalProfile
+  PROPERTIES
+  VERSION ${PROJECT_VERSION}
+  SOVERSION 1
+#  PUBLIC_HEADER "${MODEL_HEADERS}"
+  )
+
+# target dependencies on other libraries (also the header onlys)
+target_link_libraries (
+  ProcessLongitudinalProfile
+  CORSIKAenvironment
+  CORSIKAunits
+  CORSIKAparticles
+  CORSIKAgeometry
+  CORSIKAsetup
+  )
+
+target_include_directories (
+  ProcessLongitudinalProfile 
+  INTERFACE 
+  $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
+  $<INSTALL_INTERFACE:include/include>
+  )
+
+install (
+  TARGETS ProcessLongitudinalProfile
+  LIBRARY DESTINATION lib
+  ARCHIVE DESTINATION lib
+#  PUBLIC_HEADER DESTINATION include/${MODEL_NAMESPACE}
+  )
+
+
+# --------------------
+# code unit testing
+# CORSIKA_ADD_TEST(testNullModel)
+#target_link_libraries (
+#  testNullModel  ProcessNullModel
+#  CORSIKAsetup
+#  CORSIKAgeometry
+#  CORSIKAunits
+#  CORSIKAthirdparty # for catch2
+#  )
+
diff --git a/Processes/LongitudinalProfile/LongitudinalProfile.cc b/Processes/LongitudinalProfile/LongitudinalProfile.cc
new file mode 100644
index 000000000..5a088b9c5
--- /dev/null
+++ b/Processes/LongitudinalProfile/LongitudinalProfile.cc
@@ -0,0 +1,69 @@
+/*
+ * (c) Copyright 2020 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/longitudinal_profile/LongitudinalProfile.h>
+
+#include <corsika/particles/ParticleProperties.h>
+
+#include <corsika/setup/SetupStack.h>
+#include <corsika/setup/SetupTrajectory.h>
+
+#include <cmath>
+#include <iomanip>
+#include <limits>
+
+using namespace corsika::setup;
+using Particle = Stack::ParticleType;
+using Track = Trajectory;
+
+using namespace corsika::process::longitudinal_profile;
+using namespace corsika::units::si;
+
+LongitudinalProfile::LongitudinalProfile(environment::ShowerAxis const& shower_axis)
+    : shower_axis_{shower_axis}
+    , profiles_{int(shower_axis.maximumX() / dX_) + 1} {}
+
+void LongitudinalProfile::Init() {}
+
+template <>
+corsika::process::EProcessReturn LongitudinalProfile::DoContinuous(Particle const& vP,
+                                                                   Track const& vTrack) {
+  auto const pid = vP.GetPID();
+
+  GrammageType const grammageStart = shower_axis_.projectedX(vTrack.GetPosition(0));
+  GrammageType const grammageEnd = shower_axis_.projectedX(vTrack.GetPosition(1));
+
+  const int binStart = std::ceil(grammageStart / dX_);
+  const int binEnd = std::floor(grammageEnd / dX_);
+
+  for (int b = binStart; b <= binEnd; ++b) {
+    if (pid == particles::Code::MuPlus) {
+      profiles_.at(b)[ProfileIndex::MuPlus]++;
+    } else if (pid == particles::Code::MuMinus) {
+      profiles_.at(b)[ProfileIndex::MuMinus]++;
+    } else if (particles::IsHadron(pid)) {
+      profiles_.at(b)[ProfileIndex::Hadron]++;
+    }
+  }
+
+  return corsika::process::EProcessReturn::eOk;
+}
+
+void LongitudinalProfile::save(std::string const& filename) {
+  std::ofstream f{filename};
+  f << "# X / g·cm¯², mu+, mu-, all hadrons" << std::endl;
+  for (size_t b = 0; b < profiles_.size(); ++b) {
+    f << std::setprecision(5) << std::setw(11) << b * (dX_ / (1_g / 1_cm / 1_cm));
+    for (auto const& N : profiles_.at(b)) {
+      f << std::setw(width_) << std::setprecision(precision_) << std::scientific << N;
+    }
+    f << std::endl;
+  }
+}
diff --git a/Processes/LongitudinalProfile/LongitudinalProfile.h b/Processes/LongitudinalProfile/LongitudinalProfile.h
new file mode 100644
index 000000000..68de7e0a0
--- /dev/null
+++ b/Processes/LongitudinalProfile/LongitudinalProfile.h
@@ -0,0 +1,60 @@
+/*
+ * (c) Copyright 2020 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_track_writer_LongitudinalProfile_h_
+#define _Processes_track_writer_LongitudinalProfile_h_
+
+#include <corsika/environment/ShowerAxis.h>
+#include <corsika/process/ContinuousProcess.h>
+#include <corsika/units/PhysicalUnits.h>
+
+#include <array>
+#include <fstream>
+#include <limits>
+#include <string>
+
+namespace corsika::process::longitudinal_profile {
+
+  class LongitudinalProfile
+      : public corsika::process::ContinuousProcess<LongitudinalProfile> {
+
+  public:
+    LongitudinalProfile(environment::ShowerAxis const&);
+
+    void Init();
+
+    template <typename Particle, typename Track>
+    corsika::process::EProcessReturn DoContinuous(Particle const&, Track const&);
+
+    template <typename Particle, typename Track>
+    corsika::units::si::LengthType MaxStepLength(Particle const&, Track const&) {
+      return units::si::meter * std::numeric_limits<double>::infinity();
+    }
+
+    void save(std::string const&);
+
+  private:
+    units::si::GrammageType const dX_ = std::invoke([]() {
+      using namespace units::si;
+      return 10_g / square(1_cm);
+    }); // profile binning
+
+    environment::ShowerAxis const& shower_axis_;
+    using ProfileEntry = std::array<uint32_t, 3>;
+    enum ProfileIndex { MuPlus = 0, MuMinus = 1, Hadron = 2 };
+    std::vector<ProfileEntry> profiles_; // longitudinal profile
+
+    static int const width_ = 14;
+    static int const precision_ = 6;
+  };
+
+} // namespace corsika::process::longitudinal_profile
+
+#endif
-- 
GitLab