diff --git a/Environment/NuclearComposition.h b/Environment/NuclearComposition.h
index dee18a4463c340698bbff372cc504d4744748375..e7f0d6d8a3827985238117dc24db28345d25463e 100644
--- a/Environment/NuclearComposition.h
+++ b/Environment/NuclearComposition.h
@@ -13,6 +13,8 @@
 #define _include_NuclearComposition_h
 
 #include <corsika/particles/ParticleProperties.h>
+#include <corsika/units/PhysicalUnits.h>
+
 #include <cassert>
 #include <numeric>
 #include <random>
@@ -76,6 +78,28 @@ namespace corsika::environment {
       }
     }
 
+    template <typename TFunction>
+    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;
+      };
+
+      if constexpr (phys::units::is_quantity_v<ResultQuantity>) {
+        return std::inner_product(
+            fComponents.cbegin(), fComponents.cend(), fNumberFractions.cbegin(),
+            ResultQuantity::zero(), // .zero() is defined for quantity types only
+            sum, 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);
+      }
+    }
+
     auto size() const { return fNumberFractions.size(); }
 
     auto const& GetFractions() const { return fNumberFractions; }