Break down pileup processing to 1 vertex per ROF

Author: mpuccio

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

@@ -86,6 +86,7 @@ struct TrackingParameters {
   float MaxChi2NDF = 30.f;
   bool UseTrackFollower = false;
   bool FindShortTracks = false;
+  bool PerPrimaryVertexProcessing = false;
 };
 
 inline int TrackingParameters::CellMinimumLevel()

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

@@ -78,13 +78,13 @@ class Tracker
 
  private:
   void initialiseTimeFrame(int& iteration);
-  void computeTracklets(int& iteration, int& iROFslice);
+  void computeTracklets(int& iteration, int& iROFslice, int& iVertex);
   void computeCells(int& iteration);
   void findCellsNeighbours(int& iteration);
   void findRoads(int& iteration);
 
   void initialiseTimeFrameHybrid(int& iteration);
-  void computeTrackletsHybrid(int& iteration, int& iROFslice);
+  void computeTrackletsHybrid(int& iteration, int& iROFslice, int& iVertex);
   void computeCellsHybrid(int& iteration);
   void findCellsNeighboursHybrid(int& iteration);
   void findRoadsHybrid(int& iteration);

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

@@ -52,7 +52,7 @@ class TrackerTraits
   virtual ~TrackerTraits() = default;
   virtual void adoptTimeFrame(TimeFrame* tf);
   virtual void initialiseTimeFrame(const int iteration);
-  virtual void computeLayerTracklets(const int iteration, int iROFslice);
+  virtual void computeLayerTracklets(const int iteration, int iROFslice, int iVertex);
   virtual void computeLayerCells(const int iteration);
   virtual void findCellsNeighbours(const int iteration);
   virtual void findRoads(const int iteration);

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

@@ -76,6 +76,7 @@ struct TrackerParamConfig : public o2::conf::ConfigurableParamHelper<TrackerPara
   int nThreads = 1;
   int nOrbitsPerIterations = 0;
   int nROFsPerIterations = 0;
+  bool perPrimaryVertexProcessing = false;
 
   O2ParamDef(TrackerParamConfig, "ITSCATrackerParam");
 };

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

@@ -49,39 +49,64 @@ void Tracker::clustersToTracks(std::function<void(std::string s)> logger, std::f
 {
   double total{0};
   mTraits->UpdateTrackingParameters(mTrkParams);
+  int maxNvertices{-1};
+  if (mTrkParams[0].PerPrimaryVertexProcessing) {
+    for (int iROF{0}; iROF < mTimeFrame->getNrof(); ++iROF) {
+      maxNvertices = std::max(maxNvertices, (int)mTimeFrame->getPrimaryVertices(iROF).size());
+    }
+  }
+
   for (int iteration = 0; iteration < (int)mTrkParams.size(); ++iteration) {
+
+    logger(fmt::format("ITS Tracking iteration {} summary:", iteration));
+    double timeTracklets{0.}, timeCells{0.}, timeNeighbours{0.}, timeRoads{0.};
+    int nTracklets{0}, nCells{0}, nNeighbours{0}, nTracks{-static_cast<int>(mTimeFrame->getNumberOfTracks())};
+
     total += evaluateTask(&Tracker::initialiseTimeFrame, "Timeframe initialisation", logger, iteration);
     int nROFsIterations = mTrkParams[iteration].nROFsPerIterations > 0 ? mTimeFrame->getNrof() / mTrkParams[iteration].nROFsPerIterations + bool(mTimeFrame->getNrof() % mTrkParams[iteration].nROFsPerIterations) : 1;
-    for (int iROFs{0}; iROFs < nROFsIterations; ++iROFs) {
-      total += evaluateTask(&Tracker::computeTracklets, "Tracklet finding", logger, iteration, iROFs);
-      logger(fmt::format("\t- Number of tracklets: {}", mTraits->getTFNumberOfTracklets()));
-      if (!mTimeFrame->checkMemory(mTrkParams[iteration].MaxMemory)) {
-        error(fmt::format("Too much memory used during trackleting in iteration {}, check the detector status and/or the selections.", iteration));
-        break;
-      }
-      float trackletsPerCluster = mTraits->getTFNumberOfClusters() > 0 ? float(mTraits->getTFNumberOfTracklets()) / mTraits->getTFNumberOfClusters() : 0.f;
-      if (trackletsPerCluster > mTrkParams[iteration].TrackletsPerClusterLimit) {
-        error(fmt::format("Too many tracklets per cluster ({}) in iteration {}, check the detector status and/or the selections. Current limit is {}", trackletsPerCluster, iteration, mTrkParams[iteration].TrackletsPerClusterLimit));
-        break;
-      }
+    int iVertex{std::min(maxNvertices, 0)};
+
+    do {
+      for (int iROFs{0}; iROFs < nROFsIterations; ++iROFs) {
+        timeTracklets += evaluateTask(
+          &Tracker::computeTracklets, "Tracklet finding", [](std::string) {}, iteration, iROFs, iVertex);
+        nTracklets += mTraits->getTFNumberOfTracklets();
+        if (!mTimeFrame->checkMemory(mTrkParams[iteration].MaxMemory)) {
+          error(fmt::format("Too much memory used during trackleting in iteration {}, check the detector status and/or the selections.", iteration));
+          break;
+        }
+        float trackletsPerCluster = mTraits->getTFNumberOfClusters() > 0 ? float(mTraits->getTFNumberOfTracklets()) / mTraits->getTFNumberOfClusters() : 0.f;
+        if (trackletsPerCluster > mTrkParams[iteration].TrackletsPerClusterLimit) {
+          error(fmt::format("Too many tracklets per cluster ({}) in iteration {}, check the detector status and/or the selections. Current limit is {}", trackletsPerCluster, iteration, mTrkParams[iteration].TrackletsPerClusterLimit));
+          break;
+        }
 
-      total += evaluateTask(&Tracker::computeCells, "Cell finding", logger, iteration);
-      logger(fmt::format("\t- Number of Cells: {}", mTraits->getTFNumberOfCells()));
-      if (!mTimeFrame->checkMemory(mTrkParams[iteration].MaxMemory)) {
-        error(fmt::format("Too much memory used during cell finding in iteration {}, check the detector status and/or the selections.", iteration));
-        break;
-      }
-      float cellsPerCluster = mTraits->getTFNumberOfClusters() > 0 ? float(mTraits->getTFNumberOfCells()) / mTraits->getTFNumberOfClusters() : 0.f;
-      if (cellsPerCluster > mTrkParams[iteration].CellsPerClusterLimit) {
-        error(fmt::format("Too many cells per cluster ({}) in iteration {}, check the detector status and/or the selections. Current limit is {}", cellsPerCluster, iteration, mTrkParams[iteration].CellsPerClusterLimit));
-        break;
-      }
+        timeCells += evaluateTask(
+          &Tracker::computeCells, "Cell finding", [](std::string) {}, iteration);
+        nCells += mTraits->getTFNumberOfCells();
+        if (!mTimeFrame->checkMemory(mTrkParams[iteration].MaxMemory)) {
+          error(fmt::format("Too much memory used during cell finding in iteration {}, check the detector status and/or the selections.", iteration));
+          break;
+        }
+        float cellsPerCluster = mTraits->getTFNumberOfClusters() > 0 ? float(mTraits->getTFNumberOfCells()) / mTraits->getTFNumberOfClusters() : 0.f;
+        if (cellsPerCluster > mTrkParams[iteration].CellsPerClusterLimit) {
+          error(fmt::format("Too many cells per cluster ({}) in iteration {}, check the detector status and/or the selections. Current limit is {}", cellsPerCluster, iteration, mTrkParams[iteration].CellsPerClusterLimit));
+          break;
+        }
 
-      total += evaluateTask(&Tracker::findCellsNeighbours, "Neighbour finding", logger, iteration);
-      logger(fmt::format("\t- Number of neighbours: {}", mTimeFrame->getNumberOfNeighbours()));
-      total += evaluateTask(&Tracker::findRoads, "Road finding", logger, iteration);
-      logger(fmt::format("\t- Number of Tracks: {}", mTimeFrame->getNumberOfTracks()));
-    }
+        timeNeighbours += evaluateTask(
+          &Tracker::findCellsNeighbours, "Neighbour finding", [](std::string) {}, iteration);
+        nNeighbours += mTimeFrame->getNumberOfNeighbours();
+        timeRoads += evaluateTask(
+          &Tracker::findRoads, "Road finding", [](std::string) {}, iteration);
+      }
+      iVertex++;
+    } while (iVertex < maxNvertices);
+    logger(fmt::format("\t- Tracklet finding: {} tracklets in {:.2f} ms", nTracklets, timeTracklets));
+    logger(fmt::format("\t- Cell finding: {} cells found in {:.2f} ms", nCells, timeCells));
+    logger(fmt::format("\t- Neighbours finding: {} neighbours found in {:.2f} ms", nNeighbours, timeNeighbours));
+    logger(fmt::format("\t- Track finding: {} tracks found in {:.2f} ms", nTracks + mTimeFrame->getNumberOfTracks(), timeRoads));
+    total += timeTracklets + timeCells + timeNeighbours + timeRoads;
     total += evaluateTask(&Tracker::extendTracks, "Extending tracks", logger, iteration);
   }
 
@@ -107,7 +132,7 @@ void Tracker::clustersToTracksHybrid(std::function<void(std::string s)> logger,
   mTraits->UpdateTrackingParameters(mTrkParams);
   for (int iteration = 0; iteration < (int)mTrkParams.size(); ++iteration) {
     total += evaluateTask(&Tracker::initialiseTimeFrameHybrid, "Hybrid Timeframe initialisation", logger, iteration);
-    total += evaluateTask(&Tracker::computeTrackletsHybrid, "Hybrid Tracklet finding", logger, iteration, iteration); // TODO: iteration argument put just for the sake of the interface, to be updated with the proper ROF slicing
+    total += evaluateTask(&Tracker::computeTrackletsHybrid, "Hybrid Tracklet finding", logger, iteration, iteration, iteration); // TODO: iteration argument put just for the sake of the interface, to be updated with the proper ROF slicing
     logger(fmt::format("\t- Number of tracklets: {}", mTraits->getTFNumberOfTracklets()));
     if (!mTimeFrame->checkMemory(mTrkParams[iteration].MaxMemory)) {
       error("Too much memory used during trackleting, check the detector status and/or the selections.");
@@ -143,9 +168,9 @@ void Tracker::initialiseTimeFrame(int& iteration)
   mTraits->initialiseTimeFrame(iteration);
 }
 
-void Tracker::computeTracklets(int& iteration, int& iROFslice)
+void Tracker::computeTracklets(int& iteration, int& iROFslice, int& iVertex)
 {
-  mTraits->computeLayerTracklets(iteration, iROFslice);
+  mTraits->computeLayerTracklets(iteration, iROFslice, iVertex);
 }
 
 void Tracker::computeCells(int& iteration)
@@ -168,7 +193,7 @@ void Tracker::initialiseTimeFrameHybrid(int& iteration)
   mTraits->initialiseTimeFrameHybrid(iteration);
 }
 
-void Tracker::computeTrackletsHybrid(int& iteration, int&)
+void Tracker::computeTrackletsHybrid(int& iteration, int&, int&)
 {
   mTraits->computeTrackletsHybrid(iteration);
 }
@@ -416,6 +441,7 @@ void Tracker::getGlobalConfiguration()
     params.CellDeltaTanLambdaSigma *= tc.deltaTanLres > 0 ? tc.deltaTanLres : 1.f;
     params.TrackletMinPt *= tc.minPt > 0 ? tc.minPt : 1.f;
     params.nROFsPerIterations = nROFsPerIterations;
+    params.PerPrimaryVertexProcessing = tc.perPrimaryVertexProcessing;
     for (int iD{0}; iD < 3; ++iD) {
       params.Diamond[iD] = tc.diamondPos[iD];
     }

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

@@ -51,7 +51,7 @@ namespace its
 
 constexpr int debugLevel{0};
 
-void TrackerTraits::computeLayerTracklets(const int iteration, int iROFslice)
+void TrackerTraits::computeLayerTracklets(const int iteration, int iROFslice, int iVertex)
 {
   TimeFrame* tf = mTimeFrame;
 
@@ -74,6 +74,8 @@ void TrackerTraits::computeLayerTracklets(const int iteration, int iROFslice)
   int endROF{mTrkParams[iteration].nROFsPerIterations > 0 ? std::min((iROFslice + 1) * mTrkParams[iteration].nROFsPerIterations + mTrkParams[iteration].DeltaROF, tf->getNrof()) : tf->getNrof()};
   for (int rof0{startROF}; rof0 < endROF; ++rof0) {
     gsl::span<const Vertex> primaryVertices = mTrkParams[iteration].UseDiamond ? diamondSpan : tf->getPrimaryVertices(rof0);
+    const int startVtx{iVertex >= 0 ? iVertex : 0};
+    const int endVtx{iVertex >= 0 ? std::min(iVertex + 1, static_cast<int>(primaryVertices.size())) : static_cast<int>(primaryVertices.size())};
     int minRof = std::max(startROF, rof0 - mTrkParams[iteration].DeltaROF);
     int maxRof = std::min(endROF - 1, rof0 + mTrkParams[iteration].DeltaROF);
 #pragma omp parallel for num_threads(mNThreads)
@@ -94,7 +96,8 @@ void TrackerTraits::computeLayerTracklets(const int iteration, int iROFslice)
         }
         const float inverseR0{1.f / currentCluster.radius};
 
-        for (auto& primaryVertex : primaryVertices) {
+        for (int iV{startVtx}; iV < endVtx; ++iV) {
+          auto& primaryVertex{primaryVertices[iV]};
           const float resolution = std::sqrt(Sq(mTrkParams[iteration].PVres) / primaryVertex.getNContributors() + Sq(tf->getPositionResolution(iLayer)));
 
           const float tanLambda{(currentCluster.zCoordinate - primaryVertex.getZ()) * inverseR0};
@@ -402,17 +405,16 @@ void TrackerTraits::findCellsNeighbours(const int iteration)
   std::ofstream off(fmt::format("cellneighs{}.txt", iteration));
 #endif
   for (int iLayer{0}; iLayer < mTrkParams[iteration].CellsPerRoad() - 1; ++iLayer) {
-
+    const int nextLayerCellsNum{static_cast<int>(mTimeFrame->getCells()[iLayer + 1].size())};
+    mTimeFrame->getCellsNeighboursLUT()[iLayer].clear();
+    mTimeFrame->getCellsNeighboursLUT()[iLayer].resize(nextLayerCellsNum, 0);
     if (mTimeFrame->getCells()[iLayer + 1].empty() ||
         mTimeFrame->getCellsLookupTable()[iLayer].empty()) {
+      mTimeFrame->getCellsNeighbours()[iLayer].clear();
       continue;
     }
 
     int layerCellsNum{static_cast<int>(mTimeFrame->getCells()[iLayer].size())};
-    const int nextLayerCellsNum{static_cast<int>(mTimeFrame->getCells()[iLayer + 1].size())};
-
-    mTimeFrame->getCellsNeighboursLUT()[iLayer].clear();
-    mTimeFrame->getCellsNeighboursLUT()[iLayer].resize(nextLayerCellsNum, 0);
     std::vector<std::pair<int, int>> cellsNeighbours;
     cellsNeighbours.reserve(nextLayerCellsNum);