From 2dd2d97a12fdf0e4cca6fe4f11f88857041f5ae9 Mon Sep 17 00:00:00 2001
From: ralfulrich <ralf.ulrich@kit.edu>
Date: Wed, 1 Dec 2021 11:35:13 +0100
Subject: [PATCH] updated examples

---
 .../detail/modules/proposal/Interaction.inl   |  8 ++---
 .../core/EnergyMomentumOperations.hpp         | 36 +++++++++++++++++++
 corsika/framework/units/quantity.hpp          | 18 +++++-----
 examples/boundary_example.cpp                 |  9 +++--
 examples/cascade_example.cpp                  |  9 +++--
 examples/cascade_proton_example.cpp           | 10 ++++--
 examples/corsika.cpp                          | 19 +++++-----
 examples/em_shower.cpp                        | 19 ++++++----
 examples/hybrid_MC.cpp                        | 19 +++++-----
 examples/mars.cpp                             | 15 +++++---
 examples/stopping_power.cpp                   |  5 ++-
 examples/vertical_EAS.cpp                     | 16 ++++++---
 12 files changed, 126 insertions(+), 57 deletions(-)
 create mode 100644 corsika/framework/core/EnergyMomentumOperations.hpp

diff --git a/corsika/detail/modules/proposal/Interaction.inl b/corsika/detail/modules/proposal/Interaction.inl
index e81f9984b..d8916cbe9 100644
--- a/corsika/detail/modules/proposal/Interaction.inl
+++ b/corsika/detail/modules/proposal/Interaction.inl
@@ -93,12 +93,10 @@ namespace corsika::proposal {
       for (auto& s : sec) {
         auto E = s.energy * 1_MeV;
         auto vecProposal = s.direction;
-        auto vec = QuantityVector(vecProposal.GetX() * E, vecProposal.GetY() * E,
-                                  vecProposal.GetZ() * E);
-        auto p = MomentumVector(labCS, vec);
+        auto dir = DirectionVector(
+            labCS, {vecProposal.GetX(), vecProposal.GetY(), vecProposal.GetZ()});
         auto sec_code = convert_from_PDG(static_cast<PDGCode>(s.type));
-        view.addSecondary(
-            std::make_tuple(sec_code, p, projectile.getPosition(), projectile.getTime()));
+        view.addSecondary(std::make_tuple(sec_code, E - get_mass(sec_code), dir));
       }
     }
     return ProcessReturn::Ok;
diff --git a/corsika/framework/core/EnergyMomentumOperations.hpp b/corsika/framework/core/EnergyMomentumOperations.hpp
new file mode 100644
index 000000000..17509810e
--- /dev/null
+++ b/corsika/framework/core/EnergyMomentumOperations.hpp
@@ -0,0 +1,36 @@
+/*
+ * (c) Copyright 2020 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/core/PhysicalUnits.hpp>
+#include <corsika/framework/core/PhysicalGeometry.hpp>
+#include <corsika/framework/core/PhysicalConstants.hpp>
+
+/**
+ * \file EnergyMomentumOperations.hpp
+ *
+ * Relativic energy momentum calculations.
+ */
+
+namespace corsika {
+
+  auto constexpr get_total_energy_sqr(HEPMomentumType const p, HEPMassType const m) {
+    return p * p + m * m;
+  }
+
+  HEPEnergyType constexpr get_total_energy(HEPMomentumType const p, HEPMassType const m) {
+    return sqrt(get_total_energy_sqr(p, m));
+  }
+
+  HEPEnergyType constexpr get_kinetic_energy(HEPMomentumType const p,
+                                             HEPMassType const m) {
+    return get_total_energy(p, m) - m;
+  }
+
+} // namespace corsika
\ No newline at end of file
diff --git a/corsika/framework/units/quantity.hpp b/corsika/framework/units/quantity.hpp
index 7b5bd3587..1ebd98015 100644
--- a/corsika/framework/units/quantity.hpp
+++ b/corsika/framework/units/quantity.hpp
@@ -454,7 +454,7 @@ namespace phys {
       // powers and roots
 
       template <int N, typename D, typename X>
-      friend detail::Power<D, N, X> nth_power(quantity<D, X> const& x);
+      friend constexpr detail::Power<D, N, X> nth_power(quantity<D, X> const& x);
 
       template <typename D, typename X>
       friend constexpr detail::Power<D, 2, X> square(quantity<D, X> const& x);
@@ -463,13 +463,13 @@ namespace phys {
       friend constexpr detail::Power<D, 3, X> cube(quantity<D, X> const& x);
 
       template <int N, typename D, typename X>
-      friend detail::Root<D, N, X> nth_root(quantity<D, X> const& x);
+      friend detail::Root<D, N, X> constexpr nth_root(quantity<D, X> const& x);
 
       template <typename D, typename X>
-      friend detail::Root<D, 2, X> sqrt(quantity<D, X> const& x);
+      friend detail::Root<D, 2, X> constexpr sqrt(quantity<D, X> const& x);
 
       template <typename D, typename X>
-      friend detail::Root<D, 3, X> cbrt(quantity<D, X> const& x);
+      friend detail::Root<D, 3, X> constexpr cbrt(quantity<D, X> const& x);
 
       // comparison
 
@@ -629,7 +629,7 @@ namespace phys {
     /// absolute value.
 
     template <typename D, typename X>
-    constexpr quantity<D, X> abs(quantity<D, X> const& x) {
+    quantity<D, X> constexpr abs(quantity<D, X> const& x) {
       return quantity<D, X>(std::abs(x.m_value));
     }
 
@@ -638,7 +638,7 @@ namespace phys {
     /// N-th power.
 
     template <int N, typename D, typename X>
-    detail::Power<D, N, X> nth_power(quantity<D, X> const& x) {
+    detail::Power<D, N, X> constexpr nth_power(quantity<D, X> const& x) {
       return detail::Power<D, N, X>(std::pow(x.m_value, X(N)));
     }
 
@@ -663,7 +663,7 @@ namespace phys {
     /// n-th root.
 
     template <int N, typename D, typename X>
-    detail::Root<D, N, X> nth_root(quantity<D, X> const& x) {
+    detail::Root<D, N, X> constexpr nth_root(quantity<D, X> const& x) {
       static_assert(detail::root<D, N, X>::all_even_multiples,
                     "root result dimensions must be integral");
 
@@ -675,7 +675,7 @@ namespace phys {
     /// square root.
 
     template <typename D, typename X>
-    detail::Root<D, 2, X> sqrt(quantity<D, X> const& x) {
+    detail::Root<D, 2, X> constexpr sqrt(quantity<D, X> const& x) {
       static_assert(detail::root<D, 2, X>::all_even_multiples,
                     "root result dimensions must be integral");
 
@@ -685,7 +685,7 @@ namespace phys {
     /// cubic root.
 
     template <typename D, typename X>
-    detail::Root<D, 3, X> cbrt(quantity<D, X> const& x) {
+    detail::Root<D, 3, X> constexpr cbrt(quantity<D, X> const& x) {
       static_assert(detail::root<D, 3, X>::all_even_multiples,
                     "root result dimensions must be integral");
 
diff --git a/examples/boundary_example.cpp b/examples/boundary_example.cpp
index 8b2876903..caf3717b9 100644
--- a/examples/boundary_example.cpp
+++ b/examples/boundary_example.cpp
@@ -6,13 +6,14 @@
  * the license.
  */
 
-#include <corsika/framework/core/Cascade.hpp>
 #include <corsika/framework/process/ProcessSequence.hpp>
 #include <corsika/framework/geometry/Sphere.hpp>
-#include <corsika/framework/core/PhysicalUnits.hpp>
 #include <corsika/framework/random/RNGManager.hpp>
 #include <corsika/framework/utility/CorsikaFenv.hpp>
+#include <corsika/framework/core/Cascade.hpp>
 #include <corsika/framework/core/Logging.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/core/EnergyMomentumOperations.hpp>
 
 #include <corsika/output/OutputManager.hpp>
 
@@ -160,7 +161,9 @@ int main() {
         beamCode, theta, phi, plab.getComponents() / 1_GeV);
     // shoot particles from inside target out
     Point pos(rootCS, 0_m, 0_m, 0_m);
-    stack.addParticle(std::make_tuple(beamCode, plab, pos, 0_ns));
+    stack.addParticle(
+        std::make_tuple(beamCode, get_kinetic_energy(plab.getNorm(), get_mass(beamCode)),
+                        plab.normalized(), pos, 0_ns));
   }
 
   // define air shower object, run simulation
diff --git a/examples/cascade_example.cpp b/examples/cascade_example.cpp
index cf01b33ca..94e6d5d59 100644
--- a/examples/cascade_example.cpp
+++ b/examples/cascade_example.cpp
@@ -6,14 +6,15 @@
  * the license.
  */
 
-#include <corsika/framework/core/Cascade.hpp>
 #include <corsika/framework/process/ProcessSequence.hpp>
 #include <corsika/framework/process/SwitchProcessSequence.hpp>
+#include <corsika/framework/core/Cascade.hpp>
 #include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/core/Logging.hpp>
+#include <corsika/framework/core/EnergyMomentumOperations.hpp>
 #include <corsika/framework/random/RNGManager.hpp>
 #include <corsika/framework/geometry/Sphere.hpp>
 #include <corsika/framework/utility/CorsikaFenv.hpp>
-#include <corsika/framework/core/Logging.hpp>
 
 #include <corsika/output/OutputManager.hpp>
 
@@ -125,7 +126,9 @@ int main() {
     cout << "input particle: " << beamCode << endl;
     cout << "input angles: theta=" << theta << " phi=" << phi << endl;
     cout << "input momentum: " << plab.getComponents() / 1_GeV << endl;
-    stack.addParticle(std::make_tuple(beamCode, plab, injectionPos, 0_ns));
+    stack.addParticle(
+        std::make_tuple(beamCode, get_kinetic_energy(plab.getNorm(), get_mass(beamCode)),
+                        plab.normalized(), injectionPos, 0_ns));
   }
 
   // setup processes, decays and interactions
diff --git a/examples/cascade_proton_example.cpp b/examples/cascade_proton_example.cpp
index fda0a3af6..ab898b2bd 100644
--- a/examples/cascade_proton_example.cpp
+++ b/examples/cascade_proton_example.cpp
@@ -6,13 +6,15 @@
  * the license.
  */
 
-#include <corsika/framework/core/Cascade.hpp>
 #include <corsika/framework/process/ProcessSequence.hpp>
+#include <corsika/framework/process/SwitchProcessSequence.hpp>
+#include <corsika/framework/core/Cascade.hpp>
 #include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/core/Logging.hpp>
+#include <corsika/framework/core/EnergyMomentumOperations.hpp>
 #include <corsika/framework/random/RNGManager.hpp>
 #include <corsika/framework/geometry/Sphere.hpp>
 #include <corsika/framework/utility/CorsikaFenv.hpp>
-#include <corsika/framework/core/Logging.hpp>
 
 #include <corsika/output/OutputManager.hpp>
 
@@ -109,7 +111,9 @@ int main() {
     cout << "input particle: " << beamCode << endl;
     cout << "input angles: theta=" << theta << " phi=" << phi << endl;
     cout << "input momentum: " << plab.getComponents() / 1_GeV << endl;
-    stack.addParticle(std::make_tuple(beamCode, plab, injectionPos, 0_ns));
+    stack.addParticle(
+        std::make_tuple(beamCode, get_kinetic_energy(plab.getNorm(), get_mass(beamCode)),
+                        plab.normalized(), injectionPos, 0_ns));
   }
 
   // setup processes, decays and interactions
diff --git a/examples/corsika.cpp b/examples/corsika.cpp
index 7df6dcbde..c55584fff 100644
--- a/examples/corsika.cpp
+++ b/examples/corsika.cpp
@@ -12,18 +12,19 @@
 // to include it first...
 #include <corsika/framework/process/InteractionCounter.hpp>
 /* clang-format on */
-#include <corsika/framework/geometry/Plane.hpp>
-#include <corsika/framework/geometry/Sphere.hpp>
-#include <corsika/framework/core/Logging.hpp>
-#include <corsika/framework/utility/SaveBoostHistogram.hpp>
 #include <corsika/framework/process/ProcessSequence.hpp>
 #include <corsika/framework/process/SwitchProcessSequence.hpp>
 #include <corsika/framework/process/InteractionCounter.hpp>
-#include <corsika/framework/random/RNGManager.hpp>
+#include <corsika/framework/geometry/Plane.hpp>
+#include <corsika/framework/geometry/Sphere.hpp>
+#include <corsika/framework/geometry/PhysicalGeometry.hpp>
+#include <corsika/framework/core/Logging.hpp>
 #include <corsika/framework/core/PhysicalUnits.hpp>
-#include <corsika/framework/utility/CorsikaFenv.hpp>
 #include <corsika/framework/core/Cascade.hpp>
-#include <corsika/framework/geometry/PhysicalGeometry.hpp>
+#include <corsika/framework/core/EnergyMomentumOperations.hpp>
+#include <corsika/framework/utility/SaveBoostHistogram.hpp>
+#include <corsika/framework/utility/CorsikaFenv.hpp>
+#include <corsika/framework/random/RNGManager.hpp>
 
 #include <corsika/output/OutputManager.hpp>
 #include <corsika/output/NoOutput.hpp>
@@ -395,7 +396,9 @@ int main(int argc, char** argv) {
     stack.clear();
 
     // add the desired particle to the stack
-    stack.addParticle(std::make_tuple(beamCode, plab, injectionPos, 0_ns));
+    stack.addParticle(
+        std::make_tuple(beamCode, get_kinetic_energy(plab.getNorm(), get_mass(beamCode)),
+                        plab.normalized(), injectionPos, 0_ns));
 
     // if we want to fix the first location of the shower
     if (force_interaction) {
diff --git a/examples/em_shower.cpp b/examples/em_shower.cpp
index b01d7b608..bd219acda 100644
--- a/examples/em_shower.cpp
+++ b/examples/em_shower.cpp
@@ -6,17 +6,20 @@
  * the license.
  */
 
+#include <corsika/framework/process/ProcessSequence.hpp>
+#include <corsika/framework/process/SwitchProcessSequence.hpp>
+#include <corsika/framework/process/InteractionCounter.hpp>
 #include <corsika/framework/core/Cascade.hpp>
-#include <corsika/framework/utility/SaveBoostHistogram.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/core/Logging.hpp>
+#include <corsika/framework/core/EnergyMomentumOperations.hpp>
+#include <corsika/framework/random/RNGManager.hpp>
+#include <corsika/framework/geometry/Sphere.hpp>
 #include <corsika/framework/geometry/Plane.hpp>
 #include <corsika/framework/geometry/Sphere.hpp>
 #include <corsika/framework/geometry/PhysicalGeometry.hpp>
-#include <corsika/framework/process/ProcessSequence.hpp>
-#include <corsika/framework/random/RNGManager.hpp>
-#include <corsika/framework/core/PhysicalUnits.hpp>
 #include <corsika/framework/utility/CorsikaFenv.hpp>
-#include <corsika/framework/process/InteractionCounter.hpp>
-#include <corsika/framework/core/Logging.hpp>
+#include <corsika/framework/utility/SaveBoostHistogram.hpp>
 
 #include <corsika/output/OutputManager.hpp>
 
@@ -130,7 +133,9 @@ int main(int argc, char** argv) {
 
   std::cout << "point of injection: " << injectionPos.getCoordinates() << std::endl;
 
-  stack.addParticle(std::make_tuple(beamCode, plab, injectionPos, 0_ns));
+  stack.addParticle(
+      std::make_tuple(beamCode, get_kinetic_energy(plab.getNorm(), get_mass(beamCode)),
+                      plab.normalized(), injectionPos, 0_ns));
 
   std::cout << "shower axis length: " << (showerCore - injectionPos).getNorm() * 1.02
             << std::endl;
diff --git a/examples/hybrid_MC.cpp b/examples/hybrid_MC.cpp
index 0ff10af16..2a65854a3 100644
--- a/examples/hybrid_MC.cpp
+++ b/examples/hybrid_MC.cpp
@@ -12,18 +12,19 @@
 // to include it first...
 #include <corsika/framework/process/InteractionCounter.hpp>
 /* clang-format on */
-#include <corsika/framework/geometry/Plane.hpp>
-#include <corsika/framework/utility/SaveBoostHistogram.hpp>
-#include <corsika/framework/geometry/Sphere.hpp>
-#include <corsika/framework/core/Logging.hpp>
 #include <corsika/framework/process/ProcessSequence.hpp>
 #include <corsika/framework/process/SwitchProcessSequence.hpp>
 #include <corsika/framework/process/InteractionCounter.hpp>
-#include <corsika/framework/random/RNGManager.hpp>
-#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/geometry/Plane.hpp>
+#include <corsika/framework/geometry/Sphere.hpp>
+#include <corsika/framework/geometry/PhysicalGeometry.hpp>
+#include <corsika/framework/utility/SaveBoostHistogram.hpp>
 #include <corsika/framework/utility/CorsikaFenv.hpp>
+#include <corsika/framework/core/Logging.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
 #include <corsika/framework/core/Cascade.hpp>
-#include <corsika/framework/geometry/PhysicalGeometry.hpp>
+#include <corsika/framework/core/EnergyMomentumOperations.hpp>
+#include <corsika/framework/random/RNGManager.hpp>
 
 #include <corsika/output/OutputManager.hpp>
 
@@ -183,7 +184,9 @@ int main(int argc, char** argv) {
 
   std::cout << "point of injection: " << injectionPos.getCoordinates() << std::endl;
 
-  stack.addParticle(std::make_tuple(Code::Proton, plab, injectionPos, 0_ns));
+  stack.addParticle(std::make_tuple(
+      Code::Proton, get_kinetic_energy(plab.getNorm(), get_mass(beamCode)),
+      plab.normalized(), injectionPos, 0_ns));
 
   std::cout << "shower axis length: " << (showerCore - injectionPos).getNorm() * 1.02
             << std::endl;
diff --git a/examples/mars.cpp b/examples/mars.cpp
index afa8c71b7..ae41e2457 100644
--- a/examples/mars.cpp
+++ b/examples/mars.cpp
@@ -14,16 +14,19 @@
 /* clang-format on */
 #include <corsika/framework/geometry/Plane.hpp>
 #include <corsika/framework/geometry/Sphere.hpp>
+#include <corsika/framework/geometry/PhysicalGeometry.hpp>
+
 #include <corsika/framework/core/Logging.hpp>
+#include <corsika/framework/core/EnergyMomentumOperations.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/core/Cascade.hpp>
+
 #include <corsika/framework/utility/SaveBoostHistogram.hpp>
+#include <corsika/framework/utility/CorsikaFenv.hpp>
 #include <corsika/framework/process/ProcessSequence.hpp>
 #include <corsika/framework/process/SwitchProcessSequence.hpp>
 #include <corsika/framework/process/InteractionCounter.hpp>
 #include <corsika/framework/random/RNGManager.hpp>
-#include <corsika/framework/core/PhysicalUnits.hpp>
-#include <corsika/framework/utility/CorsikaFenv.hpp>
-#include <corsika/framework/core/Cascade.hpp>
-#include <corsika/framework/geometry/PhysicalGeometry.hpp>
 
 #include <corsika/output/OutputManager.hpp>
 #include <corsika/output/NoOutput.hpp>
@@ -418,7 +421,9 @@ int main(int argc, char** argv) {
     stack.clear();
 
     // add the desired particle to the stack
-    stack.addParticle(std::make_tuple(beamCode, plab, injectionPos, 0_ns));
+    stack.addParticle(
+        std::make_tuple(beamCode, get_kinetic_energy(plab.getNorm(), get_mass(beamCode)),
+                        plab.normalized(), injectionPos, 0_ns));
 
     // run the shower
     EAS.run();
diff --git a/examples/stopping_power.cpp b/examples/stopping_power.cpp
index e971a279d..c7762fb38 100644
--- a/examples/stopping_power.cpp
+++ b/examples/stopping_power.cpp
@@ -15,6 +15,7 @@
 #include <corsika/modules/BetheBlochPDG.hpp>
 #include <corsika/setup/SetupStack.hpp>
 #include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/core/EnergyMomentumOperations.hpp>
 #include <corsika/framework/utility/CorsikaFenv.hpp>
 
 #include <fstream>
@@ -75,7 +76,9 @@ int main() {
         momentumComponents(theta / 180. * constants::pi, phi / 180. * constants::pi, P0);
     auto plab = MomentumVector(rootCS, {px, py, pz});
 
-    stack.addParticle(std::make_tuple(beamCode, plab, injectionPos, 0_ns));
+    stack.addParticle(
+        std::make_tuple(beamCode, get_kinetic_energy(plab.getNorm(), get_mass(beamCode)),
+                        plab.normalized(), injectionPos, 0_ns));
 
     auto const p = stack.getNextParticle();
     HEPEnergyType dE = eLoss.getTotalEnergyLoss(p, 1_g / square(1_cm));
diff --git a/examples/vertical_EAS.cpp b/examples/vertical_EAS.cpp
index 6bb97291a..bad43e4e1 100644
--- a/examples/vertical_EAS.cpp
+++ b/examples/vertical_EAS.cpp
@@ -16,16 +16,20 @@
 /* clang-format on */
 #include <corsika/framework/geometry/Plane.hpp>
 #include <corsika/framework/geometry/Sphere.hpp>
+#include <corsika/framework/geometry/PhysicalGeometry.hpp>
+
 #include <corsika/framework/core/Logging.hpp>
+#include <corsika/framework/core/EnergyMomentumOperations.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/core/Cascade.hpp>
+
 #include <corsika/framework/utility/SaveBoostHistogram.hpp>
+#include <corsika/framework/utility/CorsikaFenv.hpp>
+
 #include <corsika/framework/process/ProcessSequence.hpp>
 #include <corsika/framework/process/SwitchProcessSequence.hpp>
 #include <corsika/framework/process/InteractionCounter.hpp>
 #include <corsika/framework/random/RNGManager.hpp>
-#include <corsika/framework/core/PhysicalUnits.hpp>
-#include <corsika/framework/utility/CorsikaFenv.hpp>
-#include <corsika/framework/core/Cascade.hpp>
-#include <corsika/framework/geometry/PhysicalGeometry.hpp>
 
 #include <corsika/output/OutputManager.hpp>
 #include <corsika/output/NoOutput.hpp>
@@ -251,7 +255,9 @@ int main(int argc, char** argv) {
   setup::Stack stack;
   stack.clear();
 
-  stack.addParticle(std::make_tuple(beamCode, plab, injectionPos, 0_ns));
+  stack.addParticle(
+      std::make_tuple(beamCode, get_kinetic_energy(plab.getNorm(), get_mass(beamCode)),
+                      plab.normalized(), injectionPos, 0_ns));
 
   BetheBlochPDG emContinuous(showerAxis);
 
-- 
GitLab