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