Fix in the cluster stat. extraction + clusters sigmas (#13648)

* Fix in the cluster stat. extraction

* Add clusters sigmas to the output

Author: shahor02

Detectors/GlobalTrackingWorkflow/study/CMakeLists.txt

@@ -9,7 +9,7 @@
 # granted to it by virtue of its status as an Intergovernmental Organization
 # or submit itself to any jurisdiction.
 
-# add_compile_options(-O0 -g -fPIC)
+#add_compile_options(-O0 -g -fPIC)
 
 o2_add_library(GlobalTrackingStudy
                SOURCES src/TPCTrackStudy.cxx

Detectors/GlobalTrackingWorkflow/study/include/GlobalTrackingStudy/TrackInfoExt.h

@@ -37,10 +37,13 @@ struct TrackInfoExt {
   float q2ptTPC = 0.f;
   float q2ptITSTPC = 0.f;
   float q2ptITSTPCTRD = 0.f;
-  int nClTPC = 0;
-  int nClITS = 0;
-  int pattITS = 0;
-  ClassDefNV(TrackInfoExt, 1);
+  uint16_t nClTPC = 0;
+  uint8_t pattITS = 0;
+  uint8_t nClITS = 0;
+  uint8_t rowMinTPC = 0;
+  uint8_t rowMaxTPC = 0;
+  uint8_t rowCountTPC = 0;
+  ClassDefNV(TrackInfoExt, 2);
 };
 
 } // namespace dataformats

Detectors/GlobalTrackingWorkflow/study/include/GlobalTrackingStudy/TrackMCStudyTypes.h

@@ -46,12 +46,13 @@ struct MCTrackInfo {
   int8_t parentDecID = -1;
   uint8_t minTPCRow = -1;
   uint8_t maxTPCRow = 0;
+  uint8_t nUsedPadRows = 0;
   uint8_t maxTPCRowInner = 0; // highest row in the sector containing the lowest one
   uint8_t minTPCRowSect = -1;
   uint8_t maxTPCRowSect = -1;
   int8_t nITSCl = 0;
   int8_t pattITSCl = 0;
-  ClassDefNV(MCTrackInfo, 3);
+  ClassDefNV(MCTrackInfo, 4);
 };
 
 struct RecTrack {
@@ -127,6 +128,10 @@ struct ClResTPCCont {
 
   int getNExt() const { return (below[0] > 1.) + (above[0] > 1.); }
 
+  float getClX() const { return xyz[0]; }
+  float getClY() const { return xyz[1]; }
+  float getClZ() const { return xyz[2]; }
+
   float getDY() const { return xyz[1] - getYRef(); }
   float getDZ() const { return xyz[2] - getZRef(); }
 
@@ -195,9 +200,14 @@ struct ClResTPC {
   uint8_t row = 0;
   uint8_t ncont = 0;
   uint8_t flags = 0;
+  uint8_t sigmaTimePacked;
+  uint8_t sigmaPadPacked;
   float qmax = 0;
   float qtot = 0;
   float occ = 0;
+  float occBin = 0;
+  float getSigmaPad() const { return float(sigmaPadPacked) * (1.f / 32); }
+  float getSigmaTime() const { return float(sigmaTimePacked) * (1.f / 32); }
 
   std::vector<ClResTPCCont> contTracks;
   int getNCont() const { return contTracks.size(); }

Detectors/GlobalTrackingWorkflow/study/src/TrackMCStudy.cxx

@@ -117,6 +117,7 @@ class TrackMCStudy : public Task
   std::shared_ptr<o2::base::GRPGeomRequest> mGGCCDBRequest;
   std::unique_ptr<o2::utils::TreeStreamRedirector> mDBGOut;
   std::vector<float> mTBinClOcc; ///< TPC occupancy histo: i-th entry is the integrated occupancy for ~1 orbit starting from the TB = i*mNTPCOccBinLength
+  std::vector<float> mTBinClOccHist; //< original occupancy
   std::vector<long> mIntBC;      ///< interaction global BC wrt TF start
   std::vector<float> mTPCOcc;    ///< TPC occupancy for this interaction time
   std::vector<int> mITSOcc;      //< N ITS clusters in the ROF containing collision
@@ -620,15 +621,15 @@ void TrackMCStudy::process(const o2::globaltracking::RecoContainer& recoData)
     }
   }
 
+  // collect ITS/TPC cluster info for selected MC particles
+  fillMCClusterInfo(recoData);
+
   // single tracks
   for (auto& entry : mSelMCTracks) {
     auto& trackFam = entry.second;
     (*mDBGOut) << "tracks" << "tr=" << trackFam << "\n";
   }
 
-  // collect ITS/TPC cluster info for selected MC particles
-  fillMCClusterInfo(recoData);
-
   // decays
   std::vector<TrackFamily> decFam;
   for (int id = 0; id < mNCheckDecays; id++) {
@@ -735,6 +736,7 @@ void TrackMCStudy::fillMCClusterInfo(const o2::globaltracking::RecoContainer& re
           ncontLb++;
         }
         const auto& clus = TPCClusterIdxStruct.clusters[sector][row][icl0];
+        int tbinH = int(clus.getTime() * mNTPCOccBinLengthInv); // time bin converted to slot of the occ. histo
         clRes.contTracks.clear();
         bool doClusRes = (params.minTPCRefsToExtractClRes > 0) && (params.rejectClustersResStat <= 0. || gRandom->Rndm() < params.rejectClustersResStat);
         for (auto lbl : labels) {
@@ -749,6 +751,9 @@ void TrackMCStudy::fillMCClusterInfo(const o2::globaltracking::RecoContainer& re
           if (row > mctr.maxTPCRow) {
             mctr.maxTPCRow = row;
             mctr.maxTPCRowSect = sector;
+            mctr.nUsedPadRows++;
+          } else if (row == 0 && mctr.nUsedPadRows == 0) {
+            mctr.nUsedPadRows++;
           }
           if (row < mctr.minTPCRow) {
             mctr.minTPCRow = row;
@@ -838,8 +843,18 @@ void TrackMCStudy::fillMCClusterInfo(const o2::globaltracking::RecoContainer& re
           clRes.qtot = clus.getQtot();
           clRes.qmax = clus.getQmax();
           clRes.flags = clus.getFlags();
+          clRes.sigmaTimePacked = clus.sigmaTimePacked;
+          clRes.sigmaPadPacked = clus.sigmaPadPacked;
           clRes.ncont = ncontLb;
           clRes.sortCont();
+
+          if (tbinH < 0) {
+            tbinH = 0;
+          } else if (tbinH >= int(mTBinClOccHist.size())) {
+            tbinH = (int)mTBinClOccHist.size() - 1;
+          }
+          clRes.occBin = mTBinClOccHist[tbinH];
+
           (*mDBGOut) << "clres" << "clr=" << clRes << "\n";
         }
       }
@@ -1147,21 +1162,22 @@ void TrackMCStudy::loadTPCOccMap(const o2::globaltracking::RecoContainer& recoDa
     int nTPCBins = NHBPerTF * o2::constants::lhc::LHCMaxBunches / 8, ninteg = 0;
     int nTPCOccBins = nTPCBins * mNTPCOccBinLengthInv, sumBins = std::max(1, int(o2::constants::lhc::LHCMaxBunches / 8 * mNTPCOccBinLengthInv));
     mTBinClOcc.resize(nTPCOccBins);
-    std::vector<float> mltHistTB(nTPCOccBins);
+    mTBinClOccHist.resize(nTPCOccBins);
     float sm = 0., tb = 0.5 * mNTPCOccBinLength;
     for (int i = 0; i < nTPCOccBins; i++) {
-      mltHistTB[i] = TPCRefitter->getParam()->GetUnscaledMult(tb);
+      mTBinClOccHist[i] = TPCRefitter->getParam()->GetUnscaledMult(tb);
       tb += mNTPCOccBinLength;
     }
     for (int i = nTPCOccBins; i--;) {
-      sm += mltHistTB[i];
+      sm += mTBinClOccHist[i];
       if (i + sumBins < nTPCOccBins) {
-        sm -= mltHistTB[i + sumBins];
+        sm -= mTBinClOccHist[i + sumBins];
       }
       mTBinClOcc[i] = sm;
     }
   } else {
     mTBinClOcc.resize(1);
+    mTBinClOccHist.resize(1);
   }
 }
 

Detectors/GlobalTrackingWorkflow/study/src/TrackingStudy.cxx

@@ -248,6 +248,28 @@ void TrackingStudySpec::process(o2::globaltracking::RecoContainer& recoData)
   auto vdrit = mTPCVDriftHelper.getVDriftObject().getVDrift();
   bool tpcTrackOK = recoData.isTrackSourceLoaded(GTrackID::TPC);
 
+  auto getTPCClInfo = [&recoData](const o2::tpc::TrackTPC& trc) {
+    const auto clRefs = recoData.getTPCTracksClusterRefs();
+    std::array<int, 3> clinfo = {};
+    if (recoData.inputsTPCclusters) {
+      uint8_t clSect = 0, clRow = 0, clRowP = -1;
+      uint32_t clIdx = 0;
+      for (int ic = 0; ic < trc.getNClusterReferences(); ic++) {
+        trc.getClusterReference(clRefs, ic, clSect, clRow, clIdx);
+        if (clRow != clRowP) {
+          clinfo[2]++;
+          clRowP = clRow;
+        }
+      }
+      const auto clRefs = recoData.getTPCTracksClusterRefs();
+      trc.getClusterReference(clRefs, trc.getNClusterReferences() - 1, clSect, clRow, clIdx);
+      clinfo[0] = clRow;
+      trc.getClusterReference(clRefs, 0, clSect, clRow, clIdx);
+      clinfo[1] = clRow;
+    }
+    return clinfo;
+  };
+
   for (int iv = 0; iv < nv; iv++) {
     LOGP(debug, "processing PV {} of {}", iv, nv);
     const auto& vtref = vtxRefs[iv];
@@ -297,13 +319,15 @@ void TrackingStudySpec::process(o2::globaltracking::RecoContainer& recoData)
         GTrackID tpcTrID;
         const o2::tpc::TrackTPC* tpcTr = nullptr;
         int nclTPC = 0;
+        std::array<int, 3> tpcClInfo{};
         if (dm[DetID::TPC] && tpcTrackOK) {
           tpcTrID = recoData.getTPCContributorGID(vid);
           tpcTr = &recoData.getTPCTrack(tpcTrID);
           nclTPC = tpcTr->getNClusters();
           if (nclTPC < mMinTPCClusters) {
             continue;
           }
+          tpcClInfo = getTPCClInfo(*tpcTr);
         }
         bool ambig = vid.isAmbiguous();
         auto trc = recoData.getTrackParam(vid);
@@ -368,6 +392,9 @@ void TrackingStudySpec::process(o2::globaltracking::RecoContainer& recoData)
           if (gidRefs[GTrackID::TPC].isIndexSet()) {
             trcExt.q2ptTPC = recoData.getTrackParam(gidRefs[GTrackID::TPC]).getQ2Pt();
             trcExt.nClTPC = nclTPC;
+            trcExt.rowMinTPC = tpcClInfo[0];
+            trcExt.rowMaxTPC = tpcClInfo[1];
+            trcExt.rowCountTPC = tpcClInfo[2];
           }
           if (gidRefs[GTrackID::ITSTPC].isIndexSet()) {
             const auto& trTPCITS = recoData.getTPCITSTrack(gidRefs[GTrackID::ITSTPC]);