Optionally fill dphi-dtgl with N Sh.clusters in trackStudy

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

@@ -86,6 +86,28 @@ struct RecTrack {
   ClassDefNV(RecTrack, 1);
 };
 
+struct TrackPairInfo {
+  RecTrack tr0;
+  RecTrack tr1;
+  uint8_t nshTPC = 0;
+  uint8_t nshTPCRow = 0;
+
+  int getComb() const { return tr0.track.getSign() != tr1.track.getSign() ? 0 : (tr0.track.getSign() > 0 ? 1 : 2); }
+  float getDPhi() const
+  {
+    float dphi = tr0.track.getPhi() - tr1.track.getPhi();
+    if (dphi < -o2::constants::math::PI) {
+      dphi += o2::constants::math::TwoPI;
+    } else if (dphi > o2::constants::math::PI) {
+      dphi -= o2::constants::math::TwoPI;
+    }
+    return dphi;
+  }
+  float getDTgl() const { return tr0.track.getTgl() - tr1.track.getTgl(); }
+
+  ClassDefNV(TrackPairInfo, 1)
+};
+
 struct TrackFamily { // set of tracks related to the same MC label
   MCTrackInfo mcTrackInfo{};
   std::vector<RecTrack> recTracks{};

Detectors/GlobalTrackingWorkflow/study/src/GlobalTrackingStudyLinkDef.h

@@ -36,5 +36,7 @@
 #pragma link C++ class o2::trackstudy::ClResTPC + ;
 #pragma link C++ class o2::trackstudy::ClResTPCCont + ;
 #pragma link C++ class std::vector < o2::trackstudy::ClResTPCCont> + ;
+#pragma link C++ class o2::trackstudy::TrackPairInfo + ;
+#pragma link C++ class std::vector < o2::trackstudy::TrackPairInfo> + ;
 
 #endif

Detectors/GlobalTrackingWorkflow/study/src/TrackingStudy.cxx

@@ -33,6 +33,7 @@
 #include "DetectorsBase/GRPGeomHelper.h"
 #include "GlobalTrackingStudy/TrackingStudy.h"
 #include "GlobalTrackingStudy/TrackInfoExt.h"
+#include "GlobalTrackingStudy/TrackMCStudyTypes.h"
 #include "TPCBase/ParameterElectronics.h"
 #include "ReconstructionDataFormats/PrimaryVertex.h"
 #include "ReconstructionDataFormats/PrimaryVertexExt.h"
@@ -108,6 +109,7 @@ class TrackingStudySpec : public Task
   int mNHBPerTF = 0;
   float mNTPCOccBinLengthInv;
   bool mStoreWithITSOnly = false;
+  bool mDoPairsCorr = false;
   std::string mDCAYFormula = "0.0105 + 0.0350 / pow(x, 1.1)";
   std::string mDCAZFormula = "0.0105 + 0.0350 / pow(x, 1.1)";
   GTrackID::mask_t mTracksSrc{};
@@ -136,6 +138,7 @@ void TrackingStudySpec::init(InitContext& ic)
   mMinTPCClusters = ic.options().get<int>("min-tpc-clusters");
   mDCAYFormula = ic.options().get<std::string>("dcay-vs-pt");
   mDCAZFormula = ic.options().get<std::string>("dcaz-vs-pt");
+  mDoPairsCorr = ic.options().get<bool>("pair-correlations");
 }
 
 void TrackingStudySpec::run(ProcessingContext& pc)
@@ -245,6 +248,7 @@ void TrackingStudySpec::process(o2::globaltracking::RecoContainer& recoData)
   float tBiasITS = alpParams.roFrameBiasInBC * o2::constants::lhc::LHCBunchSpacingMUS;
   const o2::ft0::InteractionTag& ft0Params = o2::ft0::InteractionTag::Instance();
   std::vector<o2::dataformats::TrackInfoExt> trcExtVec;
+  std::vector<o2::trackstudy::TrackPairInfo> trcPairsVec;
   auto vdrit = mTPCVDriftHelper.getVDriftObject().getVDrift();
   bool tpcTrackOK = recoData.isTrackSourceLoaded(GTrackID::TPC);
 
@@ -278,6 +282,82 @@ void TrackingStudySpec::process(o2::globaltracking::RecoContainer& recoData)
     }
   };
 
+  auto getTPCPairSharing = [&recoData, this](const o2::tpc::TrackTPC& trc0, const o2::tpc::TrackTPC& trc1) {
+    const auto clRefs = recoData.getTPCTracksClusterRefs();
+    const auto shMap = recoData.clusterShMapTPC.data();
+    uint8_t nsh = 0, nshRows = 0, lastSharedRow = -1;
+    if (recoData.inputsTPCclusters) {
+      uint8_t clSect0 = 0, clRow0 = 0, clSect1 = 0, clRow1 = 0;
+      uint32_t clIdx0 = 0, clIdx1 = 0;
+      int ic1Start = 0;
+      for (int ic0 = 0; ic0 < trc0.getNClusterReferences(); ic0++) { // outside -> inside
+        trc0.getClusterReference(clRefs, ic0, clSect0, clRow0, clIdx0);
+        for (int ic1 = ic1Start; ic1 < trc1.getNClusterReferences(); ic1++) { // outside -> inside
+          trc1.getClusterReference(clRefs, ic1, clSect1, clRow1, clIdx1);
+          if (clRow1 > clRow0) {
+            ic1Start = ic1 + 1;
+            continue; // catch up ic0
+          }
+          if (clRow1 == clRow0) {
+            if (clSect0 == clSect1 && clIdx0 == clIdx1) {
+              nsh++;
+              if (lastSharedRow != clRow0) {
+                lastSharedRow = clRow0;
+                nshRows++;
+              }
+              ic1Start = ic1 + 1;
+              break; // check next ic0
+            }
+          }
+        }
+      }
+    }
+    return std::make_pair(nsh, nshRows);
+  };
+
+  auto assignRecTrack = [&recoData, this](const o2::dataformats::TrackInfoExt& src, o2::trackstudy::RecTrack& dst) {
+    dst.track = src.track;
+    dst.gid = src.gid;
+    dst.ts.setTimeStamp(src.ttime);
+    dst.ts.setTimeStampError(src.ttimeE);
+    dst.nClITS = src.nClITS;
+    dst.nClTPC = src.nClTPC;
+    dst.pattITS = src.pattITS;
+    if (src.q2ptITS == 0. && dst.nClITS > 0) {
+      dst.pattITS |= 0x1 << 7;
+    }
+    dst.lowestPadRow = src.rowMinTPC;
+    if (this->mUseMC) {
+      auto gidSet = recoData.getSingleDetectorRefs(src.gid);
+      if (recoData.getTrackMCLabel(src.gid).isFake()) {
+        dst.flags |= RecTrack::FakeGLO;
+      }
+      auto msk = src.gid.getSourceDetectorsMask();
+      if (msk[DetID::ITS]) {
+        if (gidSet[GTrackID::ITS].isSourceSet()) { // has ITS track rather than AB tracklet
+          auto lblITS = recoData.getTrackMCLabel(gidSet[GTrackID::ITS]);
+          if (lblITS.isFake()) {
+            dst.flags |= RecTrack::FakeITS;
+          }
+        } else { // AB ITS tracklet
+          if (recoData.getTrackMCLabel(gidSet[GTrackID::ITSAB]).isFake()) {
+            dst.flags |= RecTrack::FakeITS;
+          }
+        }
+        if (msk[DetID::TPC]) { // has both ITS and TPC contribution
+          if (recoData.getTrackMCLabel(gidSet[GTrackID::ITSTPC]).isFake()) {
+            dst.flags |= RecTrack::FakeITSTPC;
+          }
+        }
+      }
+      if (msk[DetID::TPC]) {
+        if (recoData.getTrackMCLabel(gidSet[GTrackID::TPC]).isFake()) {
+          dst.flags |= RecTrack::FakeTPC;
+        }
+      }
+    }
+  };
+
   for (int iv = 0; iv < nv; iv++) {
     LOGP(debug, "processing PV {} of {}", iv, nv);
     const auto& vtref = vtxRefs[iv];
@@ -309,6 +389,7 @@ void TrackingStudySpec::process(o2::globaltracking::RecoContainer& recoData)
       pve.VtxID = iv;
     }
     trcExtVec.clear();
+    trcPairsVec.clear();
     float q2ptITS, q2ptTPC, q2ptITSTPC, q2ptITSTPCTRD;
     for (int is = 0; is < GTrackID::NSources; is++) {
       DetID::mask_t dm = GTrackID::getSourceDetectorsMask(is);
@@ -444,6 +525,42 @@ void TrackingStudySpec::process(o2::globaltracking::RecoContainer& recoData)
                << "orbit=" << recoData.startIR.orbit << "tfID=" << TFCount
                << "tpcOccBef=" << tpcOccBef << "tpcOccAft=" << tpcOccAft
                << "pve=" << pveVec[iv] << "trc=" << trcExtVec << "\n";
+
+    if (mDoPairsCorr) {
+      for (int it0 = 0; it0 < (int)trcExtVec.size(); it0++) {
+        const auto& tr0 = trcExtVec[it0];
+        if (tr0.nClTPC < 1) {
+          continue;
+        }
+        for (int it1 = it0 + 1; it1 < (int)trcExtVec.size(); it1++) {
+          const auto& tr1 = trcExtVec[it1];
+          if (tr1.nClTPC < 1) {
+            continue;
+          }
+
+          if (std::abs(tr0.track.getTgl() - tr1.track.getTgl()) > 0.25) {
+            continue;
+          }
+          auto dphi = tr0.track.getPhi() - tr1.track.getPhi();
+          if (dphi < -o2::constants::math::PI) {
+            dphi += o2::constants::math::TwoPI;
+          } else if (dphi > o2::constants::math::PI) {
+            dphi -= o2::constants::math::TwoPI;
+          }
+          if (std::abs(dphi) > 0.25) {
+            continue;
+          }
+          auto& pr = trcPairsVec.emplace_back();
+          assignRecTrack(tr0, pr.tr0);
+          assignRecTrack(tr1, pr.tr1);
+          auto shinfo = getTPCPairSharing(recoData.getTPCTrack(recoData.getTPCContributorGID(tr0.gid)), recoData.getTPCTrack(recoData.getTPCContributorGID(tr1.gid)));
+          pr.nshTPC = shinfo.first;
+          pr.nshTPCRow = shinfo.second;
+        }
+      }
+    }
+    (*mDBGOut) << "pairs"
+               << "pr=" << trcPairsVec << "\n";
   }
 
   int nvtot = mMaxNeighbours < 0 ? -1 : (int)pveVec.size();
@@ -600,6 +717,7 @@ DataProcessorSpec getTrackingStudySpec(GTrackID::mask_t srcTracks, GTrackID::mas
     {"max-eta", VariantType::Float, 1.0f, {"Cut on track eta"}},
     {"min-pt", VariantType::Float, 0.1f, {"Cut on track pT"}},
     {"with-its-only", VariantType::Bool, false, {"Store tracks with ITS only"}},
+    {"pair-correlations", VariantType::Bool, false, {"Do pairs correlation"}},
     {"min-x-prop", VariantType::Float, 100.f, {"track should be propagated to this X at least"}},
   };
   o2::tpc::VDriftHelper::requestCCDBInputs(dataRequest->inputs);