ITS: optionally restrict CA to a subset of layers + support inactive layers (#15515)

* ITS: allow CA seeding on configurable layer subsets

Add a SeedingLayers mask to TrackingParameters and use it when building the CA tracking topology. Transitions and cells are built only from layers that participate in seeding, while skipped non-seeding layers do not contribute to the effective seed length or hole budget.

Move the track-seed selection predicate to TrackHelpers so CPU and GPU road finding use the same host/device logic.

* ITS: distinguish inactive layers from track holes

Add an InactiveLayerMask to TrackingParameters to represent detector layers that are absent for a given input or configuration. Inactive layers remain skippable, but they do not inflate the effective seed length or consume the ordinary hole budget.

Guard the multiple-scattering helper against zero material so inactive layers can be configured with no material budget.

* ITS: make sure that seeds are not cut by min track lenght

* ITS: implement Felix' comments

* ITS: second round of comments

Author: mpuccio

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 nonSeedingLayerMask,
                               const std::vector<float>& layerxX0Host,
                               const o2::base::Propagator* propagator,
                               const o2::base::PropagatorF::MatCorrType matCorrType,

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

@@ -331,8 +331,9 @@ void TrackerTraitsGPU<NLayers>::findRoads(const int iteration)
                                         this->mTrkParams[iteration].MaxChi2ClusterAttachment,
                                         this->mTrkParams[iteration].MaxChi2NDF,
                                         this->mTrkParams[iteration].MaxHoles,
-                                        this->mTrkParams[iteration].MinTrackLength,
+                                        this->mTrkParams[iteration].getMinSeedingClusters(),
                                         this->mTrkParams[iteration].HoleLayerMask,
+                                        this->mTrkParams[iteration].getNonSeedingLayerMask(),
                                         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 nonSeedingLayerMask,
                               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>{constants::MaxTrackSeedQ2Pt, maxChi2NDF, startLevel, maxHoles, minSeedingClusters, holeLayerMask, nonSeedingLayerMask});
   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 nonSeedingLayerMask,
                                           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 nonSeedingLayerMask,
                                            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 getNonSeedingLayerMask() 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/Constants.h

@@ -29,6 +29,7 @@ constexpr bool DoTimeBenchmarks = true;
 constexpr bool SaveTimeBenchmarks = false;
 constexpr float Tolerance = 1e-12;                  // numerical tolerance
 constexpr int ClustersPerCell = 3;                  // number of clusters for a cell
+constexpr float MaxTrackSeedQ2Pt = 1.e3f;           // maximum q/pt for track seeds
 constexpr int UnusedIndex = -1;                     // global unused flag
 constexpr float UnsetValue = -999.f;                // global unset value
 constexpr float Radl = 9.36f;                       // Radiation length of Si [cm]

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"
@@ -34,16 +35,46 @@ GPUhdi() bool isBetter(const int nClustersA, const float chi2A, const int nClust
   return (nClustersA > nClustersB) || (nClustersA == nClustersB && chi2A < chi2B);
 }
 
-GPUhdi() bool isBetter(const o2::its::TrackITS& a, const o2::its::TrackITS& b)
+GPUhdi() bool isBetter(const auto& a, const auto& b)
 {
   return isBetter(a.getNumberOfClusters(), a.getChi2(), b.getNumberOfClusters(), b.getChi2());
 }
 
 template <int NLayers>
-GPUhdi() bool isBetter(const o2::its::TrackITSInternal<NLayers>& a, const o2::its::TrackITSInternal<NLayers>& b)
-{
-  return isBetter(a.getNumberOfClusters(), a.getChi2(), b.getNumberOfClusters(), b.getChi2());
-}
+struct TrackSeedSelector {
+  float maxQ2Pt;
+  float maxChi2;
+  int maxHoles;
+  int minTrackLength;
+  LayerMask holeLayerMask;
+  LayerMask nonSeedingLayerMask;
+
+  GPUhd() TrackSeedSelector(float maxQ2Pt, float maxChi2NDF, int startLevel, int maxHoles, int minTrackLength, LayerMask holeLayerMask, LayerMask nonSeedingLayerMask)
+    : maxQ2Pt{maxQ2Pt}, maxChi2{maxChi2NDF * ((startLevel + 2) * 2 - 5)}, maxHoles{maxHoles}, minTrackLength{minTrackLength}, holeLayerMask{holeLayerMask}, nonSeedingLayerMask{nonSeedingLayerMask}
+  {
+  }
+
+  static GPUhdi() int getEffectiveTrackLength(LayerMask hitLayerMask, LayerMask excludedLayerMask)
+  {
+    if (hitLayerMask.empty()) {
+      return 0;
+    }
+    return hitLayerMask.length() - (LayerMask::span(hitLayerMask.first(), hitLayerMask.last()) & excludedLayerMask).count();
+  }
+
+  static GPUhdi() LayerMask getEffectiveHoleMask(LayerMask hitLayerMask, LayerMask excludedLayerMask)
+  {
+    return hitLayerMask.holeMask() & ~excludedLayerMask;
+  }
+
+  GPUhd() bool operator()(const TrackSeed<NLayers>& seed) const
+  {
+    const auto hitLayerMask = seed.getHitLayerMask();
+    return !(seed.getQ2Pt() > maxQ2Pt || seed.getChi2() > maxChi2) &&
+           getEffectiveTrackLength(hitLayerMask, nonSeedingLayerMask) >= minTrackLength &&
+           getEffectiveHoleMask(hitLayerMask, nonSeedingLayerMask).isAllowedHoleMask(maxHoles, holeLayerMask);
+  }
+};
 
 // Find the populated interior layer closest to the radial midpoint.
 // If no layer can be found, return constants::UnusedIndex.

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

@@ -70,6 +70,7 @@ class TrackingTopology
     const CellTopology* cells{nullptr};
     const Range* cellsByFirstLinkIndex{nullptr};
     const Id* cellsByFirstLink{nullptr};
+    Mask seedingLayerMask{0};
     Id nLinks{0};
     Id nCells{0};
     Id nCellsByFirstLink{0};
@@ -81,7 +82,7 @@ class TrackingTopology
 #ifndef GPUCA_GPUCODE
     std::string asString() const
     {
-      std::string out = fmt::format("TrackingTopology: links={} cells={}", nLinks, nCells);
+      std::string out = fmt::format("TrackingTopology: links={} cells={} seedingLayers={}", nLinks, nCells, seedingLayerMask.asString());
       out += "\n  links:";
       for (Id linkId = 0; linkId < nLinks; ++linkId) {
         const auto& t = links[linkId];
@@ -104,15 +105,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 +136,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};
         }
       }
@@ -141,6 +151,7 @@ class TrackingTopology
                 mCells.data(),
                 mCellsByFirstLinkIndex.data(),
                 mCellsByFirstLink.data(),
+                mSeedingLayerMask,
                 mNLinks,
                 mNCells,
                 mNCellsByFirstLink};
@@ -155,6 +166,7 @@ class TrackingTopology
                 deviceCells,
                 deviceCellsByFirstLinkIndex,
                 deviceCellsByFirstLink,
+                mSeedingLayerMask,
                 mNLinks,
                 mNCells,
                 mNCellsByFirstLink};
@@ -202,9 +214,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());
   }
 }