ITS: hole-tracking

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

Author: f3sch

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

@@ -16,10 +16,14 @@
 #ifndef TRACKINGITSU_INCLUDE_CACELL_H_
 #define TRACKINGITSU_INCLUDE_CACELL_H_
 
+#include <cstdint>
+
 #include "ITStracking/Constants.h"
+#include "ITStracking/MathUtils.h"
 #include "DataFormatsITS/TimeEstBC.h"
 #include "ReconstructionDataFormats/Track.h"
 #include "GPUCommonDef.h"
+#include "GPUCommonMath.h"
 
 namespace o2::its
 {
@@ -28,7 +32,9 @@ template <int NClusters>
 class SeedBase : public o2::track::TrackParCovF
 {
  public:
-  GPUhd() int getInnerLayer() const { return getUserField(); }
+  GPUhd() uint16_t getHitLayerMask() const { return static_cast<uint16_t>(getUserField()); }
+  GPUhd() void setHitLayerMask(uint16_t mask) { setUserField(mask); }
+  GPUhd() int getInnerLayer() const { return math_utils::firstSetBit(getHitLayerMask()); }
   GPUhd() int getFirstTrackletIndex() const { return mTracklets[0]; };
   GPUhd() void setFirstTrackletIndex(int trkl) { mTracklets[0] = trkl; };
   GPUhd() int getSecondTrackletIndex() const { return mTracklets[1]; };
@@ -40,6 +46,14 @@ class SeedBase : public o2::track::TrackParCovF
   GPUhd() int* getLevelPtr() { return &mLevel; }
   GPUhd() auto& getTimeStamp() noexcept { return mTime; }
   GPUhd() const auto& getTimeStamp() const noexcept { return mTime; }
+  GPUhdi() static int layerSlot(uint16_t mask, int layer)
+  {
+    if ((mask & (uint16_t(1) << layer)) == 0) {
+      return constants::UnusedIndex;
+    }
+    const uint32_t lowerLayers = (uint32_t(1) << layer) - 1;
+    return static_cast<int>(o2::gpu::GPUCommonMath::Popcount(static_cast<uint32_t>(mask) & lowerLayers));
+  }
 
  protected:
   GPUhdDefault() SeedBase() = default;
@@ -71,9 +85,13 @@ class CellSeed final : public SeedBase<constants::ClustersPerCell>
  public:
   GPUhdDefault() CellSeed() = default;
   GPUhd() CellSeed(int innerL, int cl0, int cl1, int cl2, int trkl0, int trkl1, const o2::track::TrackParCovF& tpc, float chi2, const TimeEstBC& time)
+    : CellSeed(static_cast<uint16_t>((uint16_t(1) << innerL) | (uint16_t(1) << (innerL + 1)) | (uint16_t(1) << (innerL + 2))), cl0, cl1, cl2, trkl0, trkl1, tpc, chi2, time)
+  {
+  }
+  GPUhd() CellSeed(uint16_t hitLayerMask, int cl0, int cl1, int cl2, int trkl0, int trkl1, const o2::track::TrackParCovF& tpc, float chi2, const TimeEstBC& time)
     : Base(tpc, chi2, 1, time)
   {
-    setUserField(innerL);
+    setHitLayerMask(hitLayerMask);
     auto& clusters = this->clustersRaw();
     clusters[0] = cl0;
     clusters[1] = cl1;
@@ -92,12 +110,12 @@ class CellSeed final : public SeedBase<constants::ClustersPerCell>
   GPUhd() int getThirdClusterIndex() const { return this->clustersRaw()[2]; };
   GPUhd() auto& getClusters() { return this->clustersRaw(); }
   GPUhd() const auto& getClusters() const { return this->clustersRaw(); }
-  /// getCluster takes an ABSOLUTE layer index and returns UnusedIndex if the
-  /// layer is outside the 3 stored slots (innerL, innerL+1, innerL+2).
+  /// getCluster takes an ABSOLUTE layer index. Compact cluster slots are
+  /// mapped to absolute layers by set-bit order in the hit-layer mask.
   GPUhd() int getCluster(int layer) const
   {
-    const int rel = layer - getInnerLayer();
-    return (rel >= 0 && rel < constants::ClustersPerCell) ? this->clustersRaw()[rel] : constants::UnusedIndex;
+    const int slot = Base::layerSlot(getHitLayerMask(), layer);
+    return (slot >= 0 && slot < constants::ClustersPerCell) ? this->clustersRaw()[slot] : constants::UnusedIndex;
   }
 };
 
@@ -114,29 +132,45 @@ class TrackSeed final : public SeedBase<NLayers>
   GPUhd() TrackSeed(const CellSeed& cs)
     : Base(static_cast<const o2::track::TrackParCovF&>(cs), cs.getChi2(), cs.getLevel(), cs.getTimeStamp())
   {
-    this->setUserField(cs.getInnerLayer());
+    this->setHitLayerMask(cs.getHitLayerMask());
     this->setFirstTrackletIndex(cs.getFirstTrackletIndex());
     this->setSecondTrackletIndex(cs.getSecondTrackletIndex());
-    const int innerL = cs.getInnerLayer();
     auto& clusters = this->clustersRaw();
-    clusters[innerL + 0] = cs.getFirstClusterIndex();
-    clusters[innerL + 1] = cs.getSecondClusterIndex();
-    clusters[innerL + 2] = cs.getThirdClusterIndex();
+    int slot = 0;
+    const uint16_t hitMask = cs.getHitLayerMask();
+    for (int layer = 0; layer < NLayers; ++layer) {
+      if (hitMask & (uint16_t(1) << layer)) {
+        clusters[layer] = cs.getClusters()[slot++];
+      }
+    }
   }
   GPUhdDefault() TrackSeed(const TrackSeed&) = default;
   GPUhdDefault() ~TrackSeed() = default;
   GPUhdDefault() TrackSeed(TrackSeed&&) = default;
   GPUhdDefault() TrackSeed& operator=(const TrackSeed&) = default;
   GPUhdDefault() TrackSeed& operator=(TrackSeed&&) = default;
 
-  /// Three-cluster view of the original cell — note: innerL (UserField) is not
-  /// updated when processNeighbours extends the cluster list leftward.
-  GPUhd() int getFirstClusterIndex() const { return this->clustersRaw()[this->getUserField()]; }
-  GPUhd() int getSecondClusterIndex() const { return this->clustersRaw()[this->getUserField() + 1]; }
-  GPUhd() int getThirdClusterIndex() const { return this->clustersRaw()[this->getUserField() + 2]; }
+  GPUhd() int getFirstClusterIndex() const { return getClusterBySlot(0); }
+  GPUhd() int getSecondClusterIndex() const { return getClusterBySlot(1); }
+  GPUhd() int getThirdClusterIndex() const { return getClusterBySlot(2); }
   GPUhd() auto& getClusters() { return this->clustersRaw(); }
   GPUhd() const auto& getClusters() const { return this->clustersRaw(); }
   GPUhd() int getCluster(int layer) const { return this->clustersRaw()[layer]; }
+
+ private:
+  GPUhd() int getClusterBySlot(int requestedSlot) const
+  {
+    int slot = 0;
+    const uint16_t hitMask = this->getHitLayerMask();
+    for (int layer = 0; layer < NLayers; ++layer) {
+      if (hitMask & (uint16_t(1) << layer)) {
+        if (slot++ == requestedSlot) {
+          return this->clustersRaw()[layer];
+        }
+      }
+    }
+    return constants::UnusedIndex;
+  }
 };
 
 } // namespace o2::its

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

@@ -68,6 +68,8 @@ struct TrackingParameters {
   bool AllowSharingFirstCluster = false;
   int ClusterSharing = 0;
   int MinTrackLength = 7;
+  int MaxHoles = 0;
+  uint16_t HoleLayerMask = 0;
   float NSigmaCut = 5;
   float PVres = 1.e-2f;
   /// Trackleting cuts

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

@@ -16,6 +16,8 @@
 #ifndef O2_ITS_TRACKING_MATHUTILS_H_
 #define O2_ITS_TRACKING_MATHUTILS_H_
 
+#include <cstdint>
+
 #include "CommonConstants/MathConstants.h"
 #include "ITStracking/Constants.h"
 #include "MathUtils/Utils.h"
@@ -115,6 +117,16 @@ GPUhdi() constexpr float SqDiff(float x, float y)
   return Sq(x - y);
 }
 
+GPUhdi() int firstSetBit(uint32_t mask)
+{
+  return mask ? static_cast<int>(o2::gpu::GPUCommonMath::Ctz(mask)) : constants::UnusedIndex;
+}
+
+GPUhdi() int lastSetBit(uint32_t mask)
+{
+  return mask ? 31 - static_cast<int>(o2::gpu::GPUCommonMath::Clz(mask)) : constants::UnusedIndex;
+}
+
 GPUhdi() float MSangle(float mass, float p, float xX0)
 {
   float beta = p / o2::gpu::CAMath::Hypot(mass, p);

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

@@ -32,6 +32,7 @@
 #include "ITStracking/ExternalAllocator.h"
 #include "ITStracking/BoundedAllocator.h"
 #include "ITStracking/ROFLookupTables.h"
+#include "ITStracking/TrackingTopology.h"
 #include "SimulationDataFormat/MCCompLabel.h"
 #include "SimulationDataFormat/MCTruthContainer.h"
 
@@ -66,6 +67,7 @@ struct TimeFrame {
   using ROFOverlapTableN = ROFOverlapTable<NLayers>;
   using ROFVertexLookupTableN = ROFVertexLookupTable<NLayers>;
   using ROFMaskTableN = ROFMaskTable<NLayers>;
+  using TrackingTopologyN = TrackingTopology<NLayers>;
   using TrackSeedN = TrackSeed<NLayers>;
   friend class gpu::TimeFrameGPU<NLayers>;
 
@@ -112,10 +114,8 @@ struct TimeFrame {
   auto& getMaxRs() { return mMaxR; }
   float getMinR(int layer) const { return mMinR[layer]; }
   float getMaxR(int layer) const { return mMaxR[layer]; }
-  float getMSangle(int layer) const { return mMSangles[layer]; }
-  auto& getMSangles() { return mMSangles; }
-  float getPhiCut(int layer) const { return mPhiCuts[layer]; }
-  auto& getPhiCuts() { return mPhiCuts; }
+  float getTransitionPhiCut(int transitionId) const { return mTransitionPhiCuts[transitionId]; }
+  float getTransitionMSAngle(int transitionId) const { return mTransitionMSAngles[transitionId]; }
   float getPositionResolution(int layer) const { return mPositionResolution[layer]; }
   auto& getPositionResolutions() { return mPositionResolution; }
 
@@ -135,6 +135,8 @@ struct TimeFrame {
   const auto& getIndexTableUtils() const { return mIndexTableUtils; }
   const auto& getROFOverlapTable() const { return mROFOverlapTable; }
   const auto& getROFOverlapTableView() const { return mROFOverlapTableView; }
+  const auto& getTrackingTopologies() const { return mTrackingTopologies; }
+  const auto& getTrackingTopologyView() const { return mTrackingTopologyView; }
   void setROFOverlapTable(ROFOverlapTableN table)
   {
     mROFOverlapTable = std::move(table);
@@ -177,7 +179,10 @@ struct TimeFrame {
   auto& getCellsLabel(int layer) { return mCellLabels[layer]; }
 
   bool hasMCinformation() const { return mClusterLabels[0] != nullptr; }
-  void initialise(const TrackingParameters& trkParam, const int maxLayers = NLayers);
+  void initVertexingTopology(const TrackingParameters& trkParam);
+  void initDefaultTrackingTopology(const TrackingParameters& trkParam, const int maxLayers = NLayers);
+  void initTrackingTopologies(gsl::span<const TrackingParameters> trkParams, const int maxLayers = NLayers);
+  void initialise(const TrackingParameters& trkParam, const int maxLayers = NLayers, const int iteration = constants::UnusedIndex);
 
   bool isClusterUsed(int layer, int clusterId) const { return mUsedClusters[layer][clusterId]; }
   void markUsedCluster(int layer, int clusterId) { mUsedClusters[layer][clusterId] = true; }
@@ -193,6 +198,7 @@ struct TimeFrame {
 
   auto& getCellsLookupTable() { return mCellsLookupTable; }
   auto& getCellsNeighbours() { return mCellsNeighbours; }
+  auto& getCellsNeighboursTopology() { return mCellsNeighboursTopology; }
   auto& getCellsNeighboursLUT() { return mCellsNeighboursLUT; }
   auto& getTracks() { return mTracks; }
   auto& getTracksLabel() { return mTracksLabel; }
@@ -273,6 +279,7 @@ struct TimeFrame {
   bounded_vector<TrackITSExt> mTracks;
   bounded_vector<MCCompLabel> mTracksLabel;
   std::vector<bounded_vector<int>> mCellsNeighbours;
+  std::vector<bounded_vector<int>> mCellsNeighboursTopology;
   std::vector<bounded_vector<int>> mCellsLookupTable;
 
   const o2::base::PropagatorImpl<float>* mPropagatorDevice = nullptr; // Needed only for GPU
@@ -292,8 +299,8 @@ struct TimeFrame {
   bool isBeamPositionOverridden = false;
   std::array<float, NLayers> mMinR;
   std::array<float, NLayers> mMaxR;
-  bounded_vector<float> mMSangles;
-  bounded_vector<float> mPhiCuts;
+  bounded_vector<float> mTransitionPhiCuts;
+  bounded_vector<float> mTransitionMSAngles;
   bounded_vector<float> mPositionResolution;
   std::array<bounded_vector<uint8_t>, NLayers> mClusterSize;
 
@@ -319,6 +326,10 @@ struct TimeFrame {
   IndexTableUtilsN mIndexTableUtils;
   ROFOverlapTableN mROFOverlapTable;
   ROFOverlapTableN::View mROFOverlapTableView;
+  TrackingTopologyN mVertexingTopology;
+  TrackingTopologyN mDefaultTrackingTopology;
+  std::vector<TrackingTopologyN> mTrackingTopologies;
+  TrackingTopologyN::View mTrackingTopologyView;
   ROFVertexLookupTableN mROFVertexLookupTable;
   ROFVertexLookupTableN::View mROFVertexLookupTableView;
   ROFMaskTableN mMultiplicityCutMask;

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

@@ -45,15 +45,15 @@ class TrackerTraits
 
   virtual ~TrackerTraits() = default;
   virtual void adoptTimeFrame(TimeFrame<NLayers>* tf) { mTimeFrame = tf; }
-  virtual void initialiseTimeFrame(const int iteration) { mTimeFrame->initialise(mTrkParams[iteration], mTrkParams[iteration].NLayers); }
+  virtual void initialiseTimeFrame(const int iteration) { mTimeFrame->initialise(mTrkParams[iteration], mTrkParams[iteration].NLayers, iteration); }
 
   virtual void computeLayerTracklets(const int iteration, int iVertex);
   virtual void computeLayerCells(const int iteration);
   virtual void findCellsNeighbours(const int iteration);
   virtual void findRoads(const int iteration);
 
   template <typename InputSeed>
-  void processNeighbours(int iteration, int iLayer, int iLevel, const bounded_vector<InputSeed>& currentCellSeed, const bounded_vector<int>& currentCellId, bounded_vector<TrackSeedN>& updatedCellSeed, bounded_vector<int>& updatedCellId);
+  void processNeighbours(int iteration, int defaultCellTopologyId, int iLevel, const bounded_vector<InputSeed>& currentCellSeed, const bounded_vector<int>& currentCellId, const bounded_vector<int>& currentCellTopologyId, bounded_vector<TrackSeedN>& updatedCellSeed, bounded_vector<int>& updatedCellId, bounded_vector<int>& updatedCellTopologyId);
 
   void acceptTracks(int iteration, bounded_vector<TrackITSExt>& tracks, bounded_vector<bounded_vector<int>>& firstClusters, bounded_vector<bounded_vector<int>>& sharedFirstClusters);
   void markTracks(int iteration, bounded_vector<bounded_vector<int>>& sharedFirstClusters);

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

@@ -71,6 +71,8 @@ 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[constants::MaxIter] = {};                                        // minimum track length at each iteration, used only if >0, otherwise use code defaults
   uint8_t startLayerMask[constants::MaxIter] = {};                                     // mask of start layer for this iteration (if >0)
+  int maxHolesIter[constants::MaxIter] = {};                                           // maximum number of missing internal layers allowed in the CA topology for each iteration
+  uint16_t holeLayerMaskIter[constants::MaxIter] = {};                                 // layers that may be skipped by the CA topology for each iteration
   float minPtIterLgt[constants::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