Share the scaled TPC distortions between digitizers through a file

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

@@ -141,6 +141,12 @@ class Digitizer
   /// in case of scaled distortions, the distortions can be recalculated to ensure consistent distortions and corrections
   void recalculateDistortions();
 
+  /// Use distortions that are already scaled to the instantaneous luminosity
+  void setScaledSCDistortions(std::unique_ptr<SC> spaceCharge);
+
+  /// \return the space-charge object holding the distortions, or nullptr
+  const SC* getSCDistortions() const { return mSpaceCharge.get(); }
+
  private:
   DigitContainer mDigitContainer;      ///< Container for the Digits
   std::unique_ptr<SC> mSpaceCharge;    ///< Handler of full distortions (static + IR dependant)

Detectors/TPC/simulation/src/Digitizer.cxx

@@ -235,8 +235,21 @@ void Digitizer::setStartTime(double time)
   mDigitContainer.setStartTime(starttimebin);
 }
 
+void Digitizer::setScaledSCDistortions(std::unique_ptr<SC> spaceCharge)
+{
+  mUseSCDistortions = true;
+  mSpaceCharge = std::move(spaceCharge);
+  mSpaceCharge->initAfterReadingFromFile();
+  mSpaceCharge->printMetaData();
+  mUseScaledDistortions = true;
+}
+
 void Digitizer::setLumiScaleFactor()
 {
+  if (!mSpaceCharge || !mSpaceChargeDer) {
+    LOGP(info, "Average or derivative distortions not set");
+    return;
+  }
   mLumiScaleFactor = (CorrMapParam::Instance().lumiInst - mSpaceCharge->getMeanLumi()) / mSpaceChargeDer->getMeanLumi();
   LOGP(info, "Setting Lumi scale factor: lumiInst: {}  lumi mean: {} lumi mean derivative: {} lumi scale factor: {}", CorrMapParam::Instance().lumiInst, mSpaceCharge->getMeanLumi(), mSpaceChargeDer->getMeanLumi(), mLumiScaleFactor);
 }
@@ -274,6 +287,9 @@ void Digitizer::recalculateDistortions()
     }
     // set new lumi of avg map
     mSpaceCharge->setMeanLumi(CorrMapParam::Instance().lumiInst);
+    // the corrections were only needed to compute the scaled distortions
+    mSpaceCharge->releaseGlobalCorrections();
+    mSpaceChargeDer->releaseGlobalCorrections();
   } else {
     LOGP(info, "Inst. lumi {} is same as mean lumi {}. Skip recalculation of distortions", CorrMapParam::Instance().lumiInst, mSpaceCharge->getMeanLumi());
   }

Detectors/TPC/spacecharge/include/TPCSpaceCharge/SpaceCharge.h

@@ -638,6 +638,9 @@ class SpaceCharge
   /// \param side side of the TPC
   void setDistortionLookupTables(const DataContainer& distdZ, const DataContainer& distdR, const DataContainer& distdRPhi, const Side side);
 
+  /// release the memory of the global corrections of both sides
+  void releaseGlobalCorrections();
+
   /// set the density, potential, electric fields, local distortions/corrections, global distortions/corrections from a file. Missing objects in the file are ignored.
   /// \param file output file where the electrical fields will be written to
   /// \param side of the TPC

Detectors/TPC/spacecharge/src/SpaceCharge.cxx

@@ -2126,6 +2126,16 @@ void SpaceCharge<DataT>::init()
   }
 }
 
+template <typename DataT>
+void SpaceCharge<DataT>::releaseGlobalCorrections()
+{
+  for (int iside = 0; iside < FNSIDES; ++iside) {
+    mGlobalCorrdR[iside] = DataContainer();
+    mGlobalCorrdZ[iside] = DataContainer();
+    mGlobalCorrdRPhi[iside] = DataContainer();
+  }
+}
+
 template <typename DataT>
 void SpaceCharge<DataT>::setDistortionLookupTables(const DataContainer& distdZ, const DataContainer& distdR, const DataContainer& distdRPhi, const Side side)
 {

Steer/DigitizerWorkflow/CMakeLists.txt

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

Steer/DigitizerWorkflow/src/TPCDigitizerSpec.cxx

@@ -15,6 +15,7 @@
 #include <cassert>
 #include "Framework/RootSerializationSupport.h"
 #include "TPCDigitizerSpec.h"
+#include "TPCDistortionMapFile.h"
 #include "Framework/ControlService.h"
 #include "Framework/ParallelContext.h"
 #include "Framework/ConfigParamRegistry.h"
@@ -127,6 +128,13 @@ class TPCDPLDigitizerTask : public BaseDPLDigitizer
     mUseCalibrationsFromCCDB = ic.options().get<bool>("TPCuseCCDB");
     mMeanLumiDistortions = ic.options().get<float>("meanLumiDistortions");
     mMeanLumiDistortionsDerivative = ic.options().get<float>("meanLumiDistortionsDerivative");
+    if (const auto file = ic.options().get<std::string>("tpc-distortion-file"); !file.empty()) {
+      if (mDistortionType == 2 && mRecalcDistortions) {
+        mDistortionMapFile = std::make_unique<TPCDistortionMapFile>(file);
+      } else {
+        LOG(info) << "TPC: --tpc-distortion-file is only used for distortion type 2 with recalculation";
+      }
+    }
 
     LOG(info) << "TPC calibrations from CCDB: " << mUseCalibrationsFromCCDB;
 
@@ -230,13 +238,20 @@ class TPCDPLDigitizerTask : public BaseDPLDigitizer
     pc.inputs().get<ParameterGas*>("gasparam");
     pc.inputs().get<ParameterGEM*>("gemparam");
     if (mDistortionType) {
-      pc.inputs().get<SC*>("tpcdistortions");
-      if (mDistortionType == 2) {
-        pc.inputs().get<SC*>("tpcdistortionsderiv");
-        mDigitizer.setLumiScaleFactor();
-        if (mRecalcDistortions) {
-          mDigitizer.recalculateDistortions();
+      auto computeDistortions = [&]() {
+        pc.inputs().get<SC*>("tpcdistortions");
+        if (mDistortionType == 2) {
+          pc.inputs().get<SC*>("tpcdistortionsderiv");
+          mDigitizer.setLumiScaleFactor();
+          if (mRecalcDistortions) {
+            mDigitizer.recalculateDistortions();
+          }
         }
+      };
+      if (mDistortionMapFile) {
+        mDistortionMapFile->prepare(pc, mDigitizer, mMeanLumiDistortions, mMeanLumiDistortionsDerivative, computeDistortions);
+      } else {
+        computeDistortions();
       }
     }
 
@@ -495,6 +510,7 @@ class TPCDPLDigitizerTask : public BaseDPLDigitizer
   int mDistortionType = 0;
   float mMeanLumiDistortions = -1;
   float mMeanLumiDistortionsDerivative = -1;
+  std::unique_ptr<TPCDistortionMapFile> mDistortionMapFile; ///< file sharing the scaled distortions between digitizers
   bool mRecalcDistortions = false;
 };
 
@@ -536,6 +552,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-distortion-file", VariantType::String, "", {"File in which the digitizers of one job share the scaled distortions (distortion type 2 with recalculation), computed once; keep it in the job directory"}},
     }};
 }
 

Steer/DigitizerWorkflow/src/TPCDistortionMapFile.cxx

@@ -0,0 +1,169 @@
+// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
+// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
+// All rights not expressly granted are reserved.
+//
+// This software is distributed under the terms of the GNU General Public
+// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
+//
+// In applying this license CERN does not waive the privileges and immunities
+// granted to it by virtue of its status as an Intergovernmental Organization
+// or submit itself to any jurisdiction.
+
+#include "TPCDistortionMapFile.h"
+#include "Framework/InputRecord.h"
+#include "Framework/Logger.h"
+#include "Framework/ProcessingContext.h"
+#include "TPCCalibration/CorrMapParam.h"
+#include "TPCSimulation/Digitizer.h"
+#include "TPCSpaceCharge/SpaceCharge.h"
+#include <Compression.h>
+#include <TFile.h>
+#include <TNamed.h>
+#include <TStopwatch.h>
+#include <fcntl.h>
+#include <sys/file.h>
+#include <unistd.h>
+#include <filesystem>
+#include <memory>
+#include <fmt/format.h>
+
+using namespace o2::tpc;
+using SC = o2::tpc::SpaceCharge<float>;
+
+namespace
+{
+constexpr const char* ObjectName = "ccdb_object";
+constexpr const char* KeyName = "tpc_distortion_key";
+constexpr int KeyFormat = 1; ///< version of the file content
+
+/// lock on a file, shared or exclusive, released at destruction
+class FileLock
+{
+ public:
+  FileLock(const std::string& path, int operation)
+  {
+    mFD = ::open(path.c_str(), O_CREAT | O_RDWR, 0644);
+    if (mFD >= 0 && ::flock(mFD, operation) != 0) {
+      ::close(mFD);
+      mFD = -1;
+    }
+  }
+  ~FileLock()
+  {
+    if (mFD >= 0) {
+      ::flock(mFD, LOCK_UN);
+      ::close(mFD);
+    }
+  }
+  bool ok() const { return mFD >= 0; }
+  FileLock(const FileLock&) = delete;
+  FileLock& operator=(const FileLock&) = delete;
+
+ private:
+  int mFD = -1;
+};
+
+std::string getETag(o2::framework::ProcessingContext& pc, const char* binding)
+{
+  const auto& headers = pc.inputs().get<o2::framework::CCDBMetadataExtractor>(binding);
+  const auto it = headers.find("ETag");
+  return it == headers.end() ? std::string{} : it->second;
+}
+
+std::unique_ptr<SC> readMap(const std::string& path, const std::string& key, bool warn)
+{
+  if (!std::filesystem::exists(path)) {
+    return nullptr;
+  }
+  TFile file(path.c_str(), "READ");
+  std::unique_ptr<TNamed> storedKey(file.IsZombie() ? nullptr : file.Get<TNamed>(KeyName));
+  if (!storedKey || key != storedKey->GetTitle()) {
+    if (warn) {
+      LOGP(warning, "TPC: {} was written for other distortion inputs and will be replaced", path);
+    }
+    return nullptr;
+  }
+  return std::unique_ptr<SC>(file.Get<SC>(ObjectName));
+}
+
+bool writeMap(const std::string& path, const std::string& key, const SC& sc)
+{
+  // written under another name and renamed, so that readers never see a partial file
+  const auto tmpPath = fmt::format("{}.{}.tmp", path, ::getpid());
+  bool ok = false;
+  {
+    TFile file(tmpPath.c_str(), "RECREATE", "", ROOT::RCompressionSetting::EDefaults::kUseAnalysis);
+    if (!file.IsZombie()) {
+      TNamed named(KeyName, key.c_str());
+      ok = file.WriteObjectAny(&sc, "o2::tpc::SpaceCharge<float>", ObjectName) > 0 && file.WriteObject(&named, KeyName) > 0;
+      file.Close();
+      ok = ok && !file.TestBit(TFile::kWriteError);
+    }
+  }
+  std::error_code ec;
+  if (ok) {
+    std::filesystem::rename(tmpPath, path, ec);
+  } else {
+    std::filesystem::remove(tmpPath, ec);
+  }
+  return ok && !ec;
+}
+} // namespace
+
+void TPCDistortionMapFile::prepare(o2::framework::ProcessingContext& pc, Digitizer& digitizer, float meanLumi, float meanLumiDerivative,
+                                   const std::function<void()>& compute)
+{
+  const auto etagDist = getETag(pc, "tpcdistortions");
+  const auto etagDeriv = getETag(pc, "tpcdistortionsderiv");
+  if (etagDist.empty() || etagDeriv.empty()) {
+    LOGP(warning, "TPC: no ETag for the distortion maps, the scaled distortions are computed and not shared through {}", mPath);
+    compute();
+    return;
+  }
+  const auto key = fmt::format("format {} dist {} deriv {} lumiInst {} meanLumi {} meanLumiDerivative {}", KeyFormat, etagDist, etagDeriv,
+                               CorrMapParam::Instance().lumiInst, meanLumi, meanLumiDerivative);
+  if (key == mKey) {
+    return; // the digitizer holds them already
+  }
+
+  const auto lockPath = mPath + ".lock";
+  TStopwatch timer;
+  auto read = [&](bool warn) {
+    auto sc = readMap(mPath, key, warn);
+    if (sc) {
+      digitizer.setScaledSCDistortions(std::move(sc));
+      LOGP(info, "TPC: scaled distortions read from {} in {:.1f} s", mPath, timer.RealTime());
+      mKey = key;
+    }
+    return mKey == key;
+  };
+  {
+    // digitizers that find the file read it concurrently
+    FileLock lock(lockPath, LOCK_SH);
+    if (!lock.ok()) {
+      LOGP(warning, "TPC: cannot lock {}, the scaled distortions are computed and not shared", lockPath);
+      compute();
+      return;
+    }
+    if (read(false)) {
+      return;
+    }
+  }
+  // one digitizer computes and writes, the others wait and then read
+  FileLock lock(lockPath, LOCK_EX);
+  if (!lock.ok()) {
+    LOGP(warning, "TPC: cannot lock {}, the scaled distortions are computed and not shared", lockPath);
+    compute();
+    return;
+  }
+  if (!read(true)) {
+    compute();
+    timer.Start();
+    if (writeMap(mPath, key, *digitizer.getSCDistortions())) {
+      LOGP(info, "TPC: scaled distortions written to {} in {:.1f} s", mPath, timer.RealTime());
+      mKey = key;
+    } else {
+      LOGP(warning, "TPC: cannot write the scaled distortions to {}", mPath);
+    }
+  }
+}

Steer/DigitizerWorkflow/src/TPCDistortionMapFile.h

@@ -0,0 +1,52 @@
+// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
+// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
+// All rights not expressly granted are reserved.
+//
+// This software is distributed under the terms of the GNU General Public
+// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
+//
+// In applying this license CERN does not waive the privileges and immunities
+// granted to it by virtue of its status as an Intergovernmental Organization
+// or submit itself to any jurisdiction.
+
+#ifndef STEER_DIGITIZERWORKFLOW_TPCDISTORTIONMAPFILE_H_
+#define STEER_DIGITIZERWORKFLOW_TPCDISTORTIONMAPFILE_H_
+
+#include <functional>
+#include <string>
+
+namespace o2
+{
+namespace framework
+{
+class ProcessingContext;
+}
+namespace tpc
+{
+class Digitizer;
+
+/// File with the TPC distortions scaled to the instantaneous luminosity, shared by the digitizers of a job
+/// so that the scaled distortions are computed only once.
+///
+/// The file is keyed on the ETags of the CCDB distortion maps, TPCCorrMap.lumiInst and the mean-luminosity
+/// overrides. A file written for other inputs is replaced. A lock file serialises the digitizers, so that
+/// one computes and writes while the others wait and then read.
+class TPCDistortionMapFile
+{
+ public:
+  explicit TPCDistortionMapFile(std::string path) : mPath(std::move(path)) {}
+
+  /// Give the digitizer the scaled distortions, from the file or, if it has none for these inputs,
+  /// by calling compute, which must leave them in the digitizer; they are then written to the file
+  void prepare(framework::ProcessingContext& pc, Digitizer& digitizer, float meanLumi, float meanLumiDerivative,
+               const std::function<void()>& compute);
+
+ private:
+  std::string mPath; ///< path of the file
+  std::string mKey;  ///< key of the distortions the digitizer holds
+};
+
+} // namespace tpc
+} // namespace o2
+
+#endif