Adapt writing of aggregated TPC residuals for EPN

- Write metafile for TPC residuals
- Make output directory for residuals configurable
- Add statistics information for each TF and sort residuals in time

Author: martenole

Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/include/TPCInterpolationWorkflow/TPCResidualAggregatorSpec.h

@@ -25,6 +25,8 @@
 #include "Framework/ControlService.h"
 #include "Framework/WorkflowSpec.h"
 #include "DetectorsBase/GRPGeomHelper.h"
+#include "Framework/RawDeviceService.h"
+#include <fairmq/Device.h>
 
 using namespace o2::framework;
 
@@ -48,11 +50,24 @@ class ResidualAggregatorDevice : public o2::framework::Task
     }
     auto updateInterval = ic.options().get<uint32_t>("updateInterval");
     auto delay = ic.options().get<uint32_t>("max-delay");
+
+    std::string outputDir = o2::utils::Str::rectifyDirectory(ic.options().get<std::string>("output-dir"));
+    std::string metaFileDir = ic.options().get<std::string>("meta-output-dir");
+    bool storeMetaFile = false;
+    if (metaFileDir != "/dev/null") {
+      metaFileDir = o2::utils::Str::rectifyDirectory(metaFileDir);
+      storeMetaFile = true;
+    }
     mAggregator = std::make_unique<o2::tpc::ResidualAggregator>(minEnt);
+    mAggregator->setOutputDir(outputDir);
+    if (storeMetaFile) {
+      mAggregator->setMetaFileOutputDir(metaFileDir);
+    }
     // TODO mAggregator should get an option to set the binning externally (expose TrackResiduals::setBinning methods to user? as command line option?)
     mAggregator->setSlotLength(slotLength);
     mAggregator->setMaxSlotsDelay(delay);
     mAggregator->setCheckIntervalInfiniteSlot(updateInterval);
+    mAggregator->initOutput();
   }
 
   void finaliseCCDB(o2::framework::ConcreteDataMatcher& matcher, void* obj) final
@@ -65,6 +80,20 @@ class ResidualAggregatorDevice : public o2::framework::Task
     o2::base::GRPGeomHelper::instance().checkUpdates(pc);
     auto data = pc.inputs().get<gsl::span<o2::tpc::TrackResiduals::UnbinnedResid>>("input");
     o2::base::TFIDInfoHelper::fillTFIDInfo(pc, mAggregator->getCurrentTFInfo());
+    if (!isLHCPeriodSet) {
+      // read the LHC period information only once
+      const std::string NAStr = "NA";
+      auto LHCPeriodStr = pc.services().get<RawDeviceService>().device()->fConfig->GetProperty<std::string>("LHCPeriod", NAStr);
+      if (LHCPeriodStr == NAStr) {
+        const char* months[12] = {"JAN", "FEB", "MAR", "APR", "MAY", "JUN", "JUL", "AUG", "SEP", "OCT", "NOV", "DEC"};
+        time_t now = time(nullptr);
+        auto ltm = gmtime(&now);
+        LHCPeriodStr = months[ltm->tm_mon];
+        LOG(warning) << "LHCPeriod is not available, using current month " << LHCPeriodStr;
+      }
+      mAggregator->setLHCPeriod(LHCPeriodStr);
+      isLHCPeriodSet = true;
+    }
     LOG(debug) << "Processing TF " << mAggregator->getCurrentTFInfo().tfCounter << " with " << data.size() << " unbinned residuals";
     mAggregator->process(data);
   }
@@ -79,6 +108,7 @@ class ResidualAggregatorDevice : public o2::framework::Task
  private:
   std::unique_ptr<o2::tpc::ResidualAggregator> mAggregator;
   std::shared_ptr<o2::base::GRPGeomRequest> mCCDBRequest;
+  bool isLHCPeriodSet{false};
 };
 
 } // namespace calibration
@@ -105,7 +135,9 @@ DataProcessorSpec getTPCResidualAggregatorSpec()
       {"tf-per-slot", VariantType::UInt32, 6'000u, {"number of TFs per calibration time slot (put 0 for infinite slot length)"}},
       {"updateInterval", VariantType::UInt32, 6'000u, {"update interval in number of TFs in case slot length is infinite"}},
       {"max-delay", VariantType::UInt32, 10u, {"number of slots in past to consider"}},
-      {"min-entries", VariantType::Int, 0, {"minimum number of entries on average per voxel"}}}};
+      {"min-entries", VariantType::Int, 0, {"minimum number of entries on average per voxel"}},
+      {"output-dir", VariantType::String, "none", {"Output directory for residuals, must exist"}},
+      {"meta-output-dir", VariantType::String, "/dev/null", {"Residuals metadata output directory, must exist (if not /dev/null)"}}}};
 }
 
 } // namespace framework

Detectors/TPC/calibration/SpacePoints/CMakeLists.txt

@@ -15,6 +15,7 @@ o2_add_library(SpacePoints
                        src/TrackInterpolation.cxx
                        src/ResidualAggregator.cxx
                PUBLIC_LINK_LIBRARIES O2::DataFormatsTPC
+                                     O2::CommonUtils
                                      O2::TPCBase
                                      O2::TRDBase
                                      O2::TPCReconstruction

Detectors/TPC/calibration/SpacePoints/include/SpacePoints/ResidualAggregator.h

@@ -20,6 +20,7 @@
 #include "DetectorsCalibration/TimeSlot.h"
 #include "DataFormatsTPC/Defs.h"
 #include "SpacePoints/TrackResiduals.h"
+#include "CommonUtils/StringUtils.h"
 
 #include <vector>
 #include <array>
@@ -41,7 +42,7 @@ struct ResidualsContainer {
   ResidualsContainer& operator=(const ResidualsContainer& src) = delete;
   ~ResidualsContainer();
 
-  void init(const TrackResiduals* residualsEngine);
+  void init(const TrackResiduals* residualsEngine, std::string outputDir);
   void fillStatisticsBranches();
   uint64_t getNEntries() const { return nResidualsTotal; }
 
@@ -54,13 +55,18 @@ struct ResidualsContainer {
   std::array<std::vector<TrackResiduals::LocalResid>*, SECTORSPERSIDE * SIDES> residualsPtr{};
   std::array<std::vector<TrackResiduals::VoxStats>, SECTORSPERSIDE * SIDES> stats{}; ///< voxel statistics sent to the aggregator
   std::array<std::vector<TrackResiduals::VoxStats>*, SECTORSPERSIDE * SIDES> statsPtr{};
+  uint32_t runNumber;                                                        ///< run number (required for meta data file)
+  std::vector<uint32_t> tfOrbits, *tfOrbitsPtr{&tfOrbits};                   ///< first TF orbit
+  std::vector<uint32_t> sumOfResiduals, *sumOfResidualsPtr{&sumOfResiduals}; ///< sum of residuals for each TF
 
   std::string fileName{"o2tpc_residuals"};
   std::string treeNameResiduals{"resid"};
   std::string treeNameStats{"stats"};
+  std::string treeNameRecords{"records"};
   std::unique_ptr<TFile> fileOut{nullptr};
   std::unique_ptr<TTree> treeOutResiduals{nullptr};
   std::unique_ptr<TTree> treeOutStats{nullptr};
+  std::unique_ptr<TTree> treeOutRecords{nullptr};
 
   uint64_t nResidualsTotal{0};
 
@@ -75,13 +81,25 @@ class ResidualAggregator final : public o2::calibration::TimeSlotCalibration<Tra
   ResidualAggregator(size_t nMin = 1'000) : mMinEntries(nMin) { mTrackResiduals.init(); }
   ~ResidualAggregator() final;
 
+  void setOutputDir(std::string dir) { mOutputDir = dir.empty() ? o2::utils::Str::rectifyDirectory("./") : dir; }
+  void setMetaFileOutputDir(std::string dir)
+  {
+    mMetaOutputDir = dir;
+    mStoreMetaData = true;
+  }
+  void setLHCPeriod(std::string period) { mLHCPeriod = period; }
+
   bool hasEnoughData(const Slot& slot) const final;
   void initOutput() final;
   void finalizeSlot(Slot& slot) final;
   Slot& emplaceNewSlot(bool front, TFType tStart, TFType tEnd) final;
 
  private:
   TrackResiduals mTrackResiduals; ///< providing the functionality for voxel binning of the residuals
+  std::string mOutputDir{"./"};   ///< the directory where the output of residuals is stored
+  std::string mMetaOutputDir{"none"}; ///< the directory where the meta data file is stored
+  std::string mLHCPeriod{""};         ///< the LHC period to be put into the meta file
+  bool mStoreMetaData{false};         ///< flag, whether meta file is supposed to be stored
   size_t mMinEntries;             ///< the minimum number of residuals required for the map creation (per voxel)
 
   ClassDefOverride(ResidualAggregator, 1);

Detectors/TPC/calibration/SpacePoints/src/ResidualAggregator.cxx

@@ -15,8 +15,11 @@
 
 #include "SpacePoints/ResidualAggregator.h"
 #include "SpacePoints/SpacePointsCalibParam.h"
+#include "DetectorsCommonDataFormats/FileMetaData.h"
 
 #include <filesystem>
+#include <numeric>
+#include <algorithm>
 
 using namespace o2::tpc;
 
@@ -51,20 +54,24 @@ ResidualsContainer::ResidualsContainer(ResidualsContainer&& rhs)
   treeOutStats = std::move(rhs.treeOutStats);
   for (int iSec = 0; iSec < SECTORSPERSIDE * SIDES; ++iSec) {
     residuals[iSec] = std::move(rhs.residuals[iSec]);
-    residualsPtr[iSec] = std::move(rhs.residualsPtr[iSec]);
     stats[iSec] = std::move(rhs.stats[iSec]);
-    statsPtr[iSec] = std::move(rhs.statsPtr[iSec]);
   }
+  runNumber = rhs.runNumber;
+  tfOrbits = std::move(rhs.tfOrbits);
+  sumOfResiduals = std::move(rhs.sumOfResiduals);
 }
 
-void ResidualsContainer::init(const TrackResiduals* residualsEngine)
+void ResidualsContainer::init(const TrackResiduals* residualsEngine, std::string outputDir)
 {
   trackResiduals = residualsEngine;
   fileName += std::to_string(std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::system_clock::now().time_since_epoch()).count());
   fileName += ".root";
-  fileOut = std::make_unique<TFile>(fileName.c_str(), "recreate");
+  std::string fileNameTmp = outputDir + fileName;
+  fileNameTmp += ".part"; // to prevent premature external usage of the file use temporary name
+  fileOut = std::make_unique<TFile>(fileNameTmp.c_str(), "recreate");
   treeOutResiduals = std::make_unique<TTree>(treeNameResiduals.c_str(), "TPC binned residuals");
   treeOutStats = std::make_unique<TTree>(treeNameStats.c_str(), "Voxel statistics mean position and nEntries");
+  treeOutRecords = std::make_unique<TTree>(treeNameRecords.c_str(), "Statistics per TF slot");
   for (int iSec = 0; iSec < SECTORSPERSIDE * SIDES; ++iSec) {
     residualsPtr[iSec] = &residuals[iSec];
     statsPtr[iSec] = &stats[iSec];
@@ -81,6 +88,8 @@ void ResidualsContainer::init(const TrackResiduals* residualsEngine)
     treeOutResiduals->Branch(Form("sec%d", iSec), &residualsPtr[iSec]);
     treeOutStats->Branch(Form("sec%d", iSec), &statsPtr[iSec]);
   }
+  treeOutRecords->Branch("firstTForbit", &tfOrbitsPtr);
+  treeOutRecords->Branch("sumOfResiduals", &sumOfResidualsPtr);
 }
 
 void ResidualsContainer::fillStatisticsBranches()
@@ -89,13 +98,15 @@ void ResidualsContainer::fillStatisticsBranches()
   // remains empty and we keep the statistics in memory in the vectors
   // (since their size anyway does not change)
   treeOutStats->Fill();
+  treeOutRecords->Fill();
 }
 
 void ResidualsContainer::fill(const o2::dataformats::TFIDInfo& ti, const gsl::span<const TrackResiduals::UnbinnedResid> data)
 {
   // receives large vector of unbinned residuals and fills the sector-wise vectors
   // with binned residuals and statistics
   LOG(debug) << "Filling ResidualsContainer with vector of size " << data.size();
+  uint32_t nResidualsInTF = 0;
   for (const auto& residIn : data) {
     int sec = residIn.sec;
     auto& residVecOut = residuals[sec];
@@ -117,11 +128,15 @@ void ResidualsContainer::fill(const o2::dataformats::TFIDInfo& ti, const gsl::sp
     stat.meanPos[TrackResiduals::VoxF] = (stat.meanPos[TrackResiduals::VoxF] * oldEntries + yPos / param::RowX[residIn.row]) * norm;
     stat.meanPos[TrackResiduals::VoxZ] = (stat.meanPos[TrackResiduals::VoxZ] * oldEntries + zPos / param::RowX[residIn.row]) * norm;
     ++nResidualsTotal;
+    ++nResidualsInTF;
   }
   treeOutResiduals->Fill();
   for (auto& residVecOut : residuals) {
     residVecOut.clear();
   }
+  runNumber = ti.runNumber;
+  tfOrbits.push_back(ti.firstTForbit);
+  sumOfResiduals.push_back(nResidualsInTF);
 }
 
 void ResidualsContainer::merge(ResidualsContainer* prev)
@@ -158,6 +173,13 @@ void ResidualsContainer::merge(ResidualsContainer* prev)
   treeOutResiduals = std::move(prev->treeOutResiduals);
 
   nResidualsTotal += prev->nResidualsTotal;
+
+  // append the current vector to the vector of the previous container and afterwards swap them,
+  // since the vector of the previous container will be deleted
+  prev->tfOrbits.insert(prev->tfOrbits.end(), tfOrbits.begin(), tfOrbits.end());
+  std::swap(prev->tfOrbits, tfOrbits);
+  prev->sumOfResiduals.insert(prev->sumOfResiduals.end(), sumOfResiduals.begin(), sumOfResiduals.end());
+  std::swap(prev->sumOfResiduals, sumOfResiduals);
 }
 
 void ResidualsContainer::print()
@@ -195,15 +217,41 @@ void ResidualAggregator::finalizeSlot(Slot& slot)
   cont->treeOutResiduals.reset();
   cont->treeOutStats->Write();
   cont->treeOutStats.reset();
+  cont->treeOutRecords->Write();
+  cont->treeOutRecords.reset();
   cont->fileOut->Close();
   cont->fileOut.reset();
+  std::filesystem::rename(o2::utils::Str::concat_string(mOutputDir, cont->fileName, ".part"), mOutputDir + cont->fileName);
+  if (mStoreMetaData) {
+    o2::dataformats::FileMetaData fileMetaData; // object with information for meta data file
+    fileMetaData.fillFileData(mOutputDir + cont->fileName);
+    fileMetaData.run = cont->runNumber;
+    fileMetaData.LHCPeriod = mLHCPeriod;
+    fileMetaData.type = "calib";
+    fileMetaData.priority = "high";
+    auto metaFileNameTmp = fmt::format("{}{}.tmp", mMetaOutputDir, cont->fileName);
+    auto metaFileName = fmt::format("{}{}.done", mMetaOutputDir, cont->fileName);
+    try {
+      std::ofstream metaFileOut(metaFileNameTmp);
+      metaFileOut << fileMetaData;
+      metaFileOut << "TFOrbits: ";
+      for (size_t i = 0; i < cont->tfOrbits.size(); i++) {
+        metaFileOut << fmt::format("{}{}", i ? ", " : "", cont->tfOrbits[i]);
+      }
+      metaFileOut << '\n';
+      metaFileOut.close();
+      std::filesystem::rename(metaFileNameTmp, metaFileName);
+    } catch (std::exception const& e) {
+      LOG(error) << "Failed to store residuals meta data file " << metaFileName << ", reason: " << e.what();
+    }
+  }
 }
 
 Slot& ResidualAggregator::emplaceNewSlot(bool front, TFType tStart, TFType tEnd)
 {
   auto& cont = getSlots();
   auto& slot = front ? cont.emplace_front(tStart, tEnd) : cont.emplace_back(tStart, tEnd);
   slot.setContainer(std::make_unique<ResidualsContainer>());
-  slot.getContainer()->init(&mTrackResiduals);
+  slot.getContainer()->init(&mTrackResiduals, mOutputDir);
   return slot;
 }