diff --git a/Environment/NuclearComposition.h b/Environment/NuclearComposition.h index e7f0d6d8a3827985238117dc24db28345d25463e..a3ebfe8f819d5c51ba91a01df9a266db1167d40c 100644 --- a/Environment/NuclearComposition.h +++ b/Environment/NuclearComposition.h @@ -16,6 +16,7 @@ #include <corsika/units/PhysicalUnits.h> #include <cassert> +#include <functional> #include <numeric> #include <random> #include <stdexcept> @@ -82,7 +83,6 @@ namespace corsika::environment { auto WeightedSum(TFunction func) const { using ResultQuantity = decltype(func(*fComponents.cbegin())); - auto const sum = [](auto x, auto y) { return x + y; }; auto const prod = [&](auto const compID, auto const fraction) { return func(compID) * fraction; }; @@ -91,12 +91,12 @@ namespace corsika::environment { return std::inner_product( fComponents.cbegin(), fComponents.cend(), fNumberFractions.cbegin(), ResultQuantity::zero(), // .zero() is defined for quantity types only - sum, prod); + std::plus<ResultQuantity>(), prod); } else { return std::inner_product( fComponents.cbegin(), fComponents.cend(), fNumberFractions.cbegin(), ResultQuantity(0), // in other cases we have to use a bare 0 - sum, prod); + std::plus<ResultQuantity>(), prod); } } diff --git a/Processes/Pythia/Interaction.cc b/Processes/Pythia/Interaction.cc index ae8104350086f5cfda94b174305d0f0366faa2ce..f4272b165c842f80487c112ea8e106248dc21c0b 100644 --- a/Processes/Pythia/Interaction.cc +++ b/Processes/Pythia/Interaction.cc @@ -215,26 +215,12 @@ namespace corsika::process::pythia { const auto mediumComposition = currentNode->GetModelProperties().GetNuclearComposition(); // determine average interaction length - // weighted sum - int i = -1; - si::CrossSectionType weightedProdCrossSection = 0_mb; - // get weights of components from environment/medium - const auto& w = mediumComposition.GetFractions(); - // loop over components in medium - for (auto const targetId : mediumComposition.GetComponents()) { - i++; - cout << "Interaction: get interaction length for target: " << targetId << endl; - - auto const [productionCrossSection, elaCrossSection] = - GetCrossSection(corsikaBeamId, targetId, ECoM); - [[maybe_unused]] const auto& dummy_elaCrossSection = elaCrossSection; - - cout << "Interaction: IntLength: pythia return (mb): " - << productionCrossSection / 1_mb << endl - << "Interaction: IntLength: weight : " << w[i] << endl; - - weightedProdCrossSection += w[i] * productionCrossSection; - } + + auto const weightedProdCrossSection = + mediumComposition.WeightedSum([=](auto vTargetID) { + return std::get<0>(this->GetCrossSection(corsikaBeamId, vTargetID, ECoM)); + }); + cout << "Interaction: IntLength: weighted CrossSection (mb): " << weightedProdCrossSection / 1_mb << endl << "Interaction: IntLength: average mass number: "