moved SampleTarget to NuclearComposition, removed temporary weight vector

Environment/IMediumModel.h

@@ -17,7 +17,6 @@
 #include <corsika/geometry/Point.h>
 #include <corsika/geometry/Trajectory.h>
 #include <corsika/units/PhysicalUnits.h>
-#include <random>
 
 namespace corsika::environment {
 
@@ -39,28 +38,6 @@ namespace corsika::environment {
         corsika::units::si::GrammageType) const = 0;
 
     virtual NuclearComposition const& GetNuclearComposition() const = 0;
-
-    template <class TRNG>
-    corsika::particles::Code SampleTarget(
-        std::vector<corsika::units::si::CrossSectionType> const& sigma,
-        TRNG& randomStream) const {
-      using namespace corsika::units::si;
-
-      auto const& nuclComp = GetNuclearComposition();
-      auto const& fractions = nuclComp.GetFractions();
-      assert(sigma.size() == fractions.size());
-
-      std::vector<float> weights(fractions.size());
-
-      for (size_t i = 0; i < fractions.size(); ++i) {
-        std::cout << "HomogeneousMedium: fraction: " << fractions[i] << std::endl;
-        weights[i] = fractions[i] * sigma[i].magnitude();
-      }
-
-      std::discrete_distribution channelDist(weights.cbegin(), weights.cend());
-      const int iChannel = channelDist(randomStream);
-      return nuclComp.GetComponents()[iChannel];
-    }
   };
 
 } // namespace corsika::environment

Environment/NuclearComposition.h

@@ -15,6 +15,7 @@
 #include <corsika/particles/ParticleProperties.h>
 #include <cassert>
 #include <numeric>
+#include <random>
 #include <stdexcept>
 #include <vector>
 
@@ -26,6 +27,35 @@ namespace corsika::environment {
 
     double const fAvgMassNumber;
 
+    template <class AConstIterator, class BConstIterator>
+    class WeightProviderIterator {
+      AConstIterator fAIter;
+      BConstIterator fBIter;
+
+    public:
+      using value_type = double;
+      using iterator_category = std::input_iterator_tag;
+      using pointer = double*;
+      using reference = double&;
+      using difference_type = ptrdiff_t;
+
+      WeightProviderIterator(AConstIterator a, BConstIterator b)
+          : fAIter(a)
+          , fBIter(b) {}
+
+      double operator*() const { return ((*fAIter) * (*fBIter)).magnitude(); }
+
+      WeightProviderIterator& operator++() { // prefix ++
+        ++fAIter;
+        ++fBIter;
+        return *this;
+      }
+
+      auto operator==(WeightProviderIterator other) { return fAIter == other.fAIter; }
+
+      auto operator!=(WeightProviderIterator other) { return !(*this == other); }
+    };
+
   public:
     NuclearComposition(std::vector<corsika::particles::Code> pComponents,
                        std::vector<float> pFractions)
@@ -51,6 +81,26 @@ namespace corsika::environment {
     auto const& GetFractions() const { return fNumberFractions; }
     auto const& GetComponents() const { return fComponents; }
     auto const GetAverageMassNumber() const { return fAvgMassNumber; }
+
+    template <class TRNG>
+    corsika::particles::Code SampleTarget(
+        std::vector<corsika::units::si::CrossSectionType> const& sigma,
+        TRNG& randomStream) const {
+      using namespace corsika::units::si;
+
+      assert(sigma.size() == fNumberFractions.size());
+
+      std::discrete_distribution channelDist(
+          WeightProviderIterator<decltype(fNumberFractions.begin()),
+                                 decltype(sigma.begin())>(fNumberFractions.begin(),
+                                                          sigma.begin()),
+          WeightProviderIterator<decltype(fNumberFractions.begin()),
+                                 decltype(sigma.end())>(fNumberFractions.end(),
+                                                        sigma.end()));
+
+      auto const iChannel = channelDist(randomStream);
+      return fComponents[iChannel];
+    }
   };
 
 } // namespace corsika::environment

Environment/testEnvironment.cc

@@ -17,6 +17,8 @@
 #include <corsika/environment/VolumeTreeNode.h>
 #include <corsika/particles/ParticleProperties.h>
 #include <catch2/catch.hpp>
+#include <random>
+#include <vector>
 
 using namespace corsika::geometry;
 using namespace corsika::environment;
@@ -28,3 +30,33 @@ TEST_CASE("HomogeneousMedium") {
       std::vector<float>{1.f});
   HomogeneousMedium<IMediumModel> const medium(19.2_g / cube(1_cm), protonComposition);
 }
+
+TEST_CASE("NuclearComposition") {
+  NuclearComposition const composition(
+      std::vector<corsika::particles::Code>{corsika::particles::Code::Proton,
+                                            corsika::particles::Code::Neutron},
+      std::vector<float>{2.f / 3.f, 1.f / 3.f});
+  SECTION("SampleTarget") {
+    std::vector<CrossSectionType> crossSections{50_mbarn, 100_mbarn};
+
+    std::mt19937 rng;
+
+    int proton{0}, neutron{0};
+
+    for (int i = 0; i < 1'000'000; ++i) {
+      corsika::particles::Code p = composition.SampleTarget(crossSections, rng);
+      switch (p) {
+        case corsika::particles::Code::Proton:
+          proton++;
+          break;
+        case corsika::particles::Code::Neutron:
+          neutron++;
+          break;
+        default:
+          throw std::runtime_error("");
+      }
+    }
+
+    REQUIRE(static_cast<double>(proton) / neutron == Approx(1).epsilon(1e-2));
+  }
+}

Processes/HadronicElasticModel/HadronicElasticModel.h

@@ -156,8 +156,8 @@ namespace corsika::process::HadronicElasticModel {
         cross_section_of_components[i] = CrossSection(s);
       }
 
-      const auto targetCode = currentNode->GetModelProperties().SampleTarget(
-          cross_section_of_components, fRNG);
+      const auto targetCode =
+          mediumComposition.SampleTarget(cross_section_of_components, fRNG);
 
       auto const targetMass = corsika::particles::GetMass(targetCode);
 

Processes/Sibyll/Interaction.h

@@ -264,8 +264,8 @@ namespace corsika::process::sibyll {
           cross_section_of_components[i] = sigProd;
         }
 
-        const auto targetCode = currentNode->GetModelProperties().SampleTarget(
-            cross_section_of_components, fRNG);
+        const auto targetCode =
+            mediumComposition.SampleTarget(cross_section_of_components, fRNG);
         cout << "Interaction: target selected: " << targetCode << endl;
         /*
           FOR NOW: allow nuclei with A<18 or protons only.