added process to check particle masses and shift on-shell if necessary

Documentation/Examples/CMakeLists.txt

@@ -99,6 +99,7 @@ if (Pythia8_FOUND)
     ProcessTrackWriter
     ProcessTrackingLine
     ProcessParticleCut
+    ProcessOnShellCheck
     ProcessStackInspector
     CORSIKAprocesses
     CORSIKAcascade

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;

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)

Processes/OnShellCheck/CMakeLists.txt

+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
+  )

Processes/OnShellCheck/OnShellCheck.cc

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

Processes/OnShellCheck/OnShellCheck.h

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

Processes/OnShellCheck/ParticleCut.h

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

Processes/OnShellCheck/testOnShellCheck.cc

+/*
+ * (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));
+    }
+  }
+}