Add a multithreaded TPC digitizer device

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

@@ -104,7 +104,8 @@ class DigitContainer
   std::unique_ptr<DigitTime::PrevDigitInfoArray> mPrevDigArr; ///< Keep track of ToT and ion tail cumul from last time bin
   o2::utils::DebugStreamer mStreamer;                         ///< Debug streamer
   const CalPad* mPadParams[3] = {nullptr, nullptr, nullptr};  //!< Ion tail fraction, ion tail slope and common mode k-values per pad
-  const CalDet<bool>* mDeadMap = nullptr;                     //!< Dead channel map
+  const CalDet<bool>* mDeadMap = nullptr;                     //!< Dead channel map, points to mDeadMapCopy if set
+  std::unique_ptr<CalDet<bool>> mDeadMapCopy;                 //!< Own copy of the dead channel map, which the CDBInterface rebuilds on each access
   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

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

@@ -143,15 +143,28 @@ class Digitizer
   /// in case of scaled distortions, the distortions can be recalculated to ensure consistent distortions and corrections
   void recalculateDistortions();
 
+  /// Use the distortion maps and scaling of another digitizer, which keeps ownership of the maps
+  void setSCDistortionsFrom(const Digitizer& other);
+
   /// 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
+  /// so that the digits of a sector do not depend on which sectors were processed before.
+  /// All sectors share the ring contents and differ in the start positions, as with the singleton rings.
   /// \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; }
 
+  /// Use copies of the sector-seeding random rings of another digitizer, which fills them if needed
+  void setSectorSeedingFrom(Digitizer& other);
+
+  /// Whether flush() also flushes the process-wide debug streamer
+  void setFlushDebugStreamer(bool flush) { mFlushDebugStreamer = flush; }
+
  private:
+  /// Fill the random rings of the sector seeding, if not done yet
+  void initSectorSeeding();
+
   ElectronTransport& electronTransport();
   GEMAmplification& gemAmplification();
   SAMPAProcessing& sampaProcessing();
@@ -162,6 +175,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
+  bool mFlushDebugStreamer = true;                       //!< Whether flush() flushes the debug streamer
 
   DigitContainer mDigitContainer;        ///< Container for the Digits
   std::vector<float> mSignalArray;       //!< Workspace for the shaped signal of one electron

Detectors/TPC/simulation/src/DigitContainer.cxx

@@ -14,6 +14,7 @@
 /// \author Andi Mathis, TU München, andreas.mathis@ph.tum.de
 
 #include "TPCSimulation/DigitContainer.h"
+#include <atomic>
 #include <memory>
 #include <type_traits>
 #include <fairlogger/Logger.h>
@@ -60,15 +61,6 @@ void DigitContainer::fillOutputContainer(std::vector<Digit>& output,
     mPrevDigArr = std::make_unique<DigitTime::PrevDigitInfoArray>();
   }
 
-  static bool reportedSettings = false;
-  if (!reportedSettings) {
-    reportSettings();
-    if (mDeadMap) {
-      LOGP(info, "Using dead map with {} masked pads", mDeadMap->getSum<int>());
-    }
-    reportedSettings = true;
-  }
-
   for (auto& time : mTimeBins) {
     /// the time bins between the last event and the timing of this event are uncorrelated and can be written out
     /// OR the readout is triggered (i.e. not continuous) and we can dump everything in any case, as long it is within one drift time interval
@@ -161,12 +153,25 @@ void DigitContainer::updateCalibrations()
   }
 
   // dead channel map
-  mDeadMap = eleParam.applyDeadMap ? &cdb.getDeadChannelMap() : nullptr;
+  mDeadMap = nullptr;
+  mDeadMapCopy.reset();
+  if (eleParam.applyDeadMap) {
+    mDeadMapCopy = std::make_unique<CalDet<bool>>(cdb.getDeadChannelMap());
+    mDeadMap = mDeadMapCopy.get();
+  }
 
   // common mode correction in the FEE, needed to resolve DigiMode Auto
   mCMCEnabled = (eleParam.DigiMode == DigitzationMode::Auto) && cdb.getFEEConfig().isCMCEnabled();
 
   mCalibrationsUpdated = true;
+
+  static std::atomic<bool> reportedSettings{false};
+  if (!reportedSettings.exchange(true)) {
+    reportSettings();
+    if (mDeadMap) {
+      LOGP(info, "Using dead map with {} masked pads", mDeadMap->getSum<int>());
+    }
+  }
 }
 
 void DigitContainer::reportSettings()

Detectors/TPC/simulation/src/Digitizer.cxx

@@ -31,6 +31,7 @@
 
 #include <fairlogger/Logger.h>
 #include <TRandom3.h>
+#include <mutex>
 
 ClassImp(o2::tpc::Digitizer);
 
@@ -42,6 +43,13 @@ Digitizer::Digitizer() = default;
 
 namespace
 {
+/// serialises the calls that read calibrations through the CDBInterface, which is not thread safe
+std::mutex& calibrationMutex()
+{
+  static std::mutex mutex;
+  return mutex;
+}
+
 // splitmix64 finaliser, to derive well separated seeds from the base seed
 ULong64_t mixSeed(ULong64_t x)
 {
@@ -80,17 +88,36 @@ void Digitizer::setUseSectorSeeding(ULong64_t seed)
   mSeed = seed;
 }
 
+void Digitizer::initSectorSeeding()
+{
+  if (mUseSectorSeeding && !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());
+  }
+}
+
+void Digitizer::setSectorSeedingFrom(Digitizer& other)
+{
+  other.initSectorSeeding();
+  mUseSectorSeeding = other.mUseSectorSeeding;
+  mSeed = other.mSeed;
+  if (other.mGEMAmplification) {
+    mGEMAmplification = std::make_unique<GEMAmplification>(*other.mGEMAmplification);
+    mElectronTransport = std::make_unique<ElectronTransport>(*other.mElectronTransport);
+    mSAMPAProcessing = std::make_unique<SAMPAProcessing>(*other.mSAMPAProcessing);
+    mDigitContainer.setSAMPAProcessing(mSAMPAProcessing.get());
+  }
+}
+
 void Digitizer::init()
 {
+  std::lock_guard<std::mutex> lock(calibrationMutex());
   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());
-    }
+    initSectorSeeding();
     const ULong64_t sectorSeed = mixSeed(mSeed ^ mixSeed(static_cast<ULong64_t>(mSector) + 1));
     TRandom3 rng(toTRandom3Seed(mixSeed(sectorSeed ^ mixSeed(mSeedingTimeframe))));
     mGEMAmplification->setRandomRingPositions(rng);
@@ -230,7 +257,7 @@ void Digitizer::flush(std::vector<o2::tpc::Digit>& digits,
   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) {
+  if (mFlushDebugStreamer && ((finalFlush && mIsContinuous) || (!mIsContinuous)) && mSpaceCharge) {
     o2::utils::DebugStreamer::instance()->flush();
   }
 }
@@ -289,12 +316,25 @@ void Digitizer::setStartTime(double time)
   // 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 = this->sampaProcessing();
-  sampaProcessing.updateParameters(mVDrift);
+  {
+    std::lock_guard<std::mutex> lock(calibrationMutex());
+    sampaProcessing.updateParameters(mVDrift);
+  }
   const auto timediff = time - mOutputDigitTimeOffset;
   const auto starttimebin = sampaProcessing.getTimeBinFromTime(timediff);
   mDigitContainer.setStartTime(starttimebin);
 }
 
+void Digitizer::setSCDistortionsFrom(const Digitizer& other)
+{
+  mSpaceChargeOwned.reset();
+  mSpaceCharge = other.mSpaceCharge;
+  mSpaceChargeDer = other.mSpaceChargeDer;
+  mUseSCDistortions = other.mUseSCDistortions;
+  mLumiScaleFactor = other.mLumiScaleFactor;
+  mUseScaledDistortions = other.mUseScaledDistortions;
+}
+
 void Digitizer::setLumiScaleFactor()
 {
   mLumiScaleFactor = (CorrMapParam::Instance().lumiInst - mSpaceCharge->getMeanLumi()) / mSpaceChargeDer->getMeanLumi();

Steer/DigitizerWorkflow/CMakeLists.txt

@@ -25,6 +25,8 @@ o2_add_executable(digitizer-workflow
                           src/SimReaderSpec.cxx
                           src/SimpleDigitizerWorkflow.cxx
                           src/TPCDigitizerSpec.cxx
+                          src/TPCSectorDigitization.cxx
+                          src/TPCDigitizerMTSpec.cxx
                           src/ZDCDigitizerSpec.cxx
                           src/TOFDigitizerSpec.cxx
                           $<$<BOOL:${ENABLE_UPGRADES}>:src/ITS3DigitizerSpec.cxx>

Steer/DigitizerWorkflow/src/SimpleDigitizerWorkflow.cxx

@@ -186,6 +186,7 @@ void customize(std::vector<o2::framework::ConfigParamSpec>& workflowOptions)
   // Option to write TPC digits internaly, without forwarding to a special writer instance.
   // This is useful in GRID productions with small available memory.
   workflowOptions.push_back(ConfigParamSpec{"tpc-chunked-writer", o2::framework::VariantType::Bool, false, {"Write independent TPC digit chunks as soon as they can be flushed."}});
+  workflowOptions.push_back(ConfigParamSpec{"tpc-digi-threads", o2::framework::VariantType::Int, 0, {"Digitize TPC in a single device with this many threads instead of lanes (0 = use lanes); requires --tpc-chunked-writer"}});
   workflowOptions.push_back(ConfigParamSpec{"tpc-distortion-type", o2::framework::VariantType::Int, 0, {"Simulate distortions in the TPC (0=no distortions, 1=distortions without scaling, 2=distortions with CTP scaling)"}});
 
   std::string simhelp("Comma separated list of simulation prefixes (for background, signal productions)");
@@ -624,8 +625,16 @@ WorkflowSpec defineDataProcessing(ConfigContext const& configcontext)
 
     auto internalwrite = configcontext.options().get<bool>("tpc-chunked-writer");
     auto distortionType = configcontext.options().get<int>("tpc-distortion-type");
-    WorkflowSpec tpcPipelines = o2::tpc::getTPCDigitizerSpec(lanes, tpcsectors, mctruth, internalwrite, distortionType);
-    specs.insert(specs.end(), tpcPipelines.begin(), tpcPipelines.end());
+    auto digiThreads = configcontext.options().get<int>("tpc-digi-threads");
+    if (digiThreads > 0) {
+      if (!internalwrite) {
+        LOG(fatal) << "--tpc-digi-threads requires --tpc-chunked-writer";
+      }
+      specs.emplace_back(o2::tpc::getTPCDigitizerMTSpec(digiThreads, tpcsectors, distortionType));
+    } else {
+      WorkflowSpec tpcPipelines = o2::tpc::getTPCDigitizerSpec(lanes, tpcsectors, mctruth, internalwrite, distortionType);
+      specs.insert(specs.end(), tpcPipelines.begin(), tpcPipelines.end());
+    }
 
     if (configcontext.options().get<std::string>("tpc-reco-type").empty() == false) {
       throw std::runtime_error("option 'tpc-reco-type' is deprecated, please connect workflows on the command line by pipe");