ITS: optionally restrict CA to a subset of layers + support inactive layers

Detectors/ITSMFT/ITS/tracking/GPU/ITStrackingGPU/TrackingKernels.h

@@ -193,8 +193,9 @@ void processNeighboursHandler(const int startLevel,
                               const float MaxChi2ClusterAttachment,
                               const float maxChi2NDF,
                               const int maxHoles,
-                              const int minTrackLength,
+                              const int minSeedingClusters,
                               const LayerMask holeLayerMask,
+                              const LayerMask nonCountingLayerMask,
                               const std::vector<float>& layerxX0Host,
                               const o2::base::Propagator* propagator,
                               const o2::base::PropagatorF::MatCorrType matCorrType,

Detectors/ITSMFT/ITS/tracking/GPU/cuda/TrackerTraitsGPU.cxx

@@ -308,6 +308,8 @@ void TrackerTraitsGPU<NLayers>::findRoads(const int iteration)
   const bool extendTop = this->mTrkParams[iteration].PassFlags[IterationStep::TrackFollowerTop];
   const bool extendBot = this->mTrkParams[iteration].PassFlags[IterationStep::TrackFollowerBot];
   const bool extendTracks = extendTop || extendBot;
+  const auto notSeedingLayers = this->mTrkParams[iteration].getNotSeedingLayerMask();
+  const auto minSeedingClusters = this->mTrkParams[iteration].getMinSeedingClusters();
   for (int startLevel{this->mTrkParams[iteration].CellsPerRoad()}; startLevel >= this->mTrkParams[iteration].CellMinimumLevel(); --startLevel) {
     bounded_vector<TrackSeed<NLayers>> trackSeeds(this->getMemoryPool().get());
     for (int startCellTopologyId{0}; startCellTopologyId < hostTopology.nCells; ++startCellTopologyId) {
@@ -331,8 +333,9 @@ void TrackerTraitsGPU<NLayers>::findRoads(const int iteration)
                                         this->mTrkParams[iteration].MaxChi2ClusterAttachment,
                                         this->mTrkParams[iteration].MaxChi2NDF,
                                         this->mTrkParams[iteration].MaxHoles,
-                                        this->mTrkParams[iteration].MinTrackLength,
+                                        minSeedingClusters,
                                         this->mTrkParams[iteration].HoleLayerMask,
+                                        notSeedingLayers,
                                         this->mTrkParams[iteration].LayerxX0,
                                         mTimeFrameGPU->getDevicePropagator(),
                                         this->mTrkParams[iteration].CorrType,

Detectors/ITSMFT/ITS/tracking/GPU/cuda/TrackingKernels.cu

@@ -32,8 +32,8 @@
 #include "ITStracking/Tracklet.h"
 #include "ITStracking/Cluster.h"
 #include "ITStracking/Cell.h"
-#include "ITStracking/TrackFollower.h"
 #include "ITStracking/TrackHelpers.h"
+#include "ITStracking/TrackFollower.h"
 #include "DataFormatsITS/TrackITS.h"
 #include "ITStrackingGPU/TrackingKernels.h"
 #include "ITStrackingGPU/Utils.h"
@@ -86,23 +86,6 @@ struct is_valid_pair {
   }
 };
 
-template <int NLayers>
-struct seed_selector {
-  float mMaxQ2Pt;
-  float mMaxChi2;
-  int mMaxHoles;
-  int mMinTrackLength;
-  LayerMask mHoleLayerMask;
-
-  GPUhd() seed_selector(float maxQ2Pt, float maxChi2, int maxHoles, int minTrackLength, LayerMask holeLayerMask) : mMaxQ2Pt(maxQ2Pt), mMaxChi2(maxChi2), mMaxHoles(maxHoles), mMinTrackLength(minTrackLength), mHoleLayerMask(holeLayerMask) {}
-  GPUhd() bool operator()(const TrackSeed<NLayers>& seed) const
-  {
-    return !(seed.getQ2Pt() > mMaxQ2Pt || seed.getChi2() > mMaxChi2) &&
-           seed.getHitLayerMask().length() >= mMinTrackLength &&
-           seed.getHitLayerMask().isAllowed(mMaxHoles, mHoleLayerMask);
-  }
-};
-
 struct compare_track_chi2 {
   GPUhd() bool operator()(const TrackITSExt& a, const TrackITSExt& b) const
   {
@@ -975,8 +958,9 @@ void processNeighboursHandler(const int startLevel,
                               const float maxChi2ClusterAttachment,
                               const float maxChi2NDF,
                               const int maxHoles,
-                              const int minTrackLength,
+                              const int minSeedingClusters,
                               const LayerMask holeLayerMask,
+                              const LayerMask nonCountingLayerMask,
                               const std::vector<float>& layerxX0Host,
                               const o2::base::Propagator* propagator,
                               const o2::base::PropagatorF::MatCorrType matCorrType,
@@ -1110,7 +1094,7 @@ void processNeighboursHandler(const int startLevel,
   }
   GPUChkErrS(cudaStreamSynchronize(gpu::Stream::DefaultStream));
   thrust::device_vector<TrackSeed<NLayers>, gpu::TypedAllocator<TrackSeed<NLayers>>> outSeeds(updatedCellSeed.size(), allocTrackSeed);
-  auto end = thrust::copy_if(nosync_policy, updatedCellSeed.begin(), updatedCellSeed.end(), outSeeds.begin(), gpu::seed_selector<NLayers>(1.e3, maxChi2NDF * ((startLevel + 2) * 2 - 5), maxHoles, minTrackLength, holeLayerMask));
+  auto end = thrust::copy_if(nosync_policy, updatedCellSeed.begin(), updatedCellSeed.end(), outSeeds.begin(), track::TrackSeedSelector<NLayers>{1.e3f, maxChi2NDF * ((startLevel + 2) * 2 - 5), maxHoles, minSeedingClusters, holeLayerMask, nonCountingLayerMask});
   auto s{end - outSeeds.begin()};
   seedsHost.reserve(seedsHost.size() + s);
   thrust::copy(outSeeds.begin(), outSeeds.begin() + s, std::back_inserter(seedsHost));
@@ -1380,8 +1364,9 @@ template void processNeighboursHandler<7>(const int startLevel,
                                           const float maxChi2ClusterAttachment,
                                           const float maxChi2NDF,
                                           const int maxHoles,
-                                          const int minTrackLength,
+                                          const int minSeedingClusters,
                                           const LayerMask holeLayerMask,
+                                          const LayerMask nonCountingLayerMask,
                                           const std::vector<float>& layerxX0Host,
                                           const o2::base::Propagator* propagator,
                                           const o2::base::PropagatorF::MatCorrType matCorrType,
@@ -1580,8 +1565,9 @@ template void processNeighboursHandler<11>(const int startLevel,
                                            const float maxChi2ClusterAttachment,
                                            const float maxChi2NDF,
                                            const int maxHoles,
-                                           const int minTrackLength,
+                                           const int minSeedingClusters,
                                            const LayerMask holeLayerMask,
+                                           const LayerMask nonCountingLayerMask,
                                            const std::vector<float>& layerxX0Host,
                                            const o2::base::Propagator* propagator,
                                            const o2::base::PropagatorF::MatCorrType matCorrType,

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

@@ -46,15 +46,42 @@ enum class IterationStep : uint16_t {
 using IterationSteps = o2::utils::EnumFlags<IterationStep>;
 
 struct TrackingParameters {
-  int CellMinimumLevel() const noexcept
+  LayerMask getActiveLayerMask() const noexcept
+  {
+    return LayerMask::span(0, NLayers - 1) & ~InactiveLayerMask;
+  }
+
+  LayerMask getSeedingLayerMask() const noexcept
+  {
+    const auto activeLayers = getActiveLayerMask();
+    return SeedingLayers.empty() ? activeLayers : (SeedingLayers & activeLayers);
+  }
+
+  LayerMask getNotSeedingLayerMask() const noexcept
+  {
+    return ~getSeedingLayerMask();
+  }
+
+  int getNSeedingLayers() const noexcept
+  {
+    return getSeedingLayerMask().count();
+  }
+
+  int getMinSeedingClusters() const noexcept
   {
     const int minClusters = MinTrackLength - (MaxHoles > 0 ? MaxHoles : 0);
-    const int effectiveMinClusters = minClusters > constants::ClustersPerCell ? minClusters : constants::ClustersPerCell;
-    return effectiveMinClusters - constants::ClustersPerCell + 1;
+    const int minClustersWithCells = minClusters > constants::ClustersPerCell ? minClusters : constants::ClustersPerCell;
+    const int nSeedingLayers = getNSeedingLayers();
+    return minClustersWithCells < nSeedingLayers ? minClustersWithCells : nSeedingLayers;
+  }
+
+  int CellMinimumLevel() const noexcept
+  {
+    return getMinSeedingClusters() - constants::ClustersPerCell + 1;
   }
-  int NeighboursPerRoad() const noexcept { return NLayers - 3; }
-  int CellsPerRoad() const noexcept { return NLayers - 2; }
-  int TrackletsPerRoad() const noexcept { return NLayers - 1; }
+  int NeighboursPerRoad() const noexcept { return getNSeedingLayers() - 3; }
+  int CellsPerRoad() const noexcept { return getNSeedingLayers() - 2; }
+  int TrackletsPerRoad() const noexcept { return getNSeedingLayers() - 1; }
   std::string asString() const;
 
   IterationSteps PassFlags{IterationStep::FirstPass, IterationStep::RebuildClusterLUT};
@@ -76,6 +103,8 @@ struct TrackingParameters {
   int MinTrackLength = 7;
   int MaxHoles = 0;
   LayerMask HoleLayerMask = 0;
+  LayerMask InactiveLayerMask = 0;
+  LayerMask SeedingLayers = 0;
   float NSigmaCut = 5;
   float PVres = 1.e-2f;
   /// Trackleting cuts

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

@@ -122,6 +122,9 @@ GPUhdi() constexpr float SqDiff(float x, float y)
 
 GPUhdi() float MSangle(float mass, float p, float xX0)
 {
+  if (xX0 <= 0.f) {
+    return 0.f;
+  }
   float beta = p / o2::gpu::CAMath::Hypot(mass, p);
   return 0.0136f * o2::gpu::CAMath::Sqrt(xX0) * (1.f + 0.038f * o2::gpu::CAMath::Log(xX0)) / (beta * p);
 }

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

@@ -21,6 +21,7 @@
 #include "ITStracking/Cell.h"
 #include "ITStracking/Cluster.h"
 #include "ITStracking/Constants.h"
+#include "ITStracking/LayerMask.h"
 #include "ITStracking/MathUtils.h"
 #include "ITStracking/TrackITSInternal.h"
 #include "DetectorsBase/Propagator.h"
@@ -45,6 +46,37 @@ GPUhdi() bool isBetter(const o2::its::TrackITSInternal<NLayers>& a, const o2::it
   return isBetter(a.getNumberOfClusters(), a.getChi2(), b.getNumberOfClusters(), b.getChi2());
 }
 
+GPUhdi() int getEffectiveTrackLength(LayerMask hitLayerMask, LayerMask nonCountingLayerMask)
+{
+  if (hitLayerMask.empty()) {
+    return 0;
+  }
+  return hitLayerMask.length() - (LayerMask::span(hitLayerMask.first(), hitLayerMask.last()) & nonCountingLayerMask).count();
+}
+
+GPUhdi() LayerMask getEffectiveHoleMask(LayerMask hitLayerMask, LayerMask nonCountingLayerMask)
+{
+  return hitLayerMask.holeMask() & ~nonCountingLayerMask;
+}
+
+template <int NLayers>
+struct TrackSeedSelector {
+  float maxQ2Pt;
+  float maxChi2;
+  int maxHoles;
+  int minTrackLength;
+  LayerMask holeLayerMask;
+  LayerMask nonCountingLayerMask;
+
+  GPUhd() bool operator()(const TrackSeed<NLayers>& seed) const
+  {
+    const auto hitLayerMask = seed.getHitLayerMask();
+    return !(seed.getQ2Pt() > maxQ2Pt || seed.getChi2() > maxChi2) &&
+           getEffectiveTrackLength(hitLayerMask, nonCountingLayerMask) >= minTrackLength &&
+           getEffectiveHoleMask(hitLayerMask, nonCountingLayerMask).isAllowedHoleMask(maxHoles, holeLayerMask);
+  }
+};
+
 // Find the populated interior layer closest to the radial midpoint.
 // If no layer can be found, return constants::UnusedIndex.
 // Should minimize the sagitta bias.

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

@@ -104,15 +104,24 @@ class TrackingTopology
 #endif
   };
 
-  void init(int maxLayers, int maxHoles, Mask holeLayerMask)
+  void init(int maxLayers, int maxHoles, Mask holeLayerMask, Mask seedingLayerMask = 0)
   {
     clear();
     mMaxLayers = o2::gpu::CAMath::Max(0, o2::gpu::CAMath::Min(maxLayers, NLayers));
     mMaxHoles = o2::gpu::CAMath::Max(maxHoles, 0);
     mHoleLayerMask = holeLayerMask;
+    mSeedingLayerMask = seedingLayerMask.empty() ? Mask::span(0, mMaxLayers - 1) : (seedingLayerMask & Mask::span(0, mMaxLayers - 1));
+#ifndef GPUCA_GPUCODE
+    if (mSeedingLayerMask.count() < constants::ClustersPerCell) {
+      LOGP(fatal, "Tracking topology has {} seeding layers, but at least {} are required to build CA cells", mSeedingLayerMask.count(), constants::ClustersPerCell);
+    }
+#endif
     for (int fromLayer = 0; fromLayer < mMaxLayers; ++fromLayer) {
+      if (!mSeedingLayerMask.has(fromLayer)) {
+        continue;
+      }
       for (int toLayer = fromLayer + 1; toLayer < mMaxLayers; ++toLayer) {
-        if (Mask::skipped(fromLayer, toLayer).isAllowedHoleMask(mMaxHoles, mHoleLayerMask)) {
+        if (mSeedingLayerMask.has(toLayer) && isAllowedSeedingLink(fromLayer, toLayer)) {
           mLinks[mNLinks++] = LayerLink{static_cast<Id>(fromLayer), static_cast<Id>(toLayer)};
         }
       }
@@ -126,7 +135,7 @@ class TrackingTopology
           continue;
         }
         const Mask hitMask{first.fromLayer, first.toLayer, second.toLayer};
-        if (hitMask.isAllowed(mMaxHoles, mHoleLayerMask)) {
+        if ((hitMask.holeMask() & mSeedingLayerMask).isAllowedHoleMask(mMaxHoles, mHoleLayerMask)) {
           mCells[mNCells++] = CellTopology{firstId, secondId, hitMask};
         }
       }
@@ -202,9 +211,15 @@ class TrackingTopology
     mNCellsByFirstLink = offset;
   }
 
+  bool isAllowedSeedingLink(int fromLayer, int toLayer) const noexcept
+  {
+    return (Mask::skipped(fromLayer, toLayer) & mSeedingLayerMask).isAllowedHoleMask(mMaxHoles, mHoleLayerMask);
+  }
+
   int mMaxLayers{0};
   int mMaxHoles{0};
   Mask mHoleLayerMask{0};
+  Mask mSeedingLayerMask{0};
   Id mNLinks{0};
   Id mNCells{0};
   Id mNCellsByFirstLink{0};

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

@@ -59,6 +59,12 @@ std::string TrackingParameters::asString() const
   if (MaxHoles) {
     str += std::format(" MaxHoles:{} HoleMask:{}", MaxHoles, HoleLayerMask.asString());
   }
+  if (!InactiveLayerMask.empty()) {
+    str += std::format(" InactiveMask:{}", InactiveLayerMask.asString());
+  }
+  if (!SeedingLayers.empty()) {
+    str += std::format(" SeedingLayers:{}", SeedingLayers.asString());
+  }
   if (PassFlags[IterationStep::TrackFollowerTop] || PassFlags[IterationStep::TrackFollowerBot]) {
     const bool top = PassFlags[IterationStep::TrackFollowerTop], bot = PassFlags[IterationStep::TrackFollowerBot];
     str += std::format(" TrackFollower:{} NSigmaZ/Phi:{:.2f}/{:.2f}",

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

@@ -249,16 +249,25 @@ void TimeFrame<NLayers>::initVertexingTopology(const TrackingParameters& trkPara
 template <int NLayers>
 void TimeFrame<NLayers>::initDefaultTrackingTopology(const TrackingParameters& trkParam, const int maxLayers)
 {
-  mDefaultTrackingTopology.init(maxLayers, trkParam.MaxHoles, trkParam.HoleLayerMask);
+  if (maxLayers < trkParam.NLayers) {
+    LOGP(fatal, "Default tracking topology limited to {} layers, but the tracking parameters expect {}", maxLayers, trkParam.NLayers);
+  }
+  mDefaultTrackingTopology.init(trkParam.NLayers, trkParam.MaxHoles, trkParam.HoleLayerMask, trkParam.getSeedingLayerMask());
 }
 
 template <int NLayers>
 void TimeFrame<NLayers>::initTrackerTopologies(gsl::span<const TrackingParameters> trkParams, const int maxLayers)
 {
   mTrackerTopologies.resize(trkParams.size());
   for (size_t iteration = 0; iteration < trkParams.size(); ++iteration) {
-    const int iterationMaxLayers = std::min(maxLayers, trkParams[iteration].NLayers);
-    mTrackerTopologies[iteration].init(iterationMaxLayers, trkParams[iteration].MaxHoles, trkParams[iteration].HoleLayerMask);
+    if (maxLayers < trkParams[iteration].NLayers) {
+      LOGP(fatal, "Iteration {}: tracking topology limited to {} layers, but the tracking parameters expect {}", iteration, maxLayers, trkParams[iteration].NLayers);
+    }
+    const int nActiveLayers = trkParams[iteration].getActiveLayerMask().count();
+    if (trkParams[iteration].MinTrackLength > nActiveLayers) {
+      LOGP(fatal, "Iteration {}: MinTrackLength {} cannot be satisfied with {} active layers", iteration, trkParams[iteration].MinTrackLength, nActiveLayers);
+    }
+    mTrackerTopologies[iteration].init(trkParams[iteration].NLayers, trkParams[iteration].MaxHoles, trkParams[iteration].HoleLayerMask, trkParams[iteration].getSeedingLayerMask());
   }
 }
 

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

@@ -759,13 +759,11 @@ void TrackerTraits<NLayers>::findRoads(const int iteration)
     unsortedClusters[iLayer] = mTimeFrame->getUnsortedClusters()[iLayer].data();
   }
   const auto topology = mTimeFrame->getTrackingTopologyView();
+  const auto notSeedingLayers = mTrkParams[iteration].getNotSeedingLayerMask();
+  const auto minSeedingClusters = mTrkParams[iteration].getMinSeedingClusters();
   for (int startLevel{mTrkParams[iteration].CellsPerRoad()}; startLevel >= mTrkParams[iteration].CellMinimumLevel(); --startLevel) {
 
-    auto seedFilter = [&](const auto& seed) {
-      return seed.getHitLayerMask().isAllowed(mTrkParams[iteration].MaxHoles, mTrkParams[iteration].HoleLayerMask) &&
-             seed.getHitLayerMask().length() >= mTrkParams[iteration].MinTrackLength &&
-             seed.getQ2Pt() <= 1.e3 && seed.getChi2() <= mTrkParams[iteration].MaxChi2NDF * ((startLevel + 2) * 2 - 5);
-    };
+    const track::TrackSeedSelector<NLayers> seedFilter{1.e3f, mTrkParams[iteration].MaxChi2NDF * ((startLevel + 2) * 2 - 5), mTrkParams[iteration].MaxHoles, minSeedingClusters, mTrkParams[iteration].HoleLayerMask, notSeedingLayers};
 
     bounded_vector<TrackSeedN> trackSeeds(mMemoryPool.get());
     for (int startCellTopologyId{0}; startCellTopologyId < topology.nCells; ++startCellTopologyId) {
@@ -855,14 +853,18 @@ void TrackerTraits<NLayers>::acceptTracks(int iteration,
   auto& trks = mTimeFrame->getTracks();
   trks.reserve(trks.size() + tracks.size());
   const float smallestROFHalf = mTimeFrame->getROFOverlapTableView().getClockLayer().mROFLength * 0.5f;
+  const int minTrackLength = mTrkParams[iteration].MinTrackLength;
+  const LayerMask inactiveLayerMask = mTrkParams[iteration].InactiveLayerMask;
   for (auto& track : tracks) {
     int nShared = 0;
     bool isFirstShared{false};
     int firstLayer{-1}, firstCluster{-1};
+    LayerMask hitLayerMask{0};
     for (int iLayer{0}; iLayer < mTrkParams[iteration].NLayers; ++iLayer) {
       if (track.getClusterIndex(iLayer) == constants::UnusedIndex) {
         continue;
       }
+      hitLayerMask.set(iLayer);
       bool isShared = mTimeFrame->isClusterUsed(iLayer, track.getClusterIndex(iLayer));
       nShared += int(isShared);
       if (firstLayer < 0) {
@@ -872,6 +874,11 @@ void TrackerTraits<NLayers>::acceptTracks(int iteration,
       }
     }
 
+    /// seeds are selected with a length cut relaxed to the seeding layers: enforce the full minimum length before accepting the final track
+    if (track::getEffectiveTrackLength(hitLayerMask, inactiveLayerMask) < minTrackLength) {
+      continue;
+    }
+
     /// do not account for the first cluster in the shared clusters number if it is allowed
     if (nShared - int(isFirstShared && mTrkParams[iteration].AllowSharingFirstCluster) > mTrkParams[iteration].SharedMaxClusters) {
       continue;