TPC: adding dEdx space-charge correction (#13557)

* TPC: adding dEdx space-charge correction

- adding check for shared clusters
- adding cluster occupancy

* Adding missing include

Author: matthias-kleiner

Detectors/TPC/calibration/CMakeLists.txt

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

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

@@ -24,6 +24,7 @@
 #include "CalibdEdxContainer.h"
 #include "CorrectionMapsHelper.h"
 #include "CommonUtils/TreeStreamRedirector.h"
+#include "TPCCalibration/CorrectdEdxDistortions.h"
 
 #include <vector>
 
@@ -55,6 +56,7 @@ enum class CorrectionFlags : unsigned short {
   GainFull = 1 << 2,       ///< flag for full gain map from calibration container
   GainResidual = 1 << 3,   ///< flag for residuals gain map from calibration container
   dEdxResidual = 1 << 4,   ///< flag for residual dEdx correction
+  dEdxSC = 1 << 5,         ///< flag for space-charge dEdx correction
 };
 
 enum class ClusterFlags : unsigned short {
@@ -64,6 +66,7 @@ enum class ClusterFlags : unsigned short {
   ExcludeEdgeCl = 1 << 2,           ///< flag to exclude sector edge clusters in dEdx calculation
   ExcludeSubthresholdCl = 1 << 3,   ///< flag to exclude subthreshold clusters in dEdx calculation
   ExcludeSectorBoundaries = 1 << 4, ///< flag to exclude sector boundary clusters in subthreshold cluster treatment
+  ExcludeSharedCl = 1 << 5,         ///< flag to exclude clusters shared between tracks
 };
 
 inline CorrectionFlags operator&(CorrectionFlags a, CorrectionFlags b) { return static_cast<CorrectionFlags>(static_cast<unsigned short>(a) & static_cast<unsigned short>(b)); }
@@ -111,6 +114,9 @@ class CalculatedEdx
   /// set the debug streamer
   void setStreamer(const char* debugRootFile) { mStreamer = std::make_unique<o2::utils::TreeStreamRedirector>(debugRootFile, "recreate"); };
 
+  /// set the debug streamer of the space-charge dedx correction
+  void setSCStreamer(const char* debugRootFile = "debug_sc_corrections.root") { mSCdEdxCorrection.setStreamer(debugRootFile); }
+
   /// \return returns magnetic field in kG
   float getFieldNominalGPUBz() { return mFieldNominalGPUBz; }
 
@@ -160,7 +166,8 @@ class CalculatedEdx
 
   /// load calibration objects from CCDB
   /// \param runNumberOrTimeStamp run number or time stamp
-  void loadCalibsFromCCDB(long runNumberOrTimeStamp);
+  /// \param isMC set if dEdx space-charge corrections will be loaded for MC or real data
+  void loadCalibsFromCCDB(long runNumberOrTimeStamp, const bool isMC = false);
 
   /// load calibration objects from local CCDB folder
   /// \param localCCDBFolder local CCDB folder
@@ -208,6 +215,15 @@ class CalculatedEdx
   /// \param object name of the object to load
   void setPropagatorFromFile(const char* folder, const char* file, const char* object);
 
+  /// \param lumi set luminosity for space-charge correction map scaling
+  void setLumi(const float lumi) { mSCdEdxCorrection.setLumi(lumi); }
+
+  /// \return returns space-charge dedx correctin
+  auto& getSCCorrection() { return mSCdEdxCorrection; }
+
+  /// \return returns cluster occupancy for given cluster
+  unsigned int getOccupancy(const o2::tpc::ClusterNative& cl) const;
+
  private:
   std::vector<TrackTPC>* mTracks{nullptr};                       ///< vector containing the tpc tracks which will be processed
   std::vector<TPCClRefElem>* mTPCTrackClIdxVecInput{nullptr};    ///< input vector with TPC tracks cluster indicies
@@ -228,6 +244,8 @@ class CalculatedEdx
 
   std::array<std::vector<float>, 5> mChargeTotROC;
   std::array<std::vector<float>, 5> mChargeMaxROC;
+
+  CorrectdEdxDistortions mSCdEdxCorrection; ///< for space-charge correction of dE/dx
 };
 
 } // namespace o2::tpc

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

@@ -0,0 +1,101 @@
+// 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.
+
+///
+/// @file   CorrectdEdxDistortions.h
+/// @author Matthias Kleiner, mkleiner@ikf.uni-frankfurt.de
+///
+
+#ifndef ALICEO2_TPC_CORRECTDEDXDISTORTIONS_H
+#define ALICEO2_TPC_CORRECTDEDXDISTORTIONS_H
+
+#include "GPUTPCGeometry.h"
+#include <memory>
+
+// forward declare classes
+namespace o2::gpu
+{
+class TPCFastTransform;
+}
+
+namespace o2::utils
+{
+class TreeStreamRedirector;
+}
+
+namespace o2::tpc
+{
+
+class CorrectdEdxDistortions
+{
+ public:
+  /// Default constructor
+  CorrectdEdxDistortions();
+
+  /// Default destructor
+  ~CorrectdEdxDistortions();
+
+  /// setting the space-charge corrections from local file
+  /// \param scAvgFile average space-charge correction map
+  /// \param scDerFile derivative space-charge correction map
+  /// \param lumi luminosity which is used for scaling the correction map
+  void setSCCorrFromFile(const char* scAvgFile, const char* scDerFile, const float lumi = -1);
+
+  /// setting the space-charge corrections
+  /// \param avg average space-charge correction map
+  /// \param der derivative space-charge correction map
+  /// \param lumi luminosity which is used for scaling the correction map
+  void setCorrectionMaps(o2::gpu::TPCFastTransform* avg, o2::gpu::TPCFastTransform* der, const float lumi);
+
+  /// setting the space-charge corrections
+  /// \param avg average space-charge correction map
+  /// \param der derivative space-charge correction map
+  void setCorrectionMaps(o2::gpu::TPCFastTransform* avg, o2::gpu::TPCFastTransform* der);
+
+  /// \param lumi luminosity which is used for scaling the correction map
+  void setLumi(float lumi);
+
+  /// enable the debug streamer
+  void setStreamer(const char* debugRootFile = "debug_sc_corrections.root");
+
+  /// \return getting the correction value for the cluster charge
+  /// \param time true time information of the cluster
+  /// \param sector TPC sector of the cluster
+  /// \param padrow global padrow
+  /// \param pad pad number in the padrow
+  float getCorrection(const float time, unsigned char sector, unsigned char padrow, int pad) const;
+
+  /// set minimum allowed lx correction at lower pad
+  void setMinLX0(const float minLX) { mLX0Min = minLX; }
+
+  /// set minimum allowed lx correction upper lower pad
+  void setMinLX1(const float minLX) { mLX1Min = minLX; }
+
+  /// set minimum allowed correction value
+  void setMinCorr(const float minCorr) { mScCorrMin = minCorr; }
+
+  /// set maximum allowed correction value
+  void setMaxCorr(const float maxCorr) { mScCorrMax = maxCorr; }
+
+ private:
+  float mScaleDer = 0;
+  std::unique_ptr<o2::gpu::TPCFastTransform> mCorrAvg{nullptr};        ///< correction map for average distortions
+  std::unique_ptr<o2::gpu::TPCFastTransform> mCorrDer{nullptr};        ///< correction map for distortions
+  o2::gpu::GPUTPCGeometry mTPCGeometry;                                ///< geometry information of TPC
+  std::unique_ptr<o2::utils::TreeStreamRedirector> mStreamer{nullptr}; ///< debug streamer
+  float mLX0Min{82};                                                   ///< minimum allowed lx position after correction at position of the bottom pad
+  float mLX1Min{82};                                                   ///< minimum allowed lx position after correction at position of the top pad
+  float mScCorrMin{0.7};                                               ///< minimum allowed correction values
+  float mScCorrMax{1.6};                                               ///< maximum allowed correction values
+};
+} // namespace o2::tpc
+
+#endif

Detectors/TPC/calibration/src/CalculatedEdx.cxx

@@ -26,6 +26,7 @@
 #include "CalibdEdxTrackTopologyPol.h"
 #include "DataFormatsParameters/GRPMagField.h"
 #include "GPUO2InterfaceUtils.h"
+#include "GPUTPCGMMergedTrackHit.h"
 
 using namespace o2::tpc;
 
@@ -111,9 +112,12 @@ void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, fl
   std::vector<float> gainResidualVector;
   std::vector<float> residualCorrTotVector;
   std::vector<float> residualCorrMaxVector;
+  std::vector<float> scCorrVector;
 
   std::vector<o2::tpc::TrackTPC> trackVector;
   std::vector<o2::tpc::ClusterNative> clVector;
+  std::vector<unsigned int> occupancyVector;
+  std::vector<bool> isClusterShared;
 
   if (mDebug) {
     excludeClVector.reserve(nClusters);
@@ -134,6 +138,9 @@ void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, fl
     residualCorrMaxVector.reserve(nClusters);
     trackVector.reserve(nClusters);
     clVector.reserve(nClusters);
+    scCorrVector.reserve(nClusters);
+    occupancyVector.reserve(nClusters);
+    isClusterShared.reserve(nClusters);
   }
 
   // for missing clusters
@@ -154,6 +161,10 @@ void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, fl
     // set sectorIndex, rowIndex, clusterIndexNumb
     track.getClusterReference(*mTPCTrackClIdxVecInput, iCl, sectorIndex, rowIndex, clusterIndexNumb);
 
+    // check if the cluster is shared
+    const unsigned int absoluteIndex = mClusterIndex->clusterOffset[sectorIndex][rowIndex] + clusterIndexNumb;
+    const bool isShared = mRefit ? (mTPCRefitterShMap[absoluteIndex] & GPUCA_NAMESPACE::gpu::GPUTPCGMMergedTrackHit::flagShared) : 0;
+
     // get region, pad, stack and stack ID
     const int region = Mapper::REGION[rowIndex];
     const unsigned char pad = std::clamp(static_cast<unsigned int>(cl.getPad() + 0.5f), static_cast<unsigned int>(0), Mapper::PADSPERROW[region][Mapper::getLocalRowFromGlobalRow(rowIndex)] - 1); // the left side of the pad is defined at e.g. 3.5 and the right side at 4.5
@@ -179,6 +190,9 @@ void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, fl
     if (((clusterMask & ClusterFlags::ExcludeEdgeCl) == ClusterFlags::ExcludeEdgeCl) && ((flagsCl & ClusterNative::flagEdge) == ClusterNative::flagEdge)) {
       excludeCl += 0b100; // 4 for edge cluster
     }
+    if (((clusterMask & ClusterFlags::ExcludeSharedCl) == ClusterFlags::ExcludeSharedCl) && isShared) {
+      excludeCl += 0b10000; // for shared cluster
+    }
 
     // get the x position of the track
     const float xPosition = Mapper::instance().getPadCentre(PadPos(rowIndex, 0)).X();
@@ -214,6 +228,8 @@ void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, fl
         offsPadVector.emplace_back(offsPad);
         trackVector.emplace_back(track);
         clVector.emplace_back(cl);
+        occupancyVector.emplace_back(getOccupancy(cl));
+        isClusterShared.emplace_back(isShared);
 
         topologyCorrVector.emplace_back(-999.f);
         topologyCorrTotVector.emplace_back(-999.f);
@@ -222,6 +238,7 @@ void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, fl
         gainResidualVector.emplace_back(-999.f);
         residualCorrTotVector.emplace_back(-999.f);
         residualCorrMaxVector.emplace_back(-999.f);
+        scCorrVector.emplace_back(-999.f);
       }
       // to avoid counting the skipped cluster as a subthreshold cluster
       rowIndexOld = rowIndex;
@@ -325,6 +342,19 @@ void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, fl
       minChargeMax = chargeMax;
     };
 
+    // space-charge dEdx corrections
+    const float time = cl.getTime() - track.getTime0(); // ToDo: get correct time from ITS-TPC track if possible
+    float scCorr = 1.0f;
+    if ((correctionMask & CorrectionFlags::dEdxSC) == CorrectionFlags::dEdxSC) {
+      scCorr = mSCdEdxCorrection.getCorrection(time, sectorIndex, rowIndex, pad);
+      if (scCorr > 0) {
+        chargeTot /= scCorr;
+      };
+      if (corrMax > 0) {
+        chargeMax /= scCorr;
+      };
+    }
+
     if (stack == GEMstack::IROCgem) {
       mChargeTotROC[0].emplace_back(chargeTot);
       mChargeMaxROC[0].emplace_back(chargeMax);
@@ -359,6 +389,8 @@ void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, fl
       offsPadVector.emplace_back(offsPad);
       trackVector.emplace_back(track);
       clVector.emplace_back(cl);
+      occupancyVector.emplace_back(getOccupancy(cl));
+      isClusterShared.emplace_back(isShared);
 
       topologyCorrVector.emplace_back(effectiveLength);
       topologyCorrTotVector.emplace_back(effectiveLengthTot);
@@ -367,6 +399,7 @@ void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, fl
       gainResidualVector.emplace_back(gainResidual);
       residualCorrTotVector.emplace_back(corrTot);
       residualCorrMaxVector.emplace_back(corrMax);
+      scCorrVector.emplace_back(scCorr);
     };
   }
 
@@ -394,6 +427,10 @@ void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, fl
     output.NHitsOROC3 = nClsROC[3];
   }
 
+  // copy corrected cluster charges
+  auto chargeTotVector = mChargeTotROC[4];
+  auto chargeMaxVector = mChargeMaxROC[4];
+
   // calculate dEdx
   output.dEdxTotIROC = getTruncMean(mChargeTotROC[0], low, high);
   output.dEdxTotOROC1 = getTruncMean(mChargeTotROC[1], low, high);
@@ -428,6 +465,7 @@ void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, fl
                  << "gainResidualVector=" << gainResidualVector
                  << "residualCorrTotVector=" << residualCorrTotVector
                  << "residualCorrMaxVector=" << residualCorrMaxVector
+                 << "scCorrVector=" << scCorrVector
                  << "localXVector=" << localXVector
                  << "localYVector=" << localYVector
                  << "offsPadVector=" << offsPadVector
@@ -436,6 +474,10 @@ void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, fl
                  << "minChargeTot=" << minChargeTot
                  << "minChargeMax=" << minChargeMax
                  << "output=" << output
+                 << "occupancy=" << occupancyVector
+                 << "chargeTotVector=" << chargeTotVector
+                 << "chargeMaxVector=" << chargeMaxVector
+                 << "isClusterShared=" << isClusterShared
                  << "\n";
   }
 }
@@ -492,7 +534,7 @@ float CalculatedEdx::getTrackTopologyCorrectionPol(const o2::tpc::TrackTPC& trac
   return effectiveLength;
 }
 
-void CalculatedEdx::loadCalibsFromCCDB(long runNumberOrTimeStamp)
+void CalculatedEdx::loadCalibsFromCCDB(long runNumberOrTimeStamp, const bool isMC)
 {
   // setup CCDB manager
   auto& cm = o2::ccdb::BasicCCDBManager::instance();
@@ -542,6 +584,15 @@ void CalculatedEdx::loadCalibsFromCCDB(long runNumberOrTimeStamp)
   auto propagator = o2::base::Propagator::Instance();
   const o2::base::MatLayerCylSet* matLut = o2::base::MatLayerCylSet::rectifyPtrFromFile(cm.get<o2::base::MatLayerCylSet>("GLO/Param/MatLUT"));
   propagator->setMatLUT(matLut);
+
+  // load sc correction maps
+  auto avgMap = isMC ? cm.getForTimeStamp<o2::gpu::TPCFastTransform>(o2::tpc::CDBTypeMap.at(o2::tpc::CDBType::CalCorrMapMC), tRun) : cm.getForTimeStamp<o2::gpu::TPCFastTransform>(o2::tpc::CDBTypeMap.at(o2::tpc::CDBType::CalCorrMap), tRun);
+  avgMap->rectifyAfterReadingFromFile();
+
+  auto derMap = isMC ? cm.getForTimeStamp<o2::gpu::TPCFastTransform>(o2::tpc::CDBTypeMap.at(o2::tpc::CDBType::CalCorrDerivMapMC), tRun) : cm.getForTimeStamp<o2::gpu::TPCFastTransform>(o2::tpc::CDBTypeMap.at(o2::tpc::CDBType::CalCorrDerivMap), tRun);
+  derMap->rectifyAfterReadingFromFile();
+
+  mSCdEdxCorrection.setCorrectionMaps(avgMap, derMap);
 }
 
 void CalculatedEdx::loadCalibsFromLocalCCDBFolder(const char* localCCDBFolder)
@@ -624,4 +675,12 @@ void CalculatedEdx::setPropagatorFromFile(const char* folder, const char* file,
     o2::base::MatLayerCylSet* matLut = o2::base::MatLayerCylSet::rectifyPtrFromFile((o2::base::MatLayerCylSet*)matLutFile->Get(object));
     propagator->setMatLUT(matLut);
   }
-}
+}
+
+unsigned int CalculatedEdx::getOccupancy(const o2::tpc::ClusterNative& cl) const
+{
+  const int nTimeBinsPerOccupBin = 16;
+  const int iBinOcc = cl.getTime() / nTimeBinsPerOccupBin + 2;
+  const unsigned int occupancy = mTPCRefitterOccMap.empty() ? -1 : mTPCRefitterOccMap[iBinOcc];
+  return occupancy;
+}