IAP GITLAB

Skip to content
Snippets Groups Projects
Commit 4c838dec authored by Maximilian Reininghaus's avatar Maximilian Reininghaus :vulcan:
Browse files

use WeightedSum in GetInteractionLength

parent 8a1c6daa
No related branches found
No related tags found
No related merge requests found
......@@ -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);
}
}
......
......@@ -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: "
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment