Tracker.cxx

Go to the documentation of this file.

// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
// All rights not expressly granted are reserved.
//
// This software is distributed under the terms of the GNU General Public
// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
//
// In applying this license CERN does not waive the privileges and immunities
// granted to it by virtue of its status as an Intergovernmental Organization
// or submit itself to any jurisdiction.
 
#include "ITStracking/Tracker.h"
 
#include "ITStracking/Cell.h"
#include "ITStracking/Constants.h"
#include "ITStracking/IndexTableUtils.h"
#include "ITStracking/Smoother.h"
#include "ITStracking/Tracklet.h"
#include "ITStracking/TrackerTraits.h"
#include "ITStracking/TrackingConfigParam.h"
 
#include "ReconstructionDataFormats/Track.h"
#include <cassert>
#include <format>
#include <cstdlib>
#include <string>
#include <climits>
 
namespace o2
{
namespace its
{
using o2::its::constants::GB;
 
Tracker::Tracker(o2::its::TrackerTraits* traits) : mTraits(traits)
{
  mTrkParams.resize(1);
}
 
void Tracker::clustersToTracks(LogFunc logger, LogFunc error)
{
  LogFunc evalLog = [](const std::string&) {};
 
  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());
    }
  }
 
  bool dropTF = false;
  for (int iteration = 0; iteration < (int)mTrkParams.size(); ++iteration) {
    if (iteration == 3 && mTrkParams[0].DoUPCIteration) {
      mTimeFrame->swapMasks();
    }
    double timeTracklets{0.}, timeCells{0.}, timeNeighbours{0.}, timeRoads{0.};
    int nTracklets{0}, nCells{0}, nNeighbours{0}, nTracks{-static_cast<int>(mTimeFrame->getNumberOfTracks())};
    int nROFsIterations = mTrkParams[iteration].nROFsPerIterations > 0 ? mTimeFrame->getNrof() / mTrkParams[iteration].nROFsPerIterations + bool(mTimeFrame->getNrof() % mTrkParams[iteration].nROFsPerIterations) : 1;
    int iVertex{std::min(maxNvertices, 0)};
    logger(std::format("==== ITS {} Tracking iteration {} summary ====", mTraits->getName(), iteration));
 
    total += evaluateTask(&Tracker::initialiseTimeFrame, "Timeframe initialisation", logger, iteration);
    do {
      for (int iROFs{0}; iROFs < nROFsIterations; ++iROFs) {
        timeTracklets += evaluateTask(&Tracker::computeTracklets, "Tracklet finding", evalLog, iteration, iROFs, iVertex);
        nTracklets += mTraits->getTFNumberOfTracklets();
        if (!mTimeFrame->checkMemory(mTrkParams[iteration].MaxMemory)) {
          mTimeFrame->printSliceInfo(iROFs, mTrkParams[iteration].nROFsPerIterations);
          error(std::format("Too much memory used during trackleting in iteration {} in ROF span {}-{}: {:.2f} GB. Current limit is {:.2f} GB, check the detector status and/or the selections.",
                            iteration, iROFs, iROFs + mTrkParams[iteration].nROFsPerIterations, mTimeFrame->getArtefactsMemory() / GB, mTrkParams[iteration].MaxMemory / GB));
          if (mTrkParams[iteration].DropTFUponFailure) {
            dropTF = true;
          }
          break;
        }
        float trackletsPerCluster = mTraits->getTFNumberOfClusters() > 0 ? float(mTraits->getTFNumberOfTracklets()) / mTraits->getTFNumberOfClusters() : 0.f;
        if (trackletsPerCluster > mTrkParams[iteration].TrackletsPerClusterLimit) {
          error(std::format("Too many tracklets per cluster ({}) in iteration {} in ROF span {}-{}:, check the detector status and/or the selections. Current limit is {}",
                            trackletsPerCluster, iteration, iROFs, iROFs + mTrkParams[iteration].nROFsPerIterations, mTrkParams[iteration].TrackletsPerClusterLimit));
          break;
        }
 
        timeCells += evaluateTask(&Tracker::computeCells, "Cell finding", evalLog, iteration);
        nCells += mTraits->getTFNumberOfCells();
        if (!mTimeFrame->checkMemory(mTrkParams[iteration].MaxMemory)) {
          mTimeFrame->printSliceInfo(iROFs, mTrkParams[iteration].nROFsPerIterations);
          error(std::format("Too much memory used during cell finding in iteration {} in ROF span {}-{}: {:.2f} GB. Current limit is {:.2f} GB, check the detector status and/or the selections.",
                            iteration, iROFs, iROFs + mTrkParams[iteration].nROFsPerIterations, mTimeFrame->getArtefactsMemory() / GB, mTrkParams[iteration].MaxMemory / GB));
          if (mTrkParams[iteration].DropTFUponFailure) {
            dropTF = true;
          }
          break;
        }
        float cellsPerCluster = mTraits->getTFNumberOfClusters() > 0 ? float(mTraits->getTFNumberOfCells()) / mTraits->getTFNumberOfClusters() : 0.f;
        if (cellsPerCluster > mTrkParams[iteration].CellsPerClusterLimit) {
          error(std::format("Too many cells per cluster ({}) in iteration {} in ROF span {}-{}, check the detector status and/or the selections. Current limit is {}",
                            cellsPerCluster, iteration, iROFs, iROFs + mTrkParams[iteration].nROFsPerIterations, mTrkParams[iteration].CellsPerClusterLimit));
          break;
        }
 
        timeNeighbours += evaluateTask(&Tracker::findCellsNeighbours, "Neighbour finding", evalLog, iteration);
        nNeighbours += mTimeFrame->getNumberOfNeighbours();
        timeRoads += evaluateTask(&Tracker::findRoads, "Road finding", evalLog, iteration);
      }
      iVertex++;
    } while (iVertex < maxNvertices && !dropTF);
    logger(std::format(" - Tracklet finding: {} tracklets found in {:.2f} ms", nTracklets, timeTracklets));
    logger(std::format(" - Cell finding: {} cells found in {:.2f} ms", nCells, timeCells));
    logger(std::format(" - Neighbours finding: {} neighbours found in {:.2f} ms", nNeighbours, timeNeighbours));
    logger(std::format(" - Track finding: {} tracks found in {:.2f} ms", nTracks + mTimeFrame->getNumberOfTracks(), timeRoads));
    total += timeTracklets + timeCells + timeNeighbours + timeRoads;
    if (mTraits->supportsExtendTracks() && mTrkParams[iteration].UseTrackFollower && !dropTF) {
      int nExtendedTracks{-mTimeFrame->mNExtendedTracks}, nExtendedClusters{-mTimeFrame->mNExtendedUsedClusters};
      auto timeExtending = evaluateTask(&Tracker::extendTracks, "Extending tracks", [](const std::string&) {}, iteration);
      total += timeExtending;
      logger(std::format(" - Extending Tracks: {} extended tracks using {} clusters found in {:.2f} ms", nExtendedTracks + mTimeFrame->mNExtendedTracks, nExtendedClusters + mTimeFrame->mNExtendedUsedClusters, timeExtending));
    }
    if (dropTF) {
      error("...Dropping Timeframe...");
      mTimeFrame->dropTracks();
      ++mNumberOfDroppedTFs;
      return;
    }
  }
 
  if (mTraits->supportsFindShortPrimaries() && mTrkParams[0].FindShortTracks) {
    auto nTracksB = mTimeFrame->getNumberOfTracks();
    total += evaluateTask(&Tracker::findShortPrimaries, "Short primaries finding", logger);
    auto nTracksA = mTimeFrame->getNumberOfTracks();
    logger(std::format("  `-> found {} additional tracks", nTracksA - nTracksB));
  }
 
  if constexpr (constants::DoTimeBenchmarks) {
    logger(std::format("=== TimeFrame {} processing completed in: {:.2f} ms using {} thread(s) ===", mTimeFrameCounter, total, mTraits->getNThreads()));
  }
 
  if (mTimeFrame->hasMCinformation()) {
    computeTracksMClabels();
  }
  rectifyClusterIndices();
  ++mTimeFrameCounter;
  mTotalTime += total;
}
 
void Tracker::initialiseTimeFrame(int& iteration)
{
  mTraits->initialiseTimeFrame(iteration);
}
 
void Tracker::computeTracklets(int& iteration, int& iROFslice, int& iVertex)
{
  mTraits->computeLayerTracklets(iteration, iROFslice, iVertex);
}
 
void Tracker::computeCells(int& iteration)
{
  mTraits->computeLayerCells(iteration);
}
 
void Tracker::findCellsNeighbours(int& iteration)
{
  mTraits->findCellsNeighbours(iteration);
}
 
void Tracker::findRoads(int& iteration)
{
  mTraits->findRoads(iteration);
}
 
void Tracker::extendTracks(int& iteration)
{
  mTraits->extendTracks(iteration);
}
 
void Tracker::findShortPrimaries()
{
  mTraits->findShortPrimaries();
}
 
void Tracker::computeRoadsMClabels()
{
  if (!mTimeFrame->hasMCinformation()) {
    return;
  }
 
  mTimeFrame->initialiseRoadLabels();
 
  int roadsNum{static_cast<int>(mTimeFrame->getRoads().size())};
 
  for (int iRoad{0}; iRoad < roadsNum; ++iRoad) {
 
    Road<5>& currentRoad{mTimeFrame->getRoads()[iRoad]};
    std::vector<std::pair<MCCompLabel, size_t>> occurrences;
    bool isFakeRoad{false};
    bool isFirstRoadCell{true};
 
    for (int iCell{0}; iCell < mTrkParams[0].CellsPerRoad(); ++iCell) {
      const int currentCellIndex{currentRoad[iCell]};
 
      if (currentCellIndex == constants::its::UnusedIndex) {
        if (isFirstRoadCell) {
          continue;
        } else {
          break;
        }
      }
 
      const CellSeed& currentCell{mTimeFrame->getCells()[iCell][currentCellIndex]};
 
      if (isFirstRoadCell) {
 
        const int cl0index{mTimeFrame->getClusters()[iCell][currentCell.getFirstClusterIndex()].clusterId};
        auto cl0labs{mTimeFrame->getClusterLabels(iCell, cl0index)};
        bool found{false};
        for (size_t iOcc{0}; iOcc < occurrences.size(); ++iOcc) {
          std::pair<o2::MCCompLabel, size_t>& occurrence = occurrences[iOcc];
          for (auto& label : cl0labs) {
            if (label == occurrence.first) {
              ++occurrence.second;
              found = true;
              // break; // uncomment to stop to the first hit
            }
          }
        }
        if (!found) {
          for (auto& label : cl0labs) {
            occurrences.emplace_back(label, 1);
          }
        }
 
        const int cl1index{mTimeFrame->getClusters()[iCell + 1][currentCell.getSecondClusterIndex()].clusterId};
 
        const auto& cl1labs{mTimeFrame->getClusterLabels(iCell + 1, cl1index)};
        found = false;
        for (size_t iOcc{0}; iOcc < occurrences.size(); ++iOcc) {
          std::pair<o2::MCCompLabel, size_t>& occurrence = occurrences[iOcc];
          for (auto& label : cl1labs) {
            if (label == occurrence.first) {
              ++occurrence.second;
              found = true;
              // break; // uncomment to stop to the first hit
            }
          }
        }
        if (!found) {
          for (auto& label : cl1labs) {
            occurrences.emplace_back(label, 1);
          }
        }
 
        isFirstRoadCell = false;
      }
 
      const int cl2index{mTimeFrame->getClusters()[iCell + 2][currentCell.getThirdClusterIndex()].clusterId};
      const auto& cl2labs{mTimeFrame->getClusterLabels(iCell + 2, cl2index)};
      bool found{false};
      for (size_t iOcc{0}; iOcc < occurrences.size(); ++iOcc) {
        std::pair<o2::MCCompLabel, size_t>& occurrence = occurrences[iOcc];
        for (auto& label : cl2labs) {
          if (label == occurrence.first) {
            ++occurrence.second;
            found = true;
            // break; // uncomment to stop to the first hit
          }
        }
      }
      if (!found) {
        for (auto& label : cl2labs) {
          occurrences.emplace_back(label, 1);
        }
      }
    }
 
    std::sort(occurrences.begin(), occurrences.end(), [](auto e1, auto e2) {
      return e1.second > e2.second;
    });
 
    auto maxOccurrencesValue = occurrences[0].first;
    mTimeFrame->setRoadLabel(iRoad, maxOccurrencesValue.getRawValue(), isFakeRoad);
  }
}
 
void Tracker::computeTracksMClabels()
{
  for (int iROF{0}; iROF < mTimeFrame->getNrof(); ++iROF) {
    for (auto& track : mTimeFrame->getTracks(iROF)) {
      std::vector<std::pair<MCCompLabel, size_t>> occurrences;
      occurrences.clear();
 
      for (int iCluster = 0; iCluster < TrackITSExt::MaxClusters; ++iCluster) {
        const int index = track.getClusterIndex(iCluster);
        if (index == constants::its::UnusedIndex) {
          continue;
        }
        auto labels = mTimeFrame->getClusterLabels(iCluster, index);
        bool found{false};
        for (size_t iOcc{0}; iOcc < occurrences.size(); ++iOcc) {
          std::pair<o2::MCCompLabel, size_t>& occurrence = occurrences[iOcc];
          for (auto& label : labels) {
            if (label == occurrence.first) {
              ++occurrence.second;
              found = true;
              // break; // uncomment to stop to the first hit
            }
          }
        }
        if (!found) {
          for (auto& label : labels) {
            occurrences.emplace_back(label, 1);
          }
        }
      }
      std::sort(std::begin(occurrences), std::end(occurrences), [](auto e1, auto e2) {
        return e1.second > e2.second;
      });
 
      auto maxOccurrencesValue = occurrences[0].first;
      uint32_t pattern = track.getPattern();
      // set fake clusters pattern
      for (int ic{TrackITSExt::MaxClusters}; ic--;) {
        auto clid = track.getClusterIndex(ic);
        if (clid != constants::its::UnusedIndex) {
          auto labelsSpan = mTimeFrame->getClusterLabels(ic, clid);
          for (auto& currentLabel : labelsSpan) {
            if (currentLabel == maxOccurrencesValue) {
              pattern |= 0x1 << (16 + ic); // set bit if correct
              break;
            }
          }
        }
      }
      track.setPattern(pattern);
      if (occurrences[0].second < track.getNumberOfClusters()) {
        maxOccurrencesValue.setFakeFlag();
      }
      mTimeFrame->getTracksLabel(iROF).emplace_back(maxOccurrencesValue);
    }
  }
}
 
void Tracker::rectifyClusterIndices()
{
  for (int iROF{0}; iROF < mTimeFrame->getNrof(); ++iROF) {
    for (auto& track : mTimeFrame->getTracks(iROF)) {
      for (int iCluster = 0; iCluster < TrackITSExt::MaxClusters; ++iCluster) {
        const int index = track.getClusterIndex(iCluster);
        if (index != constants::its::UnusedIndex) {
          track.setExternalClusterIndex(iCluster, mTimeFrame->getClusterExternalIndex(iCluster, index));
        }
      }
    }
  }
}
 
void Tracker::getGlobalConfiguration()
{
  auto& tc = o2::its::TrackerParamConfig::Instance();
  if (tc.useMatCorrTGeo) {
    mTraits->setCorrType(o2::base::PropagatorImpl<float>::MatCorrType::USEMatCorrTGeo);
  } else if (tc.useFastMaterial) {
    mTraits->setCorrType(o2::base::PropagatorImpl<float>::MatCorrType::USEMatCorrNONE);
  } else {
    mTraits->setCorrType(o2::base::PropagatorImpl<float>::MatCorrType::USEMatCorrLUT);
  }
  setNThreads(tc.nThreads);
  int nROFsPerIterations = tc.nROFsPerIterations > 0 ? tc.nROFsPerIterations : -1;
  if (tc.nOrbitsPerIterations > 0) {
  }
  for (auto& params : mTrkParams) {
    if (params.NLayers == 7) {
      for (int i{0}; i < 7; ++i) {
        params.SystErrorY2[i] = tc.sysErrY2[i] > 0 ? tc.sysErrY2[i] : params.SystErrorY2[i];
        params.SystErrorZ2[i] = tc.sysErrZ2[i] > 0 ? tc.sysErrZ2[i] : params.SystErrorZ2[i];
      }
    }
    params.DeltaROF = tc.deltaRof;
    params.DoUPCIteration = tc.doUPCIteration;
    params.MaxChi2ClusterAttachment = tc.maxChi2ClusterAttachment > 0 ? tc.maxChi2ClusterAttachment : params.MaxChi2ClusterAttachment;
    params.MaxChi2NDF = tc.maxChi2NDF > 0 ? tc.maxChi2NDF : params.MaxChi2NDF;
    params.PhiBins = tc.LUTbinsPhi > 0 ? tc.LUTbinsPhi : params.PhiBins;
    params.ZBins = tc.LUTbinsZ > 0 ? tc.LUTbinsZ : params.ZBins;
    params.PVres = tc.pvRes > 0 ? tc.pvRes : params.PVres;
    params.NSigmaCut *= tc.nSigmaCut > 0 ? tc.nSigmaCut : 1.f;
    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;
    params.SaveTimeBenchmarks = tc.saveTimeBenchmarks;
    params.FataliseUponFailure = tc.fataliseUponFailure;
    params.DropTFUponFailure = tc.dropTFUponFailure;
    for (int iD{0}; iD < 3; ++iD) {
      params.Diamond[iD] = tc.diamondPos[iD];
    }
    params.UseDiamond = tc.useDiamond;
    if (tc.maxMemory) {
      params.MaxMemory = tc.maxMemory;
    }
    if (tc.useTrackFollower > 0) {
      params.UseTrackFollower = true;
      // Bit 0: Allow for mixing of top&bot extension --> implies Bits 1&2 set
      // Bit 1: Allow for top extension
      // Bit 2: Allow for bot extension
      params.UseTrackFollowerMix = ((tc.useTrackFollower & (1 << 0)) != 0);
      params.UseTrackFollowerTop = ((tc.useTrackFollower & (1 << 1)) != 0);
      params.UseTrackFollowerBot = ((tc.useTrackFollower & (1 << 2)) != 0);
      params.TrackFollowerNSigmaCutZ = tc.trackFollowerNSigmaZ;
      params.TrackFollowerNSigmaCutPhi = tc.trackFollowerNSigmaPhi;
    }
    if (tc.cellsPerClusterLimit >= 0) {
      params.CellsPerClusterLimit = tc.cellsPerClusterLimit;
    }
    if (tc.trackletsPerClusterLimit >= 0) {
      params.TrackletsPerClusterLimit = tc.trackletsPerClusterLimit;
    }
    if (tc.findShortTracks >= 0) {
      params.FindShortTracks = tc.findShortTracks;
    }
  }
}
 
void Tracker::adoptTimeFrame(TimeFrame& tf)
{
  mTimeFrame = &tf;
  mTraits->adoptTimeFrame(&tf);
}
 
void Tracker::setBz(float bz)
{
  mTraits->setBz(bz);
}
 
void Tracker::setCorrType(const o2::base::PropagatorImpl<float>::MatCorrType type)
{
  mTraits->setCorrType(type);
}
 
bool Tracker::isMatLUT() const
{
  return mTraits->isMatLUT();
}
 
void Tracker::setNThreads(int n)
{
  mTraits->setNThreads(n);
}
 
int Tracker::getNThreads() const
{
  return mTraits->getNThreads();
}
 
void Tracker::printSummary() const
{
  LOGP(info, "Tracker summary: Processed {} TFs (dropped {}) in TOT={:.2f} s, AVG/TF={:.2f} s", mTimeFrameCounter, mNumberOfDroppedTFs, mTotalTime * 1.e-3, mTotalTime * 1.e-3 / ((mTimeFrameCounter > 0) ? (double)mTimeFrameCounter : -1.0));
}
 
} // namespace its
} // namespace o2