Fix ITS/MFT clusterization for some complex shapes, O2-6424.

shahor02 · shahor02 · commit ecd2aca52cdd · 2026-01-14T23:43:39.000+01:00
diff --git a/Detectors/ITSMFT/common/reconstruction/include/ITSMFTReconstruction/Clusterer.h b/Detectors/ITSMFT/common/reconstruction/include/ITSMFTReconstruction/Clusterer.h
@@ -121,6 +121,10 @@ class Clusterer
   };
 
   struct ClustererThread {
+    struct PreCluster {
+      int head = 0; // index of precluster head in the pixels
+      int index = 0;
+    };
     int id = -1;
     Clusterer* parent = nullptr; // parent clusterer
     // buffers for entries in preClusterIndices in 2 columns, to avoid boundary checks, we reserve
@@ -132,12 +136,11 @@ class Clusterer
     // pixels[].first is the index of the next pixel of the same precluster in the pixels
     // pixels[].second is the index of the referred pixel in the ChipPixelData (element of mChips)
     std::vector<std::pair<int, uint32_t>> pixels;
-    std::vector<int> preClusterHeads; // index of precluster head in the pixels
-    std::vector<int> preClusterIndices;
     uint16_t currCol = 0xffff;               ///< Column being processed
     bool noLeftCol = true;                   ///< flag that there is no column on the left to check
     std::array<Label, MaxLabels> labelsBuff; //! temporary buffer for building cluster labels
     std::vector<PixelData> pixArrBuff;       //! temporary buffer for pattern calc.
+    std::vector<PreCluster> preClusters;     //! preclusters info
     //
     /// temporary storage for the thread output
     CompClusCont compClusters;
@@ -154,7 +157,7 @@ class Clusterer
     ///< add cluster at row (entry ip in the ChipPixeData) to the precluster with given index
     void expandPreCluster(uint32_t ip, uint16_t row, int preClusIndex)
     {
-      auto& firstIndex = preClusterHeads[preClusterIndices[preClusIndex]];
+      auto& firstIndex = preClusters[preClusters[preClusIndex].index].head;
       pixels.emplace_back(firstIndex, ip);
       firstIndex = pixels.size() - 1;
       curr[row] = preClusIndex;
@@ -163,11 +166,10 @@ class Clusterer
     ///< add new precluster at given row of current column for the fired pixel with index ip in the ChipPixelData
     void addNewPrecluster(uint32_t ip, uint16_t row)
     {
-      preClusterHeads.push_back(pixels.size());
+      int lastIndex = preClusters.size();
+      preClusters.emplace_back(pixels.size(), lastIndex);
       // new head does not point yet (-1) on other pixels, store just the entry of the pixel in the ChipPixelData
       pixels.emplace_back(-1, ip);
-      int lastIndex = preClusterIndices.size();
-      preClusterIndices.push_back(lastIndex);
       curr[row] = lastIndex; // store index of the new precluster in the current column buffer
     }
 
diff --git a/Detectors/ITSMFT/common/reconstruction/src/Clusterer.cxx b/Detectors/ITSMFT/common/reconstruction/src/Clusterer.cxx
@@ -133,15 +133,17 @@ void Clusterer::process(int nThreads, PixelReader& reader, CompClusCont* compClu
           if (stat.firstChip == chid) {
             thrStatIdx[ith]++;
             chid += stat.nChips; // next chip to look
-            const auto clbeg = mThreads[ith]->compClusters.begin() + stat.firstClus;
-            auto szold = compClus->size();
-            compClus->insert(compClus->end(), clbeg, clbeg + stat.nClus);
-            if (patterns) {
-              const auto ptbeg = mThreads[ith]->patterns.begin() + stat.firstPatt;
-              patterns->insert(patterns->end(), ptbeg, ptbeg + stat.nPatt);
-            }
-            if (labelsCl) {
-              labelsCl->mergeAtBack(mThreads[ith]->labels, stat.firstClus, stat.nClus);
+            if (stat.nClus > 0) {
+              const auto clbeg = mThreads[ith]->compClusters.begin() + stat.firstClus;
+              auto szold = compClus->size();
+              compClus->insert(compClus->end(), clbeg, clbeg + stat.nClus);
+              if (patterns) {
+                const auto ptbeg = mThreads[ith]->patterns.begin() + stat.firstPatt;
+                patterns->insert(patterns->end(), ptbeg, ptbeg + stat.nPatt);
+              }
+              if (labelsCl) {
+                labelsCl->mergeAtBack(mThreads[ith]->labels, stat.firstClus, stat.nClus);
+              }
             }
           }
         }
@@ -214,14 +216,22 @@ void Clusterer::ClustererThread::finishChip(ChipPixelData* curChipData, CompClus
                                             PatternCont* patternsPtr, const ConstMCTruth* labelsDigPtr, MCTruth* labelsClusPtr)
 {
   const auto& pixData = curChipData->getData();
-  for (int i1 = 0; i1 < preClusterHeads.size(); ++i1) {
-    auto ci = preClusterIndices[i1];
+  int nPreclusters = preClusters.size();
+  // account for the eventual reindexing of preClusters: Id2 might have been reindexed to Id1, which later was reindexed to Id0
+  for (int i = 1; i < nPreclusters; i++) {
+    if (preClusters[i].index != i) { // reindexing is always done towards smallest index
+      preClusters[i].index = preClusters[preClusters[i].index].index;
+    }
+  }
+  for (int i1 = 0; i1 < nPreclusters; ++i1) {
+    auto& preCluster = preClusters[i1];
+    auto ci = preCluster.index;
     if (ci < 0) {
       continue;
     }
     BBox bbox(curChipData->getChipID());
     int nlab = 0;
-    int next = preClusterHeads[i1];
+    int next = preCluster.head;
     pixArrBuff.clear();
     while (next >= 0) {
       const auto& pixEntry = pixels[next];
@@ -237,12 +247,13 @@ void Clusterer::ClustererThread::finishChip(ChipPixelData* curChipData, CompClus
       }
       next = pixEntry.first;
     }
-    preClusterIndices[i1] = -1;
-    for (int i2 = i1 + 1; i2 < preClusterHeads.size(); ++i2) {
-      if (preClusterIndices[i2] != ci) {
+    preCluster.index = -1;
+    for (int i2 = i1 + 1; i2 < nPreclusters; ++i2) {
+      auto& preCluster2 = preClusters[i2];
+      if (preCluster2.index != ci) {
         continue;
       }
-      next = preClusterHeads[i2];
+      next = preCluster2.head;
       while (next >= 0) {
         const auto& pixEntry = pixels[next];
         const auto pix = pixData[pixEntry.second]; // PixelData
@@ -257,7 +268,7 @@ void Clusterer::ClustererThread::finishChip(ChipPixelData* curChipData, CompClus
         }
         next = pixEntry.first;
       }
-      preClusterIndices[i2] = -1;
+      preCluster2.index = -1;
     }
     if (bbox.isAcceptableSize()) {
       parent->streamCluster(pixArrBuff, &labelsBuff, bbox, parent->mPattIdConverter, compClusPtr, patternsPtr, labelsClusPtr, nlab);
@@ -344,18 +355,15 @@ void Clusterer::ClustererThread::initChip(const ChipPixelData* curChipData, uint
   prev = column1 + 1;
   curr = column2 + 1;
   resetColumn(curr);
-
   pixels.clear();
-  preClusterHeads.clear();
-  preClusterIndices.clear();
+  preClusters.clear();
   auto pix = curChipData->getData()[first];
   currCol = pix.getCol();
   curr[pix.getRowDirect()] = 0; // can use getRowDirect since the pixel is not masked
   // start the first pre-cluster
-  preClusterHeads.push_back(0);
-  preClusterIndices.push_back(0);
+  preClusters.emplace_back();
   pixels.emplace_back(-1, first); // id of current pixel
-  noLeftCol = true;               // flag that there is no column on the left to check yet
+  noLeftCol = true;
 }
 
 //__________________________________________________
@@ -378,39 +386,58 @@ void Clusterer::ClustererThread::updateChip(const ChipPixelData* curChipData, ui
     currCol = pix.getCol();
   }
 
-  Bool_t orphan = true;
-
   if (noLeftCol) { // check only the row above
     if (curr[row - 1] >= 0) {
       expandPreCluster(ip, row, curr[row - 1]); // attach to the precluster of the previous row
-      return;
+    } else {
+      addNewPrecluster(ip, row); // start new precluster
     }
   } else {
+    // row above should be always checked
+    int nnb = 0, lowestIndex = curr[row - 1], lowestNb = 0, *nbrCol[4], nbrRow[4];
+    if (lowestIndex >= 0) {
+      nbrCol[nnb] = curr;
+      nbrRow[nnb++] = row - 1;
+    } else {
+      lowestIndex = 0x7ffff;
+      lowestNb = -1;
+    }
 #ifdef _ALLOW_DIAGONAL_ALPIDE_CLUSTERS_
-    int neighbours[]{curr[row - 1], prev[row], prev[row + 1], prev[row - 1]};
-#else
-    int neighbours[]{curr[row - 1], prev[row]};
-#endif
-    for (auto pci : neighbours) {
-      if (pci < 0) {
-        continue;
+    for (int i : {-1, 0, 1}) {
+      auto v = prev[row + i];
+      if (v >= 0) {
+        nbrCol[nnb] = prev;
+        nbrRow[nnb] = row + i;
+        if (v < lowestIndex) {
+          lowestIndex = v;
+          lowestNb = nnb;
+        }
+        nnb++;
       }
-      if (orphan) {
-        expandPreCluster(ip, row, pci); // attach to the adjascent precluster
-        orphan = false;
-        continue;
+    }
+#else
+    if (prev[row] >= 0) {
+      nbrCol[nnb] = prev;
+      nbrRow[nnb] = row;
+      if (prev[row] < lowestIndex) {
+        lowestIndex = v;
+        lowestNb = nnb;
       }
-      // reassign precluster index to smallest one
-      if (preClusterIndices[pci] < preClusterIndices[curr[row]]) {
-        preClusterIndices[curr[row]] = preClusterIndices[pci];
-      } else {
-        preClusterIndices[pci] = preClusterIndices[curr[row]];
+      nnb++;
+    }
+#endif
+    if (!nnb) {                  // no neighbours, create new precluster
+      addNewPrecluster(ip, row); // start new precluster
+    } else {
+      expandPreCluster(ip, row, lowestIndex); // attach to the adjascent precluster with smallest index
+      if (nnb > 1) {
+        for (int inb = 0; inb < nnb; inb++) { // reassign precluster index to smallest one, replicating updated values to columns caches
+          auto& prevIndex = (nbrCol[inb])[nbrRow[inb]];
+          prevIndex = preClusters[prevIndex].index = lowestIndex;
+        }
       }
     }
   }
-  if (orphan) {
-    addNewPrecluster(ip, row); // start new precluster
-  }
 }
 
 //__________________________________________________
