feat: Random number channel for Examples

Examples/Algorithms/Digitization/src/DigitizationAlgorithm.cpp

@@ -14,24 +14,20 @@
 #include "Acts/Geometry/TrackingGeometry.hpp"
 #include "Acts/Surfaces/Surface.hpp"
 #include "Acts/Utilities/BinUtility.hpp"
-#include "Acts/Utilities/Result.hpp"
 #include "ActsExamples/Digitization/ModuleClusters.hpp"
 #include "ActsExamples/EventData/GeometryContainers.hpp"
 #include "ActsExamples/EventData/Index.hpp"
-#include "ActsExamples/EventData/SimHit.hpp"
+#include "ActsExamples/EventData/SimParticle.hpp"
 #include "ActsExamples/Framework/AlgorithmContext.hpp"
-#include "ActsExamples/Utilities/GroupBy.hpp"
+#include "ActsExamples/Framework/RandomNumbers.hpp"
 #include "ActsExamples/Utilities/Range.hpp"
 #include "ActsFatras/EventData/Barcode.hpp"
 #include "ActsFatras/EventData/Hit.hpp"
 
 #include <algorithm>
 #include <array>
-#include <cmath>
-#include <cstdint>
 #include <limits>
 #include <ostream>
-#include <set>
 #include <stdexcept>
 #include <string>
 #include <utility>
@@ -155,8 +151,7 @@ ActsExamples::ProcessCode ActsExamples::DigitizationAlgorithm::execute(
   measurementParticlesMap.reserve(simHits.size());
   measurementSimHitsMap.reserve(simHits.size());
 
-  // Setup random number generator
-  auto rng = m_cfg.randomNumbers->spawnGenerator(ctx);
+  RandomSeed eventSeed = m_cfg.randomNumbers->generateSeed(ctx);
 
   // Some statistics
   std::size_t skippedHits = 0;
@@ -165,6 +160,10 @@ ActsExamples::ProcessCode ActsExamples::DigitizationAlgorithm::execute(
   // Thus we need to store the cell data from the simulation.
   CellsMap cellsMap;
 
+  // Count the number of hits per particle and module. Cleared before each
+  // module.
+  std::map<SimBarcode, std::uint32_t> particleOnModuleHitCount;
+
   ACTS_DEBUG("Starting loop over modules ...");
   for (const auto& simHitsGroup : groupByModule(simHits)) {
     // Manual pair unpacking instead of using
@@ -194,6 +193,8 @@ ActsExamples::ProcessCode ActsExamples::DigitizationAlgorithm::execute(
       ACTS_VERBOSE("Digitizer found for module " << moduleGeoId);
     }
 
+    particleOnModuleHitCount.clear();
+
     // Run the digitizer. Iterate over the hits for this surface inside the
     // visitor so we do not need to lookup the variant object per-hit.
     std::visit(
@@ -206,10 +207,23 @@ ActsExamples::ProcessCode ActsExamples::DigitizationAlgorithm::execute(
             const auto& simHit = *h;
             const auto simHitIdx = simHits.index_of(h);
 
+            const auto hitIndex =
+                particleOnModuleHitCount[simHit.particleId()]++;
+
+            if (hitIndex > 0) {
+              ACTS_DEBUG("more than one hit for particle "
+                         << simHit.particleId() << " on module " << moduleGeoId
+                         << " - still forwarding it to the digitizer");
+            }
+
+            RandomSeed hitSeed = eventSeed + moduleGeoId.value() +
+                                 simHit.particleId().value() + hitIndex;
+            RandomEngine rng(hitSeed);
+
             DigitizedParameters dParameters;
 
             if (simHit.depositedEnergy() < m_cfg.minEnergyDeposit) {
-              ACTS_VERBOSE("Skip hit because energy deposit to small");
+              ACTS_VERBOSE("Skip hit because energy deposit too small");
               continue;
             }
 

Examples/Framework/include/ActsExamples/Framework/RandomNumbers.hpp

@@ -6,28 +6,38 @@
 // License, v. 2.0. If a copy of the MPL was not distributed with this
 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
 
-//
-//  RandomNumbers.hpp
-//  ActsExamples
-//
-//  Created by Andreas Salzburger on 17/05/16.
-//
-//
-
 #pragma once
 
-#include "ActsExamples/Framework/AlgorithmContext.hpp"
-
 #include <cstdint>
 #include <random>
 
 namespace ActsExamples {
 struct AlgorithmContext;
 
 /// The random number generator used in the framework.
-using RandomEngine = std::mt19937;  ///< Mersenne Twister
+using RandomEngine = std::mt19937_64;  ///< Mersenne Twister
+
+/// The seed type used in the framework.
+using RandomSeed = std::uint64_t;
+
+class RandomNumbers;
+
+class RandomNumberChannel {
+ public:
+  RandomNumberChannel(const RandomNumbers& randomNumbers, RandomSeed seed);
+
+  RandomSeed seed() const;
+
+  RandomEngine createEngine() const;
 
-/// Provide event and algorithm specific random number generator.s
+  RandomNumberChannel createSubChannel(RandomSeed seed) const;
+
+ private:
+  const RandomNumbers* m_randomNumbers = nullptr;
+  RandomSeed m_seed = 0;
+};
+
+/// Provide an event and algorithm specific random number generator.
 ///
 /// This provides local random number generators, allowing for
 /// thread-safe, lock-free, and reproducible random number generation across
@@ -41,11 +51,16 @@ using RandomEngine = std::mt19937;  ///< Mersenne Twister
 class RandomNumbers {
  public:
   struct Config {
-    std::uint64_t seed = 1234567890u;  ///< random seed
+    RandomSeed seed = 1234567890u;  ///< random seed
   };
 
   explicit RandomNumbers(const Config& cfg);
 
+  RandomNumberChannel createChannel() const;
+
+  RandomNumberChannel createAlgorithmEventChannel(
+      const AlgorithmContext& context) const;
+
   /// Spawn an algorithm-local random number generator. To avoid inefficiencies
   /// and multiple uses of a given RNG seed, this should only be done once per
   /// Algorithm invocation, after what the generator object should be reused.
@@ -55,11 +70,11 @@ class RandomNumbers {
   /// @param context is the AlgorithmContext of the host algorithm
   RandomEngine spawnGenerator(const AlgorithmContext& context) const;
 
-  /// Generate a event and algorithm specific seed value.
+  /// Generate an event and algorithm specific seed value.
   ///
   /// This should only be used in special cases e.g. where a custom
   /// random engine is used and `spawnGenerator` can not be used.
-  std::uint64_t generateSeed(const AlgorithmContext& context) const;
+  RandomSeed generateSeed(const AlgorithmContext& context) const;
 
  private:
   Config m_cfg;

Examples/Framework/src/Framework/RandomNumbers.cpp

@@ -6,26 +6,49 @@
 // License, v. 2.0. If a copy of the MPL was not distributed with this
 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
 
-//
-//  RandomNumbers.cpp
-//  ActsExamples
-//
-//  Created by Andreas Salzburger on 17/05/16.
-//
-//
-
 #include "ActsExamples/Framework/RandomNumbers.hpp"
 
 #include "ActsExamples/Framework/AlgorithmContext.hpp"
 
-ActsExamples::RandomNumbers::RandomNumbers(const Config& cfg) : m_cfg(cfg) {}
+#include <boost/functional/hash.hpp>
+
+namespace ActsExamples {
+
+RandomNumberChannel::RandomNumberChannel(const RandomNumbers& randomNumbers,
+                                         RandomSeed seed)
+    : m_randomNumbers(&randomNumbers), m_seed(seed) {}
+
+RandomSeed RandomNumberChannel::seed() const {
+  return m_seed;
+}
+
+RandomEngine RandomNumberChannel::createEngine() const {
+  return RandomEngine(m_seed);
+}
 
-ActsExamples::RandomEngine ActsExamples::RandomNumbers::spawnGenerator(
+RandomNumberChannel RandomNumberChannel::createSubChannel(
+    RandomSeed seed) const {
+  return RandomNumberChannel(*m_randomNumbers, m_seed + seed);
+}
+
+RandomNumbers::RandomNumbers(const Config& cfg) : m_cfg(cfg) {}
+
+RandomNumberChannel RandomNumbers::createChannel() const {
+  return RandomNumberChannel(*this, m_cfg.seed);
+}
+
+RandomNumberChannel RandomNumbers::createAlgorithmEventChannel(
     const AlgorithmContext& context) const {
-  return RandomEngine(generateSeed(context));
+  return RandomNumberChannel(*this, generateSeed(context));
 }
 
-std::uint64_t ActsExamples::RandomNumbers::generateSeed(
+RandomEngine RandomNumbers::spawnGenerator(
     const AlgorithmContext& context) const {
+  return RandomEngine(generateSeed(context));
+}
+
+RandomSeed RandomNumbers::generateSeed(const AlgorithmContext& context) const {
   return m_cfg.seed + context.eventNumber;
 }
+
+}  // namespace ActsExamples

Tests/UnitTests/Examples/Io/Json/JsonDigitizationConfigTests.cpp

@@ -141,7 +141,7 @@ BOOST_AUTO_TEST_CASE(GaussianSmearing) {
     if (s.forcePositive[i]) {
       ref = std::abs(ref);
     }
-    CHECK_CLOSE_REL(par[i], ref, 0.15);
+    CHECK_CLOSE_ABS(par[i], ref, 0.15);
   }
 }