From d0b193d2755b40645a844af4d8d53b0632c59318 Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Wed, 3 Jun 2020 18:52:58 +0100
Subject: [PATCH 01/31] shift energies of particles from sibyll so that mass is
 consistent, implemented after boost to lab frame

---
 Processes/Pythia/Decay.cc       |  4 ++--
 Processes/Sibyll/Interaction.cc | 38 +++++++++++++++++++++++++++++----
 2 files changed, 36 insertions(+), 6 deletions(-)

diff --git a/Processes/Pythia/Decay.cc b/Processes/Pythia/Decay.cc
index 27890d28b..a17678d18 100644
--- a/Processes/Pythia/Decay.cc
+++ b/Processes/Pythia/Decay.cc
@@ -57,7 +57,7 @@ namespace corsika::process::pythia {
     fPythia.readString("Next:numberShowEvent = 0");
 
     fPythia.readString("Print:quiet = on");
-    fPythia.readString("Check:particleData = 1");
+    fPythia.readString("Check:particleData = 0");
 
     /*
        switching off event check in pythia is needed to allow decays that are off-shell
@@ -65,7 +65,7 @@ namespace corsika::process::pythia {
        the consistency of particle masses between event generators is an unsolved issues
     */
     cout << "Pythia::Init: switching off event checking in pythia.." << endl;
-    fPythia.readString("Check:event = 0");
+    fPythia.readString("Check:event = 1");
 
     fPythia.readString("ProcessLevel:all = off");
     fPythia.readString("ProcessLevel:resonanceDecays = off");
diff --git a/Processes/Sibyll/Interaction.cc b/Processes/Sibyll/Interaction.cc
index 916ea2ec0..e7d610dc5 100644
--- a/Processes/Sibyll/Interaction.cc
+++ b/Processes/Sibyll/Interaction.cc
@@ -316,20 +316,50 @@ namespace corsika::process::sibyll {
           HEPEnergyType const eCoM = psib.GetEnergy();
           auto const Plab = boost.fromCoM(FourVector(eCoM, pCoM));
 
+          auto const pid = process::sibyll::ConvertFromSibyll(psib.GetPID());
+          // check if on-shell in corsika
+          auto const m_kinetic = Plab.GetNorm();
+          auto const m_corsika = particles::GetMass(pid);
+          auto const e_corsika = Plab.GetTimeLikeComponent();
+          auto const m_sibyll = corsika::process::sibyll::GetSibyllMass(pid);
+          auto const m_err = abs(m_kinetic - m_corsika) / m_corsika;
+          if (m_err > 1.e-5) {
+            const HEPEnergyType e_shift_corsika = sqrt(
+                Plab.GetSpaceLikeComponents().GetSquaredNorm() + m_corsika * m_corsika);
+            auto const e_shift_relative = (e_shift_corsika / e_corsika - 1) * 100;
+            // warn about percent level shifts in particle energy
+            if (abs(e_shift_relative) > 1) {
+              std::cout << "Sibyll::Interaction: shifted particle energy by "
+                        << e_shift_relative << " %" << std::endl;
+
+              std::cout << "shift particle mass for " << pid << std::endl
+                        << "corsika mass (GeV): " << m_corsika / 1_GeV << std::endl
+                        << "kinetic mass (GeV): " << m_kinetic / 1_GeV << std::endl
+                        << "sibyll  mass (GeV): " << m_sibyll / 1_GeV << std::endl
+                        << "(m_kin-m_cor)/en: " << m_err << std::endl;
+            }
+            Plab.GetTimeLikeComponent() = e_shift_corsika;
+          }
+
           // add to corsika stack
           auto pnew = vP.AddSecondary(
               tuple<particles::Code, units::si::HEPEnergyType, stack::MomentumVector,
                     geometry::Point, units::si::TimeType>{
-                  process::sibyll::ConvertFromSibyll(psib.GetPID()),
-                  Plab.GetTimeLikeComponent(), Plab.GetSpaceLikeComponents(), pOrig,
+                  pid, Plab.GetTimeLikeComponent(), Plab.GetSpaceLikeComponents(), pOrig,
                   tOrig});
 
           Plab_final += pnew.GetMomentum();
           Elab_final += pnew.GetEnergy();
           Ecm_final += psib.GetEnergy();
         }
-        cout << "conservation (all GeV): Ecm_final=" << Ecm_final / 1_GeV << endl
-             << "Elab_final=" << Elab_final / 1_GeV
+        cout << "conservation (all GeV):" << endl
+             << "Ecm_initial=" << Ecm / 1_GeV << " Ecm_final=" << Ecm_final / 1_GeV
+             << endl
+             << "Elab_initial=" << eProjectileLab / 1_GeV
+             << " Elab_final=" << Elab_final / 1_GeV
+             << " diff (%)=" << (Elab_final / eProjectileLab / get_nwounded() - 1) * 100
+             << endl
+             << "Plab_initial=" << (pProjectileLab / 1_GeV).GetComponents()
              << ", Plab_final=" << (Plab_final / 1_GeV).GetComponents() << endl;
       }
     }
-- 
GitLab


From 08e9579daf6138ef3b549f6805ef2142f1603d8c Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Wed, 3 Jun 2020 23:17:39 +0100
Subject: [PATCH 02/31] shift particle masses in sibyll interface

---
 Processes/Sibyll/Interaction.cc | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Processes/Sibyll/Interaction.cc b/Processes/Sibyll/Interaction.cc
index e7d610dc5..0a4c7d59d 100644
--- a/Processes/Sibyll/Interaction.cc
+++ b/Processes/Sibyll/Interaction.cc
@@ -323,7 +323,7 @@ namespace corsika::process::sibyll {
           auto const e_corsika = Plab.GetTimeLikeComponent();
           auto const m_sibyll = corsika::process::sibyll::GetSibyllMass(pid);
           auto const m_err = abs(m_kinetic - m_corsika) / m_corsika;
-          if (m_err > 1.e-5) {
+          if (m_err > 1.e-5 && false) {
             const HEPEnergyType e_shift_corsika = sqrt(
                 Plab.GetSpaceLikeComponents().GetSquaredNorm() + m_corsika * m_corsika);
             auto const e_shift_relative = (e_shift_corsika / e_corsika - 1) * 100;
-- 
GitLab


From 20bf7cd5e4a85886ea254fda6ad39624047878ee Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Wed, 3 Jun 2020 23:18:40 +0100
Subject: [PATCH 03/31] added process to check particle masses and shift
 on-shell if necessary

---
 Documentation/Examples/CMakeLists.txt      |  1 +
 Documentation/Examples/vertical_EAS.cc     |  5 +-
 Processes/CMakeLists.txt                   |  3 +-
 Processes/OnShellCheck/CMakeLists.txt      | 67 ++++++++++++++++
 Processes/OnShellCheck/OnShellCheck.cc     | 71 +++++++++++++++++
 Processes/OnShellCheck/OnShellCheck.h      | 39 ++++++++++
 Processes/OnShellCheck/ParticleCut.h       | 55 +++++++++++++
 Processes/OnShellCheck/testOnShellCheck.cc | 91 ++++++++++++++++++++++
 8 files changed, 330 insertions(+), 2 deletions(-)
 create mode 100644 Processes/OnShellCheck/CMakeLists.txt
 create mode 100644 Processes/OnShellCheck/OnShellCheck.cc
 create mode 100644 Processes/OnShellCheck/OnShellCheck.h
 create mode 100644 Processes/OnShellCheck/ParticleCut.h
 create mode 100644 Processes/OnShellCheck/testOnShellCheck.cc

diff --git a/Documentation/Examples/CMakeLists.txt b/Documentation/Examples/CMakeLists.txt
index ee95e2e89..38a5ab9a2 100644
--- a/Documentation/Examples/CMakeLists.txt
+++ b/Documentation/Examples/CMakeLists.txt
@@ -99,6 +99,7 @@ if (Pythia8_FOUND)
     ProcessTrackWriter
     ProcessTrackingLine
     ProcessParticleCut
+    ProcessOnShellCheck
     ProcessStackInspector
     CORSIKAprocesses
     CORSIKAcascade
diff --git a/Documentation/Examples/vertical_EAS.cc b/Documentation/Examples/vertical_EAS.cc
index 2cdb101b6..e716e3e66 100644
--- a/Documentation/Examples/vertical_EAS.cc
+++ b/Documentation/Examples/vertical_EAS.cc
@@ -21,6 +21,7 @@
 #include <corsika/process/interaction_counter/InteractionCounter.h>
 #include <corsika/process/observation_plane/ObservationPlane.h>
 #include <corsika/process/particle_cut/ParticleCut.h>
+#include <corsika/process/on_shell_check/OnShellCheck.h>
 #include <corsika/process/pythia/Decay.h>
 #include <corsika/process/sibyll/Decay.h>
 #include <corsika/process/sibyll/Interaction.h>
@@ -175,6 +176,8 @@ int main(int argc, char** argv) {
 
   process::particle_cut::ParticleCut cut(100_GeV);
 
+  process::on_shell_check::OnShellCheck reset_particle_mass(1.e-2,1.e-2);
+  
   process::energy_loss::EnergyLoss eLoss(showerAxis);
 
   Plane const obsPlane(Point(rootCS, 0_m, 0_m, observationHeight),
@@ -189,7 +192,7 @@ int main(int argc, char** argv) {
   auto sibyllSequence = sibyllNucCounted << sibyllCounted;
   process::switch_process::SwitchProcess switchProcess(urqmd, sibyllSequence, 55_GeV);
   auto decaySequence = decayPythia << decaySibyll;
-  auto sequence = switchProcess << decaySequence << eLoss << cut << observationLevel;
+  auto sequence = switchProcess << reset_particle_mass << decaySequence << eLoss << cut << observationLevel;
 
   // define air shower object, run simulation
   tracking_line::TrackingLine tracking;
diff --git a/Processes/CMakeLists.txt b/Processes/CMakeLists.txt
index 6e93793ec..5813f03dc 100644
--- a/Processes/CMakeLists.txt
+++ b/Processes/CMakeLists.txt
@@ -20,7 +20,7 @@ add_subdirectory (StackInspector)
 # secondaries process
 # cuts, thinning, etc.
 add_subdirectory (ParticleCut)
-
+add_subdirectory (OnShellCheck)
 # meta-processes
 add_subdirectory (InteractionCounter)
 add_subdirectory (SwitchProcess)
@@ -38,3 +38,4 @@ add_dependencies(CORSIKAprocesses ProcessTrackingLine)
 add_dependencies(CORSIKAprocesses ProcessEnergyLoss)
 add_dependencies(CORSIKAprocesses ProcessUrQMD)
 add_dependencies(CORSIKAprocesses ProcessParticleCut)
+add_dependencies(CORSIKAprocesses ProcessOnShellCheck)
diff --git a/Processes/OnShellCheck/CMakeLists.txt b/Processes/OnShellCheck/CMakeLists.txt
new file mode 100644
index 000000000..0d2191ad3
--- /dev/null
+++ b/Processes/OnShellCheck/CMakeLists.txt
@@ -0,0 +1,67 @@
+set (
+  MODEL_SOURCES
+  OnShellCheck.cc
+)
+
+set (
+  MODEL_HEADERS
+  OnShellCheck.h
+  )
+
+set (
+  MODEL_NAMESPACE
+  corsika/process/on_shell_check
+  )
+
+add_library (ProcessOnShellCheck STATIC ${MODEL_SOURCES})
+CORSIKA_COPY_HEADERS_TO_NAMESPACE (ProcessOnShellCheck ${MODEL_NAMESPACE} ${MODEL_HEADERS})
+
+set_target_properties (
+  ProcessOnShellCheck
+  PROPERTIES
+  VERSION ${PROJECT_VERSION}
+  SOVERSION 1
+#  PUBLIC_HEADER "${MODEL_HEADERS}"
+  )
+
+# target dependencies on other libraries (also the header onlys)
+target_link_libraries (
+  ProcessOnShellCheck
+  CORSIKAunits
+  CORSIKAparticles
+  CORSIKAprocesssequence
+  CORSIKAsetup
+  )
+
+target_include_directories (
+  ProcessOnShellCheck 
+  INTERFACE 
+  $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
+  $<INSTALL_INTERFACE:include/include>
+  )
+
+install (
+  TARGETS ProcessOnShellCheck
+  LIBRARY DESTINATION lib
+  ARCHIVE DESTINATION lib
+#  PUBLIC_HEADER DESTINATION include/${MODEL_NAMESPACE}
+  )
+
+# --------------------
+# code unit testing
+CORSIKA_ADD_TEST(testOnShellCheck SOURCES
+  testOnShellCheck.cc
+  ${MODEL_HEADERS}
+)
+
+target_link_libraries (
+  testOnShellCheck
+  ProcessOnShellCheck
+  CORSIKAunits
+  CORSIKAstackinterface
+  CORSIKAprocesssequence
+  CORSIKAsetup
+  CORSIKAgeometry
+  CORSIKAenvironment
+  CORSIKAtesting
+  )
diff --git a/Processes/OnShellCheck/OnShellCheck.cc b/Processes/OnShellCheck/OnShellCheck.cc
new file mode 100644
index 000000000..8abfc2011
--- /dev/null
+++ b/Processes/OnShellCheck/OnShellCheck.cc
@@ -0,0 +1,71 @@
+
+/*
+ * (c) Copyright 2018 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/on_shell_check/OnShellCheck.h>
+#include <corsika/geometry/FourVector.h>
+
+using namespace std;
+
+using namespace corsika;
+using namespace corsika::process;
+using namespace corsika::units::si;
+using namespace corsika::particles;
+using namespace corsika::setup;
+
+namespace corsika::process {
+  namespace on_shell_check {
+
+    void OnShellCheck::Init(){
+      std::cout << "OnShellCheck: mass tolerance is set to " << mass_tolerance_ * 100
+                << "%" << endl
+                << "              energy tolerance is set to " << energy_tolerance_ * 100
+                << "%" << std::endl;
+    }
+
+    EProcessReturn OnShellCheck::DoSecondaries(corsika::setup::StackView& vS) {
+      for (auto& p : vS) {
+        auto const pid = p.GetPID();
+	//if(pid==particles::Code::Gamma || particles::IsNeutrino(pid) || particles::IsNucleus(pid)) continue;
+	if(!particles::IsHadron(pid)) continue;
+        auto const e_original = p.GetEnergy();
+        auto const p_original = p.GetMomentum();
+        auto const Plab = corsika::geometry::FourVector(e_original, p_original);
+        auto const m_kinetic = Plab.GetNorm();
+        auto const m_corsika = particles::GetMass(pid);
+        auto const m_err = abs(m_kinetic - m_corsika) / m_corsika;
+        if (m_err > mass_tolerance_) {
+          const HEPEnergyType e_shifted =
+              sqrt(p_original.GetSquaredNorm() + m_corsika * m_corsika);
+          auto const e_shift_relative = (e_shifted / e_original - 1);
+	  /* 
+	     For now we warn if the necessary shift is larger than 1%.
+	     we could promote this to an error.
+	   */
+	  if (abs(e_shift_relative) > energy_tolerance_) {
+	    std::cout << "OnShellCheck::DoSecondaries: warning! shifted particle energy by "
+		      << e_shift_relative*100 << " %" << std::endl;
+	  }
+          std::cout << "OnShellCheck::DoSecondaries: shift particle mass for " << pid
+                    << std::endl
+                    << std::setw(20) << std::setfill(' ')
+                    << "corsika mass (GeV): " << m_corsika / 1_GeV << std::endl
+                    << "kinetic mass (GeV): " << m_kinetic / 1_GeV << std::endl
+                    << "(m_kin-m_cor)/en: " << m_err << std::endl
+                    << "mass tolerance: " << mass_tolerance_ << std::endl;
+          // reset energy
+          p.SetEnergy(e_shifted);
+        } else
+	  std::cout << "OnShellCheck::DoSecondaries: particle mass for " << pid << " OK" << std::endl;
+      }
+      return EProcessReturn::eOk;
+    }
+  }
+} // namespace on_shell_check
diff --git a/Processes/OnShellCheck/OnShellCheck.h b/Processes/OnShellCheck/OnShellCheck.h
new file mode 100644
index 000000000..f9a047324
--- /dev/null
+++ b/Processes/OnShellCheck/OnShellCheck.h
@@ -0,0 +1,39 @@
+/*
+ * (c) Copyright 2018 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 _corsika_process_on_shell_check_OnShellCheck_h_
+#define _corsika_process_on_shell_check_OnShellCheck_h_
+
+#include <corsika/particles/ParticleProperties.h>
+#include <corsika/process/SecondariesProcess.h>
+#include <corsika/setup/SetupStack.h>
+#include <corsika/units/PhysicalUnits.h>
+
+namespace corsika::process {
+  namespace on_shell_check {
+    class OnShellCheck : public process::SecondariesProcess<OnShellCheck> {
+
+    public:
+      OnShellCheck(const double vMassTolerance, const double vEnergyTolerance)
+          : mass_tolerance_(vMassTolerance)
+          , energy_tolerance_(vEnergyTolerance) {}
+
+      EProcessReturn DoSecondaries(corsika::setup::StackView&);
+
+      void Init();
+
+    private:
+      double mass_tolerance_;
+      double energy_tolerance_;
+    };
+  } // namespace on_shell_check
+} // namespace corsika::process
+
+#endif
diff --git a/Processes/OnShellCheck/ParticleCut.h b/Processes/OnShellCheck/ParticleCut.h
new file mode 100644
index 000000000..739a1919e
--- /dev/null
+++ b/Processes/OnShellCheck/ParticleCut.h
@@ -0,0 +1,55 @@
+/*
+ * (c) Copyright 2018 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 _corsika_process_particle_cut_ParticleCut_h_
+#define _corsika_process_particle_cut_ParticleCut_h_
+
+#include <corsika/particles/ParticleProperties.h>
+#include <corsika/process/SecondariesProcess.h>
+#include <corsika/setup/SetupStack.h>
+#include <corsika/units/PhysicalUnits.h>
+
+namespace corsika::process {
+  namespace particle_cut {
+    class ParticleCut : public process::SecondariesProcess<ParticleCut> {
+
+      units::si::HEPEnergyType const fECut;
+
+      units::si::HEPEnergyType fEnergy = 0 * units::si::electronvolt;
+      units::si::HEPEnergyType fEmEnergy = 0 * units::si::electronvolt;
+      unsigned int fEmCount = 0;
+      units::si::HEPEnergyType fInvEnergy = 0 * units::si::electronvolt;
+      unsigned int fInvCount = 0;
+
+    public:
+      ParticleCut(const units::si::HEPEnergyType vCut)
+          : fECut(vCut) {}
+
+      bool ParticleIsInvisible(particles::Code) const;
+      EProcessReturn DoSecondaries(corsika::setup::StackView&);
+
+      template <typename TParticle>
+      bool ParticleIsBelowEnergyCut(TParticle const&) const;
+
+      bool ParticleIsEmParticle(particles::Code) const;
+
+      void Init();
+      void ShowResults();
+
+      units::si::HEPEnergyType GetInvEnergy() const { return fInvEnergy; }
+      units::si::HEPEnergyType GetCutEnergy() const { return fEnergy; }
+      units::si::HEPEnergyType GetEmEnergy() const { return fEmEnergy; }
+      unsigned int GetNumberEmParticles() const { return fEmCount; }
+      unsigned int GetNumberInvParticles() const { return fInvCount; }
+    };
+  } // namespace particle_cut
+} // namespace corsika::process
+
+#endif
diff --git a/Processes/OnShellCheck/testOnShellCheck.cc b/Processes/OnShellCheck/testOnShellCheck.cc
new file mode 100644
index 000000000..4337446a9
--- /dev/null
+++ b/Processes/OnShellCheck/testOnShellCheck.cc
@@ -0,0 +1,91 @@
+/*
+ * (c) Copyright 2018 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/on_shell_check/OnShellCheck.h>
+
+#include <corsika/environment/Environment.h>
+#include <corsika/geometry/Point.h>
+#include <corsika/geometry/RootCoordinateSystem.h>
+#include <corsika/geometry/Vector.h>
+#include <corsika/geometry/FourVector.h>
+#include <corsika/units/PhysicalUnits.h>
+#include <corsika/utl/CorsikaFenv.h>
+
+#include <corsika/setup/SetupStack.h>
+
+#include <catch2/catch.hpp>
+
+using namespace corsika;
+using namespace corsika::process::on_shell_check;
+using namespace corsika::units;
+using namespace corsika::units::si;
+
+TEST_CASE("OnShellCheck", "[processes]") {
+  feenableexcept(FE_INVALID);
+  using EnvType = environment::Environment<setup::IEnvironmentModel>;
+  EnvType env;
+  const geometry::CoordinateSystem& rootCS = env.GetCoordinateSystem();
+
+  // setup empty particle stack
+  setup::Stack stack;
+  stack.Clear();
+  // two energies
+  const HEPEnergyType E = 10_GeV;
+  // list of arbitrary particles
+  std::array<particles::Code, 2> particleList = {
+      particles::Code::PiPlus,   particles::Code::PiMinus};
+
+  std::array<double, 2> mass_shifts = {1.1, 1.001};
+
+  SECTION("check particle masses") {
+
+    OnShellCheck check(1.e-2, 0.01);
+
+    check.Init();
+    
+    // add primary particle to stack
+    auto particle = stack.AddParticle(
+        std::tuple<particles::Code, units::si::HEPEnergyType,
+                   corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
+            particles::Code::Proton, E,
+            corsika::stack::MomentumVector(rootCS, {0_GeV, 0_GeV, 0_GeV}),
+            geometry::Point(rootCS, 0_m, 0_m, 0_m), 0_ns});
+    // view on secondary particles
+    corsika::stack::SecondaryView view(particle);
+    // ref. to primary particle through the secondary view.
+    // only this way the secondary view is populated
+    auto projectile = view.GetProjectile();
+    // add secondaries, all with energies above the threshold
+    // only cut is by species
+    int count = -1;
+    for (auto proType : particleList) {
+      count++;
+      const auto pz = sqrt((E - particles::GetMass(proType)*mass_shifts[count]) *
+                           (E + particles::GetMass(proType)*mass_shifts[count]));
+      auto const momentum = corsika::stack::MomentumVector(rootCS, {0_GeV, 0_GeV, pz});
+      projectile.AddSecondary(std::tuple<particles::Code, units::si::HEPEnergyType,
+                                         corsika::stack::MomentumVector, geometry::Point,
+                                         units::si::TimeType>{
+          proType, E, momentum, geometry::Point(rootCS, 0_m, 0_m, 0_m), 0_ns});
+
+      check.DoSecondaries(view);
+    }
+    int i = -1;
+    for ( auto& p : view) {
+      i++;
+      auto const Plab = corsika::geometry::FourVector(p.GetEnergy(), p.GetMomentum());
+      auto const m_kinetic = Plab.GetNorm();
+      if(i==0)
+	REQUIRE(m_kinetic / particles::PiPlus::GetMass() == Approx(1));
+      else
+	REQUIRE_FALSE(m_kinetic / particles::PiMinus::GetMass() == Approx(1));
+    }
+  }
+}
-- 
GitLab


From aec9a1c7ec2eddf0cd945896fdf79bbca19ac9a3 Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Thu, 4 Jun 2020 08:55:22 +0100
Subject: [PATCH 04/31] output format etc

---
 Processes/OnShellCheck/OnShellCheck.cc | 21 ++++++++++++---------
 1 file changed, 12 insertions(+), 9 deletions(-)

diff --git a/Processes/OnShellCheck/OnShellCheck.cc b/Processes/OnShellCheck/OnShellCheck.cc
index 8abfc2011..55a389421 100644
--- a/Processes/OnShellCheck/OnShellCheck.cc
+++ b/Processes/OnShellCheck/OnShellCheck.cc
@@ -40,8 +40,8 @@ namespace corsika::process {
         auto const Plab = corsika::geometry::FourVector(e_original, p_original);
         auto const m_kinetic = Plab.GetNorm();
         auto const m_corsika = particles::GetMass(pid);
-        auto const m_err = abs(m_kinetic - m_corsika) / m_corsika;
-        if (m_err > mass_tolerance_) {
+        auto const m_err_abs = abs(m_kinetic - m_corsika);
+        if (m_err_abs >= mass_tolerance_ * m_corsika) {
           const HEPEnergyType e_shifted =
               sqrt(p_original.GetSquaredNorm() + m_corsika * m_corsika);
           auto const e_shift_relative = (e_shifted / e_original - 1);
@@ -50,20 +50,23 @@ namespace corsika::process {
 	     we could promote this to an error.
 	   */
 	  if (abs(e_shift_relative) > energy_tolerance_) {
-	    std::cout << "OnShellCheck::DoSecondaries: warning! shifted particle energy by "
-		      << e_shift_relative*100 << " %" << std::endl;
+	    std::cout << "OnShellCheck: warning! shifted particle energy by "
+		      << e_shift_relative*100 << " %" << std::endl;	    
 	  }
-          std::cout << "OnShellCheck::DoSecondaries: shift particle mass for " << pid
+          std::cout << "OnShellCheck: shift particle mass for " << pid
                     << std::endl
-                    << std::setw(20) << std::setfill(' ')
+                    << std::setw(35) << std::setfill(' ')
                     << "corsika mass (GeV): " << m_corsika / 1_GeV << std::endl
+		    << std::setw(35) << std::setfill(' ')
                     << "kinetic mass (GeV): " << m_kinetic / 1_GeV << std::endl
-                    << "(m_kin-m_cor)/en: " << m_err << std::endl
-                    << "mass tolerance: " << mass_tolerance_ << std::endl;
+		    << std::setw(35) << std::setfill(' ')
+                    << "m_kin-m_cor (GeV): " << m_err_abs / 1_GeV << std::endl
+		    << std::setw(35) << std::setfill(' ')
+                    << "mass tolerance (GeV): " << ( m_corsika * mass_tolerance_ ) / 1_GeV << std::endl;
           // reset energy
           p.SetEnergy(e_shifted);
         } else
-	  std::cout << "OnShellCheck::DoSecondaries: particle mass for " << pid << " OK" << std::endl;
+          std::cout << "OnShellCheck: particle mass for " << pid << " OK" << std::endl;
       }
       return EProcessReturn::eOk;
     }
-- 
GitLab


From b91272823d7213e0d69ae05e9023e567edc2ca7d Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Thu, 4 Jun 2020 08:55:39 +0100
Subject: [PATCH 05/31] fixed test

---
 Processes/OnShellCheck/testOnShellCheck.cc | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/Processes/OnShellCheck/testOnShellCheck.cc b/Processes/OnShellCheck/testOnShellCheck.cc
index 4337446a9..62539e36a 100644
--- a/Processes/OnShellCheck/testOnShellCheck.cc
+++ b/Processes/OnShellCheck/testOnShellCheck.cc
@@ -73,10 +73,9 @@ TEST_CASE("OnShellCheck", "[processes]") {
       projectile.AddSecondary(std::tuple<particles::Code, units::si::HEPEnergyType,
                                          corsika::stack::MomentumVector, geometry::Point,
                                          units::si::TimeType>{
-          proType, E, momentum, geometry::Point(rootCS, 0_m, 0_m, 0_m), 0_ns});
-
-      check.DoSecondaries(view);
+          proType, E, momentum, geometry::Point(rootCS, 0_m, 0_m, 0_m), 0_ns});      
     }
+    check.DoSecondaries(view);
     int i = -1;
     for ( auto& p : view) {
       i++;
-- 
GitLab


From eed49b529c93e4a39c8d75e84f3cda29b9cd2595 Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Thu, 4 Jun 2020 09:05:51 +0100
Subject: [PATCH 06/31] set tolerance in vertical eas

---
 Documentation/Examples/vertical_EAS.cc | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/Documentation/Examples/vertical_EAS.cc b/Documentation/Examples/vertical_EAS.cc
index e716e3e66..16400cf58 100644
--- a/Documentation/Examples/vertical_EAS.cc
+++ b/Documentation/Examples/vertical_EAS.cc
@@ -176,8 +176,8 @@ int main(int argc, char** argv) {
 
   process::particle_cut::ParticleCut cut(100_GeV);
 
-  process::on_shell_check::OnShellCheck reset_particle_mass(1.e-2,1.e-2);
-  
+  process::on_shell_check::OnShellCheck reset_particle_mass(1.e-3, 1.e-2);
+
   process::energy_loss::EnergyLoss eLoss(showerAxis);
 
   Plane const obsPlane(Point(rootCS, 0_m, 0_m, observationHeight),
-- 
GitLab


From be4c20a9dfff4996b6fc8e9eae6250d3b775788b Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Thu, 4 Jun 2020 09:08:11 +0100
Subject: [PATCH 07/31] removed check from sibyll interface

---
 Processes/Sibyll/Interaction.cc | 23 -----------------------
 1 file changed, 23 deletions(-)

diff --git a/Processes/Sibyll/Interaction.cc b/Processes/Sibyll/Interaction.cc
index 0a4c7d59d..1a6e65598 100644
--- a/Processes/Sibyll/Interaction.cc
+++ b/Processes/Sibyll/Interaction.cc
@@ -317,29 +317,6 @@ namespace corsika::process::sibyll {
           auto const Plab = boost.fromCoM(FourVector(eCoM, pCoM));
 
           auto const pid = process::sibyll::ConvertFromSibyll(psib.GetPID());
-          // check if on-shell in corsika
-          auto const m_kinetic = Plab.GetNorm();
-          auto const m_corsika = particles::GetMass(pid);
-          auto const e_corsika = Plab.GetTimeLikeComponent();
-          auto const m_sibyll = corsika::process::sibyll::GetSibyllMass(pid);
-          auto const m_err = abs(m_kinetic - m_corsika) / m_corsika;
-          if (m_err > 1.e-5 && false) {
-            const HEPEnergyType e_shift_corsika = sqrt(
-                Plab.GetSpaceLikeComponents().GetSquaredNorm() + m_corsika * m_corsika);
-            auto const e_shift_relative = (e_shift_corsika / e_corsika - 1) * 100;
-            // warn about percent level shifts in particle energy
-            if (abs(e_shift_relative) > 1) {
-              std::cout << "Sibyll::Interaction: shifted particle energy by "
-                        << e_shift_relative << " %" << std::endl;
-
-              std::cout << "shift particle mass for " << pid << std::endl
-                        << "corsika mass (GeV): " << m_corsika / 1_GeV << std::endl
-                        << "kinetic mass (GeV): " << m_kinetic / 1_GeV << std::endl
-                        << "sibyll  mass (GeV): " << m_sibyll / 1_GeV << std::endl
-                        << "(m_kin-m_cor)/en: " << m_err << std::endl;
-            }
-            Plab.GetTimeLikeComponent() = e_shift_corsika;
-          }
 
           // add to corsika stack
           auto pnew = vP.AddSecondary(
-- 
GitLab


From abac61ac0b92dc338f7738c1659877fe5979217d Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Thu, 4 Jun 2020 09:42:24 +0100
Subject: [PATCH 08/31] clang

---
 Documentation/Examples/vertical_EAS.cc     |  5 +--
 Processes/OnShellCheck/OnShellCheck.cc     | 41 +++++++++++-----------
 Processes/OnShellCheck/testOnShellCheck.cc | 22 ++++++------
 3 files changed, 35 insertions(+), 33 deletions(-)

diff --git a/Documentation/Examples/vertical_EAS.cc b/Documentation/Examples/vertical_EAS.cc
index 16400cf58..b0de5727b 100644
--- a/Documentation/Examples/vertical_EAS.cc
+++ b/Documentation/Examples/vertical_EAS.cc
@@ -20,8 +20,8 @@
 #include <corsika/process/energy_loss/EnergyLoss.h>
 #include <corsika/process/interaction_counter/InteractionCounter.h>
 #include <corsika/process/observation_plane/ObservationPlane.h>
-#include <corsika/process/particle_cut/ParticleCut.h>
 #include <corsika/process/on_shell_check/OnShellCheck.h>
+#include <corsika/process/particle_cut/ParticleCut.h>
 #include <corsika/process/pythia/Decay.h>
 #include <corsika/process/sibyll/Decay.h>
 #include <corsika/process/sibyll/Interaction.h>
@@ -192,7 +192,8 @@ int main(int argc, char** argv) {
   auto sibyllSequence = sibyllNucCounted << sibyllCounted;
   process::switch_process::SwitchProcess switchProcess(urqmd, sibyllSequence, 55_GeV);
   auto decaySequence = decayPythia << decaySibyll;
-  auto sequence = switchProcess << reset_particle_mass << decaySequence << eLoss << cut << observationLevel;
+  auto sequence = switchProcess << reset_particle_mass << decaySequence << eLoss << cut
+                                << observationLevel;
 
   // define air shower object, run simulation
   tracking_line::TrackingLine tracking;
diff --git a/Processes/OnShellCheck/OnShellCheck.cc b/Processes/OnShellCheck/OnShellCheck.cc
index 55a389421..beee3a70e 100644
--- a/Processes/OnShellCheck/OnShellCheck.cc
+++ b/Processes/OnShellCheck/OnShellCheck.cc
@@ -9,8 +9,8 @@
  * the license.
  */
 
-#include <corsika/process/on_shell_check/OnShellCheck.h>
 #include <corsika/geometry/FourVector.h>
+#include <corsika/process/on_shell_check/OnShellCheck.h>
 
 using namespace std;
 
@@ -23,7 +23,7 @@ using namespace corsika::setup;
 namespace corsika::process {
   namespace on_shell_check {
 
-    void OnShellCheck::Init(){
+    void OnShellCheck::Init() {
       std::cout << "OnShellCheck: mass tolerance is set to " << mass_tolerance_ * 100
                 << "%" << endl
                 << "              energy tolerance is set to " << energy_tolerance_ * 100
@@ -33,8 +33,9 @@ namespace corsika::process {
     EProcessReturn OnShellCheck::DoSecondaries(corsika::setup::StackView& vS) {
       for (auto& p : vS) {
         auto const pid = p.GetPID();
-	//if(pid==particles::Code::Gamma || particles::IsNeutrino(pid) || particles::IsNucleus(pid)) continue;
-	if(!particles::IsHadron(pid)) continue;
+        // if(pid==particles::Code::Gamma || particles::IsNeutrino(pid) ||
+        // particles::IsNucleus(pid)) continue;
+        if (!particles::IsHadron(pid)) continue;
         auto const e_original = p.GetEnergy();
         auto const p_original = p.GetMomentum();
         auto const Plab = corsika::geometry::FourVector(e_original, p_original);
@@ -45,24 +46,24 @@ namespace corsika::process {
           const HEPEnergyType e_shifted =
               sqrt(p_original.GetSquaredNorm() + m_corsika * m_corsika);
           auto const e_shift_relative = (e_shifted / e_original - 1);
-	  /* 
-	     For now we warn if the necessary shift is larger than 1%.
-	     we could promote this to an error.
-	   */
-	  if (abs(e_shift_relative) > energy_tolerance_) {
-	    std::cout << "OnShellCheck: warning! shifted particle energy by "
-		      << e_shift_relative*100 << " %" << std::endl;	    
-	  }
-          std::cout << "OnShellCheck: shift particle mass for " << pid
-                    << std::endl
+          /*
+             For now we warn if the necessary shift is larger than 1%.
+             we could promote this to an error.
+           */
+          if (abs(e_shift_relative) > energy_tolerance_) {
+            std::cout << "OnShellCheck: warning! shifted particle energy by "
+                      << e_shift_relative * 100 << " %" << std::endl;
+          }
+          std::cout << "OnShellCheck: shift particle mass for " << pid << std::endl
                     << std::setw(35) << std::setfill(' ')
                     << "corsika mass (GeV): " << m_corsika / 1_GeV << std::endl
-		    << std::setw(35) << std::setfill(' ')
+                    << std::setw(35) << std::setfill(' ')
                     << "kinetic mass (GeV): " << m_kinetic / 1_GeV << std::endl
-		    << std::setw(35) << std::setfill(' ')
+                    << std::setw(35) << std::setfill(' ')
                     << "m_kin-m_cor (GeV): " << m_err_abs / 1_GeV << std::endl
-		    << std::setw(35) << std::setfill(' ')
-                    << "mass tolerance (GeV): " << ( m_corsika * mass_tolerance_ ) / 1_GeV << std::endl;
+                    << std::setw(35) << std::setfill(' ')
+                    << "mass tolerance (GeV): " << (m_corsika * mass_tolerance_) / 1_GeV
+                    << std::endl;
           // reset energy
           p.SetEnergy(e_shifted);
         } else
@@ -70,5 +71,5 @@ namespace corsika::process {
       }
       return EProcessReturn::eOk;
     }
-  }
-} // namespace on_shell_check
+  } // namespace on_shell_check
+} // namespace corsika::process
diff --git a/Processes/OnShellCheck/testOnShellCheck.cc b/Processes/OnShellCheck/testOnShellCheck.cc
index 62539e36a..bfe938305 100644
--- a/Processes/OnShellCheck/testOnShellCheck.cc
+++ b/Processes/OnShellCheck/testOnShellCheck.cc
@@ -11,10 +11,10 @@
 #include <corsika/process/on_shell_check/OnShellCheck.h>
 
 #include <corsika/environment/Environment.h>
+#include <corsika/geometry/FourVector.h>
 #include <corsika/geometry/Point.h>
 #include <corsika/geometry/RootCoordinateSystem.h>
 #include <corsika/geometry/Vector.h>
-#include <corsika/geometry/FourVector.h>
 #include <corsika/units/PhysicalUnits.h>
 #include <corsika/utl/CorsikaFenv.h>
 
@@ -39,8 +39,8 @@ TEST_CASE("OnShellCheck", "[processes]") {
   // two energies
   const HEPEnergyType E = 10_GeV;
   // list of arbitrary particles
-  std::array<particles::Code, 2> particleList = {
-      particles::Code::PiPlus,   particles::Code::PiMinus};
+  std::array<particles::Code, 2> particleList = {particles::Code::PiPlus,
+                                                 particles::Code::PiMinus};
 
   std::array<double, 2> mass_shifts = {1.1, 1.001};
 
@@ -49,7 +49,7 @@ TEST_CASE("OnShellCheck", "[processes]") {
     OnShellCheck check(1.e-2, 0.01);
 
     check.Init();
-    
+
     // add primary particle to stack
     auto particle = stack.AddParticle(
         std::tuple<particles::Code, units::si::HEPEnergyType,
@@ -67,24 +67,24 @@ TEST_CASE("OnShellCheck", "[processes]") {
     int count = -1;
     for (auto proType : particleList) {
       count++;
-      const auto pz = sqrt((E - particles::GetMass(proType)*mass_shifts[count]) *
-                           (E + particles::GetMass(proType)*mass_shifts[count]));
+      const auto pz = sqrt((E - particles::GetMass(proType) * mass_shifts[count]) *
+                           (E + particles::GetMass(proType) * mass_shifts[count]));
       auto const momentum = corsika::stack::MomentumVector(rootCS, {0_GeV, 0_GeV, pz});
       projectile.AddSecondary(std::tuple<particles::Code, units::si::HEPEnergyType,
                                          corsika::stack::MomentumVector, geometry::Point,
                                          units::si::TimeType>{
-          proType, E, momentum, geometry::Point(rootCS, 0_m, 0_m, 0_m), 0_ns});      
+          proType, E, momentum, geometry::Point(rootCS, 0_m, 0_m, 0_m), 0_ns});
     }
     check.DoSecondaries(view);
     int i = -1;
-    for ( auto& p : view) {
+    for (auto& p : view) {
       i++;
       auto const Plab = corsika::geometry::FourVector(p.GetEnergy(), p.GetMomentum());
       auto const m_kinetic = Plab.GetNorm();
-      if(i==0)
-	REQUIRE(m_kinetic / particles::PiPlus::GetMass() == Approx(1));
+      if (i == 0)
+        REQUIRE(m_kinetic / particles::PiPlus::GetMass() == Approx(1));
       else
-	REQUIRE_FALSE(m_kinetic / particles::PiMinus::GetMass() == Approx(1));
+        REQUIRE_FALSE(m_kinetic / particles::PiMinus::GetMass() == Approx(1));
     }
   }
 }
-- 
GitLab


From dc5e96efac0bf6c0f7ca3f5ff74462288c43848a Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Thu, 4 Jun 2020 10:50:54 +0100
Subject: [PATCH 09/31] removed obsolete file

---
 Processes/OnShellCheck/ParticleCut.h | 55 ----------------------------
 1 file changed, 55 deletions(-)
 delete mode 100644 Processes/OnShellCheck/ParticleCut.h

diff --git a/Processes/OnShellCheck/ParticleCut.h b/Processes/OnShellCheck/ParticleCut.h
deleted file mode 100644
index 739a1919e..000000000
--- a/Processes/OnShellCheck/ParticleCut.h
+++ /dev/null
@@ -1,55 +0,0 @@
-/*
- * (c) Copyright 2018 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 _corsika_process_particle_cut_ParticleCut_h_
-#define _corsika_process_particle_cut_ParticleCut_h_
-
-#include <corsika/particles/ParticleProperties.h>
-#include <corsika/process/SecondariesProcess.h>
-#include <corsika/setup/SetupStack.h>
-#include <corsika/units/PhysicalUnits.h>
-
-namespace corsika::process {
-  namespace particle_cut {
-    class ParticleCut : public process::SecondariesProcess<ParticleCut> {
-
-      units::si::HEPEnergyType const fECut;
-
-      units::si::HEPEnergyType fEnergy = 0 * units::si::electronvolt;
-      units::si::HEPEnergyType fEmEnergy = 0 * units::si::electronvolt;
-      unsigned int fEmCount = 0;
-      units::si::HEPEnergyType fInvEnergy = 0 * units::si::electronvolt;
-      unsigned int fInvCount = 0;
-
-    public:
-      ParticleCut(const units::si::HEPEnergyType vCut)
-          : fECut(vCut) {}
-
-      bool ParticleIsInvisible(particles::Code) const;
-      EProcessReturn DoSecondaries(corsika::setup::StackView&);
-
-      template <typename TParticle>
-      bool ParticleIsBelowEnergyCut(TParticle const&) const;
-
-      bool ParticleIsEmParticle(particles::Code) const;
-
-      void Init();
-      void ShowResults();
-
-      units::si::HEPEnergyType GetInvEnergy() const { return fInvEnergy; }
-      units::si::HEPEnergyType GetCutEnergy() const { return fEnergy; }
-      units::si::HEPEnergyType GetEmEnergy() const { return fEmEnergy; }
-      unsigned int GetNumberEmParticles() const { return fEmCount; }
-      unsigned int GetNumberInvParticles() const { return fInvCount; }
-    };
-  } // namespace particle_cut
-} // namespace corsika::process
-
-#endif
-- 
GitLab


From 7b3620e984d114b6d9754dfa14b9e1e987dc2308 Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Thu, 4 Jun 2020 11:25:23 +0100
Subject: [PATCH 10/31] avoid test for nuclei

---
 Processes/OnShellCheck/OnShellCheck.cc     |  2 +-
 Processes/OnShellCheck/testOnShellCheck.cc | 12 ++++++++----
 2 files changed, 9 insertions(+), 5 deletions(-)

diff --git a/Processes/OnShellCheck/OnShellCheck.cc b/Processes/OnShellCheck/OnShellCheck.cc
index beee3a70e..a5662e434 100644
--- a/Processes/OnShellCheck/OnShellCheck.cc
+++ b/Processes/OnShellCheck/OnShellCheck.cc
@@ -35,7 +35,7 @@ namespace corsika::process {
         auto const pid = p.GetPID();
         // if(pid==particles::Code::Gamma || particles::IsNeutrino(pid) ||
         // particles::IsNucleus(pid)) continue;
-        if (!particles::IsHadron(pid)) continue;
+        if (!particles::IsHadron(pid) || particles::IsNucleus(pid)) continue;
         auto const e_original = p.GetEnergy();
         auto const p_original = p.GetMomentum();
         auto const Plab = corsika::geometry::FourVector(e_original, p_original);
diff --git a/Processes/OnShellCheck/testOnShellCheck.cc b/Processes/OnShellCheck/testOnShellCheck.cc
index bfe938305..08e6cc323 100644
--- a/Processes/OnShellCheck/testOnShellCheck.cc
+++ b/Processes/OnShellCheck/testOnShellCheck.cc
@@ -39,10 +39,13 @@ TEST_CASE("OnShellCheck", "[processes]") {
   // two energies
   const HEPEnergyType E = 10_GeV;
   // list of arbitrary particles
-  std::array<particles::Code, 2> particleList = {particles::Code::PiPlus,
-                                                 particles::Code::PiMinus};
+  std::array<particles::Code, 4> particleList = {
+      particles::Code::PiPlus,
+      particles::Code::PiMinus,
+      particles::Code::Helium,
+      particles::Code::Gamma};
 
-  std::array<double, 2> mass_shifts = {1.1, 1.001};
+  std::array<double, 4> mass_shifts = {1.1, 1.001, 1.0, 1.0};
 
   SECTION("check particle masses") {
 
@@ -83,8 +86,9 @@ TEST_CASE("OnShellCheck", "[processes]") {
       auto const m_kinetic = Plab.GetNorm();
       if (i == 0)
         REQUIRE(m_kinetic / particles::PiPlus::GetMass() == Approx(1));
-      else
+      else if (i == 1)
         REQUIRE_FALSE(m_kinetic / particles::PiMinus::GetMass() == Approx(1));
+	
     }
   }
 }
-- 
GitLab


From 2153a4f3d0bd3a9d6de8e4e2160f0a75027f9cf1 Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Thu, 4 Jun 2020 11:26:02 +0100
Subject: [PATCH 11/31] clang

---
 Processes/OnShellCheck/testOnShellCheck.cc | 7 ++-----
 1 file changed, 2 insertions(+), 5 deletions(-)

diff --git a/Processes/OnShellCheck/testOnShellCheck.cc b/Processes/OnShellCheck/testOnShellCheck.cc
index 08e6cc323..bd366d9b5 100644
--- a/Processes/OnShellCheck/testOnShellCheck.cc
+++ b/Processes/OnShellCheck/testOnShellCheck.cc
@@ -40,9 +40,7 @@ TEST_CASE("OnShellCheck", "[processes]") {
   const HEPEnergyType E = 10_GeV;
   // list of arbitrary particles
   std::array<particles::Code, 4> particleList = {
-      particles::Code::PiPlus,
-      particles::Code::PiMinus,
-      particles::Code::Helium,
+      particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::Helium,
       particles::Code::Gamma};
 
   std::array<double, 4> mass_shifts = {1.1, 1.001, 1.0, 1.0};
@@ -61,7 +59,7 @@ TEST_CASE("OnShellCheck", "[processes]") {
             corsika::stack::MomentumVector(rootCS, {0_GeV, 0_GeV, 0_GeV}),
             geometry::Point(rootCS, 0_m, 0_m, 0_m), 0_ns});
     // view on secondary particles
-    corsika::stack::SecondaryView view(particle);
+    corsika::stack::SecondaryView view(part./ icle);
     // ref. to primary particle through the secondary view.
     // only this way the secondary view is populated
     auto projectile = view.GetProjectile();
@@ -88,7 +86,6 @@ TEST_CASE("OnShellCheck", "[processes]") {
         REQUIRE(m_kinetic / particles::PiPlus::GetMass() == Approx(1));
       else if (i == 1)
         REQUIRE_FALSE(m_kinetic / particles::PiMinus::GetMass() == Approx(1));
-	
     }
   }
 }
-- 
GitLab


From 58bffbed56ab37d3e0538cd3f55522d27ac26802 Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Thu, 4 Jun 2020 11:27:20 +0100
Subject: [PATCH 12/31] there you go

---
 Processes/OnShellCheck/testOnShellCheck.cc | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Processes/OnShellCheck/testOnShellCheck.cc b/Processes/OnShellCheck/testOnShellCheck.cc
index bd366d9b5..9f50e0b37 100644
--- a/Processes/OnShellCheck/testOnShellCheck.cc
+++ b/Processes/OnShellCheck/testOnShellCheck.cc
@@ -59,7 +59,7 @@ TEST_CASE("OnShellCheck", "[processes]") {
             corsika::stack::MomentumVector(rootCS, {0_GeV, 0_GeV, 0_GeV}),
             geometry::Point(rootCS, 0_m, 0_m, 0_m), 0_ns});
     // view on secondary particles
-    corsika::stack::SecondaryView view(part./ icle);
+    corsika::stack::SecondaryView view(particle);
     // ref. to primary particle through the secondary view.
     // only this way the secondary view is populated
     auto projectile = view.GetProjectile();
-- 
GitLab


From d6c6be6c9a8867b6ecf3ffe56ad96c7423452b58 Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Fri, 12 Jun 2020 13:06:04 +0100
Subject: [PATCH 13/31] added summary

---
 Processes/OnShellCheck/OnShellCheck.cc |  3 +++
 Processes/OnShellCheck/OnShellCheck.h  | 12 ++++++++++++
 2 files changed, 15 insertions(+)

diff --git a/Processes/OnShellCheck/OnShellCheck.cc b/Processes/OnShellCheck/OnShellCheck.cc
index a5662e434..133af10c6 100644
--- a/Processes/OnShellCheck/OnShellCheck.cc
+++ b/Processes/OnShellCheck/OnShellCheck.cc
@@ -46,6 +46,9 @@ namespace corsika::process {
           const HEPEnergyType e_shifted =
               sqrt(p_original.GetSquaredNorm() + m_corsika * m_corsika);
           auto const e_shift_relative = (e_shifted / e_original - 1);
+          count_ = count_ + 1;
+          average_shift_ += abs(e_shift_relative);
+          if (abs(e_shift_relative) > max_shift_) max_shift_ = abs(e_shift_relative);
           /*
              For now we warn if the necessary shift is larger than 1%.
              we could promote this to an error.
diff --git a/Processes/OnShellCheck/OnShellCheck.h b/Processes/OnShellCheck/OnShellCheck.h
index f9a047324..f1c588be0 100644
--- a/Processes/OnShellCheck/OnShellCheck.h
+++ b/Processes/OnShellCheck/OnShellCheck.h
@@ -19,12 +19,24 @@
 namespace corsika::process {
   namespace on_shell_check {
     class OnShellCheck : public process::SecondariesProcess<OnShellCheck> {
+      double average_shift_ = 0;
+      double max_shift_ = 0;
+      double count_ = 0;
 
     public:
       OnShellCheck(const double vMassTolerance, const double vEnergyTolerance)
           : mass_tolerance_(vMassTolerance)
           , energy_tolerance_(vEnergyTolerance) {}
 
+      ~OnShellCheck() {
+        std::cout << "OnShellCheck: summary" << std::endl
+                  << " particles shifted: " << int(count_) << std::endl;
+        if (count_)
+          std::cout << " average energy shift (%): " << average_shift_ / count_ * 100.
+                    << std::endl
+                    << " max. energy shift (%): " << max_shift_ * 100. << std::endl;
+      };
+
       EProcessReturn DoSecondaries(corsika::setup::StackView&);
 
       void Init();
-- 
GitLab


From 067fe0914c0d5b0f726d2467a4804efdfe6cdaa4 Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Fri, 12 Jun 2020 15:04:40 +0100
Subject: [PATCH 14/31] added error for extreme cases

---
 Documentation/Examples/vertical_EAS.cc     |  2 +-
 Processes/OnShellCheck/OnShellCheck.cc     | 31 +++++++++++-----------
 Processes/OnShellCheck/OnShellCheck.h      |  7 +++--
 Processes/OnShellCheck/testOnShellCheck.cc |  2 +-
 4 files changed, 23 insertions(+), 19 deletions(-)

diff --git a/Documentation/Examples/vertical_EAS.cc b/Documentation/Examples/vertical_EAS.cc
index b0de5727b..417c2f2f2 100644
--- a/Documentation/Examples/vertical_EAS.cc
+++ b/Documentation/Examples/vertical_EAS.cc
@@ -176,7 +176,7 @@ int main(int argc, char** argv) {
 
   process::particle_cut::ParticleCut cut(100_GeV);
 
-  process::on_shell_check::OnShellCheck reset_particle_mass(1.e-3, 1.e-2);
+  process::on_shell_check::OnShellCheck reset_particle_mass(1.e-3, 1.e-1, false);
 
   process::energy_loss::EnergyLoss eLoss(showerAxis);
 
diff --git a/Processes/OnShellCheck/OnShellCheck.cc b/Processes/OnShellCheck/OnShellCheck.cc
index 133af10c6..f47f94d25 100644
--- a/Processes/OnShellCheck/OnShellCheck.cc
+++ b/Processes/OnShellCheck/OnShellCheck.cc
@@ -1,4 +1,3 @@
-
 /*
  * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
  *
@@ -33,8 +32,6 @@ namespace corsika::process {
     EProcessReturn OnShellCheck::DoSecondaries(corsika::setup::StackView& vS) {
       for (auto& p : vS) {
         auto const pid = p.GetPID();
-        // if(pid==particles::Code::Gamma || particles::IsNeutrino(pid) ||
-        // particles::IsNucleus(pid)) continue;
         if (!particles::IsHadron(pid) || particles::IsNucleus(pid)) continue;
         auto const e_original = p.GetEnergy();
         auto const p_original = p.GetMomentum();
@@ -49,24 +46,28 @@ namespace corsika::process {
           count_ = count_ + 1;
           average_shift_ += abs(e_shift_relative);
           if (abs(e_shift_relative) > max_shift_) max_shift_ = abs(e_shift_relative);
-          /*
-             For now we warn if the necessary shift is larger than 1%.
-             we could promote this to an error.
-           */
-          if (abs(e_shift_relative) > energy_tolerance_) {
-            std::cout << "OnShellCheck: warning! shifted particle energy by "
-                      << e_shift_relative * 100 << " %" << std::endl;
-          }
           std::cout << "OnShellCheck: shift particle mass for " << pid << std::endl
-                    << std::setw(35) << std::setfill(' ')
+                    << std::setw(40) << std::setfill(' ')
                     << "corsika mass (GeV): " << m_corsika / 1_GeV << std::endl
-                    << std::setw(35) << std::setfill(' ')
+                    << std::setw(40) << std::setfill(' ')
                     << "kinetic mass (GeV): " << m_kinetic / 1_GeV << std::endl
-                    << std::setw(35) << std::setfill(' ')
+                    << std::setw(40) << std::setfill(' ')
                     << "m_kin-m_cor (GeV): " << m_err_abs / 1_GeV << std::endl
-                    << std::setw(35) << std::setfill(' ')
+                    << std::setw(40) << std::setfill(' ')
                     << "mass tolerance (GeV): " << (m_corsika * mass_tolerance_) / 1_GeV
                     << std::endl;
+          /*
+            For now we warn if the necessary shift is larger than 1%.
+            we could promote this to an error.
+          */
+          if (abs(e_shift_relative) > energy_tolerance_) {
+            std::cout << "OnShellCheck: warning! shifted particle energy by "
+                      << e_shift_relative * 100 << " %" << std::endl;
+            if (throw_error_)
+              throw std::runtime_error(
+                  "OnShellCheck: error! shifted energy by large amount!");
+          }
+
           // reset energy
           p.SetEnergy(e_shifted);
         } else
diff --git a/Processes/OnShellCheck/OnShellCheck.h b/Processes/OnShellCheck/OnShellCheck.h
index f1c588be0..d51d5e8f7 100644
--- a/Processes/OnShellCheck/OnShellCheck.h
+++ b/Processes/OnShellCheck/OnShellCheck.h
@@ -24,9 +24,11 @@ namespace corsika::process {
       double count_ = 0;
 
     public:
-      OnShellCheck(const double vMassTolerance, const double vEnergyTolerance)
+      OnShellCheck(const double vMassTolerance, const double vEnergyTolerance,
+                   const bool vError)
           : mass_tolerance_(vMassTolerance)
-          , energy_tolerance_(vEnergyTolerance) {}
+          , energy_tolerance_(vEnergyTolerance)
+          , throw_error_(vError) {}
 
       ~OnShellCheck() {
         std::cout << "OnShellCheck: summary" << std::endl
@@ -44,6 +46,7 @@ namespace corsika::process {
     private:
       double mass_tolerance_;
       double energy_tolerance_;
+      bool throw_error_;
     };
   } // namespace on_shell_check
 } // namespace corsika::process
diff --git a/Processes/OnShellCheck/testOnShellCheck.cc b/Processes/OnShellCheck/testOnShellCheck.cc
index 9f50e0b37..f67272bd2 100644
--- a/Processes/OnShellCheck/testOnShellCheck.cc
+++ b/Processes/OnShellCheck/testOnShellCheck.cc
@@ -47,7 +47,7 @@ TEST_CASE("OnShellCheck", "[processes]") {
 
   SECTION("check particle masses") {
 
-    OnShellCheck check(1.e-2, 0.01);
+    OnShellCheck check(1.e-2, 0.01, false);
 
     check.Init();
 
-- 
GitLab


From 0af4a0df8e113b66d8da9461b4ccd1673672f393 Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Wed, 3 Jun 2020 18:52:58 +0100
Subject: [PATCH 15/31] shift energies of particles from sibyll so that mass is
 consistent, implemented after boost to lab frame

---
 Processes/Pythia/Decay.cc       |  4 ++--
 Processes/Sibyll/Interaction.cc | 38 +++++++++++++++++++++++++++++----
 2 files changed, 36 insertions(+), 6 deletions(-)

diff --git a/Processes/Pythia/Decay.cc b/Processes/Pythia/Decay.cc
index 27890d28b..a17678d18 100644
--- a/Processes/Pythia/Decay.cc
+++ b/Processes/Pythia/Decay.cc
@@ -57,7 +57,7 @@ namespace corsika::process::pythia {
     fPythia.readString("Next:numberShowEvent = 0");
 
     fPythia.readString("Print:quiet = on");
-    fPythia.readString("Check:particleData = 1");
+    fPythia.readString("Check:particleData = 0");
 
     /*
        switching off event check in pythia is needed to allow decays that are off-shell
@@ -65,7 +65,7 @@ namespace corsika::process::pythia {
        the consistency of particle masses between event generators is an unsolved issues
     */
     cout << "Pythia::Init: switching off event checking in pythia.." << endl;
-    fPythia.readString("Check:event = 0");
+    fPythia.readString("Check:event = 1");
 
     fPythia.readString("ProcessLevel:all = off");
     fPythia.readString("ProcessLevel:resonanceDecays = off");
diff --git a/Processes/Sibyll/Interaction.cc b/Processes/Sibyll/Interaction.cc
index 916ea2ec0..e7d610dc5 100644
--- a/Processes/Sibyll/Interaction.cc
+++ b/Processes/Sibyll/Interaction.cc
@@ -316,20 +316,50 @@ namespace corsika::process::sibyll {
           HEPEnergyType const eCoM = psib.GetEnergy();
           auto const Plab = boost.fromCoM(FourVector(eCoM, pCoM));
 
+          auto const pid = process::sibyll::ConvertFromSibyll(psib.GetPID());
+          // check if on-shell in corsika
+          auto const m_kinetic = Plab.GetNorm();
+          auto const m_corsika = particles::GetMass(pid);
+          auto const e_corsika = Plab.GetTimeLikeComponent();
+          auto const m_sibyll = corsika::process::sibyll::GetSibyllMass(pid);
+          auto const m_err = abs(m_kinetic - m_corsika) / m_corsika;
+          if (m_err > 1.e-5) {
+            const HEPEnergyType e_shift_corsika = sqrt(
+                Plab.GetSpaceLikeComponents().GetSquaredNorm() + m_corsika * m_corsika);
+            auto const e_shift_relative = (e_shift_corsika / e_corsika - 1) * 100;
+            // warn about percent level shifts in particle energy
+            if (abs(e_shift_relative) > 1) {
+              std::cout << "Sibyll::Interaction: shifted particle energy by "
+                        << e_shift_relative << " %" << std::endl;
+
+              std::cout << "shift particle mass for " << pid << std::endl
+                        << "corsika mass (GeV): " << m_corsika / 1_GeV << std::endl
+                        << "kinetic mass (GeV): " << m_kinetic / 1_GeV << std::endl
+                        << "sibyll  mass (GeV): " << m_sibyll / 1_GeV << std::endl
+                        << "(m_kin-m_cor)/en: " << m_err << std::endl;
+            }
+            Plab.GetTimeLikeComponent() = e_shift_corsika;
+          }
+
           // add to corsika stack
           auto pnew = vP.AddSecondary(
               tuple<particles::Code, units::si::HEPEnergyType, stack::MomentumVector,
                     geometry::Point, units::si::TimeType>{
-                  process::sibyll::ConvertFromSibyll(psib.GetPID()),
-                  Plab.GetTimeLikeComponent(), Plab.GetSpaceLikeComponents(), pOrig,
+                  pid, Plab.GetTimeLikeComponent(), Plab.GetSpaceLikeComponents(), pOrig,
                   tOrig});
 
           Plab_final += pnew.GetMomentum();
           Elab_final += pnew.GetEnergy();
           Ecm_final += psib.GetEnergy();
         }
-        cout << "conservation (all GeV): Ecm_final=" << Ecm_final / 1_GeV << endl
-             << "Elab_final=" << Elab_final / 1_GeV
+        cout << "conservation (all GeV):" << endl
+             << "Ecm_initial=" << Ecm / 1_GeV << " Ecm_final=" << Ecm_final / 1_GeV
+             << endl
+             << "Elab_initial=" << eProjectileLab / 1_GeV
+             << " Elab_final=" << Elab_final / 1_GeV
+             << " diff (%)=" << (Elab_final / eProjectileLab / get_nwounded() - 1) * 100
+             << endl
+             << "Plab_initial=" << (pProjectileLab / 1_GeV).GetComponents()
              << ", Plab_final=" << (Plab_final / 1_GeV).GetComponents() << endl;
       }
     }
-- 
GitLab


From d12b7c37f40da1525cf388a15e23d1c7186276f6 Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Wed, 3 Jun 2020 23:17:39 +0100
Subject: [PATCH 16/31] shift particle masses in sibyll interface

---
 Processes/Sibyll/Interaction.cc | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Processes/Sibyll/Interaction.cc b/Processes/Sibyll/Interaction.cc
index e7d610dc5..0a4c7d59d 100644
--- a/Processes/Sibyll/Interaction.cc
+++ b/Processes/Sibyll/Interaction.cc
@@ -323,7 +323,7 @@ namespace corsika::process::sibyll {
           auto const e_corsika = Plab.GetTimeLikeComponent();
           auto const m_sibyll = corsika::process::sibyll::GetSibyllMass(pid);
           auto const m_err = abs(m_kinetic - m_corsika) / m_corsika;
-          if (m_err > 1.e-5) {
+          if (m_err > 1.e-5 && false) {
             const HEPEnergyType e_shift_corsika = sqrt(
                 Plab.GetSpaceLikeComponents().GetSquaredNorm() + m_corsika * m_corsika);
             auto const e_shift_relative = (e_shift_corsika / e_corsika - 1) * 100;
-- 
GitLab


From 61df23dc8f4c7db48ee3f8c0fbeeefc478100d31 Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Wed, 3 Jun 2020 23:18:40 +0100
Subject: [PATCH 17/31] added process to check particle masses and shift
 on-shell if necessary

---
 Documentation/Examples/CMakeLists.txt      |  1 +
 Documentation/Examples/vertical_EAS.cc     |  6 +-
 Processes/CMakeLists.txt                   |  3 +-
 Processes/OnShellCheck/CMakeLists.txt      | 67 ++++++++++++++++
 Processes/OnShellCheck/OnShellCheck.cc     | 71 +++++++++++++++++
 Processes/OnShellCheck/OnShellCheck.h      | 39 ++++++++++
 Processes/OnShellCheck/ParticleCut.h       | 55 +++++++++++++
 Processes/OnShellCheck/testOnShellCheck.cc | 91 ++++++++++++++++++++++
 8 files changed, 331 insertions(+), 2 deletions(-)
 create mode 100644 Processes/OnShellCheck/CMakeLists.txt
 create mode 100644 Processes/OnShellCheck/OnShellCheck.cc
 create mode 100644 Processes/OnShellCheck/OnShellCheck.h
 create mode 100644 Processes/OnShellCheck/ParticleCut.h
 create mode 100644 Processes/OnShellCheck/testOnShellCheck.cc

diff --git a/Documentation/Examples/CMakeLists.txt b/Documentation/Examples/CMakeLists.txt
index fc7c01522..0118d7c87 100644
--- a/Documentation/Examples/CMakeLists.txt
+++ b/Documentation/Examples/CMakeLists.txt
@@ -99,6 +99,7 @@ if (Pythia8_FOUND)
     ProcessTrackWriter
     ProcessTrackingLine
     ProcessParticleCut
+    ProcessOnShellCheck
     ProcessStackInspector
     ProcessLongitudinalProfile
     CORSIKAprocesses
diff --git a/Documentation/Examples/vertical_EAS.cc b/Documentation/Examples/vertical_EAS.cc
index a81414dc1..974398dfe 100644
--- a/Documentation/Examples/vertical_EAS.cc
+++ b/Documentation/Examples/vertical_EAS.cc
@@ -23,6 +23,7 @@
 #include <corsika/process/longitudinal_profile/LongitudinalProfile.h>
 #include <corsika/process/observation_plane/ObservationPlane.h>
 #include <corsika/process/particle_cut/ParticleCut.h>
+#include <corsika/process/on_shell_check/OnShellCheck.h>
 #include <corsika/process/pythia/Decay.h>
 #include <corsika/process/sibyll/Decay.h>
 #include <corsika/process/sibyll/Interaction.h>
@@ -183,6 +184,8 @@ int main(int argc, char** argv) {
 
   process::particle_cut::ParticleCut cut{60_GeV};
 
+  process::on_shell_check::OnShellCheck reset_particle_mass(1.e-2,1.e-2);
+  
   process::energy_loss::EnergyLoss eLoss(showerAxis);
   process::longitudinal_profile::LongitudinalProfile longprof{showerAxis};
 
@@ -199,7 +202,8 @@ int main(int argc, char** argv) {
   process::switch_process::SwitchProcess switchProcess(urqmdCounted, sibyllSequence,
                                                        55_GeV);
   auto decaySequence = decayPythia << decaySibyll;
-  auto sequence = switchProcess << decaySequence << longprof << eLoss << cut
+
+  auto sequence = switchProcess << reset_particle_mass << decaySequence << longprof << eLoss << cut
                                 << observationLevel;
 
   // define air shower object, run simulation
diff --git a/Processes/CMakeLists.txt b/Processes/CMakeLists.txt
index 536077c5b..f27d027a1 100644
--- a/Processes/CMakeLists.txt
+++ b/Processes/CMakeLists.txt
@@ -24,7 +24,7 @@ add_subdirectory (StackInspector)
 # secondaries process
 # cuts, thinning, etc.
 add_subdirectory (ParticleCut)
-
+add_subdirectory (OnShellCheck)
 # meta-processes
 add_subdirectory (InteractionCounter)
 add_subdirectory (SwitchProcess)
@@ -42,3 +42,4 @@ add_dependencies(CORSIKAprocesses ProcessTrackingLine)
 add_dependencies(CORSIKAprocesses ProcessEnergyLoss)
 add_dependencies(CORSIKAprocesses ProcessUrQMD)
 add_dependencies(CORSIKAprocesses ProcessParticleCut)
+add_dependencies(CORSIKAprocesses ProcessOnShellCheck)
diff --git a/Processes/OnShellCheck/CMakeLists.txt b/Processes/OnShellCheck/CMakeLists.txt
new file mode 100644
index 000000000..0d2191ad3
--- /dev/null
+++ b/Processes/OnShellCheck/CMakeLists.txt
@@ -0,0 +1,67 @@
+set (
+  MODEL_SOURCES
+  OnShellCheck.cc
+)
+
+set (
+  MODEL_HEADERS
+  OnShellCheck.h
+  )
+
+set (
+  MODEL_NAMESPACE
+  corsika/process/on_shell_check
+  )
+
+add_library (ProcessOnShellCheck STATIC ${MODEL_SOURCES})
+CORSIKA_COPY_HEADERS_TO_NAMESPACE (ProcessOnShellCheck ${MODEL_NAMESPACE} ${MODEL_HEADERS})
+
+set_target_properties (
+  ProcessOnShellCheck
+  PROPERTIES
+  VERSION ${PROJECT_VERSION}
+  SOVERSION 1
+#  PUBLIC_HEADER "${MODEL_HEADERS}"
+  )
+
+# target dependencies on other libraries (also the header onlys)
+target_link_libraries (
+  ProcessOnShellCheck
+  CORSIKAunits
+  CORSIKAparticles
+  CORSIKAprocesssequence
+  CORSIKAsetup
+  )
+
+target_include_directories (
+  ProcessOnShellCheck 
+  INTERFACE 
+  $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
+  $<INSTALL_INTERFACE:include/include>
+  )
+
+install (
+  TARGETS ProcessOnShellCheck
+  LIBRARY DESTINATION lib
+  ARCHIVE DESTINATION lib
+#  PUBLIC_HEADER DESTINATION include/${MODEL_NAMESPACE}
+  )
+
+# --------------------
+# code unit testing
+CORSIKA_ADD_TEST(testOnShellCheck SOURCES
+  testOnShellCheck.cc
+  ${MODEL_HEADERS}
+)
+
+target_link_libraries (
+  testOnShellCheck
+  ProcessOnShellCheck
+  CORSIKAunits
+  CORSIKAstackinterface
+  CORSIKAprocesssequence
+  CORSIKAsetup
+  CORSIKAgeometry
+  CORSIKAenvironment
+  CORSIKAtesting
+  )
diff --git a/Processes/OnShellCheck/OnShellCheck.cc b/Processes/OnShellCheck/OnShellCheck.cc
new file mode 100644
index 000000000..8abfc2011
--- /dev/null
+++ b/Processes/OnShellCheck/OnShellCheck.cc
@@ -0,0 +1,71 @@
+
+/*
+ * (c) Copyright 2018 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/on_shell_check/OnShellCheck.h>
+#include <corsika/geometry/FourVector.h>
+
+using namespace std;
+
+using namespace corsika;
+using namespace corsika::process;
+using namespace corsika::units::si;
+using namespace corsika::particles;
+using namespace corsika::setup;
+
+namespace corsika::process {
+  namespace on_shell_check {
+
+    void OnShellCheck::Init(){
+      std::cout << "OnShellCheck: mass tolerance is set to " << mass_tolerance_ * 100
+                << "%" << endl
+                << "              energy tolerance is set to " << energy_tolerance_ * 100
+                << "%" << std::endl;
+    }
+
+    EProcessReturn OnShellCheck::DoSecondaries(corsika::setup::StackView& vS) {
+      for (auto& p : vS) {
+        auto const pid = p.GetPID();
+	//if(pid==particles::Code::Gamma || particles::IsNeutrino(pid) || particles::IsNucleus(pid)) continue;
+	if(!particles::IsHadron(pid)) continue;
+        auto const e_original = p.GetEnergy();
+        auto const p_original = p.GetMomentum();
+        auto const Plab = corsika::geometry::FourVector(e_original, p_original);
+        auto const m_kinetic = Plab.GetNorm();
+        auto const m_corsika = particles::GetMass(pid);
+        auto const m_err = abs(m_kinetic - m_corsika) / m_corsika;
+        if (m_err > mass_tolerance_) {
+          const HEPEnergyType e_shifted =
+              sqrt(p_original.GetSquaredNorm() + m_corsika * m_corsika);
+          auto const e_shift_relative = (e_shifted / e_original - 1);
+	  /* 
+	     For now we warn if the necessary shift is larger than 1%.
+	     we could promote this to an error.
+	   */
+	  if (abs(e_shift_relative) > energy_tolerance_) {
+	    std::cout << "OnShellCheck::DoSecondaries: warning! shifted particle energy by "
+		      << e_shift_relative*100 << " %" << std::endl;
+	  }
+          std::cout << "OnShellCheck::DoSecondaries: shift particle mass for " << pid
+                    << std::endl
+                    << std::setw(20) << std::setfill(' ')
+                    << "corsika mass (GeV): " << m_corsika / 1_GeV << std::endl
+                    << "kinetic mass (GeV): " << m_kinetic / 1_GeV << std::endl
+                    << "(m_kin-m_cor)/en: " << m_err << std::endl
+                    << "mass tolerance: " << mass_tolerance_ << std::endl;
+          // reset energy
+          p.SetEnergy(e_shifted);
+        } else
+	  std::cout << "OnShellCheck::DoSecondaries: particle mass for " << pid << " OK" << std::endl;
+      }
+      return EProcessReturn::eOk;
+    }
+  }
+} // namespace on_shell_check
diff --git a/Processes/OnShellCheck/OnShellCheck.h b/Processes/OnShellCheck/OnShellCheck.h
new file mode 100644
index 000000000..f9a047324
--- /dev/null
+++ b/Processes/OnShellCheck/OnShellCheck.h
@@ -0,0 +1,39 @@
+/*
+ * (c) Copyright 2018 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 _corsika_process_on_shell_check_OnShellCheck_h_
+#define _corsika_process_on_shell_check_OnShellCheck_h_
+
+#include <corsika/particles/ParticleProperties.h>
+#include <corsika/process/SecondariesProcess.h>
+#include <corsika/setup/SetupStack.h>
+#include <corsika/units/PhysicalUnits.h>
+
+namespace corsika::process {
+  namespace on_shell_check {
+    class OnShellCheck : public process::SecondariesProcess<OnShellCheck> {
+
+    public:
+      OnShellCheck(const double vMassTolerance, const double vEnergyTolerance)
+          : mass_tolerance_(vMassTolerance)
+          , energy_tolerance_(vEnergyTolerance) {}
+
+      EProcessReturn DoSecondaries(corsika::setup::StackView&);
+
+      void Init();
+
+    private:
+      double mass_tolerance_;
+      double energy_tolerance_;
+    };
+  } // namespace on_shell_check
+} // namespace corsika::process
+
+#endif
diff --git a/Processes/OnShellCheck/ParticleCut.h b/Processes/OnShellCheck/ParticleCut.h
new file mode 100644
index 000000000..739a1919e
--- /dev/null
+++ b/Processes/OnShellCheck/ParticleCut.h
@@ -0,0 +1,55 @@
+/*
+ * (c) Copyright 2018 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 _corsika_process_particle_cut_ParticleCut_h_
+#define _corsika_process_particle_cut_ParticleCut_h_
+
+#include <corsika/particles/ParticleProperties.h>
+#include <corsika/process/SecondariesProcess.h>
+#include <corsika/setup/SetupStack.h>
+#include <corsika/units/PhysicalUnits.h>
+
+namespace corsika::process {
+  namespace particle_cut {
+    class ParticleCut : public process::SecondariesProcess<ParticleCut> {
+
+      units::si::HEPEnergyType const fECut;
+
+      units::si::HEPEnergyType fEnergy = 0 * units::si::electronvolt;
+      units::si::HEPEnergyType fEmEnergy = 0 * units::si::electronvolt;
+      unsigned int fEmCount = 0;
+      units::si::HEPEnergyType fInvEnergy = 0 * units::si::electronvolt;
+      unsigned int fInvCount = 0;
+
+    public:
+      ParticleCut(const units::si::HEPEnergyType vCut)
+          : fECut(vCut) {}
+
+      bool ParticleIsInvisible(particles::Code) const;
+      EProcessReturn DoSecondaries(corsika::setup::StackView&);
+
+      template <typename TParticle>
+      bool ParticleIsBelowEnergyCut(TParticle const&) const;
+
+      bool ParticleIsEmParticle(particles::Code) const;
+
+      void Init();
+      void ShowResults();
+
+      units::si::HEPEnergyType GetInvEnergy() const { return fInvEnergy; }
+      units::si::HEPEnergyType GetCutEnergy() const { return fEnergy; }
+      units::si::HEPEnergyType GetEmEnergy() const { return fEmEnergy; }
+      unsigned int GetNumberEmParticles() const { return fEmCount; }
+      unsigned int GetNumberInvParticles() const { return fInvCount; }
+    };
+  } // namespace particle_cut
+} // namespace corsika::process
+
+#endif
diff --git a/Processes/OnShellCheck/testOnShellCheck.cc b/Processes/OnShellCheck/testOnShellCheck.cc
new file mode 100644
index 000000000..4337446a9
--- /dev/null
+++ b/Processes/OnShellCheck/testOnShellCheck.cc
@@ -0,0 +1,91 @@
+/*
+ * (c) Copyright 2018 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/on_shell_check/OnShellCheck.h>
+
+#include <corsika/environment/Environment.h>
+#include <corsika/geometry/Point.h>
+#include <corsika/geometry/RootCoordinateSystem.h>
+#include <corsika/geometry/Vector.h>
+#include <corsika/geometry/FourVector.h>
+#include <corsika/units/PhysicalUnits.h>
+#include <corsika/utl/CorsikaFenv.h>
+
+#include <corsika/setup/SetupStack.h>
+
+#include <catch2/catch.hpp>
+
+using namespace corsika;
+using namespace corsika::process::on_shell_check;
+using namespace corsika::units;
+using namespace corsika::units::si;
+
+TEST_CASE("OnShellCheck", "[processes]") {
+  feenableexcept(FE_INVALID);
+  using EnvType = environment::Environment<setup::IEnvironmentModel>;
+  EnvType env;
+  const geometry::CoordinateSystem& rootCS = env.GetCoordinateSystem();
+
+  // setup empty particle stack
+  setup::Stack stack;
+  stack.Clear();
+  // two energies
+  const HEPEnergyType E = 10_GeV;
+  // list of arbitrary particles
+  std::array<particles::Code, 2> particleList = {
+      particles::Code::PiPlus,   particles::Code::PiMinus};
+
+  std::array<double, 2> mass_shifts = {1.1, 1.001};
+
+  SECTION("check particle masses") {
+
+    OnShellCheck check(1.e-2, 0.01);
+
+    check.Init();
+    
+    // add primary particle to stack
+    auto particle = stack.AddParticle(
+        std::tuple<particles::Code, units::si::HEPEnergyType,
+                   corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
+            particles::Code::Proton, E,
+            corsika::stack::MomentumVector(rootCS, {0_GeV, 0_GeV, 0_GeV}),
+            geometry::Point(rootCS, 0_m, 0_m, 0_m), 0_ns});
+    // view on secondary particles
+    corsika::stack::SecondaryView view(particle);
+    // ref. to primary particle through the secondary view.
+    // only this way the secondary view is populated
+    auto projectile = view.GetProjectile();
+    // add secondaries, all with energies above the threshold
+    // only cut is by species
+    int count = -1;
+    for (auto proType : particleList) {
+      count++;
+      const auto pz = sqrt((E - particles::GetMass(proType)*mass_shifts[count]) *
+                           (E + particles::GetMass(proType)*mass_shifts[count]));
+      auto const momentum = corsika::stack::MomentumVector(rootCS, {0_GeV, 0_GeV, pz});
+      projectile.AddSecondary(std::tuple<particles::Code, units::si::HEPEnergyType,
+                                         corsika::stack::MomentumVector, geometry::Point,
+                                         units::si::TimeType>{
+          proType, E, momentum, geometry::Point(rootCS, 0_m, 0_m, 0_m), 0_ns});
+
+      check.DoSecondaries(view);
+    }
+    int i = -1;
+    for ( auto& p : view) {
+      i++;
+      auto const Plab = corsika::geometry::FourVector(p.GetEnergy(), p.GetMomentum());
+      auto const m_kinetic = Plab.GetNorm();
+      if(i==0)
+	REQUIRE(m_kinetic / particles::PiPlus::GetMass() == Approx(1));
+      else
+	REQUIRE_FALSE(m_kinetic / particles::PiMinus::GetMass() == Approx(1));
+    }
+  }
+}
-- 
GitLab


From 8a8051b2da5e1799103c92e229b62e12c61d6978 Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Thu, 4 Jun 2020 08:55:22 +0100
Subject: [PATCH 18/31] output format etc

---
 Processes/OnShellCheck/OnShellCheck.cc | 21 ++++++++++++---------
 1 file changed, 12 insertions(+), 9 deletions(-)

diff --git a/Processes/OnShellCheck/OnShellCheck.cc b/Processes/OnShellCheck/OnShellCheck.cc
index 8abfc2011..55a389421 100644
--- a/Processes/OnShellCheck/OnShellCheck.cc
+++ b/Processes/OnShellCheck/OnShellCheck.cc
@@ -40,8 +40,8 @@ namespace corsika::process {
         auto const Plab = corsika::geometry::FourVector(e_original, p_original);
         auto const m_kinetic = Plab.GetNorm();
         auto const m_corsika = particles::GetMass(pid);
-        auto const m_err = abs(m_kinetic - m_corsika) / m_corsika;
-        if (m_err > mass_tolerance_) {
+        auto const m_err_abs = abs(m_kinetic - m_corsika);
+        if (m_err_abs >= mass_tolerance_ * m_corsika) {
           const HEPEnergyType e_shifted =
               sqrt(p_original.GetSquaredNorm() + m_corsika * m_corsika);
           auto const e_shift_relative = (e_shifted / e_original - 1);
@@ -50,20 +50,23 @@ namespace corsika::process {
 	     we could promote this to an error.
 	   */
 	  if (abs(e_shift_relative) > energy_tolerance_) {
-	    std::cout << "OnShellCheck::DoSecondaries: warning! shifted particle energy by "
-		      << e_shift_relative*100 << " %" << std::endl;
+	    std::cout << "OnShellCheck: warning! shifted particle energy by "
+		      << e_shift_relative*100 << " %" << std::endl;	    
 	  }
-          std::cout << "OnShellCheck::DoSecondaries: shift particle mass for " << pid
+          std::cout << "OnShellCheck: shift particle mass for " << pid
                     << std::endl
-                    << std::setw(20) << std::setfill(' ')
+                    << std::setw(35) << std::setfill(' ')
                     << "corsika mass (GeV): " << m_corsika / 1_GeV << std::endl
+		    << std::setw(35) << std::setfill(' ')
                     << "kinetic mass (GeV): " << m_kinetic / 1_GeV << std::endl
-                    << "(m_kin-m_cor)/en: " << m_err << std::endl
-                    << "mass tolerance: " << mass_tolerance_ << std::endl;
+		    << std::setw(35) << std::setfill(' ')
+                    << "m_kin-m_cor (GeV): " << m_err_abs / 1_GeV << std::endl
+		    << std::setw(35) << std::setfill(' ')
+                    << "mass tolerance (GeV): " << ( m_corsika * mass_tolerance_ ) / 1_GeV << std::endl;
           // reset energy
           p.SetEnergy(e_shifted);
         } else
-	  std::cout << "OnShellCheck::DoSecondaries: particle mass for " << pid << " OK" << std::endl;
+          std::cout << "OnShellCheck: particle mass for " << pid << " OK" << std::endl;
       }
       return EProcessReturn::eOk;
     }
-- 
GitLab


From 6830f3febb6b49c2bc76d7a91b98e64bd091e992 Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Thu, 4 Jun 2020 08:55:39 +0100
Subject: [PATCH 19/31] fixed test

---
 Processes/OnShellCheck/testOnShellCheck.cc | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/Processes/OnShellCheck/testOnShellCheck.cc b/Processes/OnShellCheck/testOnShellCheck.cc
index 4337446a9..62539e36a 100644
--- a/Processes/OnShellCheck/testOnShellCheck.cc
+++ b/Processes/OnShellCheck/testOnShellCheck.cc
@@ -73,10 +73,9 @@ TEST_CASE("OnShellCheck", "[processes]") {
       projectile.AddSecondary(std::tuple<particles::Code, units::si::HEPEnergyType,
                                          corsika::stack::MomentumVector, geometry::Point,
                                          units::si::TimeType>{
-          proType, E, momentum, geometry::Point(rootCS, 0_m, 0_m, 0_m), 0_ns});
-
-      check.DoSecondaries(view);
+          proType, E, momentum, geometry::Point(rootCS, 0_m, 0_m, 0_m), 0_ns});      
     }
+    check.DoSecondaries(view);
     int i = -1;
     for ( auto& p : view) {
       i++;
-- 
GitLab


From 859037596a7dde7a74a44e85c5ae251d2ec8fb1a Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Thu, 4 Jun 2020 09:05:51 +0100
Subject: [PATCH 20/31] set tolerance in vertical eas

---
 Documentation/Examples/vertical_EAS.cc | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/Documentation/Examples/vertical_EAS.cc b/Documentation/Examples/vertical_EAS.cc
index 974398dfe..7c70e6997 100644
--- a/Documentation/Examples/vertical_EAS.cc
+++ b/Documentation/Examples/vertical_EAS.cc
@@ -184,8 +184,8 @@ int main(int argc, char** argv) {
 
   process::particle_cut::ParticleCut cut{60_GeV};
 
-  process::on_shell_check::OnShellCheck reset_particle_mass(1.e-2,1.e-2);
-  
+  process::on_shell_check::OnShellCheck reset_particle_mass(1.e-3, 1.e-2);
+
   process::energy_loss::EnergyLoss eLoss(showerAxis);
   process::longitudinal_profile::LongitudinalProfile longprof{showerAxis};
 
-- 
GitLab


From 392d158e4160430e0a322d75a6987e594b2238d0 Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Thu, 4 Jun 2020 09:08:11 +0100
Subject: [PATCH 21/31] removed check from sibyll interface

---
 Processes/Sibyll/Interaction.cc | 23 -----------------------
 1 file changed, 23 deletions(-)

diff --git a/Processes/Sibyll/Interaction.cc b/Processes/Sibyll/Interaction.cc
index 0a4c7d59d..1a6e65598 100644
--- a/Processes/Sibyll/Interaction.cc
+++ b/Processes/Sibyll/Interaction.cc
@@ -317,29 +317,6 @@ namespace corsika::process::sibyll {
           auto const Plab = boost.fromCoM(FourVector(eCoM, pCoM));
 
           auto const pid = process::sibyll::ConvertFromSibyll(psib.GetPID());
-          // check if on-shell in corsika
-          auto const m_kinetic = Plab.GetNorm();
-          auto const m_corsika = particles::GetMass(pid);
-          auto const e_corsika = Plab.GetTimeLikeComponent();
-          auto const m_sibyll = corsika::process::sibyll::GetSibyllMass(pid);
-          auto const m_err = abs(m_kinetic - m_corsika) / m_corsika;
-          if (m_err > 1.e-5 && false) {
-            const HEPEnergyType e_shift_corsika = sqrt(
-                Plab.GetSpaceLikeComponents().GetSquaredNorm() + m_corsika * m_corsika);
-            auto const e_shift_relative = (e_shift_corsika / e_corsika - 1) * 100;
-            // warn about percent level shifts in particle energy
-            if (abs(e_shift_relative) > 1) {
-              std::cout << "Sibyll::Interaction: shifted particle energy by "
-                        << e_shift_relative << " %" << std::endl;
-
-              std::cout << "shift particle mass for " << pid << std::endl
-                        << "corsika mass (GeV): " << m_corsika / 1_GeV << std::endl
-                        << "kinetic mass (GeV): " << m_kinetic / 1_GeV << std::endl
-                        << "sibyll  mass (GeV): " << m_sibyll / 1_GeV << std::endl
-                        << "(m_kin-m_cor)/en: " << m_err << std::endl;
-            }
-            Plab.GetTimeLikeComponent() = e_shift_corsika;
-          }
 
           // add to corsika stack
           auto pnew = vP.AddSecondary(
-- 
GitLab


From 31742a07ce51613222ef06f65d9e829329271c2e Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Thu, 4 Jun 2020 09:42:24 +0100
Subject: [PATCH 22/31] clang

---
 Documentation/Examples/vertical_EAS.cc     |  2 +-
 Processes/OnShellCheck/OnShellCheck.cc     | 41 +++++++++++-----------
 Processes/OnShellCheck/testOnShellCheck.cc | 22 ++++++------
 3 files changed, 33 insertions(+), 32 deletions(-)

diff --git a/Documentation/Examples/vertical_EAS.cc b/Documentation/Examples/vertical_EAS.cc
index 7c70e6997..563e802cc 100644
--- a/Documentation/Examples/vertical_EAS.cc
+++ b/Documentation/Examples/vertical_EAS.cc
@@ -22,8 +22,8 @@
 #include <corsika/process/interaction_counter/InteractionCounter.h>
 #include <corsika/process/longitudinal_profile/LongitudinalProfile.h>
 #include <corsika/process/observation_plane/ObservationPlane.h>
-#include <corsika/process/particle_cut/ParticleCut.h>
 #include <corsika/process/on_shell_check/OnShellCheck.h>
+#include <corsika/process/particle_cut/ParticleCut.h>
 #include <corsika/process/pythia/Decay.h>
 #include <corsika/process/sibyll/Decay.h>
 #include <corsika/process/sibyll/Interaction.h>
diff --git a/Processes/OnShellCheck/OnShellCheck.cc b/Processes/OnShellCheck/OnShellCheck.cc
index 55a389421..beee3a70e 100644
--- a/Processes/OnShellCheck/OnShellCheck.cc
+++ b/Processes/OnShellCheck/OnShellCheck.cc
@@ -9,8 +9,8 @@
  * the license.
  */
 
-#include <corsika/process/on_shell_check/OnShellCheck.h>
 #include <corsika/geometry/FourVector.h>
+#include <corsika/process/on_shell_check/OnShellCheck.h>
 
 using namespace std;
 
@@ -23,7 +23,7 @@ using namespace corsika::setup;
 namespace corsika::process {
   namespace on_shell_check {
 
-    void OnShellCheck::Init(){
+    void OnShellCheck::Init() {
       std::cout << "OnShellCheck: mass tolerance is set to " << mass_tolerance_ * 100
                 << "%" << endl
                 << "              energy tolerance is set to " << energy_tolerance_ * 100
@@ -33,8 +33,9 @@ namespace corsika::process {
     EProcessReturn OnShellCheck::DoSecondaries(corsika::setup::StackView& vS) {
       for (auto& p : vS) {
         auto const pid = p.GetPID();
-	//if(pid==particles::Code::Gamma || particles::IsNeutrino(pid) || particles::IsNucleus(pid)) continue;
-	if(!particles::IsHadron(pid)) continue;
+        // if(pid==particles::Code::Gamma || particles::IsNeutrino(pid) ||
+        // particles::IsNucleus(pid)) continue;
+        if (!particles::IsHadron(pid)) continue;
         auto const e_original = p.GetEnergy();
         auto const p_original = p.GetMomentum();
         auto const Plab = corsika::geometry::FourVector(e_original, p_original);
@@ -45,24 +46,24 @@ namespace corsika::process {
           const HEPEnergyType e_shifted =
               sqrt(p_original.GetSquaredNorm() + m_corsika * m_corsika);
           auto const e_shift_relative = (e_shifted / e_original - 1);
-	  /* 
-	     For now we warn if the necessary shift is larger than 1%.
-	     we could promote this to an error.
-	   */
-	  if (abs(e_shift_relative) > energy_tolerance_) {
-	    std::cout << "OnShellCheck: warning! shifted particle energy by "
-		      << e_shift_relative*100 << " %" << std::endl;	    
-	  }
-          std::cout << "OnShellCheck: shift particle mass for " << pid
-                    << std::endl
+          /*
+             For now we warn if the necessary shift is larger than 1%.
+             we could promote this to an error.
+           */
+          if (abs(e_shift_relative) > energy_tolerance_) {
+            std::cout << "OnShellCheck: warning! shifted particle energy by "
+                      << e_shift_relative * 100 << " %" << std::endl;
+          }
+          std::cout << "OnShellCheck: shift particle mass for " << pid << std::endl
                     << std::setw(35) << std::setfill(' ')
                     << "corsika mass (GeV): " << m_corsika / 1_GeV << std::endl
-		    << std::setw(35) << std::setfill(' ')
+                    << std::setw(35) << std::setfill(' ')
                     << "kinetic mass (GeV): " << m_kinetic / 1_GeV << std::endl
-		    << std::setw(35) << std::setfill(' ')
+                    << std::setw(35) << std::setfill(' ')
                     << "m_kin-m_cor (GeV): " << m_err_abs / 1_GeV << std::endl
-		    << std::setw(35) << std::setfill(' ')
-                    << "mass tolerance (GeV): " << ( m_corsika * mass_tolerance_ ) / 1_GeV << std::endl;
+                    << std::setw(35) << std::setfill(' ')
+                    << "mass tolerance (GeV): " << (m_corsika * mass_tolerance_) / 1_GeV
+                    << std::endl;
           // reset energy
           p.SetEnergy(e_shifted);
         } else
@@ -70,5 +71,5 @@ namespace corsika::process {
       }
       return EProcessReturn::eOk;
     }
-  }
-} // namespace on_shell_check
+  } // namespace on_shell_check
+} // namespace corsika::process
diff --git a/Processes/OnShellCheck/testOnShellCheck.cc b/Processes/OnShellCheck/testOnShellCheck.cc
index 62539e36a..bfe938305 100644
--- a/Processes/OnShellCheck/testOnShellCheck.cc
+++ b/Processes/OnShellCheck/testOnShellCheck.cc
@@ -11,10 +11,10 @@
 #include <corsika/process/on_shell_check/OnShellCheck.h>
 
 #include <corsika/environment/Environment.h>
+#include <corsika/geometry/FourVector.h>
 #include <corsika/geometry/Point.h>
 #include <corsika/geometry/RootCoordinateSystem.h>
 #include <corsika/geometry/Vector.h>
-#include <corsika/geometry/FourVector.h>
 #include <corsika/units/PhysicalUnits.h>
 #include <corsika/utl/CorsikaFenv.h>
 
@@ -39,8 +39,8 @@ TEST_CASE("OnShellCheck", "[processes]") {
   // two energies
   const HEPEnergyType E = 10_GeV;
   // list of arbitrary particles
-  std::array<particles::Code, 2> particleList = {
-      particles::Code::PiPlus,   particles::Code::PiMinus};
+  std::array<particles::Code, 2> particleList = {particles::Code::PiPlus,
+                                                 particles::Code::PiMinus};
 
   std::array<double, 2> mass_shifts = {1.1, 1.001};
 
@@ -49,7 +49,7 @@ TEST_CASE("OnShellCheck", "[processes]") {
     OnShellCheck check(1.e-2, 0.01);
 
     check.Init();
-    
+
     // add primary particle to stack
     auto particle = stack.AddParticle(
         std::tuple<particles::Code, units::si::HEPEnergyType,
@@ -67,24 +67,24 @@ TEST_CASE("OnShellCheck", "[processes]") {
     int count = -1;
     for (auto proType : particleList) {
       count++;
-      const auto pz = sqrt((E - particles::GetMass(proType)*mass_shifts[count]) *
-                           (E + particles::GetMass(proType)*mass_shifts[count]));
+      const auto pz = sqrt((E - particles::GetMass(proType) * mass_shifts[count]) *
+                           (E + particles::GetMass(proType) * mass_shifts[count]));
       auto const momentum = corsika::stack::MomentumVector(rootCS, {0_GeV, 0_GeV, pz});
       projectile.AddSecondary(std::tuple<particles::Code, units::si::HEPEnergyType,
                                          corsika::stack::MomentumVector, geometry::Point,
                                          units::si::TimeType>{
-          proType, E, momentum, geometry::Point(rootCS, 0_m, 0_m, 0_m), 0_ns});      
+          proType, E, momentum, geometry::Point(rootCS, 0_m, 0_m, 0_m), 0_ns});
     }
     check.DoSecondaries(view);
     int i = -1;
-    for ( auto& p : view) {
+    for (auto& p : view) {
       i++;
       auto const Plab = corsika::geometry::FourVector(p.GetEnergy(), p.GetMomentum());
       auto const m_kinetic = Plab.GetNorm();
-      if(i==0)
-	REQUIRE(m_kinetic / particles::PiPlus::GetMass() == Approx(1));
+      if (i == 0)
+        REQUIRE(m_kinetic / particles::PiPlus::GetMass() == Approx(1));
       else
-	REQUIRE_FALSE(m_kinetic / particles::PiMinus::GetMass() == Approx(1));
+        REQUIRE_FALSE(m_kinetic / particles::PiMinus::GetMass() == Approx(1));
     }
   }
 }
-- 
GitLab


From e31665ddac855e3b79bf4979727c9a11971d92ef Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Thu, 4 Jun 2020 10:50:54 +0100
Subject: [PATCH 23/31] removed obsolete file

---
 Processes/OnShellCheck/ParticleCut.h | 55 ----------------------------
 1 file changed, 55 deletions(-)
 delete mode 100644 Processes/OnShellCheck/ParticleCut.h

diff --git a/Processes/OnShellCheck/ParticleCut.h b/Processes/OnShellCheck/ParticleCut.h
deleted file mode 100644
index 739a1919e..000000000
--- a/Processes/OnShellCheck/ParticleCut.h
+++ /dev/null
@@ -1,55 +0,0 @@
-/*
- * (c) Copyright 2018 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 _corsika_process_particle_cut_ParticleCut_h_
-#define _corsika_process_particle_cut_ParticleCut_h_
-
-#include <corsika/particles/ParticleProperties.h>
-#include <corsika/process/SecondariesProcess.h>
-#include <corsika/setup/SetupStack.h>
-#include <corsika/units/PhysicalUnits.h>
-
-namespace corsika::process {
-  namespace particle_cut {
-    class ParticleCut : public process::SecondariesProcess<ParticleCut> {
-
-      units::si::HEPEnergyType const fECut;
-
-      units::si::HEPEnergyType fEnergy = 0 * units::si::electronvolt;
-      units::si::HEPEnergyType fEmEnergy = 0 * units::si::electronvolt;
-      unsigned int fEmCount = 0;
-      units::si::HEPEnergyType fInvEnergy = 0 * units::si::electronvolt;
-      unsigned int fInvCount = 0;
-
-    public:
-      ParticleCut(const units::si::HEPEnergyType vCut)
-          : fECut(vCut) {}
-
-      bool ParticleIsInvisible(particles::Code) const;
-      EProcessReturn DoSecondaries(corsika::setup::StackView&);
-
-      template <typename TParticle>
-      bool ParticleIsBelowEnergyCut(TParticle const&) const;
-
-      bool ParticleIsEmParticle(particles::Code) const;
-
-      void Init();
-      void ShowResults();
-
-      units::si::HEPEnergyType GetInvEnergy() const { return fInvEnergy; }
-      units::si::HEPEnergyType GetCutEnergy() const { return fEnergy; }
-      units::si::HEPEnergyType GetEmEnergy() const { return fEmEnergy; }
-      unsigned int GetNumberEmParticles() const { return fEmCount; }
-      unsigned int GetNumberInvParticles() const { return fInvCount; }
-    };
-  } // namespace particle_cut
-} // namespace corsika::process
-
-#endif
-- 
GitLab


From c39c73152f4598148b88089f1d010e81bb3cd1bc Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Thu, 4 Jun 2020 11:25:23 +0100
Subject: [PATCH 24/31] avoid test for nuclei

---
 Processes/OnShellCheck/OnShellCheck.cc     |  2 +-
 Processes/OnShellCheck/testOnShellCheck.cc | 12 ++++++++----
 2 files changed, 9 insertions(+), 5 deletions(-)

diff --git a/Processes/OnShellCheck/OnShellCheck.cc b/Processes/OnShellCheck/OnShellCheck.cc
index beee3a70e..a5662e434 100644
--- a/Processes/OnShellCheck/OnShellCheck.cc
+++ b/Processes/OnShellCheck/OnShellCheck.cc
@@ -35,7 +35,7 @@ namespace corsika::process {
         auto const pid = p.GetPID();
         // if(pid==particles::Code::Gamma || particles::IsNeutrino(pid) ||
         // particles::IsNucleus(pid)) continue;
-        if (!particles::IsHadron(pid)) continue;
+        if (!particles::IsHadron(pid) || particles::IsNucleus(pid)) continue;
         auto const e_original = p.GetEnergy();
         auto const p_original = p.GetMomentum();
         auto const Plab = corsika::geometry::FourVector(e_original, p_original);
diff --git a/Processes/OnShellCheck/testOnShellCheck.cc b/Processes/OnShellCheck/testOnShellCheck.cc
index bfe938305..08e6cc323 100644
--- a/Processes/OnShellCheck/testOnShellCheck.cc
+++ b/Processes/OnShellCheck/testOnShellCheck.cc
@@ -39,10 +39,13 @@ TEST_CASE("OnShellCheck", "[processes]") {
   // two energies
   const HEPEnergyType E = 10_GeV;
   // list of arbitrary particles
-  std::array<particles::Code, 2> particleList = {particles::Code::PiPlus,
-                                                 particles::Code::PiMinus};
+  std::array<particles::Code, 4> particleList = {
+      particles::Code::PiPlus,
+      particles::Code::PiMinus,
+      particles::Code::Helium,
+      particles::Code::Gamma};
 
-  std::array<double, 2> mass_shifts = {1.1, 1.001};
+  std::array<double, 4> mass_shifts = {1.1, 1.001, 1.0, 1.0};
 
   SECTION("check particle masses") {
 
@@ -83,8 +86,9 @@ TEST_CASE("OnShellCheck", "[processes]") {
       auto const m_kinetic = Plab.GetNorm();
       if (i == 0)
         REQUIRE(m_kinetic / particles::PiPlus::GetMass() == Approx(1));
-      else
+      else if (i == 1)
         REQUIRE_FALSE(m_kinetic / particles::PiMinus::GetMass() == Approx(1));
+	
     }
   }
 }
-- 
GitLab


From 5fa15a5d385311432caafaecf58694574397c8e4 Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Thu, 4 Jun 2020 11:26:02 +0100
Subject: [PATCH 25/31] clang

---
 Processes/OnShellCheck/testOnShellCheck.cc | 7 ++-----
 1 file changed, 2 insertions(+), 5 deletions(-)

diff --git a/Processes/OnShellCheck/testOnShellCheck.cc b/Processes/OnShellCheck/testOnShellCheck.cc
index 08e6cc323..bd366d9b5 100644
--- a/Processes/OnShellCheck/testOnShellCheck.cc
+++ b/Processes/OnShellCheck/testOnShellCheck.cc
@@ -40,9 +40,7 @@ TEST_CASE("OnShellCheck", "[processes]") {
   const HEPEnergyType E = 10_GeV;
   // list of arbitrary particles
   std::array<particles::Code, 4> particleList = {
-      particles::Code::PiPlus,
-      particles::Code::PiMinus,
-      particles::Code::Helium,
+      particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::Helium,
       particles::Code::Gamma};
 
   std::array<double, 4> mass_shifts = {1.1, 1.001, 1.0, 1.0};
@@ -61,7 +59,7 @@ TEST_CASE("OnShellCheck", "[processes]") {
             corsika::stack::MomentumVector(rootCS, {0_GeV, 0_GeV, 0_GeV}),
             geometry::Point(rootCS, 0_m, 0_m, 0_m), 0_ns});
     // view on secondary particles
-    corsika::stack::SecondaryView view(particle);
+    corsika::stack::SecondaryView view(part./ icle);
     // ref. to primary particle through the secondary view.
     // only this way the secondary view is populated
     auto projectile = view.GetProjectile();
@@ -88,7 +86,6 @@ TEST_CASE("OnShellCheck", "[processes]") {
         REQUIRE(m_kinetic / particles::PiPlus::GetMass() == Approx(1));
       else if (i == 1)
         REQUIRE_FALSE(m_kinetic / particles::PiMinus::GetMass() == Approx(1));
-	
     }
   }
 }
-- 
GitLab


From 288a94a0bf1618e8dfc2d2703bae24d7a2a5f287 Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Thu, 4 Jun 2020 11:27:20 +0100
Subject: [PATCH 26/31] there you go

---
 Processes/OnShellCheck/testOnShellCheck.cc | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Processes/OnShellCheck/testOnShellCheck.cc b/Processes/OnShellCheck/testOnShellCheck.cc
index bd366d9b5..9f50e0b37 100644
--- a/Processes/OnShellCheck/testOnShellCheck.cc
+++ b/Processes/OnShellCheck/testOnShellCheck.cc
@@ -59,7 +59,7 @@ TEST_CASE("OnShellCheck", "[processes]") {
             corsika::stack::MomentumVector(rootCS, {0_GeV, 0_GeV, 0_GeV}),
             geometry::Point(rootCS, 0_m, 0_m, 0_m), 0_ns});
     // view on secondary particles
-    corsika::stack::SecondaryView view(part./ icle);
+    corsika::stack::SecondaryView view(particle);
     // ref. to primary particle through the secondary view.
     // only this way the secondary view is populated
     auto projectile = view.GetProjectile();
-- 
GitLab


From cd6fbba275868017990549f60ee689a6302e5693 Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Fri, 12 Jun 2020 13:06:04 +0100
Subject: [PATCH 27/31] added summary

---
 Processes/OnShellCheck/OnShellCheck.cc |  3 +++
 Processes/OnShellCheck/OnShellCheck.h  | 12 ++++++++++++
 2 files changed, 15 insertions(+)

diff --git a/Processes/OnShellCheck/OnShellCheck.cc b/Processes/OnShellCheck/OnShellCheck.cc
index a5662e434..133af10c6 100644
--- a/Processes/OnShellCheck/OnShellCheck.cc
+++ b/Processes/OnShellCheck/OnShellCheck.cc
@@ -46,6 +46,9 @@ namespace corsika::process {
           const HEPEnergyType e_shifted =
               sqrt(p_original.GetSquaredNorm() + m_corsika * m_corsika);
           auto const e_shift_relative = (e_shifted / e_original - 1);
+          count_ = count_ + 1;
+          average_shift_ += abs(e_shift_relative);
+          if (abs(e_shift_relative) > max_shift_) max_shift_ = abs(e_shift_relative);
           /*
              For now we warn if the necessary shift is larger than 1%.
              we could promote this to an error.
diff --git a/Processes/OnShellCheck/OnShellCheck.h b/Processes/OnShellCheck/OnShellCheck.h
index f9a047324..f1c588be0 100644
--- a/Processes/OnShellCheck/OnShellCheck.h
+++ b/Processes/OnShellCheck/OnShellCheck.h
@@ -19,12 +19,24 @@
 namespace corsika::process {
   namespace on_shell_check {
     class OnShellCheck : public process::SecondariesProcess<OnShellCheck> {
+      double average_shift_ = 0;
+      double max_shift_ = 0;
+      double count_ = 0;
 
     public:
       OnShellCheck(const double vMassTolerance, const double vEnergyTolerance)
           : mass_tolerance_(vMassTolerance)
           , energy_tolerance_(vEnergyTolerance) {}
 
+      ~OnShellCheck() {
+        std::cout << "OnShellCheck: summary" << std::endl
+                  << " particles shifted: " << int(count_) << std::endl;
+        if (count_)
+          std::cout << " average energy shift (%): " << average_shift_ / count_ * 100.
+                    << std::endl
+                    << " max. energy shift (%): " << max_shift_ * 100. << std::endl;
+      };
+
       EProcessReturn DoSecondaries(corsika::setup::StackView&);
 
       void Init();
-- 
GitLab


From 1034449036d421433167f0da30a0dbe57a7a0cb1 Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Fri, 12 Jun 2020 15:04:40 +0100
Subject: [PATCH 28/31] added error for extreme cases

---
 Documentation/Examples/vertical_EAS.cc     |  2 +-
 Processes/OnShellCheck/OnShellCheck.cc     | 31 +++++++++++-----------
 Processes/OnShellCheck/OnShellCheck.h      |  7 +++--
 Processes/OnShellCheck/testOnShellCheck.cc |  2 +-
 4 files changed, 23 insertions(+), 19 deletions(-)

diff --git a/Documentation/Examples/vertical_EAS.cc b/Documentation/Examples/vertical_EAS.cc
index 563e802cc..5925841bf 100644
--- a/Documentation/Examples/vertical_EAS.cc
+++ b/Documentation/Examples/vertical_EAS.cc
@@ -184,7 +184,7 @@ int main(int argc, char** argv) {
 
   process::particle_cut::ParticleCut cut{60_GeV};
 
-  process::on_shell_check::OnShellCheck reset_particle_mass(1.e-3, 1.e-2);
+  process::on_shell_check::OnShellCheck reset_particle_mass(1.e-3, 1.e-1, false);
 
   process::energy_loss::EnergyLoss eLoss(showerAxis);
   process::longitudinal_profile::LongitudinalProfile longprof{showerAxis};
diff --git a/Processes/OnShellCheck/OnShellCheck.cc b/Processes/OnShellCheck/OnShellCheck.cc
index 133af10c6..f47f94d25 100644
--- a/Processes/OnShellCheck/OnShellCheck.cc
+++ b/Processes/OnShellCheck/OnShellCheck.cc
@@ -1,4 +1,3 @@
-
 /*
  * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
  *
@@ -33,8 +32,6 @@ namespace corsika::process {
     EProcessReturn OnShellCheck::DoSecondaries(corsika::setup::StackView& vS) {
       for (auto& p : vS) {
         auto const pid = p.GetPID();
-        // if(pid==particles::Code::Gamma || particles::IsNeutrino(pid) ||
-        // particles::IsNucleus(pid)) continue;
         if (!particles::IsHadron(pid) || particles::IsNucleus(pid)) continue;
         auto const e_original = p.GetEnergy();
         auto const p_original = p.GetMomentum();
@@ -49,24 +46,28 @@ namespace corsika::process {
           count_ = count_ + 1;
           average_shift_ += abs(e_shift_relative);
           if (abs(e_shift_relative) > max_shift_) max_shift_ = abs(e_shift_relative);
-          /*
-             For now we warn if the necessary shift is larger than 1%.
-             we could promote this to an error.
-           */
-          if (abs(e_shift_relative) > energy_tolerance_) {
-            std::cout << "OnShellCheck: warning! shifted particle energy by "
-                      << e_shift_relative * 100 << " %" << std::endl;
-          }
           std::cout << "OnShellCheck: shift particle mass for " << pid << std::endl
-                    << std::setw(35) << std::setfill(' ')
+                    << std::setw(40) << std::setfill(' ')
                     << "corsika mass (GeV): " << m_corsika / 1_GeV << std::endl
-                    << std::setw(35) << std::setfill(' ')
+                    << std::setw(40) << std::setfill(' ')
                     << "kinetic mass (GeV): " << m_kinetic / 1_GeV << std::endl
-                    << std::setw(35) << std::setfill(' ')
+                    << std::setw(40) << std::setfill(' ')
                     << "m_kin-m_cor (GeV): " << m_err_abs / 1_GeV << std::endl
-                    << std::setw(35) << std::setfill(' ')
+                    << std::setw(40) << std::setfill(' ')
                     << "mass tolerance (GeV): " << (m_corsika * mass_tolerance_) / 1_GeV
                     << std::endl;
+          /*
+            For now we warn if the necessary shift is larger than 1%.
+            we could promote this to an error.
+          */
+          if (abs(e_shift_relative) > energy_tolerance_) {
+            std::cout << "OnShellCheck: warning! shifted particle energy by "
+                      << e_shift_relative * 100 << " %" << std::endl;
+            if (throw_error_)
+              throw std::runtime_error(
+                  "OnShellCheck: error! shifted energy by large amount!");
+          }
+
           // reset energy
           p.SetEnergy(e_shifted);
         } else
diff --git a/Processes/OnShellCheck/OnShellCheck.h b/Processes/OnShellCheck/OnShellCheck.h
index f1c588be0..d51d5e8f7 100644
--- a/Processes/OnShellCheck/OnShellCheck.h
+++ b/Processes/OnShellCheck/OnShellCheck.h
@@ -24,9 +24,11 @@ namespace corsika::process {
       double count_ = 0;
 
     public:
-      OnShellCheck(const double vMassTolerance, const double vEnergyTolerance)
+      OnShellCheck(const double vMassTolerance, const double vEnergyTolerance,
+                   const bool vError)
           : mass_tolerance_(vMassTolerance)
-          , energy_tolerance_(vEnergyTolerance) {}
+          , energy_tolerance_(vEnergyTolerance)
+          , throw_error_(vError) {}
 
       ~OnShellCheck() {
         std::cout << "OnShellCheck: summary" << std::endl
@@ -44,6 +46,7 @@ namespace corsika::process {
     private:
       double mass_tolerance_;
       double energy_tolerance_;
+      bool throw_error_;
     };
   } // namespace on_shell_check
 } // namespace corsika::process
diff --git a/Processes/OnShellCheck/testOnShellCheck.cc b/Processes/OnShellCheck/testOnShellCheck.cc
index 9f50e0b37..f67272bd2 100644
--- a/Processes/OnShellCheck/testOnShellCheck.cc
+++ b/Processes/OnShellCheck/testOnShellCheck.cc
@@ -47,7 +47,7 @@ TEST_CASE("OnShellCheck", "[processes]") {
 
   SECTION("check particle masses") {
 
-    OnShellCheck check(1.e-2, 0.01);
+    OnShellCheck check(1.e-2, 0.01, false);
 
     check.Init();
 
-- 
GitLab


From 3e846b3d05b54143ddb06b1a9579392bb4728614 Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Wed, 15 Jul 2020 15:30:33 +0100
Subject: [PATCH 29/31] fixed merge

---
 Documentation/Examples/vertical_EAS.cc | 2 --
 1 file changed, 2 deletions(-)

diff --git a/Documentation/Examples/vertical_EAS.cc b/Documentation/Examples/vertical_EAS.cc
index b79022bcf..5925841bf 100644
--- a/Documentation/Examples/vertical_EAS.cc
+++ b/Documentation/Examples/vertical_EAS.cc
@@ -186,8 +186,6 @@ int main(int argc, char** argv) {
 
   process::on_shell_check::OnShellCheck reset_particle_mass(1.e-3, 1.e-1, false);
 
-  process::on_shell_check::OnShellCheck reset_particle_mass(1.e-3, 1.e-1, false);
-
   process::energy_loss::EnergyLoss eLoss(showerAxis);
   process::longitudinal_profile::LongitudinalProfile longprof{showerAxis};
 
-- 
GitLab


From b20b3efee92f0eb3ca7183aa964e0acb8fa98653 Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Wed, 15 Jul 2020 15:33:14 +0100
Subject: [PATCH 30/31] clang

---
 Documentation/Examples/vertical_EAS.cc | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/Documentation/Examples/vertical_EAS.cc b/Documentation/Examples/vertical_EAS.cc
index 5925841bf..c4fd63e69 100644
--- a/Documentation/Examples/vertical_EAS.cc
+++ b/Documentation/Examples/vertical_EAS.cc
@@ -203,8 +203,8 @@ int main(int argc, char** argv) {
                                                        55_GeV);
   auto decaySequence = decayPythia << decaySibyll;
 
-  auto sequence = switchProcess << reset_particle_mass << decaySequence << longprof << eLoss << cut
-                                << observationLevel;
+  auto sequence = switchProcess << reset_particle_mass << decaySequence << longprof
+                                << eLoss << cut << observationLevel;
 
   // define air shower object, run simulation
   tracking_line::TrackingLine tracking;
-- 
GitLab


From 1f6e53e85c069bbda492e27f74c8cfecfd97ca22 Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Thu, 16 Jul 2020 12:15:42 +0100
Subject: [PATCH 31/31] replace REQUIRE with CHECK

---
 Processes/OnShellCheck/testOnShellCheck.cc | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/Processes/OnShellCheck/testOnShellCheck.cc b/Processes/OnShellCheck/testOnShellCheck.cc
index f67272bd2..ca2c56455 100644
--- a/Processes/OnShellCheck/testOnShellCheck.cc
+++ b/Processes/OnShellCheck/testOnShellCheck.cc
@@ -83,9 +83,9 @@ TEST_CASE("OnShellCheck", "[processes]") {
       auto const Plab = corsika::geometry::FourVector(p.GetEnergy(), p.GetMomentum());
       auto const m_kinetic = Plab.GetNorm();
       if (i == 0)
-        REQUIRE(m_kinetic / particles::PiPlus::GetMass() == Approx(1));
+        CHECK(m_kinetic / particles::PiPlus::GetMass() == Approx(1));
       else if (i == 1)
-        REQUIRE_FALSE(m_kinetic / particles::PiMinus::GetMass() == Approx(1));
+        CHECK_FALSE(m_kinetic / particles::PiMinus::GetMass() == Approx(1));
     }
   }
 }
-- 
GitLab