diff --git a/corsika/detail/modules/proposal/Interaction.inl b/corsika/detail/modules/proposal/Interaction.inl
index e81f9984b7b880ddd208556b087df0d0302b3d59..d8916cbe976e77b8865d7d0edaacc6bfc0135e0b 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 0000000000000000000000000000000000000000..17509810e777adb91ed8fd4c878834b23f42b198
--- /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 7b5bd35877b755fe1d751965975c148d42a8f20e..1ebd98015c757bfe451b838f89eea5995d22688e 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 8b28769038504483db2cc3b516d62f162d2dbd10..caf3717b940525ecb19023614d049f81e01a557d 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 cf01b33ca499cf8d1df3d61bcc142f0aee921cdf..94e6d5d594f354da86616d56b4dd9dcec26b4e31 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 fda0a3af6062b3011ee073c30e7fd2438349ae91..ab898b2bdac4fcf70c175a1175421372fa2b1735 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 7df6dcbde54c0d40190da9a377c166d9472c588b..c55584fffe0918b9d0f1d698e1511c255295973d 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 b01d7b608b85d661d7fdce6004a1ef9efcd3cd4b..bd219acda6732cc4d258b12f67c8d01969122188 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 0ff10af164ff51fc3db6cde5bbb7324dea410d24..2a65854a324b2abc897eb98974094bdd83c1d567 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 afa8c71b70b8e2526785312acc38e1c7b684484c..ae41e24571b61e4d63872a691f93099adca1c897 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 e971a279d27edd940098b228d335efdca9d27594..c7762fb38c127a6a1e6c7b7a347f7bfb1a5bd864 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 6bb97291aa2b2d7d3106fcf68a6b887aaf99c7d9..bad43e4e1506d831959072dadf87ff89d312dc46 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);