IAP GITLAB

Commit 5842aff8 authored by Ralf Ulrich's avatar Ralf Ulrich

Merge branch '288-follow-up-from-resolve-handling-of-off-shell-particles' into 'master'

Resolve "Follow-up from "Resolve "handling of off-shell particles"""

Closes #288

See merge request !308
parents 18f99dec ea5d5d40
Pipeline #3437 passed with stages
in 23 minutes and 55 seconds
......@@ -10,6 +10,7 @@
#include <corsika/framework/geometry/FourVector.hpp>
#include <corsika/modules/OnShellCheck.hpp>
#include <corsika/framework/core/Logging.hpp>
namespace corsika {
......@@ -18,26 +19,26 @@ namespace corsika {
: mass_tolerance_(vMassTolerance)
, energy_tolerance_(vEnergyTolerance)
, throw_error_(vError) {
std::cout << "OnShellCheck: mass tolerance is set to " << mass_tolerance_ * 100 << "%"
<< std::endl
<< " energy tolerance is set to " << energy_tolerance_ * 100
<< "%" << std::endl;
CORSIKA_LOGGER_DEBUG(logger_, "mass tolerance is set to {:3.2f}%",
mass_tolerance_ * 100);
CORSIKA_LOGGER_DEBUG(logger_, "energy tolerance is set to {:3.2f}%",
energy_tolerance_ * 100);
}
OnShellCheck::~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;
logger_->info(
" summary \n"
" particles shifted: {} \n"
" average energy shift (%): {} \n"
" max. energy shift (%): {} ",
int(count_), (count_ ? average_shift_ / count_ * 100 : 0), max_shift_ * 100.);
}
template <typename TView>
void OnShellCheck::doSecondaries(TView& vS) {
for (auto& p : vS) {
auto const pid = p.getPID();
if (!is_hadron(pid) || is_nucleus(pid)) continue;
if (is_nucleus(pid)) continue;
auto const e_original = p.getEnergy();
auto const p_original = p.getMomentum();
auto const Plab = FourVector(e_original, p_original);
......@@ -51,23 +52,24 @@ namespace corsika {
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;
CORSIKA_LOGGER_TRACE(
logger_,
"shift particle mass for {} \n"
"{:>45} {:7.5f} \n"
"{:>45} {:7.5f} \n"
"{:>45} {:7.5f} \n"
"{:>45} {:7.5f} \n",
pid, "corsika mass (GeV):", m_corsika / 1_GeV,
"kinetic mass (GeV): ", m_kinetic / 1_GeV,
"m_kin-m_cor (GeV): ", m_err_abs / 1_GeV,
"mass tolerance (GeV): ", (m_corsika * mass_tolerance_) / 1_GeV);
/*
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;
logger_->warn("warning! shifted particle energy by {} %",
e_shift_relative * 100);
if (throw_error_)
throw std::runtime_error(
"OnShellCheck: error! shifted energy by large amount!");
......@@ -76,7 +78,7 @@ namespace corsika {
// reset energy
p.setEnergy(e_shifted);
} else
std::cout << "OnShellCheck: particle mass for " << pid << " OK" << std::endl;
CORSIKA_LOGGER_DEBUG(logger_, "particle mass for {} OK", pid);
}
}
......
......@@ -31,7 +31,7 @@ namespace corsika {
double average_shift_ = 0;
double max_shift_ = 0;
double count_ = 0;
std::shared_ptr<spdlog::logger> logger_ = get_logger("corsika_OnShellCheck");
double mass_tolerance_;
double energy_tolerance_;
bool throw_error_;
......
......@@ -24,7 +24,7 @@ using namespace corsika;
TEST_CASE("OnShellCheck", "[processes]") {
logging::set_level(logging::level::info);
logging::set_level(logging::level::debug);
corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
feenableexcept(FE_INVALID);
......@@ -38,9 +38,10 @@ TEST_CASE("OnShellCheck", "[processes]") {
// two energies
const HEPEnergyType E = 10_GeV;
// list of arbitrary particles
std::array const particleList{Code::PiPlus, Code::PiMinus, Code::Helium, Code::Gamma};
std::array const particleList{Code::PiPlus, Code::PiMinus, Code::Helium,
Code::Gamma, Code::Electron, Code::MuPlus};
std::array const mass_shifts{1.1, 1.001, 1.0, 1.0};
std::array const mass_shifts{1.1, 1.001, 1.0, 1.0, 1.01, 1.0};
SECTION("check particle masses") {
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment