added error for extreme cases

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};

Processes/OnShellCheck/OnShellCheck.cc

-
 /*
  * (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

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

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();