TPC: Multiple fixes and extensions (#13533)

* fix completion policy

* Add possibility for a common file name extension saving canvases

* Extend FW parameter file creation

* Allow for different thresholds for each ROC Type in pedestal file
  creation
* Add common mode k-values in float precision and inverse k-factors with
  truncated precision

* Digit class extension for simple drawing

* Extend TPC refitter

* possibility to dump ITS-TPC
* Tsallis downsampling
* MB sampling
* shared flag added to native clusters
* Charge info now via native clusters
* Write full native clusters
* add occupancy info
* add cosmics refitting with TOF cluster time

* Extend time gain calibration

* Undo previous timeGain corrections to allow for residual calibration
* slot length in seconds
* slot extension for CCDB (for residual calibration)
* simple looper cut
* settable debug output name
* dumping of calibration histograms per slot
* debug output on track level
* average entries for calibration

* Update event display

- Pad vs time view for a single row
  * visualisation of hovered time bin in ADC vs. time
- cluster overlay in pad vs. time view
  * selection of cluster flags
- Pad vs. row view for a single time bin
  * visualisation of hovered time bin in ADC vs. time

DataFormats/Detectors/TPC/include/DataFormatsTPC/CalibdEdxCorrection.h

@@ -109,6 +109,12 @@ class CalibdEdxCorrection
   /// Single fit parameters averaged over all sectors for a stack type
   float getMeanParam(const GEMstack stack, ChargeType charge, uint32_t param) const;
 
+  /// Single fit parameters averaged over all sectors for a stack type
+  float getMeanEntries(ChargeType charge) const;
+
+  /// Single fit parameters averaged over all sectors for a stack type
+  float getMeanEntries(const GEMstack stack, ChargeType charge) const;
+
 #endif
 
  private:

DataFormats/Detectors/TPC/include/DataFormatsTPC/TrackCuts.h

@@ -46,13 +46,16 @@ class TrackCuts
   void setNClusMin(float NClusMin) { mNClusMin = NClusMin; }
   void setdEdxMin(float dEdxMin) { mdEdxMin = dEdxMin; }
   void setdEdxMax(float dEdxMax) { mdEdxMax = dEdxMax; }
+  /// try to remove looper cutting on (abs(z_out) - abs(z_in)) < -10), not very precise
+  void setCutLooper(bool cut) { mCutLooper = cut; }
 
  private:
-  float mPMin{0};       ///< min momentum allowed
-  float mPMax{1e10};    ///< max momentum allowed
-  float mNClusMin{0};   ///< min number of clusters in track allowed
-  float mdEdxMin{0};    ///< min dEdx
-  float mdEdxMax{1e10}; ///< max dEdx
+  float mPMin{0};         ///< min momentum allowed
+  float mPMax{1e10};      ///< max momentum allowed
+  float mNClusMin{0};     ///< min number of clusters in track allowed
+  float mdEdxMin{0};      ///< min dEdx
+  float mdEdxMax{1e10};   ///< max dEdx
+  bool mCutLooper{false}; ///< cut looper comparing zout-zin
 
   ClassDefNV(TrackCuts, 1)
 };

DataFormats/Detectors/TPC/src/CalibdEdxCorrection.cxx

@@ -135,3 +135,23 @@ float CalibdEdxCorrection::getMeanParam(const GEMstack stack, ChargeType charge,
 
   return mean / (SECTORSPERSIDE * SIDES);
 }
+
+float CalibdEdxCorrection::getMeanEntries(ChargeType charge) const
+{
+  float mean{};
+  for (int index = 0; index < FitSize / 2; ++index) {
+    mean += mEntries[index + charge * FitSize / 2];
+  }
+
+  return mean / (0.5f * FitSize);
+}
+
+float CalibdEdxCorrection::getMeanEntries(const GEMstack stack, ChargeType charge) const
+{
+  float mean{};
+  for (int index = 0; index < SECTORSPERSIDE * SIDES; ++index) {
+    mean += getEntries(StackID{index, stack}, charge);
+  }
+
+  return mean / (SECTORSPERSIDE * SIDES);
+}

DataFormats/Detectors/TPC/src/TrackCuts.cxx

@@ -49,5 +49,8 @@ bool TrackCuts::goodTrack(o2::tpc::TrackTPC const& track)
   if (dEdx < mdEdxMin) {
     return false;
   }
+  if ((std::abs(track.getOuterParam().getZ()) - std::abs(track.getZ())) < -10) {
+    return false;
+  }
   return true;
 }

Detectors/Base/include/DetectorsBase/DPLWorkflowUtils.h

@@ -25,6 +25,7 @@
 #include "Framework/CompletionPolicy.h"
 #include "Framework/CompletionPolicyHelpers.h"
 #include "Framework/DeviceSpec.h"
+#include "Framework/Task.h"
 #include "Framework/DataSpecUtils.h"
 #include <vector>
 #include <unordered_map>

Detectors/TPC/base/include/TPCBase/CRUCalibHelpers.h

@@ -20,6 +20,7 @@
 #include <gsl/span>
 #include <filesystem>
 #include <type_traits>
+#include <vector>
 namespace fs = std::filesystem;
 
 #include "Rtypes.h"
@@ -246,7 +247,9 @@ o2::tpc::CalDet<float> getCalPad(const std::string_view fileName, const std::str
   return calPad;
 }
 
-std::unordered_map<std::string, CalPad> preparePedestalFiles(const CalPad& pedestals, const CalPad& noise, float sigmaNoise = 3, float minADC = 2, float pedestalOffset = 0, bool onlyFilled = false, bool maskBad = true, float noisyChannelThreshold = 1.5, float sigmaNoiseNoisyChannels = 4, float badChannelThreshold = 6, bool fixedSize = false);
+/// \param sigmaNoiseROCType can be either one value for all ROC types, or {IROC, OROC}, or {IROC, OROC1, OROC2, OROC3}
+/// \param minADCROCType can be either one value for all ROC types, or {IROC, OROC}, or {IROC, OROC1, OROC2, OROC3}
+std::unordered_map<std::string, CalPad> preparePedestalFiles(const CalPad& pedestals, const CalPad& noise, std::vector<float> sigmaNoiseROCType = {3, 3, 3, 3}, std::vector<float> minADCROCType = {2, 2, 2, 2}, float pedestalOffset = 0, bool onlyFilled = false, bool maskBad = true, float noisyChannelThreshold = 1.5, float sigmaNoiseNoisyChannels = 4, float badChannelThreshold = 6, bool fixedSize = false);
 
 } // namespace o2::tpc::cru_calib_helpers
 

Detectors/TPC/base/include/TPCBase/Utils.h

@@ -17,6 +17,7 @@
 /// \author Jens Wiechula, [email protected]
 ///
 
+#include <string>
 #include <vector>
 #include <string_view>
 
@@ -44,9 +45,9 @@ namespace utils
 const std::vector<std::string> tokenize(const std::string_view input, const std::string_view pattern);
 TH1* getBinInfoXY(int& binx, int& biny, float& bincx, float& bincy);
 void addFECInfo();
-void saveCanvases(TObjArray& arr, std::string_view outDir, std::string_view types = "png,pdf", std::string_view rootFileName = "");
-void saveCanvases(std::vector<TCanvas*>& canvases, std::string_view outDir, std::string_view types = "png,pdf", std::string_view rootFileName = "");
-void saveCanvas(TCanvas& c, std::string_view outDir, std::string_view types);
+void saveCanvases(TObjArray& arr, std::string_view outDir, std::string_view types = "png,pdf", std::string_view rootFileName = "", std::string nameAdd = "");
+void saveCanvases(std::vector<TCanvas*>& canvases, std::string_view outDir, std::string_view types = "png,pdf", std::string_view rootFileName = "", std::string nameAdd = "");
+void saveCanvas(TCanvas& c, std::string_view outDir, std::string_view types, std::string nameAdd = "");
 std::vector<CalPad*> readCalPads(const std::string_view fileName, const std::vector<std::string>& calPadNames);
 std::vector<CalPad*> readCalPads(const std::string_view fileName, const std::string_view calPadNames);
 
@@ -69,6 +70,19 @@ void mergeCalPads(std::string_view outputFileName, std::string_view inputFileNam
 /// \param treeTitle title of the tree
 TChain* buildChain(std::string_view command, std::string_view treeName, std::string_view treeTitle = "", bool checkSubDir = false);
 
+template <typename Iterator>
+std::string elementsToString(Iterator begin, Iterator end, const std::string separator = ", ")
+{
+  return std::accumulate(std::next(begin), end, std::to_string(*begin),
+                         [&separator](auto s, auto f) { return std::move(s) + separator + std::to_string(f); });
+}
+
+template <typename T>
+std::string elementsToString(const T& val, const std::string separator = ", ")
+{
+  return elementsToString(val.begin(), val.end(), separator);
+}
+
 } // namespace utils
 } // namespace o2::tpc
 

Detectors/TPC/base/src/CRUCalibHelpers.cxx

@@ -12,7 +12,9 @@
 #include <fstream>
 #include <iostream>
 #include <numeric>
+#include <vector>
 
+#include "TPCBase/Utils.h"
 #include "TROOT.h"
 
 #include "Framework/Logger.h"
@@ -102,10 +104,26 @@ void cru_calib_helpers::debugDiff(std::string_view file1, std::string_view file2
   }
 }
 
-std::unordered_map<std::string, CalPad> cru_calib_helpers::preparePedestalFiles(const CalPad& pedestals, const CalPad& noise, float sigmaNoise, float minADC, float pedestalOffset, bool onlyFilled, bool maskBad, float noisyChannelThreshold, float sigmaNoiseNoisyChannels, float badChannelThreshold, bool fixedSize)
+std::unordered_map<std::string, CalPad> cru_calib_helpers::preparePedestalFiles(const CalPad& pedestals, const CalPad& noise, std::vector<float> sigmaNoiseROCType, std::vector<float> minADCROCType, float pedestalOffset, bool onlyFilled, bool maskBad, float noisyChannelThreshold, float sigmaNoiseNoisyChannels, float badChannelThreshold, bool fixedSize)
 {
   const auto& mapper = Mapper::instance();
 
+  auto expandVector = [](std::vector<float>& vec, const std::string name) {
+    if (vec.size() == 1) {
+      vec.resize(4);
+      std::fill_n(&vec[1], 3, vec[0]);
+    } else if (vec.size() == 2) {
+      vec.resize(4);
+      std::fill_n(&vec[2], 2, vec[1]);
+    } else if (vec.size() != 4) {
+      LOGP(fatal, "{} definition must be either one value for all ROC types, or {{IROC, OROC}}, or {{IROC, OROC1, OROC2, OROC3}}", name);
+    }
+    LOGP(info, "Using {} = {{{}}}", name, utils::elementsToString(vec));
+  };
+
+  expandVector(sigmaNoiseROCType, "sigmaNoiseROCType");
+  expandVector(minADCROCType, "minADCROCType");
+
   std::unordered_map<std::string, CalPad> pedestalsThreshold;
   pedestalsThreshold["Pedestals"] = CalPad("Pedestals");
   pedestalsThreshold["ThresholdMap"] = CalPad("ThresholdMap");
@@ -147,6 +165,9 @@ std::unordered_map<std::string, CalPad> cru_calib_helpers::preparePedestalFiles(
       const int fecInPartition = fecInfo.getIndex() - partInfo.getSectorFECOffset();
       // const int dataWrapperID = fecInPartition >= fecOffset;
       // const int globalLinkID = (fecInPartition % fecOffset) + dataWrapperID * 12;
+      const int rocType = roc.isIROC() ? 0 : cru.partition() - 1;
+      const float sigmaNoise = sigmaNoiseROCType[rocType];
+      const float minADC = minADCROCType[rocType];
 
       const auto traceLength = traceLengths[ipad];
 
@@ -174,6 +195,7 @@ std::unordered_map<std::string, CalPad> cru_calib_helpers::preparePedestalFiles(
       float threshold = (noise > 0) ? std::max(sigmaNoise * noise, minADC) : 0;
       threshold = std::min(threshold, 1023.f);
       float thresholdHighNoise = (noiseCorr > noisyChannelThreshold) ? std::max(sigmaNoiseNoisyChannels * noise, minADC) : threshold;
+      thresholdHighNoise = std::min(thresholdHighNoise, 1023.f);
 
       float pedestalHighNoise = pedestal;
       if (noiseCorr > badChannelThreshold) {

Detectors/TPC/base/src/Utils.cxx

@@ -133,11 +133,11 @@ void utils::addFECInfo()
   h->SetTitle(title.data());
 }
 
-void utils::saveCanvases(TObjArray& arr, std::string_view outDir, std::string_view types, std::string_view rootFileName)
+void utils::saveCanvases(TObjArray& arr, std::string_view outDir, std::string_view types, std::string_view rootFileName, std::string nameAdd)
 {
   if (types.size()) {
     for (auto c : arr) {
-      utils::saveCanvas(*static_cast<TCanvas*>(c), outDir, types);
+      utils::saveCanvas(*static_cast<TCanvas*>(c), outDir, types, nameAdd);
     }
   }
 
@@ -148,24 +148,24 @@ void utils::saveCanvases(TObjArray& arr, std::string_view outDir, std::string_vi
   }
 }
 
-void utils::saveCanvases(std::vector<TCanvas*>& canvases, std::string_view outDir, std::string_view types, std::string_view rootFileName)
+void utils::saveCanvases(std::vector<TCanvas*>& canvases, std::string_view outDir, std::string_view types, std::string_view rootFileName, std::string nameAdd)
 {
   TObjArray arr;
   for (auto c : canvases) {
     arr.Add(c);
   }
 
-  saveCanvases(arr, outDir, types, rootFileName);
+  saveCanvases(arr, outDir, types, rootFileName, nameAdd);
 }
 
-void utils::saveCanvas(TCanvas& c, std::string_view outDir, std::string_view types)
+void utils::saveCanvas(TCanvas& c, std::string_view outDir, std::string_view types, std::string nameAdd)
 {
   if (!types.size()) {
     return;
   }
   const auto typesVec = tokenize(types, ",");
   for (const auto& type : typesVec) {
-    c.SaveAs(fmt::format("{}/{}.{}", outDir, c.GetName(), type).data());
+    c.SaveAs(fmt::format("{}/{}{}.{}", outDir, c.GetName(), nameAdd, type).data());
   }
 }
 

Detectors/TPC/calibration/CMakeLists.txt

@@ -55,6 +55,7 @@ o2_add_library(TPCCalibration
                        src/TPCScaler.cxx
                        src/CorrMapParam.cxx
                        src/TPCMShapeCorrection.cxx
+                       src/DigitAdd.cxx
                PUBLIC_LINK_LIBRARIES O2::DataFormatsTPC O2::TPCBase
                                      O2::TPCReconstruction ROOT::Minuit
                                      Microsoft.GSL::GSL
@@ -109,7 +110,8 @@ o2_target_root_dictionary(TPCCalibration
                                   include/TPCCalibration/CalculatedEdx.h
                                   include/TPCCalibration/TPCScaler.h
                                   include/TPCCalibration/CorrMapParam.h
-                                  include/TPCCalibration/TPCMShapeCorrection.h)
+                                  include/TPCCalibration/TPCMShapeCorrection.h
+                                  include/TPCCalibration/DigitAdd.h)
 
 o2_add_test_root_macro(macro/comparePedestalsAndNoise.C
                        PUBLIC_LINK_LIBRARIES O2::TPCBase

Detectors/TPC/calibration/include/TPCCalibration/CalibdEdx.h

@@ -19,14 +19,15 @@
 #include <cstddef>
 #include <gsl/span>
 #include <string_view>
-#include <array>
 
 // o2 includes
+#include "CommonDataFormat/TFIDInfo.h"
 #include "DataFormatsTPC/TrackTPC.h"
 #include "DataFormatsTPC/TrackCuts.h"
 #include "DataFormatsTPC/Defs.h"
 #include "DataFormatsTPC/CalibdEdxCorrection.h"
 #include "DetectorsBase/Propagator.h"
+#include "CommonUtils/TreeStreamRedirector.h"
 
 // boost includes
 #include <boost/histogram.hpp>
@@ -54,6 +55,8 @@ class CalibdEdx
   // Float axis to store data, without under and overflow bins.
   using FloatAxis = boost::histogram::axis::regular<float, boost::histogram::use_default, boost::histogram::use_default, boost::histogram::axis::option::none_t>;
 
+  using TFIDInfo = o2::dataformats::TFIDInfo;
+
   // Histogram axes types
   using AxesType = std::tuple<
     FloatAxis, // dEdx
@@ -66,6 +69,9 @@ class CalibdEdx
 
   using Hist = boost::histogram::histogram<AxesType>;
 
+  /// copy ctor
+  CalibdEdx(const CalibdEdx& other);
+
   /// \param angularBins number of bins for Tgl and Snp
   /// \param fitSnp enable Snp correction
   CalibdEdx(int dEdxBins = 60, float mindEdx = 20, float maxdEdx = 90, int angularBins = 36, bool fitSnp = false);
@@ -106,6 +112,13 @@ class CalibdEdx
   void fill(const gsl::span<const TrackTPC>);
   void fill(const std::vector<TrackTPC>& tracks) { fill(gsl::span(tracks)); }
 
+  void fill(const TFIDInfo& tfid, const gsl::span<const TrackTPC> tracks)
+  {
+    mTFID = tfid;
+    fill(tracks);
+  }
+  void fill(const TFIDInfo& tfid, const std::vector<TrackTPC>& tracks) { fill(tfid, gsl::span(tracks)); }
+
   /// Add counts from another container.
   void merge(const CalibdEdx* other);
 
@@ -120,6 +133,10 @@ class CalibdEdx
 
   const CalibdEdxCorrection& getCalib() const { return mCalib; }
 
+  /// calibration used during reconstruction
+  void setCalibrationInput(const CalibdEdxCorrection& calib) { mCalibIn = calib; }
+  const CalibdEdxCorrection& getCalibrationInput() const { return mCalibIn; }
+
   /// Return the number of hist entries of the gem stack with less statistics
   int minStackEntries() const;
 
@@ -136,12 +153,25 @@ class CalibdEdx
   /// Save the histograms to a TTree.
   void writeTTree(std::string_view fileName) const;
 
+  /// Enable debug output to file of the time slots calibrations outputs and dE/dx histograms
+  void enableDebugOutput(std::string_view fileName);
+
+  /// Disable debug output to file. Also writes and closes stored time slots.
+  void disableDebugOutput();
+
+  /// Write debug output to file
+  void finalizeDebugOutput() const;
+
+  /// \return if debug output is enabled
+  bool hasDebugOutput() const { return static_cast<bool>(mDebugOutputStreamer); }
+
   constexpr static float MipScale = 1.0 / 50.0; ///< Inverse of target dE/dx value for MIPs
 
   constexpr static float scaleTgl(float tgl, GEMstack rocType) { return tgl / conf_dedx_corr::TglScale[rocType]; }
   constexpr static float recoverTgl(float scaledTgl, GEMstack rocType) { return scaledTgl * conf_dedx_corr::TglScale[rocType]; }
 
  private:
+  // ATTENTION: Adjust copy constructor
   bool mFitSnp{};
   bool mApplyCuts{true};         ///< Wether or not to apply tracks cuts
   TrackCuts mCuts{0.3, 0.7, 60}; ///< MIP
@@ -151,13 +181,16 @@ class CalibdEdx
   float mFitCut = 0.2;           ///< dEdx cut value used to remove electron tracks
   float mFitLowCutFactor = 1.5;  ///< dEdx cut multiplier for the lower dE/dx range
   int mFitPasses = 3;            ///< number of fit passes used to remove electron tracks
+  TFIDInfo mTFID{};              ///< current TFID
 
-  Hist mHist;                   ///< dEdx multidimensional histogram
-  CalibdEdxCorrection mCalib{}; ///< Calibration output
+  Hist mHist;                     ///< dEdx multidimensional histogram
+  CalibdEdxCorrection mCalib{};   ///< Calibration output
+  CalibdEdxCorrection mCalibIn{}; ///< Calibration output
 
   o2::base::Propagator::MatCorrType mMatType{}; ///< material type for track propagation
 
-  ClassDefNV(CalibdEdx, 3);
+  std::unique_ptr<o2::utils::TreeStreamRedirector> mDebugOutputStreamer; ///< Debug output streamer
+  ClassDefNV(CalibdEdx, 4);
 };
 
 } // namespace o2::tpc

Detectors/TPC/calibration/include/TPCCalibration/CalibratordEdx.h

@@ -95,13 +95,23 @@ class CalibratordEdx final : public o2::calibration::TimeSlotCalibration<o2::tpc
   /// Write debug output to file
   void finalizeDebugOutput() const;
 
+  /// Track debugging
+  void setTrackDebug(bool debug) { mEnableTrackDebug = debug; }
+
+  /// Dump histograms
+  void setDumpHistograms(uint32_t flags) { mDumpHistograms = flags; }
+
+  uint32_t getDumpHistograms() const { return mDumpHistograms; }
+
  private:
   int mdEdxBins{};                              ///< Number of dEdx bins
   float mMindEdx{};                             ///< Minimum value for the dEdx histograms
   float mMaxdEdx{};                             ///< Maximum value for the dEdx histograms
   int mAngularBins{};                           ///< Number of bins for angular data, like Tgl and Snp
   bool mFitSnp{};                               ///< enable Snp correction
+  bool mEnableTrackDebug{false};                ///< enable track-by-track debugging
   int mMinEntries{};                            ///< Minimum amount of tracks in each time slot, to get enough statics
+  uint32_t mDumpHistograms{0};                  ///< dump histograms bitmask: 0x1 = as THn; 0x2 as TTree
   std::array<int, 3> mFitThreshold{};           ///< Minimum entries per stack to perform sector, 1D and 2D fit
   bool mApplyCuts{true};                        ///< Flag to enable tracks cuts
   std::tuple<float, int, float> mElectronCut{}; ///< Values passed to CalibdEdx::setElectronCut