diff --git a/Documentation/Examples/CMakeLists.txt b/Documentation/Examples/CMakeLists.txt
index fc7c01522f9ce087bb86868fac4c75356d0822c2..0118d7c87713be4f5225ec1e8de13d4a0a32c1f7 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 a81414dc1e373aa21cea0fb4fa1e44a94b9b8336..c4fd63e69fa697c5f3602ad73339d851c1e2e600 100644
--- a/Documentation/Examples/vertical_EAS.cc
+++ b/Documentation/Examples/vertical_EAS.cc
@@ -22,6 +22,7 @@
#include <corsika/process/interaction_counter/InteractionCounter.h>
#include <corsika/process/longitudinal_profile/LongitudinalProfile.h>
#include <corsika/process/observation_plane/ObservationPlane.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>
@@ -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-3, 1.e-1, false);
+
process::energy_loss::EnergyLoss eLoss(showerAxis);
process::longitudinal_profile::LongitudinalProfile longprof{showerAxis};
@@ -199,8 +202,9 @@ 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
- << observationLevel;
+
+ auto sequence = switchProcess << reset_particle_mass << decaySequence << longprof
+ << eLoss << cut << observationLevel;
// define air shower object, run simulation
tracking_line::TrackingLine tracking;
diff --git a/Processes/CMakeLists.txt b/Processes/CMakeLists.txt
index 536077c5bd0487b5b23da6e3bef7224aa281a6b2..f27d027a10fafc1aabbc5bb2790fe4b8262ca62d 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 0000000000000000000000000000000000000000..0d2191ad34fd420b322efdf39c39d024a0f8184d
--- /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 0000000000000000000000000000000000000000..f47f94d258feec91e13a87878306c47f76cc796f
--- /dev/null
+++ b/Processes/OnShellCheck/OnShellCheck.cc
@@ -0,0 +1,79 @@
+/*
+ * (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/geometry/FourVector.h>
+#include <corsika/process/on_shell_check/OnShellCheck.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 (!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);
+ auto const m_kinetic = Plab.GetNorm();
+ auto const m_corsika = particles::GetMass(pid);
+ 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);
+ count_ = count_ + 1;
+ average_shift_ += abs(e_shift_relative);
+ if (abs(e_shift_relative) > max_shift_) max_shift_ = abs(e_shift_relative);
+ std::cout << "OnShellCheck: shift particle mass for " << pid << std::endl
+ << std::setw(40) << std::setfill(' ')
+ << "corsika mass (GeV): " << m_corsika / 1_GeV << std::endl
+ << std::setw(40) << std::setfill(' ')
+ << "kinetic mass (GeV): " << m_kinetic / 1_GeV << std::endl
+ << std::setw(40) << std::setfill(' ')
+ << "m_kin-m_cor (GeV): " << m_err_abs / 1_GeV << std::endl
+ << 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
+ std::cout << "OnShellCheck: particle mass for " << pid << " OK" << std::endl;
+ }
+ return EProcessReturn::eOk;
+ }
+ } // namespace on_shell_check
+} // namespace corsika::process
diff --git a/Processes/OnShellCheck/OnShellCheck.h b/Processes/OnShellCheck/OnShellCheck.h
new file mode 100644
index 0000000000000000000000000000000000000000..d51d5e8f7d204dd37d84c72ba73d8eb66c233e4a
--- /dev/null
+++ b/Processes/OnShellCheck/OnShellCheck.h
@@ -0,0 +1,54 @@
+/*
+ * (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> {
+ double average_shift_ = 0;
+ double max_shift_ = 0;
+ double count_ = 0;
+
+ public:
+ OnShellCheck(const double vMassTolerance, const double vEnergyTolerance,
+ const bool vError)
+ : mass_tolerance_(vMassTolerance)
+ , energy_tolerance_(vEnergyTolerance)
+ , throw_error_(vError) {}
+
+ ~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();
+
+ private:
+ double mass_tolerance_;
+ double energy_tolerance_;
+ bool throw_error_;
+ };
+ } // namespace on_shell_check
+} // namespace corsika::process
+
+#endif
diff --git a/Processes/OnShellCheck/testOnShellCheck.cc b/Processes/OnShellCheck/testOnShellCheck.cc
new file mode 100644
index 0000000000000000000000000000000000000000..ca2c564553b55d078907379c7067881d0071336f
--- /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/FourVector.h>
+#include <corsika/geometry/Point.h>
+#include <corsika/geometry/RootCoordinateSystem.h>
+#include <corsika/geometry/Vector.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, 4> particleList = {
+ 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};
+
+ SECTION("check particle masses") {
+
+ OnShellCheck check(1.e-2, 0.01, false);
+
+ 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)
+ CHECK(m_kinetic / particles::PiPlus::GetMass() == Approx(1));
+ else if (i == 1)
+ CHECK_FALSE(m_kinetic / particles::PiMinus::GetMass() == Approx(1));
+ }
+ }
+}
diff --git a/Processes/Pythia/Decay.cc b/Processes/Pythia/Decay.cc
index 27890d28bfe1ce42576fc81d08550ea80daf70de..a17678d18c9932a1f6417f792ab5cb0e40946f6a 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 916ea2ec0829f6de41b4055428bac17b78f0fb62..1a6e655982d437e0d302698c4e0e5e72ba94fb01 100644
--- a/Processes/Sibyll/Interaction.cc
+++ b/Processes/Sibyll/Interaction.cc
@@ -316,20 +316,27 @@ 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());
+
// 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;
}
}