Add optional per-sector seeding to the TPC digitizer (#16085)

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.
  
  - Two shared static workspaces of the digitizer are now members, and the SAMPAProcessing used for the output is
  passed down from DigitContainer.
  - With --tpc-sector-seeding, each Digitizer owns its ElectronTransport, GEMAmplification and SAMPAProcessing, 
  fills their random rings once from DigiParams.seed, and moves the rings to positions derived from the seed, the
  sector and the timeframe.
  - RandomRing can be filled from a given generator, moved to a random position, and wraps vector reads at its 
  end.
  - DigitContainer uses one call for all digitization modes.
  - Without the option, the output is unchanged.
  
  https://its.cern.ch/jira/browse/O2-7309

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
@@ -98,6 +101,10 @@ class RandomRing
     // within this header file (to reduce memory problems during compilation).
     // The hope is that the calling user calls this with a
     // correct Vc type (Vc::float_v) in a source file.
+    // a vector must not run past the end of the ring, also after setRandomPosition()
+    if (mRingPosition + VcType::size() > mRandomNumbers.size()) {
+      mRingPosition = 0;
+    }
     const VcType value = VcType(&mRandomNumbers[mRingPosition]);
     mRingPosition += VcType::size();
     if (mRingPosition >= mRandomNumbers.size()) {
@@ -110,6 +117,10 @@ class RandomRing
   /// @return position in the ring buffer
   unsigned int getRingPosition() const { return mRingPosition; }
 
+  /// move to a random position in the ring buffer
+  /// @param [in] rng random generator to draw the position from
+  void setRandomPosition(TRandom& rng) { mRingPosition = rng.Integer(N); }
+
  private:
   // =========================================================================
   // ===| members |===========================================================
@@ -123,12 +134,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 +153,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 +177,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/DigitContainer.h

@@ -21,6 +21,7 @@
 #include "TPCBase/CRU.h"
 #include "DataFormatsTPC/Defs.h"
 #include "TPCSimulation/DigitTime.h"
+#include "TPCSimulation/SAMPAProcessing.h"
 #include "TPCBase/CalDet.h"
 #include "CommonUtils/DebugStreamer.h"
 #include "TPCBase/ParameterDetector.h"
@@ -87,6 +88,9 @@ class DigitContainer
   /// Get the size of the container for one event
   size_t size() const { return mTimeBins.size(); }
 
+  /// Use this SAMPAProcessing instead of the singleton when filling the output
+  void setSAMPAProcessing(SAMPAProcessing* sampaProcessing) { mSAMPAProcessing = sampaProcessing; }
+
   /// Fetch the calibration objects used by fillOutputContainer from the CDBInterface
   /// Must be called after the CDBInterface is configured; fillOutputContainer calls it lazily otherwise
   void updateCalibrations();
@@ -103,6 +107,8 @@ class DigitContainer
   const CalDet<bool>* mDeadMap = nullptr;                     //!< Dead channel map
   bool mCMCEnabled = false;                                   //!< Common mode correction enabled in the FEE config (DigiMode Auto)
   bool mCalibrationsUpdated = false;                          //!< Whether updateCalibrations was called
+  SAMPAProcessing* mSAMPAProcessing = nullptr;                //!< SAMPAProcessing used for the output, the singleton if not set
+  std::vector<std::pair<MCCompLabel, int>> mLabelCollector;   //!< Workspace to sort the MC labels of a digit
 
   void reportSettings();
 };

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

@@ -97,6 +97,7 @@ class DigitGlobalPad
                            const CRU& cru, TimeBin timeBin,
                            GlobalPadNumber globalPad,
                            o2::dataformats::LabelContainer<std::pair<MCCompLabel, int>, false>& labelContainer,
+                           SAMPAProcessing& sampaProcessing, std::vector<std::pair<MCCompLabel, int>>& labelCollector,
                            float commonMode = 0.f, const PrevDigitInfo& prevDigit = PrevDigitInfo(), Streamer* debugStream = nullptr, const CalDet<bool>* deadMap = nullptr);
 
  private:
@@ -210,12 +211,11 @@ inline void DigitGlobalPad::fillOutputContainer(std::vector<Digit>& output,
                                                 const CRU& cru, TimeBin timeBin,
                                                 GlobalPadNumber globalPad,
                                                 o2::dataformats::LabelContainer<std::pair<MCCompLabel, int>, false>& labels,
+                                                SAMPAProcessing& sampaProcessing, std::vector<std::pair<MCCompLabel, int>>& labelCollector,
                                                 float commonMode, const PrevDigitInfo& prevDigit, Streamer* debugStream, const CalDet<bool>* deadMap)
 {
   const Mapper& mapper = Mapper::instance();
-  SAMPAProcessing& sampaProcessing = SAMPAProcessing::instance();
   const PadPos pad = mapper.padPos(globalPad);
-  static std::vector<std::pair<MCCompLabel, int>> labelCollector; // static workspace container for sorting
 
   /// The charge accumulated on that pad is converted into ADC counts, saturation of the SAMPA is applied and a Digit
   /// is created in written out

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

@@ -78,6 +78,7 @@ class DigitTime
   template <DigitzationMode MODE>
   void fillOutputContainer(std::vector<Digit>& output, dataformats::MCTruthContainer<MCCompLabel>& mcTruth,
                            std::vector<CommonMode>& commonModeOutput, const Sector& sector, TimeBin timeBin,
+                           SAMPAProcessing& sampaProcessing, std::vector<std::pair<MCCompLabel, int>>& labelCollector,
                            PrevDigitInfoArray* prevTime = nullptr, Streamer* debugStream = nullptr,
                            const CalPad* padParams[3] = nullptr, const CalDet<bool>* deadMap = nullptr);
 
@@ -129,6 +130,7 @@ inline float DigitTime::getCommonMode(const GEMstack& gemstack) const
 template <DigitzationMode MODE>
 inline void DigitTime::fillOutputContainer(std::vector<Digit>& output, dataformats::MCTruthContainer<MCCompLabel>& mcTruth,
                                            std::vector<CommonMode>& commonModeOutput, const Sector& sector, TimeBin timeBin,
+                                           SAMPAProcessing& sampaProcessing, std::vector<std::pair<MCCompLabel, int>>& labelCollector,
                                            PrevDigitInfoArray* prevTime, Streamer* debugStream, const CalPad* padParams[3],
                                            const CalDet<bool>* deadMap)
 {
@@ -167,7 +169,7 @@ inline void DigitTime::fillOutputContainer(std::vector<Digit>& output, dataforma
         prevDigit = (*prevTime)[iPad];
       }
       const CRU cru = mapper.getCRU(sector, iPad);
-      digit.fillOutputContainer<MODE>(output, mcTruth, cru, timeBin, iPad, mLabels, getCommonMode(cru), prevDigit, debugStream, deadMap);
+      digit.fillOutputContainer<MODE>(output, mcTruth, cru, timeBin, iPad, mLabels, sampaProcessing, labelCollector, getCommonMode(cru), prevDigit, debugStream, deadMap);
     }
   }
 }

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,8 +143,28 @@ 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> mSpaceChargeOwned; //!< Distortions created by the digitizer itself
   SC* mSpaceCharge = nullptr;            //!< Handler of full distortions (static + IR dependant), owned by mSpaceChargeOwned or by the caller
   SC* mSpaceChargeDer = nullptr;         //!< Handler of reference static distortions, owned by the caller

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/DigitContainer.cxx

@@ -15,6 +15,7 @@
 
 #include "TPCSimulation/DigitContainer.h"
 #include <memory>
+#include <type_traits>
 #include <fairlogger/Logger.h>
 #include "TPCBase/Mapper.h"
 #include "TPCBaseRecSim/CDBInterface.h"
@@ -48,6 +49,9 @@ void DigitContainer::fillOutputContainer(std::vector<Digit>& output,
     updateCalibrations();
   }
 
+  // after updateCalibrations, which may reload the FEE parameters that SAMPAProcessing picks up
+  auto& sampaProcessing = mSAMPAProcessing ? *mSAMPAProcessing : SAMPAProcessing::instance();
+
   const bool needsPrevDigArray = eleParam.doIonTail || eleParam.doIonTailPerPad || eleParam.doSaturationTail;
   const bool needsEmptyTimeBins = needsPrevDigArray || eleParam.doNoiseEmptyPads;
 
@@ -89,39 +93,37 @@ void DigitContainer::fillOutputContainer(std::vector<Digit>& output,
     // fmt::print("Processing secotor: {}, time bin: {}, mFirstTimeBin: {}, dgitTime: {}\n", sector.getSector(), timeBin, mFirstTimeBin, (void*)time);
 
     if (time) {
+      // the digitization mode is a template parameter of fillOutputContainer
+      auto fill = [&](auto mode) {
+        time->fillOutputContainer<decltype(mode)::value>(output, mcTruth, commonModeOutput, sector, timeBin, sampaProcessing, mLabelCollector, mPrevDigArr.get(), debugStream, mPadParams, mDeadMap);
+      };
+      using M = DigitzationMode;
       switch (digitizationMode) {
-        case DigitzationMode::FullMode: {
-          time->fillOutputContainer<DigitzationMode::FullMode>(output, mcTruth, commonModeOutput, sector, timeBin, mPrevDigArr.get(), debugStream, mPadParams, mDeadMap);
+        case M::FullMode:
+          fill(std::integral_constant<M, M::FullMode>{});
           break;
-        }
-        case DigitzationMode::ZeroSuppression: {
-          time->fillOutputContainer<DigitzationMode::ZeroSuppression>(output, mcTruth, commonModeOutput, sector, timeBin, mPrevDigArr.get(), debugStream, mPadParams, mDeadMap);
+        case M::ZeroSuppression:
+          fill(std::integral_constant<M, M::ZeroSuppression>{});
           break;
-        }
-        case DigitzationMode::ZeroSuppressionCMCorr: {
-          time->fillOutputContainer<DigitzationMode::ZeroSuppressionCMCorr>(output, mcTruth, commonModeOutput, sector, timeBin, mPrevDigArr.get(), debugStream, mPadParams, mDeadMap);
+        case M::ZeroSuppressionCMCorr:
+          fill(std::integral_constant<M, M::ZeroSuppressionCMCorr>{});
           break;
-        }
-        case DigitzationMode::SubtractPedestal: {
-          time->fillOutputContainer<DigitzationMode::SubtractPedestal>(output, mcTruth, commonModeOutput, sector, timeBin, mPrevDigArr.get(), debugStream, mPadParams, mDeadMap);
+        case M::SubtractPedestal:
+          fill(std::integral_constant<M, M::SubtractPedestal>{});
           break;
-        }
-        case DigitzationMode::NoSaturation: {
-          time->fillOutputContainer<DigitzationMode::NoSaturation>(output, mcTruth, commonModeOutput, sector, timeBin, mPrevDigArr.get(), debugStream, mPadParams, mDeadMap);
+        case M::NoSaturation:
+          fill(std::integral_constant<M, M::NoSaturation>{});
           break;
-        }
-        case DigitzationMode::PropagateADC: {
-          time->fillOutputContainer<DigitzationMode::PropagateADC>(output, mcTruth, commonModeOutput, sector, timeBin, mPrevDigArr.get(), debugStream, mPadParams, mDeadMap);
+        case M::PropagateADC:
+          fill(std::integral_constant<M, M::PropagateADC>{});
           break;
-        }
-        case DigitzationMode::Auto: {
+        case M::Auto:
           if (mCMCEnabled) {
-            time->fillOutputContainer<DigitzationMode::ZeroSuppressionCMCorr>(output, mcTruth, commonModeOutput, sector, timeBin, mPrevDigArr.get(), debugStream, mPadParams, mDeadMap);
+            fill(std::integral_constant<M, M::ZeroSuppressionCMCorr>{});
           } else {
-            time->fillOutputContainer<DigitzationMode::ZeroSuppression>(output, mcTruth, commonModeOutput, sector, timeBin, mPrevDigArr.get(), debugStream, mPadParams, mDeadMap);
+            fill(std::integral_constant<M, M::ZeroSuppression>{});
           }
           break;
-        }
       }
     }