Account for TPC clusters non-monotonous sorting

Author: shahor02

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

@@ -69,6 +69,10 @@ struct MCTrackInfo {
   void setBit(int bit) { flags |= BitMask & (0x1 << bit); }
   void resetBit(int bit) { flags &= ~(BitMask & (0x1 << bit)); }
 
+  o2::track::TrackPar getTrackParTPC(float b, float x = 90) const;
+  float getTrackParTPCPar(int i, float b, float x = 90) const;
+  float getTrackParTPCPhiSec(float b, float x = 90) const;
+
   ClassDefNV(MCTrackInfo, 7);
 };
 
@@ -80,28 +84,33 @@ struct RecTrack {
     FakeTOF = 0x1 << 3,
     FakeITSTPC = 0x1 << 4,
     FakeITSTPCTRD = 0x1 << 5,
+    HASACSides = 0x1 << 6,
     FakeGLO = 0x1 << 7
   };
   o2::track::TrackParCov track{};
   o2::dataformats::VtxTrackIndex gid{};
   o2::dataformats::TimeStampWithError<float, float> ts{};
   o2::MCEventLabel pvLabel{};
   short pvID = -1;
+  uint8_t nClTPCShared = 0;
   uint8_t flags = 0;
   uint8_t nClITS = 0;
   uint8_t nClTPC = 0;
   uint8_t pattITS = 0;
   int8_t lowestPadRow = -1;
   int8_t padFromEdge = -1;
+  uint8_t rowMaxTPC = 0;
+  uint8_t rowCountTPC = 0;
 
   bool isFakeGLO() const { return flags & FakeGLO; }
   bool isFakeITS() const { return flags & FakeITS; }
   bool isFakeTPC() const { return flags & FakeTPC; }
   bool isFakeTRD() const { return flags & FakeTRD; }
   bool isFakeTOF() const { return flags & FakeTOF; }
   bool isFakeITSTPC() const { return flags & FakeITSTPC; }
+  bool hasACSides() const { return flags & HASACSides; }
 
-  ClassDefNV(RecTrack, 2);
+  ClassDefNV(RecTrack, 3);
 };
 
 struct TrackPairInfo {
@@ -151,6 +160,13 @@ struct TrackFamily { // set of tracks related to the same MC label
   const RecTrack& getTrackWithTPC() const { return entTPC < 0 ? dummyRecTrack : recTracks[entTPC]; }
   const RecTrack& getTrackWithITSTPC() const { return entITSTPC < 0 ? dummyRecTrack : recTracks[entITSTPC]; }
   const RecTrack& getTrackWithITSFound() const { return entITSFound < 0 ? dummyRecTrack : recTracks[entITSFound]; }
+  const RecTrack& getLongestTPCTrack() const
+  {
+    int n = getLongestTPCTrackEntry();
+    return n < 0 ? dummyRecTrack : recTracks[n];
+  }
+  int getLongestTPCTrackEntry() const;
+  int getNTPCClones() const;
   static RecTrack dummyRecTrack; //
 
   ClassDefNV(TrackFamily, 1);

Detectors/GlobalTrackingWorkflow/study/src/TrackMCStudy.cxx

@@ -280,20 +280,33 @@ void TrackMCStudy::process(const o2::globaltracking::RecoContainer& recoData)
     return patt;
   };
 
-  auto getLowestPadrow = [&recoData](const o2::tpc::TrackTPC& trc, RecTrack& tref) {
+  auto fillTPCClusterInfo = [&recoData](const o2::tpc::TrackTPC& trc, RecTrack& tref) {
     if (recoData.inputsTPCclusters) {
-      uint8_t clSect = 0, clRow = 0;
+      uint8_t clSect = 0, clRow = 0, lowestR = -1;
       uint32_t clIdx = 0;
       const auto clRefs = recoData.getTPCTracksClusterRefs();
       const auto tpcClusAcc = recoData.getTPCClusters();
-      trc.getClusterReference(clRefs, trc.getNClusterReferences() - 1, clSect, clRow, clIdx);
+      const auto shMap = recoData.clusterShMapTPC;
+      for (int ic = 0; ic < trc.getNClusterReferences(); ic++) { // outside -> inside ordering, but on the sector boundaries backward jumps are possible
+        trc.getClusterReference(clRefs, ic, clSect, clRow, clIdx);
+        if (clRow < lowestR) {
+          tref.rowCountTPC++;
+          lowestR = clRow;
+        }
+        unsigned int absoluteIndex = tpcClusAcc.clusterOffset[clSect][clRow] + clIdx;
+        if (shMap[absoluteIndex] & o2::gpu::GPUTPCGMMergedTrackHit::flagShared) {
+          tref.nClTPCShared++;
+        }
+      }
+      tref.lowestPadRow = lowestR;
       const auto& clus = tpcClusAcc.clusters[clSect][clRow][clIdx];
       int padFromEdge = int(clus.getPad()), npads = o2::gpu::GPUTPCGeometry::NPads(clRow);
       if (padFromEdge > npads / 2) {
         padFromEdge = npads - 1 - padFromEdge;
       }
       tref.padFromEdge = uint8_t(padFromEdge);
-      tref.lowestPadRow = clRow;
+      trc.getClusterReference(clRefs, 0, clSect, clRow, clIdx);
+      tref.rowMaxTPC = clRow;
     }
   };
 
@@ -557,7 +570,10 @@ void TrackMCStudy::process(const o2::globaltracking::RecoContainer& recoData)
         if (msk[DetID::TPC]) {
           const auto& trtpc = recoData.getTPCTrack(gidSet[GTrackID::TPC]);
           tref.nClTPC = trtpc.getNClusters();
-          getLowestPadrow(trtpc, tref);
+          if (trtpc.hasBothSidesClusters()) {
+            tref.flags |= RecTrack::HASACSides;
+          }
+          fillTPCClusterInfo(trtpc, tref);
           flagTPCClusters(trtpc, entry.first);
           if (trackFam.entTPC < 0) {
             trackFam.entTPC = tcnt;
@@ -748,8 +764,8 @@ void TrackMCStudy::fillMCClusterInfo(const o2::globaltracking::RecoContainer& re
   const auto& params = o2::trackstudy::TrackMCStudyConfig::Instance();
 
   ClResTPC clRes{};
-  for (uint8_t sector = 0; sector < 36; sector++) {
-    for (uint8_t row = 0; row < 152; row++) {
+  for (uint8_t row = 0; row < 152; row++) { // we need to go in increasing row, so this should be the outer loop
+    for (uint8_t sector = 0; sector < 36; sector++) {
       unsigned int offs = TPCClusterIdxStruct.clusterOffset[sector][row];
       for (unsigned int icl0 = 0; icl0 < TPCClusterIdxStruct.nClusters[sector][row]; icl0++) {
         const auto labels = TPCClMClab->getLabels(icl0 + offs);

Detectors/GlobalTrackingWorkflow/study/src/TrackMCStudyTypes.cxx

@@ -77,4 +77,66 @@ int MCTrackInfo::getHighestITSLayer() const
   return -1;
 }
 
+o2::track::TrackPar MCTrackInfo::getTrackParTPC(float b, float x) const
+{
+  o2::track::TrackPar t(track);
+  int ntri = 0;
+  while (ntri < 2) {
+    int sector0 = o2::math_utils::angle2Sector(t.getAlpha());
+    if (!t.propagateParamTo(x, b)) {
+      t.invalidate();
+      break;
+    }
+    int sector = o2::math_utils::angle2Sector(t.getPhiPos());
+    float alpha = o2::math_utils::sector2Angle(sector);
+    if (!t.rotateParam(alpha)) {
+      t.invalidate();
+      break;
+    }
+    if (sector != sector0) {
+      ntri++;
+      continue;
+    }
+    break;
+  }
+  //  printf("%s ->\n%s <-\n",track.asString().c_str(), t.asString().c_str());
+  return t;
+}
+
+float MCTrackInfo::getTrackParTPCPar(int i, float b, float x) const
+{
+  auto t = getTrackParTPC(b, x);
+  return t.isValid() ? t.getParam(i) : -999.;
+}
+
+float MCTrackInfo::getTrackParTPCPhiSec(float b, float x) const
+{
+  auto t = getTrackParTPC(b, x);
+  return t.isValid() ? std::atan2(t.getY(), t.getX()) : -999.;
+}
+
+int TrackFamily::getLongestTPCTrackEntry() const
+{
+  int n = -1, ncl = 0;
+  int ntr = recTracks.size();
+  for (int i = 0; i < ntr; i++) {
+    if (recTracks[i].nClTPC > ncl) {
+      ncl = recTracks[i].nClTPC;
+      n = i;
+    }
+  }
+  return n;
+}
+
+int TrackFamily::getNTPCClones() const
+{
+  int n = 0;
+  for (auto& t : recTracks) {
+    if (t.nClTPC > 0) {
+      n++;
+    }
+  }
+  return n;
+}
+
 } // namespace o2::trackstudy

Detectors/GlobalTrackingWorkflow/study/src/TrackingStudy.cxx

@@ -274,25 +274,25 @@ void TrackingStudySpec::process(o2::globaltracking::RecoContainer& recoData)
     const auto clRefs = recoData.getTPCTracksClusterRefs();
     const auto tpcClusAcc = recoData.getTPCClusters();
     const auto shMap = recoData.clusterShMapTPC;
+
     if (recoData.inputsTPCclusters) {
-      uint8_t clSect = 0, clRow = 0, clRowP = -1;
+      uint8_t clSect = 0, clRow = 0, lowestR = -1;
       uint32_t clIdx = 0;
-      for (int ic = 0; ic < trc.getNClusterReferences(); ic++) {
+      for (int ic = 0; ic < trc.getNClusterReferences(); ic++) { // outside -> inside ordering, but on the sector boundaries backward jumps are possible
         trc.getClusterReference(clRefs, ic, clSect, clRow, clIdx);
-        if (clRow != clRowP) {
+        if (clRow < lowestR) {
           trExt.rowCountTPC++;
-          clRowP = clRow;
+          lowestR = clRow;
         }
         unsigned int absoluteIndex = tpcClusAcc.clusterOffset[clSect][clRow] + clIdx;
         if (shMap[absoluteIndex] & o2::gpu::GPUTPCGMMergedTrackHit::flagShared) {
           trExt.nClTPCShared++;
         }
       }
-      trc.getClusterReference(clRefs, trc.getNClusterReferences() - 1, clSect, clRow, clIdx);
-      trExt.rowMinTPC = clRow;
+      trExt.rowMinTPC = lowestR;
       const auto& clus = tpcClusAcc.clusters[clSect][clRow][clIdx];
       trExt.padFromEdge = uint8_t(clus.getPad());
-      int npads = o2::gpu::GPUTPCGeometry::NPads(clRow);
+      int npads = o2::gpu::GPUTPCGeometry::NPads(lowestR);
       if (trExt.padFromEdge > npads / 2) {
         trExt.padFromEdge = npads - 1 - trExt.padFromEdge;
       }
@@ -314,9 +314,9 @@ void TrackingStudySpec::process(o2::globaltracking::RecoContainer& recoData)
       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
+      for (int ic0 = 0; ic0 < trc0.getNClusterReferences(); ic0++) { // outside -> inside, but on the sector boundaries backward jumps are possible
         trc0.getClusterReference(clRefs, ic0, clSect0, clRow0, clIdx0);
-        for (int ic1 = ic1Start; ic1 < trc1.getNClusterReferences(); ic1++) { // outside -> inside
+        for (int ic1 = ic1Start; ic1 < trc1.getNClusterReferences(); ic1++) { // outside -> inside, but on the sector boundaries backward jumps are possible
           trc1.getClusterReference(clRefs, ic1, clSect1, clRow1, clIdx1);
           if (clRow1 > clRow0) {
             ic1Start = ic1 + 1;