From 00753aa6cd74bfb70e8d008248de30fbb95a9a51 Mon Sep 17 00:00:00 2001
From: Maximilian Reininghaus <maximilian.reininghaus@tu-dortmund.de>
Date: Tue, 25 Jun 2024 15:17:00 +0200
Subject: [PATCH] proper handling&checking of units in cross-section tables
---
.../detail/modules/pythia8/Interaction.inl | 19 +++++++++++++------
.../framework/utility/CrossSectionTable.hpp | 3 +++
corsika/modules/pythia8/Interaction.hpp | 2 ++
3 files changed, 18 insertions(+), 6 deletions(-)
diff --git a/corsika/detail/modules/pythia8/Interaction.inl b/corsika/detail/modules/pythia8/Interaction.inl
index 263c2f9e2..8795cbc97 100644
--- a/corsika/detail/modules/pythia8/Interaction.inl
+++ b/corsika/detail/modules/pythia8/Interaction.inl
@@ -121,9 +121,13 @@ namespace corsika::pythia8 {
}
auto xs_it = boost::make_transform_iterator(total_xs.cbegin(), millibarn_mult);
- crossSectionTables_.insert({std::pair{projectile, target},
- CrossSectionTable<InterpolationTransforms::Log>{
- energies.cbegin(), energies.cend(), xs_it}});
+ auto e_it_beg = boost::make_transform_iterator(energies.cbegin(), GeV_mult);
+ decltype(e_it_beg) e_it_end =
+ boost::make_transform_iterator(energies.cend(), GeV_mult);
+ crossSectionTables_.insert(
+ {std::pair{projectile, target},
+ CrossSectionTable<InterpolationTransforms::Log>{
+ std::move(e_it_beg), std::move(e_it_end), std::move(xs_it)}});
}
}
}
@@ -144,9 +148,12 @@ namespace corsika::pythia8 {
"Pythia8 table corrupt (elab size = {}; xs size = {})", e_size, xs_size)};
}
- auto const xs_it = boost::make_transform_iterator(xs.cbegin(), millibarn_mult);
- return CrossSectionTable<InterpolationTransforms::Log>{energies.cbegin(),
- energies.cend(), xs_it};
+ auto xs_it = boost::make_transform_iterator(xs.cbegin(), millibarn_mult);
+ auto e_it_beg = boost::make_transform_iterator(energies.cbegin(), GeV_mult);
+ decltype(e_it_beg) e_it_end =
+ boost::make_transform_iterator(energies.cend(), GeV_mult);
+ return CrossSectionTable<InterpolationTransforms::Log>{
+ std::move(e_it_beg), std::move(e_it_end), std::move(xs_it)};
}
inline bool Interaction::isValid(Code const projectileId, Code const targetId,
diff --git a/corsika/framework/utility/CrossSectionTable.hpp b/corsika/framework/utility/CrossSectionTable.hpp
index 67515dd4a..36a28c217 100644
--- a/corsika/framework/utility/CrossSectionTable.hpp
+++ b/corsika/framework/utility/CrossSectionTable.hpp
@@ -40,6 +40,9 @@ namespace corsika {
template <typename InputIt1, typename InputIt2>
CrossSectionTable(InputIt1 energiesFirst, InputIt1 energiesLast,
InputIt2 crosssectionsFirst) {
+ static_assert(std::is_same_v<typename InputIt1::value_type, HEPEnergyType>);
+ static_assert(std::is_same_v<typename InputIt2::value_type, CrossSectionType>);
+
table_.reserve(std::distance(energiesFirst, energiesLast));
auto itE = boost::make_transform_iterator(std::move(energiesFirst), transform_);
diff --git a/corsika/modules/pythia8/Interaction.hpp b/corsika/modules/pythia8/Interaction.hpp
index 9f6b0616d..310a753b7 100644
--- a/corsika/modules/pythia8/Interaction.hpp
+++ b/corsika/modules/pythia8/Interaction.hpp
@@ -81,7 +81,9 @@ namespace corsika::pythia8 {
using key_type = std::pair<corsika::Code, corsika::Code>;
private:
+ static auto constexpr GeV_mult = [](float x) { return x * 1_GeV; };
static auto constexpr millibarn_mult = [](float x) { return x * 1_mb; };
+
static CrossSectionTable<InterpolationTransforms::Log> loadPPTable(
boost::filesystem::path const& dataPath, char const* key);
--
GitLab