From cbc41d933fab3d793042934ca380652380561edf Mon Sep 17 00:00:00 2001 From: Sandro Wenzel Date: Mon, 9 Nov 2026 01:07:29 +0100 Subject: [PATCH 1/5] Pass per-digitizer workspaces down the TPC digitization output path This replaces two shared static workspaces of the TPC digitizer by members, in preparation for a multithreaded digitizer. - The shaped-signal buffer of Digitizer::process and the label sorting buffer of DigitGlobalPad::fillOutputContainer are now members. - The SAMPAProcessing used for the output is passed down from DigitContainer, which uses the singleton unless one is set. https://its.cern.ch/jira/browse/O2-7309 Co-Authored-By: Claude Opus 5.5 --- .../include/TPCSimulation/DigitContainer.h | 6 ++++++ .../include/TPCSimulation/DigitGlobalPad.h | 4 ++-- .../include/TPCSimulation/DigitTime.h | 4 +++- .../include/TPCSimulation/Digitizer.h | 1 + .../TPC/simulation/src/DigitContainer.cxx | 19 +++++++++++-------- Detectors/TPC/simulation/src/Digitizer.cxx | 7 +++---- 6 files changed, 26 insertions(+), 15 deletions(-) diff --git a/Detectors/TPC/simulation/include/TPCSimulation/DigitContainer.h b/Detectors/TPC/simulation/include/TPCSimulation/DigitContainer.h index 472a3b2c74618..c44bababfb98e 100644 --- a/Detectors/TPC/simulation/include/TPCSimulation/DigitContainer.h +++ b/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* 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> mLabelCollector; //!< Workspace to sort the MC labels of a digit void reportSettings(); }; diff --git a/Detectors/TPC/simulation/include/TPCSimulation/DigitGlobalPad.h b/Detectors/TPC/simulation/include/TPCSimulation/DigitGlobalPad.h index 743dfc3bd6977..6c08abbfca799 100644 --- a/Detectors/TPC/simulation/include/TPCSimulation/DigitGlobalPad.h +++ b/Detectors/TPC/simulation/include/TPCSimulation/DigitGlobalPad.h @@ -97,6 +97,7 @@ class DigitGlobalPad const CRU& cru, TimeBin timeBin, GlobalPadNumber globalPad, o2::dataformats::LabelContainer, false>& labelContainer, + SAMPAProcessing& sampaProcessing, std::vector>& labelCollector, float commonMode = 0.f, const PrevDigitInfo& prevDigit = PrevDigitInfo(), Streamer* debugStream = nullptr, const CalDet* deadMap = nullptr); private: @@ -210,12 +211,11 @@ inline void DigitGlobalPad::fillOutputContainer(std::vector& output, const CRU& cru, TimeBin timeBin, GlobalPadNumber globalPad, o2::dataformats::LabelContainer, false>& labels, + SAMPAProcessing& sampaProcessing, std::vector>& labelCollector, float commonMode, const PrevDigitInfo& prevDigit, Streamer* debugStream, const CalDet* deadMap) { const Mapper& mapper = Mapper::instance(); - SAMPAProcessing& sampaProcessing = SAMPAProcessing::instance(); const PadPos pad = mapper.padPos(globalPad); - static std::vector> 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 diff --git a/Detectors/TPC/simulation/include/TPCSimulation/DigitTime.h b/Detectors/TPC/simulation/include/TPCSimulation/DigitTime.h index 9f951afc2442e..48a84b4fea004 100644 --- a/Detectors/TPC/simulation/include/TPCSimulation/DigitTime.h +++ b/Detectors/TPC/simulation/include/TPCSimulation/DigitTime.h @@ -78,6 +78,7 @@ class DigitTime template void fillOutputContainer(std::vector& output, dataformats::MCTruthContainer& mcTruth, std::vector& commonModeOutput, const Sector& sector, TimeBin timeBin, + SAMPAProcessing& sampaProcessing, std::vector>& labelCollector, PrevDigitInfoArray* prevTime = nullptr, Streamer* debugStream = nullptr, const CalPad* padParams[3] = nullptr, const CalDet* deadMap = nullptr); @@ -129,6 +130,7 @@ inline float DigitTime::getCommonMode(const GEMstack& gemstack) const template inline void DigitTime::fillOutputContainer(std::vector& output, dataformats::MCTruthContainer& mcTruth, std::vector& commonModeOutput, const Sector& sector, TimeBin timeBin, + SAMPAProcessing& sampaProcessing, std::vector>& labelCollector, PrevDigitInfoArray* prevTime, Streamer* debugStream, const CalPad* padParams[3], const CalDet* deadMap) { @@ -167,7 +169,7 @@ inline void DigitTime::fillOutputContainer(std::vector& output, dataforma prevDigit = (*prevTime)[iPad]; } const CRU cru = mapper.getCRU(sector, iPad); - digit.fillOutputContainer(output, mcTruth, cru, timeBin, iPad, mLabels, getCommonMode(cru), prevDigit, debugStream, deadMap); + digit.fillOutputContainer(output, mcTruth, cru, timeBin, iPad, mLabels, sampaProcessing, labelCollector, getCommonMode(cru), prevDigit, debugStream, deadMap); } } } diff --git a/Detectors/TPC/simulation/include/TPCSimulation/Digitizer.h b/Detectors/TPC/simulation/include/TPCSimulation/Digitizer.h index 94ce8f6a9d8b1..4db3e0401675d 100644 --- a/Detectors/TPC/simulation/include/TPCSimulation/Digitizer.h +++ b/Detectors/TPC/simulation/include/TPCSimulation/Digitizer.h @@ -143,6 +143,7 @@ class Digitizer private: DigitContainer mDigitContainer; ///< Container for the Digits + std::vector mSignalArray; //!< Workspace for the shaped signal of one electron std::unique_ptr 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 diff --git a/Detectors/TPC/simulation/src/DigitContainer.cxx b/Detectors/TPC/simulation/src/DigitContainer.cxx index 4c8ca8c2002db..b5103c9907574 100644 --- a/Detectors/TPC/simulation/src/DigitContainer.cxx +++ b/Detectors/TPC/simulation/src/DigitContainer.cxx @@ -48,6 +48,9 @@ void DigitContainer::fillOutputContainer(std::vector& 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; @@ -91,34 +94,34 @@ void DigitContainer::fillOutputContainer(std::vector& output, if (time) { switch (digitizationMode) { case DigitzationMode::FullMode: { - time->fillOutputContainer(output, mcTruth, commonModeOutput, sector, timeBin, mPrevDigArr.get(), debugStream, mPadParams, mDeadMap); + time->fillOutputContainer(output, mcTruth, commonModeOutput, sector, timeBin, sampaProcessing, mLabelCollector, mPrevDigArr.get(), debugStream, mPadParams, mDeadMap); break; } case DigitzationMode::ZeroSuppression: { - time->fillOutputContainer(output, mcTruth, commonModeOutput, sector, timeBin, mPrevDigArr.get(), debugStream, mPadParams, mDeadMap); + time->fillOutputContainer(output, mcTruth, commonModeOutput, sector, timeBin, sampaProcessing, mLabelCollector, mPrevDigArr.get(), debugStream, mPadParams, mDeadMap); break; } case DigitzationMode::ZeroSuppressionCMCorr: { - time->fillOutputContainer(output, mcTruth, commonModeOutput, sector, timeBin, mPrevDigArr.get(), debugStream, mPadParams, mDeadMap); + time->fillOutputContainer(output, mcTruth, commonModeOutput, sector, timeBin, sampaProcessing, mLabelCollector, mPrevDigArr.get(), debugStream, mPadParams, mDeadMap); break; } case DigitzationMode::SubtractPedestal: { - time->fillOutputContainer(output, mcTruth, commonModeOutput, sector, timeBin, mPrevDigArr.get(), debugStream, mPadParams, mDeadMap); + time->fillOutputContainer(output, mcTruth, commonModeOutput, sector, timeBin, sampaProcessing, mLabelCollector, mPrevDigArr.get(), debugStream, mPadParams, mDeadMap); break; } case DigitzationMode::NoSaturation: { - time->fillOutputContainer(output, mcTruth, commonModeOutput, sector, timeBin, mPrevDigArr.get(), debugStream, mPadParams, mDeadMap); + time->fillOutputContainer(output, mcTruth, commonModeOutput, sector, timeBin, sampaProcessing, mLabelCollector, mPrevDigArr.get(), debugStream, mPadParams, mDeadMap); break; } case DigitzationMode::PropagateADC: { - time->fillOutputContainer(output, mcTruth, commonModeOutput, sector, timeBin, mPrevDigArr.get(), debugStream, mPadParams, mDeadMap); + time->fillOutputContainer(output, mcTruth, commonModeOutput, sector, timeBin, sampaProcessing, mLabelCollector, mPrevDigArr.get(), debugStream, mPadParams, mDeadMap); break; } case DigitzationMode::Auto: { if (mCMCEnabled) { - time->fillOutputContainer(output, mcTruth, commonModeOutput, sector, timeBin, mPrevDigArr.get(), debugStream, mPadParams, mDeadMap); + time->fillOutputContainer(output, mcTruth, commonModeOutput, sector, timeBin, sampaProcessing, mLabelCollector, mPrevDigArr.get(), debugStream, mPadParams, mDeadMap); } else { - time->fillOutputContainer(output, mcTruth, commonModeOutput, sector, timeBin, mPrevDigArr.get(), debugStream, mPadParams, mDeadMap); + time->fillOutputContainer(output, mcTruth, commonModeOutput, sector, timeBin, sampaProcessing, mLabelCollector, mPrevDigArr.get(), debugStream, mPadParams, mDeadMap); } break; } diff --git a/Detectors/TPC/simulation/src/Digitizer.cxx b/Detectors/TPC/simulation/src/Digitizer.cxx index 295d1286bb271..c87e921abe223 100644 --- a/Detectors/TPC/simulation/src/Digitizer.cxx +++ b/Detectors/TPC/simulation/src/Digitizer.cxx @@ -65,8 +65,7 @@ void Digitizer::process(const std::vector& hits, const int nShapedPoints = eleParam.NShapedPoints; const auto amplificationMode = gemParam.AmplMode; - static std::vector signalArray; - signalArray.resize(nShapedPoints); + mSignalArray.resize(nShapedPoints); /// Reserve space in the digit container for the current event mDigitContainer.reserve(sampaProcessing.getTimeBinFromTime(mEventTime - mOutputDigitTimeOffset)); @@ -154,11 +153,11 @@ void Digitizer::process(const std::vector& hits, const GlobalPadNumber globalPad = mapper.globalPadNumber(digiPadPos.getGlobalPadPos()); const float ADCsignal = sampaProcessing.getADCvalue(static_cast(nElectronsGEM)); const MCCompLabel label(MCTrackID, eventID, sourceID, false); - sampaProcessing.getShapedSignal(ADCsignal, absoluteTime, signalArray); + sampaProcessing.getShapedSignal(ADCsignal, absoluteTime, mSignalArray); for (float i = 0; i < nShapedPoints; ++i) { const float time = absoluteTime + i * eleParam.ZbinWidth; mDigitContainer.addDigit(label, digiPadPos.getCRU(), sampaProcessing.getTimeBinFromTime(time), globalPad, - signalArray[i]); + mSignalArray[i]); } /// TODO: add ion backflow to space-charge density } From b37372a80644d583f1d9babd77ec53453ee03f2e Mon Sep 17 00:00:00 2001 From: Sandro Wenzel Date: Mon, 9 Nov 2026 10:42:58 +0100 Subject: [PATCH 2/5] 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 --- .../MathUtils/include/MathUtils/RandomRing.h | 31 ++++---- .../include/TPCSimulation/Digitizer.h | 20 ++++++ .../include/TPCSimulation/ElectronTransport.h | 9 ++- .../include/TPCSimulation/GEMAmplification.h | 9 ++- .../include/TPCSimulation/SAMPAProcessing.h | 8 ++- Detectors/TPC/simulation/src/Digitizer.cxx | 72 ++++++++++++++++--- .../TPC/simulation/src/ElectronTransport.cxx | 8 ++- .../TPC/simulation/src/GEMAmplification.cxx | 27 +++++-- .../TPC/simulation/src/SAMPAProcessing.cxx | 2 +- .../src/TPCDigitizerSpec.cxx | 13 ++++ 10 files changed, 166 insertions(+), 33 deletions(-) diff --git a/Common/MathUtils/include/MathUtils/RandomRing.h b/Common/MathUtils/include/MathUtils/RandomRing.h index 59825f0a0ce09..a9915df430a37 100644 --- a/Common/MathUtils/include/MathUtils/RandomRing.h +++ b/Common/MathUtils/include/MathUtils/RandomRing.h @@ -54,7 +54,8 @@ 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 @@ -62,11 +63,13 @@ class RandomRing /// 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 -inline RandomRing::RandomRing(const RandomType randomType) +inline RandomRing::RandomRing(const RandomType randomType, TRandom* rng) : mRandomType(randomType), mRandomNumbers() { - initialize(randomType); + initialize(randomType, rng); } //______________________________________________________________________________ @@ -142,18 +149,18 @@ inline RandomRing::RandomRing(TF1& function) //______________________________________________________________________________ template -inline void RandomRing::initialize(const RandomType randomType) +inline void RandomRing::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::initialize(const RandomType randomType) //______________________________________________________________________________ template -inline void RandomRing::initialize(TF1& function) +inline void RandomRing::initialize(TF1& function, TRandom* rng) { mRandomType = RandomType::CustomTF1; for (auto& v : mRandomNumbers) { - v = function.GetRandom(); + v = function.GetRandom(rng); } } diff --git a/Detectors/TPC/simulation/include/TPCSimulation/Digitizer.h b/Detectors/TPC/simulation/include/TPCSimulation/Digitizer.h index 4db3e0401675d..e0fc1354fda5d 100644 --- a/Detectors/TPC/simulation/include/TPCSimulation/Digitizer.h +++ b/Detectors/TPC/simulation/include/TPCSimulation/Digitizer.h @@ -31,6 +31,8 @@ namespace tpc { class DigitContainer; +class ElectronTransport; +class GEMAmplification; template 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 mElectronTransport; //!< Digitizer-owned instance with sector seeding, else the singleton is used + std::unique_ptr mGEMAmplification; //!< Digitizer-owned instance with sector seeding, else the singleton is used + std::unique_ptr 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 mSignalArray; //!< Workspace for the shaped signal of one electron std::unique_ptr mSpaceChargeOwned; //!< Distortions created by the digitizer itself diff --git a/Detectors/TPC/simulation/include/TPCSimulation/ElectronTransport.h b/Detectors/TPC/simulation/include/TPCSimulation/ElectronTransport.h index 46d2aaf625d81..06a7a58d01f53 100644 --- a/Detectors/TPC/simulation/include/TPCSimulation/ElectronTransport.h +++ b/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; diff --git a/Detectors/TPC/simulation/include/TPCSimulation/GEMAmplification.h b/Detectors/TPC/simulation/include/TPCSimulation/GEMAmplification.h index 8dbfa21febc69..c30302644a07e 100644 --- a/Detectors/TPC/simulation/include/TPCSimulation/GEMAmplification.h +++ b/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; diff --git a/Detectors/TPC/simulation/include/TPCSimulation/SAMPAProcessing.h b/Detectors/TPC/simulation/include/TPCSimulation/SAMPAProcessing.h index 3f21b6723be14..8d6e7acf9d0fa 100644 --- a/Detectors/TPC/simulation/include/TPCSimulation/SAMPAProcessing.h +++ b/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 diff --git a/Detectors/TPC/simulation/src/Digitizer.cxx b/Detectors/TPC/simulation/src/Digitizer.cxx index c87e921abe223..7bd78c97b3c20 100644 --- a/Detectors/TPC/simulation/src/Digitizer.cxx +++ b/Detectors/TPC/simulation/src/Digitizer.cxx @@ -30,6 +30,7 @@ #include "TPCCalibration/CorrMapParam.h" #include +#include 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(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(&rng); + mElectronTransport = std::make_unique(&rng); + mSAMPAProcessing = std::make_unique(&rng); + mDigitContainer.setSAMPAProcessing(mSAMPAProcessing.get()); + } + const ULong64_t sectorSeed = mixSeed(mSeed ^ mixSeed(static_cast(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& 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& digits, std::vector& 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) { @@ -232,7 +288,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); diff --git a/Detectors/TPC/simulation/src/ElectronTransport.cxx b/Detectors/TPC/simulation/src/ElectronTransport.cxx index f9e36aa642158..e9afbe81fba3e 100644 --- a/Detectors/TPC/simulation/src/ElectronTransport.cxx +++ b/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()); diff --git a/Detectors/TPC/simulation/src/GEMAmplification.cxx b/Detectors/TPC/simulation/src/GEMAmplification.cxx index 8d47464e9ef53..0dadc00ee1f95 100644 --- a/Detectors/TPC/simulation/src/GEMAmplification.cxx +++ b/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(); diff --git a/Detectors/TPC/simulation/src/SAMPAProcessing.cxx b/Detectors/TPC/simulation/src/SAMPAProcessing.cxx index 462008846fa04..bee8b02f06e8f 100644 --- a/Detectors/TPC/simulation/src/SAMPAProcessing.cxx +++ b/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(); } diff --git a/Steer/DigitizerWorkflow/src/TPCDigitizerSpec.cxx b/Steer/DigitizerWorkflow/src/TPCDigitizerSpec.cxx index 8b0df106d81cc..a5d345dceb49c 100644 --- a/Steer/DigitizerWorkflow/src/TPCDigitizerSpec.cxx +++ b/Steer/DigitizerWorkflow/src/TPCDigitizerSpec.cxx @@ -44,6 +44,7 @@ #include "TPCCalibration/VDriftHelper.h" #include "CommonDataFormat/RangeReference.h" #include "SimConfig/DigiParams.h" +#include #include #include "Framework/CCDBParamSpec.h" @@ -133,6 +134,16 @@ class TPCDPLDigitizerTask : public BaseDPLDigitizer mDigitizer.setContinuousReadout(!triggeredMode); mDigitizer.setDistortionScaleType(mDistortionType); + if (ic.options().get("tpc-sector-seeding")) { + const auto seed = o2::conf::DigiParams::Instance().seed; + if (seed == 0) { + LOG(warning) << "TPC sector seeding with DigiParams.seed=0 is not reproducible"; + } + // with seed 0, gRandom was seeded randomly in the base class + mDigitizer.setUseSectorSeeding(seed != 0 ? seed : gRandom->Integer(kMaxUInt)); + LOG(info) << "TPC: using sector seeding"; + } + // we send the GRP data once if the corresponding output channel is available // and set the flag to false after mWriteGRP = true; @@ -381,6 +392,7 @@ class TPCDPLDigitizerTask : public BaseDPLDigitizer } mDigitizer.setSector(sector); + mDigitizer.setSeedingTimeframe(dh->tfCounter); mDigitizer.init(); auto& eventParts = context->getEventParts(); @@ -536,6 +548,7 @@ o2::framework::DataProcessorSpec getTPCDigitizerSpec(int channel, bool writeGRP, {"meanLumiDistortionsDerivative", VariantType::Float, -1.f, {"override lumi of derivative distortion object if >=0"}}, {"do-not-recalculate-distortions", VariantType::Bool, false, {"Do not recalculate the distortions"}}, {"n-threads-distortions", VariantType::Int, 4, {"Number of threads used for the calculation of the distortions"}}, + {"tpc-sector-seeding", VariantType::Bool, false, {"Seed the random numbers of each sector from DigiParams.seed and the sector, so that the digits do not depend on the lane assignment"}}, }}; } From 4bb74b49ce8735d26f83038a77ff02cffbcd4874 Mon Sep 17 00:00:00 2001 From: Sandro Wenzel Date: Mon, 9 Nov 2026 10:48:43 +0100 Subject: [PATCH 3/5] Fix formatting of the TPC Digitizer members --- Detectors/TPC/simulation/include/TPCSimulation/Digitizer.h | 1 + 1 file changed, 1 insertion(+) diff --git a/Detectors/TPC/simulation/include/TPCSimulation/Digitizer.h b/Detectors/TPC/simulation/include/TPCSimulation/Digitizer.h index e0fc1354fda5d..4d2fc765116c6 100644 --- a/Detectors/TPC/simulation/include/TPCSimulation/Digitizer.h +++ b/Detectors/TPC/simulation/include/TPCSimulation/Digitizer.h @@ -162,6 +162,7 @@ class Digitizer 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 mSignalArray; //!< Workspace for the shaped signal of one electron std::unique_ptr mSpaceChargeOwned; //!< Distortions created by the digitizer itself From 549f0e3de3c3ca5d5050b091f3fa31164269581a Mon Sep 17 00:00:00 2001 From: Sandro Wenzel Date: Mon, 9 Nov 2026 17:21:39 +0100 Subject: [PATCH 4/5] Keep RandomRing vector reads inside the ring This addresses a review comment on the random ring repositioning. - setRandomPosition() can place the ring at any index. - getNextValueVc() now wraps to the start when a vector would run past the end. https://its.cern.ch/jira/browse/O2-7309 Co-Authored-By: Claude Opus 5.5 --- Common/MathUtils/include/MathUtils/RandomRing.h | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/Common/MathUtils/include/MathUtils/RandomRing.h b/Common/MathUtils/include/MathUtils/RandomRing.h index a9915df430a37..d3e61696bedf7 100644 --- a/Common/MathUtils/include/MathUtils/RandomRing.h +++ b/Common/MathUtils/include/MathUtils/RandomRing.h @@ -101,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()) { @@ -113,7 +117,7 @@ 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() + /// 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); } From 7e6c65d1adc13ac2679500592fd43e15922f830a Mon Sep 17 00:00:00 2001 From: Sandro Wenzel Date: Mon, 9 Nov 2026 17:21:39 +0100 Subject: [PATCH 5/5] Use one call for all TPC digitization modes in DigitContainer This addresses a review comment on the repeated fillOutputContainer calls. - The eight cases now pass the mode as a type to one lambda. https://its.cern.ch/jira/browse/O2-7309 Co-Authored-By: Claude Opus 5.5 --- .../TPC/simulation/src/DigitContainer.cxx | 43 +++++++++---------- 1 file changed, 21 insertions(+), 22 deletions(-) diff --git a/Detectors/TPC/simulation/src/DigitContainer.cxx b/Detectors/TPC/simulation/src/DigitContainer.cxx index b5103c9907574..cb41bf9631500 100644 --- a/Detectors/TPC/simulation/src/DigitContainer.cxx +++ b/Detectors/TPC/simulation/src/DigitContainer.cxx @@ -15,6 +15,7 @@ #include "TPCSimulation/DigitContainer.h" #include +#include #include #include "TPCBase/Mapper.h" #include "TPCBaseRecSim/CDBInterface.h" @@ -92,39 +93,37 @@ void DigitContainer::fillOutputContainer(std::vector& 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(output, mcTruth, commonModeOutput, sector, timeBin, sampaProcessing, mLabelCollector, mPrevDigArr.get(), debugStream, mPadParams, mDeadMap); + }; + using M = DigitzationMode; switch (digitizationMode) { - case DigitzationMode::FullMode: { - time->fillOutputContainer(output, mcTruth, commonModeOutput, sector, timeBin, sampaProcessing, mLabelCollector, mPrevDigArr.get(), debugStream, mPadParams, mDeadMap); + case M::FullMode: + fill(std::integral_constant{}); break; - } - case DigitzationMode::ZeroSuppression: { - time->fillOutputContainer(output, mcTruth, commonModeOutput, sector, timeBin, sampaProcessing, mLabelCollector, mPrevDigArr.get(), debugStream, mPadParams, mDeadMap); + case M::ZeroSuppression: + fill(std::integral_constant{}); break; - } - case DigitzationMode::ZeroSuppressionCMCorr: { - time->fillOutputContainer(output, mcTruth, commonModeOutput, sector, timeBin, sampaProcessing, mLabelCollector, mPrevDigArr.get(), debugStream, mPadParams, mDeadMap); + case M::ZeroSuppressionCMCorr: + fill(std::integral_constant{}); break; - } - case DigitzationMode::SubtractPedestal: { - time->fillOutputContainer(output, mcTruth, commonModeOutput, sector, timeBin, sampaProcessing, mLabelCollector, mPrevDigArr.get(), debugStream, mPadParams, mDeadMap); + case M::SubtractPedestal: + fill(std::integral_constant{}); break; - } - case DigitzationMode::NoSaturation: { - time->fillOutputContainer(output, mcTruth, commonModeOutput, sector, timeBin, sampaProcessing, mLabelCollector, mPrevDigArr.get(), debugStream, mPadParams, mDeadMap); + case M::NoSaturation: + fill(std::integral_constant{}); break; - } - case DigitzationMode::PropagateADC: { - time->fillOutputContainer(output, mcTruth, commonModeOutput, sector, timeBin, sampaProcessing, mLabelCollector, mPrevDigArr.get(), debugStream, mPadParams, mDeadMap); + case M::PropagateADC: + fill(std::integral_constant{}); break; - } - case DigitzationMode::Auto: { + case M::Auto: if (mCMCEnabled) { - time->fillOutputContainer(output, mcTruth, commonModeOutput, sector, timeBin, sampaProcessing, mLabelCollector, mPrevDigArr.get(), debugStream, mPadParams, mDeadMap); + fill(std::integral_constant{}); } else { - time->fillOutputContainer(output, mcTruth, commonModeOutput, sector, timeBin, sampaProcessing, mLabelCollector, mPrevDigArr.get(), debugStream, mPadParams, mDeadMap); + fill(std::integral_constant{}); } break; - } } }