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] 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