asdlk

Signed-off-by: Felix Schlepper <felix.schlepper@cern.ch>

Author: f3sch

Detectors/ITSMFT/ITS/tracking/include/ITStracking/Configuration.h

@@ -67,6 +67,7 @@ struct TrackingParameters {
   int ReseedIfShorter = 6; // reseed for the final fit track with the length shorter than this
   std::vector<float> MinPt = {0.f, 0.f, 0.f, 0.f};
   uint16_t StartLayerMask = 0x7F;
+  int PNChunkSize = 0;
   bool RepeatRefitOut = false;   // repeat outward refit using inward refit as a seed
   bool ShiftRefToCluster = true; // TrackFit: after update shift the linearization reference to cluster
   bool PerPrimaryVertexProcessing = false;

Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackerTraits.h

@@ -53,7 +53,7 @@ class TrackerTraits
   virtual void computeLayerCells(const int iteration);
   virtual void findCellsNeighbours(const int iteration);
   virtual void findRoads(const int iteration);
-  virtual void processNeighbours(int iLayer, int iLevel, const bounded_vector<CellSeedN>& currentCellSeed, const bounded_vector<int>& currentCellId, bounded_vector<CellSeedN>& updatedCellSeed, bounded_vector<int>& updatedCellId);
+  virtual void processNeighbours(int iLayer, int iLevel, const bounded_vector<CellSeedN>& currentCellSeed, const bounded_vector<int>& currentCellId, bounded_vector<CellSeedN>& updatedCellSeed, bounded_vector<int>& updatedCellId, int beginCell = 0, int endCell = -1);
 
   void updateTrackingParameters(const std::vector<TrackingParameters>& trkPars) { mTrkParams = trkPars; }
   TimeFrame<NLayers>* getTimeFrame() { return mTimeFrame; }

Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackingConfigParam.h

@@ -36,9 +36,9 @@ struct VertexerParamConfig : public o2::conf::ConfigurableParamHelper<VertexerPa
   float duplicateZCut = 0.7985643f;
   float finalSelectionZCut = 0.2932624f;
   float duplicateDistance2Cut = 0.0223001f;
-  float tanLambdaCut = 0.002f;        // tanLambda = deltaZ/deltaR
+  float tanLambdaCut = 0.002f; // tanLambda = deltaZ/deltaR
   float nSigmaCut = 0.0479011f;
-  float maxZPositionAllowed = 25.f;   // 4x sZ of the beam
+  float maxZPositionAllowed = 25.f; // 4x sZ of the beam
 
   // Artefacts selections
   int clusterContributorsCut = 3; // minimum number of contributors for an accepted final vertex
@@ -71,6 +71,7 @@ struct TrackerParamConfig : public o2::conf::ConfigurableParamHelper<TrackerPara
   int addTimeError[7] = {0};                                                // configure the width of the window in BC to be considered for the tracking.
   int minTrackLgtIter[MaxIter] = {};                                        // minimum track length at each iteration, used only if >0, otherwise use code defaults
   uint8_t startLayerMask[MaxIter] = {};                                     // mask of start layer for this iteration (if >0)
+  int pnChunkSize[MaxIter] = {};                                            // chunk level in input cells for processNeighbours at each iteration (if >0)
   float minPtIterLgt[MaxIter * (MaxTrackLength - MinTrackLength + 1)] = {}; // min.pT for given track length at this iteration, used only if >0, otherwise use code defaults
   float sysErrY2[7] = {0};                                                  // systematic error^2 in Y per layer
   float sysErrZ2[7] = {0};                                                  // systematic error^2 in Z per layer

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

@@ -26,6 +26,9 @@ std::string TrackingParameters::asString() const
 {
   std::string str = std::format("NZb:{} NPhB:{} PerVtx:{} DropFail:{} ClSh:{} TtklMinPt:{:.2f} MinCl:{}",
                                 ZBins, PhiBins, PerPrimaryVertexProcessing, DropTFUponFailure, ClusterSharing, TrackletMinPt, MinTrackLength);
+  if (PNChunkSize > 0) {
+    str += std::format(" PNCh:{}", PNChunkSize);
+  }
   bool first = true;
   for (int il = NLayers; il >= MinTrackLength; il--) {
     int slot = NLayers - il;
@@ -145,6 +148,7 @@ std::vector<TrackingParameters> TrackingMode::getTrackingParameters(TrackingMode
         if (tc.startLayerMask[ip] > 0) {
           trackParams[2].StartLayerMask = tc.startLayerMask[ip];
         }
+        param.PNChunkSize = tc.pnChunkSize[ip] > 0 ? 1 << tc.pnChunkSize[ip] : 0;
         if (tc.minTrackLgtIter[ip] > 0) {
           param.MinTrackLength = tc.minTrackLgtIter[ip];
         }

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

@@ -514,12 +514,22 @@ void TrackerTraits<NLayers>::findCellsNeighbours(const int iteration)
 }
 
 template <int NLayers>
-void TrackerTraits<NLayers>::processNeighbours(int iLayer, int iLevel, const bounded_vector<CellSeedN>& currentCellSeed, const bounded_vector<int>& currentCellId, bounded_vector<CellSeedN>& updatedCellSeeds, bounded_vector<int>& updatedCellsIds)
+void TrackerTraits<NLayers>::processNeighbours(int iLayer, int iLevel, const bounded_vector<CellSeedN>& currentCellSeed, const bounded_vector<int>& currentCellId, bounded_vector<CellSeedN>& updatedCellSeeds, bounded_vector<int>& updatedCellsIds, int beginCell, int endCell)
 {
-  auto propagator = o2::base::Propagator::Instance();
+  updatedCellSeeds.clear();
+  updatedCellsIds.clear();
+
+  const int totalCells = static_cast<int>(currentCellSeed.size());
+  beginCell = std::clamp(beginCell, 0, totalCells);
+  endCell = endCell < 0 ? totalCells : std::clamp(endCell, beginCell, totalCells);
+  if (beginCell >= endCell) {
+    return;
+  }
 
+  auto propagator = o2::base::Propagator::Instance();
   mTaskArena->execute([&] {
-    auto forCellNeighbours = [&](auto Tag, int iCell, int offset = 0) -> int {
+    auto forCellNeighbours = [&](auto Tag, int localCell, int offset = 0) -> int {
+      const int iCell = beginCell + localCell;
       const auto& currentCell{currentCellSeed[iCell]};
 
       if constexpr (decltype(Tag)::value != PassMode::TwoPassInsert::value) {
@@ -604,15 +614,15 @@ void TrackerTraits<NLayers>::processNeighbours(int iLayer, int iLevel, const bou
       return foundSeeds;
     };
 
-    const int nCells = static_cast<int>(currentCellSeed.size());
+    const int nCells = endCell - beginCell;
     if (mTaskArena->max_concurrency() <= 1) {
-      for (int iCell{0}; iCell < nCells; ++iCell) {
-        forCellNeighbours(PassMode::OnePass{}, iCell);
+      for (int localCell{0}; localCell < nCells; ++localCell) {
+        forCellNeighbours(PassMode::OnePass{}, localCell);
       }
     } else {
       bounded_vector<int> perCellCount(nCells + 1, 0, mMemoryPool.get());
-      tbb::parallel_for(0, nCells, [&](const int iCell) {
-        perCellCount[iCell] = forCellNeighbours(PassMode::TwoPassCount{}, iCell);
+      tbb::parallel_for(0, nCells, [&](const int localCell) {
+        perCellCount[localCell] = forCellNeighbours(PassMode::TwoPassCount{}, localCell);
       });
 
       std::exclusive_scan(perCellCount.begin(), perCellCount.end(), perCellCount.begin(), 0);
@@ -623,12 +633,12 @@ void TrackerTraits<NLayers>::processNeighbours(int iLayer, int iLevel, const bou
       updatedCellSeeds.resize(totalNeighbours);
       updatedCellsIds.resize(totalNeighbours);
 
-      tbb::parallel_for(0, nCells, [&](const int iCell) {
-        int offset = perCellCount[iCell];
-        if (offset == perCellCount[iCell + 1]) {
+      tbb::parallel_for(0, nCells, [&](const int localCell) {
+        int offset = perCellCount[localCell];
+        if (offset == perCellCount[localCell + 1]) {
           return;
         }
-        forCellNeighbours(PassMode::TwoPassInsert{}, iCell, offset);
+        forCellNeighbours(PassMode::TwoPassInsert{}, localCell, offset);
       });
     }
   });
@@ -644,119 +654,157 @@ void TrackerTraits<NLayers>::findRoads(const int iteration)
   for (int startLevel{mTrkParams[iteration].CellsPerRoad()}; startLevel >= mTrkParams[iteration].CellMinimumLevel(); --startLevel) {
 
     auto seedFilter = [&](const auto& seed) {
-      return seed.getQ2Pt() <= 1.e3 && seed.getChi2() <= mTrkParams[0].MaxChi2NDF * ((startLevel + 2) * 2 - 5);
+      return seed.getQ2Pt() <= 1.e3 && seed.getChi2() <= mTrkParams[iteration].MaxChi2NDF * ((startLevel + 2) * 2 - 5);
     };
 
-    bounded_vector<CellSeedN> trackSeeds(mMemoryPool.get());
-    for (int startLayer{mTrkParams[iteration].NeighboursPerRoad()}; startLayer >= startLevel - 1; --startLayer) {
-      if ((mTrkParams[iteration].StartLayerMask & (1 << (startLayer + 2))) == 0) {
-        continue;
-      }
-
-      bounded_vector<int> lastCellId(mMemoryPool.get()), updatedCellId(mMemoryPool.get());
-      bounded_vector<CellSeedN> lastCellSeed(mMemoryPool.get()), updatedCellSeed(mMemoryPool.get());
-
-      processNeighbours(startLayer, startLevel, mTimeFrame->getCells()[startLayer], lastCellId, updatedCellSeed, updatedCellId);
-
-      int level = startLevel;
-      for (int iLayer{startLayer - 1}; iLayer > 0 && level > 2; --iLayer) {
-        lastCellSeed.swap(updatedCellSeed);
-        lastCellId.swap(updatedCellId);
-        deepVectorClear(updatedCellSeed); /// tame the memory peaks
-        deepVectorClear(updatedCellId);   /// tame the memory peaks
-        processNeighbours(iLayer, --level, lastCellSeed, lastCellId, updatedCellSeed, updatedCellId);
-      }
-      deepVectorClear(lastCellId);   /// tame the memory peaks
-      deepVectorClear(lastCellSeed); /// tame the memory peaks
-
-      if (!updatedCellSeed.empty()) {
-        trackSeeds.reserve(trackSeeds.size() + std::count_if(updatedCellSeed.begin(), updatedCellSeed.end(), seedFilter));
-        std::copy_if(updatedCellSeed.begin(), updatedCellSeed.end(), std::back_inserter(trackSeeds), seedFilter);
-      }
-    }
-
-    if (trackSeeds.empty()) {
-      continue;
-    }
-
     bounded_vector<TrackITSExt> tracks(mMemoryPool.get());
-    mTaskArena->execute([&] {
-      auto forSeed = [&](auto Tag, int iSeed, int offset = 0) {
-        TrackITSExt temporaryTrack = seedTrackForRefit(trackSeeds[iSeed]);
-        o2::track::TrackPar linRef{temporaryTrack};
-        bool fitSuccess = fitTrack(temporaryTrack, 0, mTrkParams[0].NLayers, 1, mTrkParams[0].MaxChi2ClusterAttachment, mTrkParams[0].MaxChi2NDF, o2::constants::math::VeryBig, 0, &linRef);
-        if (!fitSuccess) {
-          return 0;
-        }
-        temporaryTrack.getParamOut() = temporaryTrack.getParamIn();
-        linRef = temporaryTrack.getParamOut(); // use refitted track as lin.reference
-        temporaryTrack.resetCovariance();
-        temporaryTrack.setCov(temporaryTrack.getQ2Pt() * temporaryTrack.getQ2Pt() * temporaryTrack.getCov()[o2::track::CovLabels::kSigQ2Pt2], o2::track::CovLabels::kSigQ2Pt2);
-        temporaryTrack.setChi2(0);
-        fitSuccess = fitTrack(temporaryTrack, mTrkParams[0].NLayers - 1, -1, -1, mTrkParams[0].MaxChi2ClusterAttachment, mTrkParams[0].MaxChi2NDF, 50.f, 0, &linRef);
-        if (!fitSuccess || temporaryTrack.getPt() < mTrkParams[iteration].MinPt[mTrkParams[iteration].NLayers - temporaryTrack.getNClusters()]) {
-          return 0;
-        }
-        if (mTrkParams[0].RepeatRefitOut) { // repeat outward refit seeding and linearizing with the stable inward fit result
-          o2::track::TrackParCov saveInw{temporaryTrack};
-          linRef = saveInw; // use refitted track as lin.reference
-          float saveChi2 = temporaryTrack.getChi2();
+    auto fitSeeds = [&](const bounded_vector<CellSeedN>& finalSeeds) {
+      bounded_vector<TrackITSExt> trackChunk(mMemoryPool.get());
+      mTaskArena->execute([&] {
+        auto forSeed = [&](auto Tag, int iSeed, int offset = 0) {
+          if (!seedFilter(finalSeeds[iSeed])) {
+            return 0;
+          }
+
+          TrackITSExt temporaryTrack = seedTrackForRefit(finalSeeds[iSeed]);
+          o2::track::TrackPar linRef{temporaryTrack};
+          bool fitSuccess = fitTrack(temporaryTrack, 0, mTrkParams[0].NLayers, 1, mTrkParams[0].MaxChi2ClusterAttachment, mTrkParams[0].MaxChi2NDF, o2::constants::math::VeryBig, 0, &linRef);
+          if (!fitSuccess) {
+            return 0;
+          }
+          temporaryTrack.getParamOut() = temporaryTrack.getParamIn();
+          linRef = temporaryTrack.getParamOut(); // use refitted track as lin.reference
           temporaryTrack.resetCovariance();
           temporaryTrack.setCov(temporaryTrack.getQ2Pt() * temporaryTrack.getQ2Pt() * temporaryTrack.getCov()[o2::track::CovLabels::kSigQ2Pt2], o2::track::CovLabels::kSigQ2Pt2);
           temporaryTrack.setChi2(0);
-          fitSuccess = fitTrack(temporaryTrack, 0, mTrkParams[0].NLayers, 1, mTrkParams[0].MaxChi2ClusterAttachment, mTrkParams[0].MaxChi2NDF, o2::constants::math::VeryBig, 0, &linRef);
-          if (!fitSuccess) {
+          fitSuccess = fitTrack(temporaryTrack, mTrkParams[0].NLayers - 1, -1, -1, mTrkParams[0].MaxChi2ClusterAttachment, mTrkParams[0].MaxChi2NDF, 50.f, 0, &linRef);
+          if (!fitSuccess || temporaryTrack.getPt() < mTrkParams[iteration].MinPt[mTrkParams[iteration].NLayers - temporaryTrack.getNClusters()]) {
             return 0;
           }
-          temporaryTrack.getParamOut() = temporaryTrack.getParamIn();
-          temporaryTrack.getParamIn() = saveInw;
-          temporaryTrack.setChi2(saveChi2);
-        }
+          if (mTrkParams[0].RepeatRefitOut) { // repeat outward refit seeding and linearizing with the stable inward fit result
+            o2::track::TrackParCov saveInw{temporaryTrack};
+            linRef = saveInw; // use refitted track as lin.reference
+            float saveChi2 = temporaryTrack.getChi2();
+            temporaryTrack.resetCovariance();
+            temporaryTrack.setCov(temporaryTrack.getQ2Pt() * temporaryTrack.getQ2Pt() * temporaryTrack.getCov()[o2::track::CovLabels::kSigQ2Pt2], o2::track::CovLabels::kSigQ2Pt2);
+            temporaryTrack.setChi2(0);
+            fitSuccess = fitTrack(temporaryTrack, 0, mTrkParams[0].NLayers, 1, mTrkParams[0].MaxChi2ClusterAttachment, mTrkParams[0].MaxChi2NDF, o2::constants::math::VeryBig, 0, &linRef);
+            if (!fitSuccess) {
+              return 0;
+            }
+            temporaryTrack.getParamOut() = temporaryTrack.getParamIn();
+            temporaryTrack.getParamIn() = saveInw;
+            temporaryTrack.setChi2(saveChi2);
+          }
 
-        if constexpr (decltype(Tag)::value == PassMode::OnePass::value) {
-          tracks.push_back(temporaryTrack);
-        } else if constexpr (decltype(Tag)::value == PassMode::TwoPassCount::value) {
-          // nothing to do
-        } else if constexpr (decltype(Tag)::value == PassMode::TwoPassInsert::value) {
-          tracks[offset] = temporaryTrack;
+          if constexpr (decltype(Tag)::value == PassMode::OnePass::value) {
+            trackChunk.push_back(temporaryTrack);
+          } else if constexpr (decltype(Tag)::value == PassMode::TwoPassCount::value) {
+            // nothing to do
+          } else if constexpr (decltype(Tag)::value == PassMode::TwoPassInsert::value) {
+            trackChunk[offset] = temporaryTrack;
+          } else {
+            static_assert(false, "Unknown mode!");
+          }
+          return 1;
+        };
+
+        const int nSeeds = static_cast<int>(finalSeeds.size());
+        if (mTaskArena->max_concurrency() <= 1) {
+          for (int iSeed{0}; iSeed < nSeeds; ++iSeed) {
+            forSeed(PassMode::OnePass{}, iSeed);
+          }
         } else {
-          static_assert(false, "Unknown mode!");
-        }
-        return 1;
-      };
+          bounded_vector<int> perSeedCount(nSeeds + 1, 0, mMemoryPool.get());
+          tbb::parallel_for(0, nSeeds, [&](const int iSeed) {
+            perSeedCount[iSeed] = forSeed(PassMode::TwoPassCount{}, iSeed);
+          });
+
+          std::exclusive_scan(perSeedCount.begin(), perSeedCount.end(), perSeedCount.begin(), 0);
+          auto totalTracks{perSeedCount.back()};
+          if (totalTracks == 0) {
+            return;
+          }
+          trackChunk.resize(totalTracks);
 
-      const int nSeeds = static_cast<int>(trackSeeds.size());
-      if (mTaskArena->max_concurrency() <= 1) {
-        for (int iSeed{0}; iSeed < nSeeds; ++iSeed) {
-          forSeed(PassMode::OnePass{}, iSeed);
+          tbb::parallel_for(0, nSeeds, [&](const int iSeed) {
+            if (perSeedCount[iSeed] == perSeedCount[iSeed + 1]) {
+              return;
+            }
+            forSeed(PassMode::TwoPassInsert{}, iSeed, perSeedCount[iSeed]);
+          });
         }
-      } else {
-        bounded_vector<int> perSeedCount(nSeeds + 1, 0, mMemoryPool.get());
-        tbb::parallel_for(0, nSeeds, [&](const int iSeed) {
-          perSeedCount[iSeed] = forSeed(PassMode::TwoPassCount{}, iSeed);
-        });
+      });
+      if (!trackChunk.empty()) {
+        tracks.reserve(tracks.size() + trackChunk.size());
+        tracks.insert(tracks.end(), std::make_move_iterator(trackChunk.begin()), std::make_move_iterator(trackChunk.end()));
+      }
+    };
 
-        std::exclusive_scan(perSeedCount.begin(), perSeedCount.end(), perSeedCount.begin(), 0);
-        auto totalTracks{perSeedCount.back()};
-        if (totalTracks == 0) {
-          return;
-        }
-        tracks.resize(totalTracks);
+    for (int startLayer{mTrkParams[iteration].NeighboursPerRoad()}; startLayer >= startLevel - 1; --startLayer) {
+      if ((mTrkParams[iteration].StartLayerMask & (1 << (startLayer + 2))) == 0) {
+        continue;
+      }
 
-        tbb::parallel_for(0, nSeeds, [&](const int iSeed) {
-          if (perSeedCount[iSeed] == perSeedCount[iSeed + 1]) {
-            return;
+      bounded_vector<int> inputCellIds(mMemoryPool.get()), outputCellIds(mMemoryPool.get());
+      bounded_vector<CellSeedN> inputCellSeeds(mMemoryPool.get()), outputCellSeeds(mMemoryPool.get());
+      const bounded_vector<CellSeedN>* currentSeeds = &mTimeFrame->getCells()[startLayer];
+      const bounded_vector<int>* currentIds = &inputCellIds;
+      int layer = startLayer;
+      int level = startLevel;
+
+      while (!currentSeeds->empty() && level > 1) {
+        if (!mTrkParams[iteration].PNChunkSize) {
+          processNeighbours(layer, level, *currentSeeds, *currentIds, outputCellSeeds, outputCellIds);
+          if (level == 2) {
+            fitSeeds(outputCellSeeds);
           }
-          forSeed(PassMode::TwoPassInsert{}, iSeed, perSeedCount[iSeed]);
-        });
+        } else {
+          for (int beginCell{0}; beginCell < static_cast<int>(currentSeeds->size()); beginCell += mTrkParams[iteration].PNChunkSize) {
+            const int endCell = std::min(beginCell + mTrkParams[iteration].PNChunkSize, static_cast<int>(currentSeeds->size()));
+            bounded_vector<CellSeedN> chunkSeeds(mMemoryPool.get());
+            bounded_vector<int> chunkIds(mMemoryPool.get());
+            processNeighbours(layer, level, *currentSeeds, *currentIds, chunkSeeds, chunkIds, beginCell, endCell);
+            if (chunkSeeds.empty()) {
+              continue;
+            }
+            if (level == 2) {
+              fitSeeds(chunkSeeds);
+              continue;
+            }
+            outputCellSeeds.reserve(outputCellSeeds.size() + chunkSeeds.size());
+            outputCellIds.reserve(outputCellIds.size() + chunkIds.size());
+            outputCellSeeds.insert(outputCellSeeds.end(), std::make_move_iterator(chunkSeeds.begin()), std::make_move_iterator(chunkSeeds.end()));
+            outputCellIds.insert(outputCellIds.end(), std::make_move_iterator(chunkIds.begin()), std::make_move_iterator(chunkIds.end()));
+          }
+        }
+        if (outputCellSeeds.empty() || level == 2) {
+          break;
+        }
+
+        --layer;
+        --level;
+        inputCellSeeds.swap(outputCellSeeds);
+        inputCellIds.swap(outputCellIds);
+        deepVectorClear(outputCellSeeds);
+        deepVectorClear(outputCellIds);
+        currentSeeds = &inputCellSeeds;
+        currentIds = &inputCellIds;
       }
+    }
 
-      deepVectorClear(trackSeeds);
+    if (tracks.empty()) {
+      continue;
+    }
+
+    // sort by quality do to disambiguation
+    mTaskArena->execute([&] {
       tbb::parallel_sort(tracks.begin(), tracks.end(), [](const auto& a, const auto& b) {
         return a.getChi2() < b.getChi2();
       });
     });
 
+    // create final tracks
     const float smallestROFHalf = mTimeFrame->getROFOverlapTableView().getClockLayer().mROFLength * 0.5f;
     for (auto& track : tracks) {
       int nShared = 0;