TPC: Adding option to make timegain calib with gausian fits (#13641)

* TPC: Adding option to make timegain calib with gausian fits

- setting tighter default momentum range for filtering MIPs to reject
kaons and electrons

- adding option to write CalibdEdx to local file by converting the boost
histogram to a ROOT histogram

- lowering minimum dEdx to 10 for better stabillity in A16

other changes:
- TPC: make calculation of dedx parallelisable

* TPC: Fixing variable names

Author: matthias-kleiner

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

@@ -91,8 +91,8 @@ class CalibdEdxCorrection
 
   void clear();
 
-  void writeToFile(std::string_view fileName) const;
-  void loadFromFile(std::string_view fileName);
+  void writeToFile(std::string_view fileName, std::string_view objName = "CalibdEdxCorrection") const;
+  void loadFromFile(std::string_view fileName, std::string_view objName = "CalibdEdxCorrection");
 
   /// \param outFileName name of the output file
   void dumpToTree(const char* outFileName = "calib_dedx.root") const;

DataFormats/Detectors/TPC/src/CalibdEdxCorrection.cxx

@@ -36,16 +36,16 @@ void CalibdEdxCorrection::clear()
   mDims = -1;
 }
 
-void CalibdEdxCorrection::writeToFile(std::string_view fileName) const
+void CalibdEdxCorrection::writeToFile(std::string_view fileName, std::string_view objName) const
 {
   std::unique_ptr<TFile> file(TFile::Open(fileName.data(), "recreate"));
-  file->WriteObject(this, "CalibdEdxCorrection");
+  file->WriteObject(this, objName.data());
 }
 
-void CalibdEdxCorrection::loadFromFile(std::string_view fileName)
+void CalibdEdxCorrection::loadFromFile(std::string_view fileName, std::string_view objName)
 {
   std::unique_ptr<TFile> file(TFile::Open(fileName.data()));
-  auto tmp = file->Get<CalibdEdxCorrection>("CalibdEdxCorrection");
+  auto tmp = file->Get<CalibdEdxCorrection>(objName.data());
   if (tmp != nullptr) {
     *this = *tmp;
   }

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

@@ -130,7 +130,7 @@ class CalculatedEdx
   float getMinChargeMaxThreshold() { return mMinChargeMaxThreshold; }
 
   /// fill missing clusters with minimum charge (method=0) or minimum charge/2 (method=1) or Landau (method=2)
-  void fillMissingClusters(int missingClusters[4], float minChargeTot, float minChargeMax, int method);
+  void fillMissingClusters(int missingClusters[4], float minChargeTot, float minChargeMax, int method, std::array<std::vector<float>, 5>& chargeTotROC, std::array<std::vector<float>, 5>& chargeMaxROC);
 
   /// get the truncated mean for the input track with the truncation range, charge type, region and corrections
   /// the cluster charge is normalized by effective length*gain, you can turn off the normalization by setting all corrections to false
@@ -242,9 +242,6 @@ class CalculatedEdx
   CalibdEdxContainer mCalibCont;                                       ///< calibration container
   std::unique_ptr<o2::utils::TreeStreamRedirector> mStreamer{nullptr}; ///< debug streamer
 
-  std::array<std::vector<float>, 5> mChargeTotROC;
-  std::array<std::vector<float>, 5> mChargeMaxROC;
-
   CorrectdEdxDistortions mSCdEdxCorrection; ///< for space-charge correction of dE/dx
 };
 

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

@@ -33,6 +33,8 @@
 #include <boost/histogram.hpp>
 #include "THn.h"
 
+class TLinearFitter;
+
 namespace o2::tpc
 {
 
@@ -74,7 +76,7 @@ class CalibdEdx
 
   /// \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);
+  CalibdEdx(int dEdxBins = 70, float mindEdx = 10, float maxdEdx = 90, int angularBins = 36, bool fitSnp = false);
 
   void setCuts(const TrackCuts& cuts) { mCuts = cuts; }
   void setApplyCuts(bool apply) { mApplyCuts = apply; }
@@ -124,7 +126,8 @@ class CalibdEdx
 
   /// Compute MIP position from dEdx histograms and save result in the correction container.
   /// To retrieve the correction call `CalibdEdx::getCalib()`
-  void finalize();
+  /// \param useGausFits make gaussian fits of dEdx vs tgl instead of fitting the mean dEdx
+  void finalize(const bool useGausFits = true);
 
   /// Return calib data histogram
   const Hist& getHist() const { return mHist; }
@@ -170,6 +173,17 @@ class CalibdEdx
   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]; }
 
+  /// dump this object to a file - the boost histogram is converted to a ROOT histogram -
+  void dumpToFile(const char* outFile, const char* outName) const;
+
+  /// read the object from a file
+  static CalibdEdx readFromFile(const char* inFile, const char* inName);
+
+  /// set lower and upper range in units of sigma which are used for the gaussian fits
+  /// \param lowerSigma low sigma range
+  /// \param upperSigma upper sigma range
+  void setSigmaFitRange(const float lowerSigma, const float upperSigma);
+
  private:
   // ATTENTION: Adjust copy constructor
   bool mFitSnp{};
@@ -182,15 +196,26 @@ class CalibdEdx
   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
+  float mSigmaUpper = 1.;        ///< mSigma*sigma_gaus used for cutting electrons in case gaussian fits are performed
+  float mSigmaLower = 1.5;       ///< mSigma*sigma_gaus used for cutting electrons in case gaussian fits are performed
 
-  Hist mHist;                     ///< dEdx multidimensional histogram
+  Hist mHist;                     ///<! dEdx multidimensional histogram
   CalibdEdxCorrection mCalib{};   ///< Calibration output
   CalibdEdxCorrection mCalibIn{}; ///< Calibration output
 
   o2::base::Propagator::MatCorrType mMatType{}; ///< material type for track propagation
 
-  std::unique_ptr<o2::utils::TreeStreamRedirector> mDebugOutputStreamer; ///< Debug output streamer
-  ClassDefNV(CalibdEdx, 4);
+  std::unique_ptr<o2::utils::TreeStreamRedirector> mDebugOutputStreamer; ///<! Debug output streamer
+
+  THnF* getTHnF() const;
+
+  /// make fits of dEdx as a function of tgl and perform the calibration fit
+  void fitHistGaus(TLinearFitter& fitter, CalibdEdxCorrection& corr, const CalibdEdxCorrection* stackMean = nullptr);
+
+  /// set boost histogram containing the data from root histogram
+  void setFromRootHist(const THnF* hist);
+
+  ClassDefNV(CalibdEdx, 5);
 };
 
 } // namespace o2::tpc

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

@@ -60,6 +60,7 @@ class CalibratordEdx final : public o2::calibration::TimeSlotCalibration<o2::tpc
   void setApplyCuts(bool apply) { mApplyCuts = apply; }
   void setElectronCut(std::tuple<float, int, float> values) { mElectronCut = values; }
   void setMaterialType(o2::base::Propagator::MatCorrType materialType) { mMatType = materialType; }
+  void setMakeGaussianFits(const bool makeGaussianFits) { mMakeGaussianFits = makeGaussianFits; }
 
   /// \brief Check if there are enough data to compute the calibration.
   /// \return false if any of the histograms has less entries than mMinEntries
@@ -117,14 +118,15 @@ class CalibratordEdx final : public o2::calibration::TimeSlotCalibration<o2::tpc
   std::tuple<float, int, float> mElectronCut{}; ///< Values passed to CalibdEdx::setElectronCut
   TrackCuts mCuts;                              ///< Cut object
   o2::base::Propagator::MatCorrType mMatType{}; ///< material type for track propagation
+  bool mMakeGaussianFits{};                     ///< fit mean of gaussian fits instead of mean dedx
 
   TFinterval mTFIntervals;     ///< start and end time frame IDs of each calibration time slots
   TimeInterval mTimeIntervals; ///< start and end times of each calibration time slots
   CalibVector mCalibs;         ///< vector of MIP positions, each element is filled in "process" when we finalize one slot (multiple can be finalized during the same "process", which is why we have a vector. Each element is to be considered the output of the device
 
   std::unique_ptr<o2::utils::TreeStreamRedirector> mDebugOutputStreamer; ///< Debug output streamer
 
-  ClassDefOverride(CalibratordEdx, 2);
+  ClassDefOverride(CalibratordEdx, 3);
 };
 
 } // namespace o2::tpc

Detectors/TPC/calibration/src/CalculatedEdx.cxx

@@ -53,7 +53,7 @@ void CalculatedEdx::setRefit(const unsigned int nHbfPerTf)
   mRefit = std::make_unique<o2::gpu::GPUO2InterfaceRefit>(mClusterIndex, &mTPCCorrMapsHelper, mFieldNominalGPUBz, mTPCTrackClIdxVecInput->data(), nHbfPerTf, mTPCRefitterShMap.data(), mTPCRefitterOccMap.data(), mTPCRefitterOccMap.size());
 }
 
-void CalculatedEdx::fillMissingClusters(int missingClusters[4], float minChargeTot, float minChargeMax, int method)
+void CalculatedEdx::fillMissingClusters(int missingClusters[4], float minChargeTot, float minChargeMax, int method, std::array<std::vector<float>, 5>& chargeTotROC, std::array<std::vector<float>, 5>& chargeMaxROC)
 {
   if (method != 0 && method != 1) {
     LOGP(info, "Unrecognized subthreshold cluster treatment. Not adding virtual charges to the track!");
@@ -65,11 +65,11 @@ void CalculatedEdx::fillMissingClusters(int missingClusters[4], float minChargeT
       float chargeTot = (method == 1) ? minChargeTot / 2.f : minChargeTot;
       float chargeMax = (method == 1) ? minChargeMax / 2.f : minChargeMax;
 
-      mChargeTotROC[roc].emplace_back(chargeTot);
-      mChargeTotROC[4].emplace_back(chargeTot);
+      chargeTotROC[roc].emplace_back(chargeTot);
+      chargeTotROC[4].emplace_back(chargeTot);
 
-      mChargeMaxROC[roc].emplace_back(chargeMax);
-      mChargeMaxROC[4].emplace_back(chargeMax);
+      chargeMaxROC[roc].emplace_back(chargeMax);
+      chargeMaxROC[4].emplace_back(chargeMax);
     }
   }
 }
@@ -82,17 +82,13 @@ void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, fl
   int nClsROC[4] = {0, 0, 0, 0};
   int nClsSubThreshROC[4] = {0, 0, 0, 0};
 
-  mChargeTotROC[0].clear();
-  mChargeTotROC[1].clear();
-  mChargeTotROC[2].clear();
-  mChargeTotROC[3].clear();
-  mChargeTotROC[4].clear();
-
-  mChargeMaxROC[0].clear();
-  mChargeMaxROC[1].clear();
-  mChargeMaxROC[2].clear();
-  mChargeMaxROC[3].clear();
-  mChargeMaxROC[4].clear();
+  const int nType = 5;
+  std::array<std::vector<float>, nType> chargeTotROC;
+  std::array<std::vector<float>, nType> chargeMaxROC;
+  for (int i = 0; i < nType; ++i) {
+    chargeTotROC[i].reserve(Mapper::PADROWS);
+    chargeMaxROC[i].reserve(Mapper::PADROWS);
+  }
 
   // debug vectors
   std::vector<int> excludeClVector;
@@ -356,25 +352,25 @@ void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, fl
     }
 
     if (stack == GEMstack::IROCgem) {
-      mChargeTotROC[0].emplace_back(chargeTot);
-      mChargeMaxROC[0].emplace_back(chargeMax);
+      chargeTotROC[0].emplace_back(chargeTot);
+      chargeMaxROC[0].emplace_back(chargeMax);
       nClsROC[0]++;
     } else if (stack == GEMstack::OROC1gem) {
-      mChargeTotROC[1].emplace_back(chargeTot);
-      mChargeMaxROC[1].emplace_back(chargeMax);
+      chargeTotROC[1].emplace_back(chargeTot);
+      chargeMaxROC[1].emplace_back(chargeMax);
       nClsROC[1]++;
     } else if (stack == GEMstack::OROC2gem) {
-      mChargeTotROC[2].emplace_back(chargeTot);
-      mChargeMaxROC[2].emplace_back(chargeMax);
+      chargeTotROC[2].emplace_back(chargeTot);
+      chargeMaxROC[2].emplace_back(chargeMax);
       nClsROC[2]++;
     } else if (stack == GEMstack::OROC3gem) {
-      mChargeTotROC[3].emplace_back(chargeTot);
-      mChargeMaxROC[3].emplace_back(chargeMax);
+      chargeTotROC[3].emplace_back(chargeTot);
+      chargeMaxROC[3].emplace_back(chargeMax);
       nClsROC[3]++;
     };
 
-    mChargeTotROC[4].emplace_back(chargeTot);
-    mChargeMaxROC[4].emplace_back(chargeMax);
+    chargeTotROC[4].emplace_back(chargeTot);
+    chargeMaxROC[4].emplace_back(chargeMax);
 
     // for debugging
     if (mDebug) {
@@ -418,7 +414,7 @@ void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, fl
 
     // fill subthreshold clusters if not excluded
     if (((clusterMask & ClusterFlags::ExcludeSubthresholdCl) == ClusterFlags::None)) {
-      fillMissingClusters(nClsSubThreshROC, minChargeTot, minChargeMax, subthresholdMethod);
+      fillMissingClusters(nClsSubThreshROC, minChargeTot, minChargeMax, subthresholdMethod, chargeTotROC, chargeMaxROC);
     }
   } else {
     output.NHitsIROC = nClsROC[0];
@@ -428,21 +424,21 @@ void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, fl
   }
 
   // copy corrected cluster charges
-  auto chargeTotVector = mChargeTotROC[4];
-  auto chargeMaxVector = mChargeMaxROC[4];
+  auto chargeTotVector = mDebug ? chargeTotROC[4] : std::vector<float>();
+  auto chargeMaxVector = mDebug ? chargeMaxROC[4] : std::vector<float>();
 
   // calculate dEdx
-  output.dEdxTotIROC = getTruncMean(mChargeTotROC[0], low, high);
-  output.dEdxTotOROC1 = getTruncMean(mChargeTotROC[1], low, high);
-  output.dEdxTotOROC2 = getTruncMean(mChargeTotROC[2], low, high);
-  output.dEdxTotOROC3 = getTruncMean(mChargeTotROC[3], low, high);
-  output.dEdxTotTPC = getTruncMean(mChargeTotROC[4], low, high);
-
-  output.dEdxMaxIROC = getTruncMean(mChargeMaxROC[0], low, high);
-  output.dEdxMaxOROC1 = getTruncMean(mChargeMaxROC[1], low, high);
-  output.dEdxMaxOROC2 = getTruncMean(mChargeMaxROC[2], low, high);
-  output.dEdxMaxOROC3 = getTruncMean(mChargeMaxROC[3], low, high);
-  output.dEdxMaxTPC = getTruncMean(mChargeMaxROC[4], low, high);
+  output.dEdxTotIROC = getTruncMean(chargeTotROC[0], low, high);
+  output.dEdxTotOROC1 = getTruncMean(chargeTotROC[1], low, high);
+  output.dEdxTotOROC2 = getTruncMean(chargeTotROC[2], low, high);
+  output.dEdxTotOROC3 = getTruncMean(chargeTotROC[3], low, high);
+  output.dEdxTotTPC = getTruncMean(chargeTotROC[4], low, high);
+
+  output.dEdxMaxIROC = getTruncMean(chargeMaxROC[0], low, high);
+  output.dEdxMaxOROC1 = getTruncMean(chargeMaxROC[1], low, high);
+  output.dEdxMaxOROC2 = getTruncMean(chargeMaxROC[2], low, high);
+  output.dEdxMaxOROC3 = getTruncMean(chargeMaxROC[3], low, high);
+  output.dEdxMaxTPC = getTruncMean(chargeMaxROC[4], low, high);
 
   // for debugging
   if (mDebug) {