Add optional per-sector seeding to the TPC digitizer

This adds an option that makes the TPC digits independent of the lane assignment, as a preparation for MT digitization.

- The random rings of ElectronTransport, GEMAmplification and SAMPAProcessing kept their positions from one sector to the next, so the digits depended on the number of TPC lanes.
- With --tpc-sector-seeding, each Digitizer owns these three objects, fills their rings once from DigiParams.seed, and moves the rings to positions derived from the seed, the sector and the timeframe.
- RandomRing can now be filled from a given generator and moved to a random position.
- Without the option, the output is unchanged.

https://its.cern.ch/jira/browse/O2-7309

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

Author: sawenzel

Common/MathUtils/include/MathUtils/RandomRing.h

@@ -54,19 +54,22 @@ class RandomRing
 
   /// constructor
   /// @param [in] randomType type of the random generator
-  RandomRing(const RandomType randomType = RandomType::Gaus);
+  /// @param [in] rng random generator to draw from, gRandom if null
+  RandomRing(const RandomType randomType = RandomType::Gaus, TRandom* rng = nullptr);
 
   /// constructor accepting TF1
   /// @param [in] function TF1 function
   RandomRing(TF1& function);
 
   /// initialisation of the random ring
   /// @param [in] randomType type of the random generator
-  void initialize(const RandomType randomType = RandomType::Gaus);
+  /// @param [in] rng random generator to draw from, gRandom if null
+  void initialize(const RandomType randomType = RandomType::Gaus, TRandom* rng = nullptr);
 
   /// initialisation of the random ring
-  /// @param [in] randomType type of the random generator
-  void initialize(TF1& function);
+  /// @param [in] function TF1 function
+  /// @param [in] rng random generator to draw from, gRandom if null
+  void initialize(TF1& function, TRandom* rng = nullptr);
 
   /// initialisation of the random ring
   /// @param [in] randomType type of the random generator
@@ -110,6 +113,10 @@ class RandomRing
   /// @return position in the ring buffer
   unsigned int getRingPosition() const { return mRingPosition; }
 
+  /// move to a random position in the ring buffer; not for rings read with getNextValueVc()
+  /// @param [in] rng random generator to draw the position from
+  void setRandomPosition(TRandom& rng) { mRingPosition = rng.Integer(N); }
+
  private:
   // =========================================================================
   // ===| members |===========================================================
@@ -123,12 +130,12 @@ class RandomRing
 
 //______________________________________________________________________________
 template <size_t N>
-inline RandomRing<N>::RandomRing(const RandomType randomType)
+inline RandomRing<N>::RandomRing(const RandomType randomType, TRandom* rng)
   : mRandomType(randomType),
     mRandomNumbers()
 
 {
-  initialize(randomType);
+  initialize(randomType, rng);
 }
 
 //______________________________________________________________________________
@@ -142,18 +149,18 @@ inline RandomRing<N>::RandomRing(TF1& function)
 
 //______________________________________________________________________________
 template <size_t N>
-inline void RandomRing<N>::initialize(const RandomType randomType)
+inline void RandomRing<N>::initialize(const RandomType randomType, TRandom* rng)
 {
-
+  TRandom* random = rng ? rng : gRandom;
   for (auto& v : mRandomNumbers) {
     // TODO: configurable mean and sigma
     switch (randomType) {
       case RandomType::Gaus: {
-        v = gRandom->Gaus(0, 1);
+        v = random->Gaus(0, 1);
         break;
       }
       case RandomType::Flat: {
-        v = gRandom->Rndm();
+        v = random->Rndm();
         break;
       }
       default: {
@@ -166,11 +173,11 @@ inline void RandomRing<N>::initialize(const RandomType randomType)
 
 //______________________________________________________________________________
 template <size_t N>
-inline void RandomRing<N>::initialize(TF1& function)
+inline void RandomRing<N>::initialize(TF1& function, TRandom* rng)
 {
   mRandomType = RandomType::CustomTF1;
   for (auto& v : mRandomNumbers) {
-    v = function.GetRandom();
+    v = function.GetRandom(rng);
   }
 }
 

Detectors/TPC/simulation/include/TPCSimulation/Digitizer.h

@@ -31,6 +31,8 @@ namespace tpc
 {
 
 class DigitContainer;
+class ElectronTransport;
+class GEMAmplification;
 
 template <class T>
 class SpaceCharge;
@@ -141,7 +143,25 @@ class Digitizer
   /// in case of scaled distortions, the distortions can be recalculated to ensure consistent distortions and corrections
   void recalculateDistortions();
 
+  /// Use digitizer-owned random rings filled from seed and positioned per sector and timeframe in init(),
+  /// so that the digits of a sector do not depend on which sectors were processed before
+  /// \param seed base seed, typically DigiParams.seed
+  void setUseSectorSeeding(ULong64_t seed);
+
+  /// Set the timeframe identifier that enters the sector seeding
+  void setSeedingTimeframe(ULong64_t timeframe) { mSeedingTimeframe = timeframe; }
+
  private:
+  ElectronTransport& electronTransport();
+  GEMAmplification& gemAmplification();
+  SAMPAProcessing& sampaProcessing();
+
+  std::unique_ptr<ElectronTransport> mElectronTransport; //!< Digitizer-owned instance with sector seeding, else the singleton is used
+  std::unique_ptr<GEMAmplification> mGEMAmplification;   //!< Digitizer-owned instance with sector seeding, else the singleton is used
+  std::unique_ptr<SAMPAProcessing> mSAMPAProcessing;     //!< Digitizer-owned instance with sector seeding, else the singleton is used
+  ULong64_t mSeed = 0;                                   //!< Base seed for sector seeding
+  ULong64_t mSeedingTimeframe = 0;                       //!< Timeframe identifier for sector seeding
+  bool mUseSectorSeeding = false;                        //!< Whether sector seeding is used
   DigitContainer mDigitContainer;      ///< Container for the Digits
   std::vector<float> mSignalArray;     //!< Workspace for the shaped signal of one electron
   std::unique_ptr<SC> mSpaceCharge;    ///< Handler of full distortions (static + IR dependant)

Detectors/TPC/simulation/include/TPCSimulation/ElectronTransport.h

@@ -41,9 +41,16 @@ class ElectronTransport
     return electronTransport;
   }
 
+  /// Constructor
+  /// \param rng random generator to fill the random rings from, gRandom if null
+  explicit ElectronTransport(TRandom* rng = nullptr);
+
   /// Destructor
   ~ElectronTransport() = default;
 
+  /// Move the random rings to new positions drawn from rng
+  void setRandomRingPositions(TRandom& rng);
+
   /// Update the OCDB parameters cached in the class. To be called once per event
   void updateParameters(float vdrift = 0);
 
@@ -72,8 +79,6 @@ class ElectronTransport
   float getDriftTime(float zPos, float signChange = 1.f) const;
 
  private:
-  ElectronTransport();
-
   /// Circular random buffer containing random values of the Gauss distribution to take into account diffusion of the
   /// electrons
   math_utils::RandomRing<> mRandomGaus;

Detectors/TPC/simulation/include/TPCSimulation/GEMAmplification.h

@@ -43,9 +43,16 @@ class GEMAmplification
     return gemAmplification;
   }
 
+  /// Constructor
+  /// \param rng random generator to fill the random rings from, gRandom if null
+  explicit GEMAmplification(TRandom* rng = nullptr);
+
   /// Destructor
   ~GEMAmplification() = default;
 
+  /// Move the random rings to new positions drawn from rng
+  void setRandomRingPositions(TRandom& rng);
+
   /// Update the OCDB parameters cached in the class. To be called once per event
   void updateParameters();
 
@@ -90,8 +97,6 @@ class GEMAmplification
   int getGEMMultiplication(int nElectrons, int GEM);
 
  private:
-  GEMAmplification();
-
   /// Circular random buffer containing random Gaus values for gain fluctuation if the number of electrons is larger
   /// (central limit theorem)
   math_utils::RandomRing<> mRandomGaus;

Detectors/TPC/simulation/include/TPCSimulation/SAMPAProcessing.h

@@ -45,9 +45,16 @@ class SAMPAProcessing
     static SAMPAProcessing sampaProcessing;
     return sampaProcessing;
   }
+  /// Constructor
+  /// \param rng random generator to fill the random ring from, gRandom if null
+  explicit SAMPAProcessing(TRandom* rng = nullptr);
+
   /// Destructor
   ~SAMPAProcessing() = default;
 
+  /// Move the random ring to a new position drawn from rng
+  void setRandomRingPositions(TRandom& rng) { mRandomNoiseRing.setRandomPosition(rng); }
+
   /// Update the OCDB parameters cached in the class. To be called once per event
   void updateParameters(float vdrift = 0);
 
@@ -134,7 +141,6 @@ class SAMPAProcessing
   float getPedestalCRU(const int sector, const int globalPadInSector) const;
 
  private:
-  SAMPAProcessing();
   const ParameterGas* mGasParam;             ///< Caching of the parameter class to avoid multiple CDB calls
   const ParameterDetector* mDetParam;        ///< Caching of the parameter class to avoid multiple CDB calls
   const ParameterElectronics* mEleParam;     ///< Caching of the parameter class to avoid multiple CDB calls

Detectors/TPC/simulation/src/Digitizer.cxx

@@ -30,6 +30,7 @@
 #include "TPCCalibration/CorrMapParam.h"
 
 #include <fairlogger/Logger.h>
+#include <TRandom3.h>
 
 ClassImp(o2::tpc::Digitizer);
 
@@ -39,13 +40,68 @@ Digitizer::~Digitizer() = default;
 
 Digitizer::Digitizer() = default;
 
+namespace
+{
+// splitmix64 finaliser, to derive well separated seeds from the base seed
+ULong64_t mixSeed(ULong64_t x)
+{
+  x += 0x9e3779b97f4a7c15ULL;
+  x = (x ^ (x >> 30)) * 0xbf58476d1ce4e5b9ULL;
+  x = (x ^ (x >> 27)) * 0x94d049bb133111ebULL;
+  return x ^ (x >> 31);
+}
+
+// TRandom3 takes 32 bits of the seed and treats 0 as a request for a random seed
+UInt_t toTRandom3Seed(ULong64_t x)
+{
+  const UInt_t seed = static_cast<UInt_t>(x ^ (x >> 32));
+  return seed ? seed : 1;
+}
+} // namespace
+
+ElectronTransport& Digitizer::electronTransport()
+{
+  return mElectronTransport ? *mElectronTransport : ElectronTransport::instance();
+}
+
+GEMAmplification& Digitizer::gemAmplification()
+{
+  return mGEMAmplification ? *mGEMAmplification : GEMAmplification::instance();
+}
+
+SAMPAProcessing& Digitizer::sampaProcessing()
+{
+  return mSAMPAProcessing ? *mSAMPAProcessing : SAMPAProcessing::instance();
+}
+
+void Digitizer::setUseSectorSeeding(ULong64_t seed)
+{
+  mUseSectorSeeding = true;
+  mSeed = seed;
+}
+
 void Digitizer::init()
 {
-  auto& gemAmplification = GEMAmplification::instance();
+  if (mUseSectorSeeding) {
+    if (!mGEMAmplification) {
+      // the rings are filled once, in the construction order of the singletons
+      TRandom3 rng(toTRandom3Seed(mixSeed(mSeed)));
+      mGEMAmplification = std::make_unique<GEMAmplification>(&rng);
+      mElectronTransport = std::make_unique<ElectronTransport>(&rng);
+      mSAMPAProcessing = std::make_unique<SAMPAProcessing>(&rng);
+      mDigitContainer.setSAMPAProcessing(mSAMPAProcessing.get());
+    }
+    const ULong64_t sectorSeed = mixSeed(mSeed ^ mixSeed(static_cast<ULong64_t>(mSector) + 1));
+    TRandom3 rng(toTRandom3Seed(mixSeed(sectorSeed ^ mixSeed(mSeedingTimeframe))));
+    mGEMAmplification->setRandomRingPositions(rng);
+    mElectronTransport->setRandomRingPositions(rng);
+    mSAMPAProcessing->setRandomRingPositions(rng);
+  }
+  auto& gemAmplification = this->gemAmplification();
   gemAmplification.updateParameters();
-  auto& electronTransport = ElectronTransport::instance();
+  auto& electronTransport = this->electronTransport();
   electronTransport.updateParameters(mVDrift);
-  auto& sampaProcessing = SAMPAProcessing::instance();
+  auto& sampaProcessing = this->sampaProcessing();
   // the calibration update may reload FEE parameters from file, so SAMPAProcessing picks up its pointers afterwards
   mDigitContainer.updateCalibrations();
   sampaProcessing.updateParameters(mVDrift);
@@ -59,9 +115,9 @@ void Digitizer::process(const std::vector<o2::tpc::HitGroup>& hits,
   auto& eleParam = ParameterElectronics::Instance();
   auto& gemParam = ParameterGEM::Instance();
 
-  auto& gemAmplification = GEMAmplification::instance();
-  auto& electronTransport = ElectronTransport::instance();
-  auto& sampaProcessing = SAMPAProcessing::instance();
+  auto& gemAmplification = this->gemAmplification();
+  auto& electronTransport = this->electronTransport();
+  auto& sampaProcessing = this->sampaProcessing();
 
   const int nShapedPoints = eleParam.NShapedPoints;
   const auto amplificationMode = gemParam.AmplMode;
@@ -171,7 +227,7 @@ void Digitizer::flush(std::vector<o2::tpc::Digit>& digits,
                       std::vector<o2::tpc::CommonMode>& commonModeOutput,
                       bool finalFlush)
 {
-  SAMPAProcessing& sampaProcessing = SAMPAProcessing::instance();
+  SAMPAProcessing& sampaProcessing = this->sampaProcessing();
   mDigitContainer.fillOutputContainer(digits, labels, commonModeOutput, mSector, sampaProcessing.getTimeBinFromTime(mEventTime - mOutputDigitTimeOffset), mIsContinuous, finalFlush);
   // flushing debug output to file
   if (((finalFlush && mIsContinuous) || (!mIsContinuous)) && mSpaceCharge) {
@@ -227,7 +283,7 @@ void Digitizer::setStartTime(double time)
   // this is setting the first timebin index for the digit container
   // note that negative times w.r.t start of timeframe/data-taking == mOutputDigitTimeOffset
   // will yield the 0-th bin (due to casting logic in sampaProcessing)
-  SAMPAProcessing& sampaProcessing = SAMPAProcessing::instance();
+  SAMPAProcessing& sampaProcessing = this->sampaProcessing();
   sampaProcessing.updateParameters(mVDrift);
   const auto timediff = time - mOutputDigitTimeOffset;
   const auto starttimebin = sampaProcessing.getTimeBinFromTime(timediff);

Detectors/TPC/simulation/src/ElectronTransport.cxx

@@ -21,11 +21,17 @@
 using namespace o2::tpc;
 using namespace o2::math_utils;
 
-ElectronTransport::ElectronTransport() : mRandomGaus(), mRandomFlat(RandomRing<>::RandomType::Flat)
+ElectronTransport::ElectronTransport(TRandom* rng) : mRandomGaus(RandomRing<>::RandomType::Gaus, rng), mRandomFlat(RandomRing<>::RandomType::Flat, rng)
 {
   updateParameters();
 }
 
+void ElectronTransport::setRandomRingPositions(TRandom& rng)
+{
+  mRandomGaus.setRandomPosition(rng);
+  mRandomFlat.setRandomPosition(rng);
+}
+
 void ElectronTransport::updateParameters(float vdrift)
 {
   mGasParam = &(ParameterGas::Instance());

Detectors/TPC/simulation/src/GEMAmplification.cxx

@@ -26,10 +26,15 @@ using namespace o2::tpc;
 using namespace o2::math_utils;
 using boost::format;
 
-GEMAmplification::GEMAmplification()
-  : mRandomGaus(),
-    mRandomFlat(RandomRing<>::RandomType::Flat),
-    mGain()
+GEMAmplification::GEMAmplification(TRandom* rng)
+  : mRandomGaus(RandomRing<>::RandomType::Gaus, rng),
+    mRandomFlat(RandomRing<>::RandomType::Flat, rng),
+    // the gain rings are refilled from the Polya distributions below; with an explicit generator they start empty
+    mGain{RandomRing<>(rng ? RandomRing<>::RandomType::CustomTF1 : RandomRing<>::RandomType::Gaus),
+          RandomRing<>(rng ? RandomRing<>::RandomType::CustomTF1 : RandomRing<>::RandomType::Gaus),
+          RandomRing<>(rng ? RandomRing<>::RandomType::CustomTF1 : RandomRing<>::RandomType::Gaus),
+          RandomRing<>(rng ? RandomRing<>::RandomType::CustomTF1 : RandomRing<>::RandomType::Gaus)},
+    mGainFullStack(rng ? RandomRing<>::RandomType::CustomTF1 : RandomRing<>::RandomType::Gaus)
 {
   updateParameters();
 
@@ -65,7 +70,7 @@ GEMAmplification::GEMAmplification()
       polyaDistribution = (o2::math_utils::CachingTF1*)outfile->Get(TString::Format("func%d", i).Data());
       // FIXME: verify that distribution corresponds to the parameters used here
     }
-    mGain[i].initialize(*polyaDistribution);
+    mGain[i].initialize(*polyaDistribution, rng);
 
     if (!cacheexists) {
       outfile->WriteTObject(polyaDistribution, TString::Format("func%d", i).Data());
@@ -88,7 +93,7 @@ GEMAmplification::GEMAmplification()
   } else {
     polyaDistribution = (o2::math_utils::CachingTF1*)outfile->Get("polyaStack");
   }
-  mGainFullStack.initialize(*polyaDistribution);
+  mGainFullStack.initialize(*polyaDistribution, rng);
 
   if (!cacheexists) {
     outfile->WriteTObject(polyaDistribution, "polyaStack");
@@ -102,6 +107,16 @@ GEMAmplification::GEMAmplification()
   LOG(info) << "TPC: GEM setup (polya) took " << watch.CpuTime();
 }
 
+void GEMAmplification::setRandomRingPositions(TRandom& rng)
+{
+  mRandomGaus.setRandomPosition(rng);
+  mRandomFlat.setRandomPosition(rng);
+  for (auto& gain : mGain) {
+    gain.setRandomPosition(rng);
+  }
+  mGainFullStack.setRandomPosition(rng);
+}
+
 void GEMAmplification::updateParameters()
 {
   auto& cdb = CDBInterface::instance();

Detectors/TPC/simulation/src/SAMPAProcessing.cxx

@@ -24,7 +24,7 @@
 
 using namespace o2::tpc;
 
-SAMPAProcessing::SAMPAProcessing() : mRandomNoiseRing()
+SAMPAProcessing::SAMPAProcessing(TRandom* rng) : mRandomNoiseRing(o2::math_utils::RandomRing<>::RandomType::Gaus, rng)
 {
   updateParameters();
 }