ITS: rework prepareClusters

Author: f3sch

Detectors/ITSMFT/ITS/tracking/src/TimeFrame.cxx

@@ -25,6 +25,7 @@
 #include "ITStracking/TrackingConfigParam.h"
 
 #include <iostream>
+#include <numeric>
 
 namespace
 {
@@ -214,34 +215,34 @@ int TimeFrame<nLayers>::loadROFrameData(gsl::span<o2::itsmft::ROFRecord> rofs,
 template <int nLayers>
 void TimeFrame<nLayers>::prepareClusters(const TrackingParameters& trkParam, const int maxLayers)
 {
+  const int numBins{trkParam.PhiBins * trkParam.ZBins};
+  const int stride{numBins + 1};
   bounded_vector<ClusterHelper> cHelper(mMemoryPool.get());
-  bounded_vector<int> clsPerBin(trkParam.PhiBins * trkParam.ZBins, 0, mMemoryPool.get());
+  bounded_vector<int> clsPerBin(numBins, 0, mMemoryPool.get());
+  bounded_vector<int> lutPerBin(numBins, 0, mMemoryPool.get());
   for (int rof{0}; rof < mNrof; ++rof) {
     if ((int)mMultiplicityCutMask.size() == mNrof && !mMultiplicityCutMask[rof]) {
       continue;
     }
-    for (int iLayer{0}; iLayer < std::min(trkParam.NLayers, maxLayers); ++iLayer) {
-      std::fill(clsPerBin.begin(), clsPerBin.end(), 0);
-      const auto unsortedClusters{getUnsortedClustersOnLayer(rof, iLayer)};
+    for (int iLayer{0}, stopLayer = std::min(trkParam.NLayers, maxLayers); iLayer < stopLayer; ++iLayer) {
+      const auto& unsortedClusters{getUnsortedClustersOnLayer(rof, iLayer)};
       const int clustersNum{static_cast<int>(unsortedClusters.size())};
+      auto* tableBase = mIndexTables[iLayer].data() + rof * stride;
 
-      deepVectorClear(cHelper);
       cHelper.resize(clustersNum);
 
       for (int iCluster{0}; iCluster < clustersNum; ++iCluster) {
-
         const Cluster& c = unsortedClusters[iCluster];
         ClusterHelper& h = cHelper[iCluster];
-        float x = c.xCoordinate - mBeamPos[0];
-        float y = c.yCoordinate - mBeamPos[1];
-        const float& z = c.zCoordinate;
+
+        const float x = c.xCoordinate - mBeamPos[0];
+        const float y = c.yCoordinate - mBeamPos[1];
+        const float z = c.zCoordinate;
+
         float phi = math_utils::computePhi(x, y);
         int zBin{mIndexTableUtils.getZBinIndex(iLayer, z)};
-        if (zBin < 0) {
-          zBin = 0;
-          mBogusClusters[iLayer]++;
-        } else if (zBin >= trkParam.ZBins) {
-          zBin = trkParam.ZBins - 1;
+        if (zBin < 0 || zBin >= trkParam.ZBins) {
+          zBin = std::clamp(zBin, 0, trkParam.ZBins - 1);
           mBogusClusters[iLayer]++;
         }
         int bin = mIndexTableUtils.getBinIndex(zBin, mIndexTableUtils.getPhiBinIndex(phi));
@@ -252,28 +253,23 @@ void TimeFrame<nLayers>::prepareClusters(const TrackingParameters& trkParam, con
         h.bin = bin;
         h.ind = clsPerBin[bin]++;
       }
-      bounded_vector<int> lutPerBin(clsPerBin.size(), 0, mMemoryPool.get());
-      lutPerBin[0] = 0;
-      for (unsigned int iB{1}; iB < lutPerBin.size(); ++iB) {
-        lutPerBin[iB] = lutPerBin[iB - 1] + clsPerBin[iB - 1];
-      }
+      std::exclusive_scan(clsPerBin.begin(), clsPerBin.end(), lutPerBin.begin(), 0);
 
       auto clusters2beSorted{getClustersOnLayer(rof, iLayer)};
       for (int iCluster{0}; iCluster < clustersNum; ++iCluster) {
         const ClusterHelper& h = cHelper[iCluster];
-
         Cluster& c = clusters2beSorted[lutPerBin[h.bin] + h.ind];
+
         c = unsortedClusters[iCluster];
         c.phi = h.phi;
         c.radius = h.r;
         c.indexTableBinIndex = h.bin;
       }
-      for (int iB{0}; iB < (int)clsPerBin.size(); ++iB) {
-        mIndexTables[iLayer][rof * (trkParam.ZBins * trkParam.PhiBins + 1) + iB] = lutPerBin[iB];
-      }
-      for (auto iB{clsPerBin.size()}; iB < (trkParam.ZBins * trkParam.PhiBins + 1); iB++) {
-        mIndexTables[iLayer][rof * (trkParam.ZBins * trkParam.PhiBins + 1) + iB] = clustersNum;
-      }
+      std::copy_n(lutPerBin.data(), clsPerBin.size(), tableBase);
+      std::fill_n(tableBase + clsPerBin.size(), stride - clsPerBin.size(), clustersNum);
+
+      std::fill(clsPerBin.begin(), clsPerBin.end(), 0);
+      cHelper.clear();
     }
   }
 }