ITS: hole-tracking

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

Author: f3sch

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

@@ -16,6 +16,8 @@
 #ifndef TRACKINGITSU_INCLUDE_CACELL_H_
 #define TRACKINGITSU_INCLUDE_CACELL_H_
 
+#include <cstdint>
+
 #include "ITStracking/Constants.h"
 #include "DataFormatsITS/TimeEstBC.h"
 #include "ReconstructionDataFormats/Track.h"
@@ -28,7 +30,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 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 +44,26 @@ 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; }
+  GPUhd() static int firstSetBit(uint16_t mask)
+  {
+    for (int layer = 0; layer < 16; ++layer) {
+      if (mask & (uint16_t(1) << layer)) {
+        return layer;
+      }
+    }
+    return constants::UnusedIndex;
+  }
+  GPUhd() static int layerSlot(uint16_t mask, int layer)
+  {
+    if ((mask & (uint16_t(1) << layer)) == 0) {
+      return constants::UnusedIndex;
+    }
+    int slot = 0;
+    for (int iLayer = 0; iLayer < layer; ++iLayer) {
+      slot += (mask >> iLayer) & uint16_t(1);
+    }
+    return slot;
+  }
 
  protected:
   GPUhdDefault() SeedBase() = default;
@@ -69,11 +93,18 @@ class CellSeed final : public SeedBase<constants::ClustersPerCell>
   using Base = SeedBase<constants::ClustersPerCell>;
 
  public:
+  struct HitLayerMaskTag {
+  };
+
   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(HitLayerMaskTag{}, 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(HitLayerMaskTag, 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 +123,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 +145,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/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>;
 
@@ -135,6 +137,8 @@ struct TimeFrame {
   const auto& getIndexTableUtils() const { return mIndexTableUtils; }
   const auto& getROFOverlapTable() const { return mROFOverlapTable; }
   const auto& getROFOverlapTableView() const { return mROFOverlapTableView; }
+  const auto& getTrackingTopology() const { return mTrackingTopology; }
+  const auto& getTrackingTopologyView() const { return mTrackingTopologyView; }
   void setROFOverlapTable(ROFOverlapTableN table)
   {
     mROFOverlapTable = std::move(table);
@@ -193,6 +197,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 +278,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
@@ -319,6 +325,8 @@ struct TimeFrame {
   IndexTableUtilsN mIndexTableUtils;
   ROFOverlapTableN mROFOverlapTable;
   ROFOverlapTableN::View mROFOverlapTableView;
+  TrackingTopologyN mTrackingTopology;
+  TrackingTopologyN::View mTrackingTopologyView;
   ROFVertexLookupTableN mROFVertexLookupTable;
   ROFVertexLookupTableN::View mROFVertexLookupTableView;
   ROFMaskTableN mMultiplicityCutMask;

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

@@ -53,7 +53,7 @@ class TrackerTraits
   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