ITS: fix time assignment

Author: f3sch

DataFormats/Detectors/ITSMFT/ITS/include/DataFormatsITS/TrackITS.h

@@ -37,7 +37,7 @@ namespace its
 class TrackITS : public o2::track::TrackParCov
 {
   enum UserBits {
-    kSharedClusters = 1 << 29
+    kSharedClusters = 1 << 28
   };
 
   using Cluster = o2::itsmft::Cluster;
@@ -93,7 +93,7 @@ class TrackITS : public o2::track::TrackParCov
 
   GPUhdi() void setChi2(float chi2) { mChi2 = chi2; }
 
-  bool isBetter(const TrackITS& best, float maxChi2) const;
+  bool isBetter(const TrackITS& best, float maxChi2 = o2::constants::math::VeryBig) const;
 
   auto& getTimeStamp() { return mTime; }
   const auto& getTimeStamp() const { return mTime; }

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

@@ -116,7 +116,7 @@ class TrackerTraits
 
   std::shared_ptr<BoundedMemoryResource> mMemoryPool;
   std::shared_ptr<tbb::task_arena> mTaskArena;
-  dbscan::DBSCAN mDBScan;
+  dbscan::DBSCAN mDBScan{};
 
  protected:
   o2::gpu::GPUChainITS* mChain = nullptr;

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

@@ -0,0 +1,91 @@
+// 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 <vector>
+#include <array>
+#include <algorithm>
+#include <limits>
+#include <cstdint>
+
+namespace o2::its::utils
+{
+
+using Bracket = std::pair<int32_t, int32_t>;
+constexpr Bracket InvalidBracket{std::numeric_limits<int32_t>::min(), std::numeric_limits<int32_t>::min()};
+
+template <size_t N>
+Bracket computeSmallestBracket(const std::array<Bracket, N>& intervals)
+{
+  // Collect valid intervals
+  std::vector<Bracket> validIntervals;
+  for (const auto& inter : intervals) {
+    if (inter != InvalidBracket) {
+      validIntervals.push_back(inter);
+    }
+  }
+
+  if (validIntervals.empty()) {
+    return InvalidBracket;
+  }
+
+  // Check for gaps: sort by start, verify each overlaps/touches the next
+  std::sort(validIntervals.begin(), validIntervals.end());
+  for (size_t i = 1; i < validIntervals.size(); ++i) {
+    if (validIntervals[i].first > validIntervals[i - 1].second) {
+      // Gap found: intervals are disconnected
+      return InvalidBracket;
+    }
+  }
+
+  // Collect unique endpoints
+  std::array<int32_t, 2 * N> ends;
+  size_t count = 0;
+  for (const auto& inter : validIntervals) {
+    ends[count++] = inter.first;
+    ends[count++] = inter.second;
+  }
+
+  std::sort(ends.begin(), ends.begin() + count);
+  auto last = std::unique(ends.begin(), ends.begin() + count);
+  count = last - ends.begin();
+
+  int32_t bestLength = std::numeric_limits<int32_t>::max();
+  Bracket bestBracket{InvalidBracket};
+
+  for (size_t i = 0; i < count; ++i) {
+    const int32_t L = ends[i];
+    for (size_t j = i + 1; j < count; ++j) {
+      const int32_t R = ends[j];
+      int32_t length = R - L;
+      if (length >= bestLength) {
+        break;
+      }
+
+      bool overlaps = true;
+      for (const auto& inter : validIntervals) {
+        if (L >= inter.second || R <= inter.first) {
+          overlaps = false;
+          break;
+        }
+      }
+
+      if (overlaps) {
+        bestLength = length;
+        bestBracket = {L, R};
+        break;
+      }
+    }
+  }
+
+  return bestBracket;
+}
+
+} // namespace o2::its::utils

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

@@ -24,7 +24,7 @@ using namespace o2::its;
 
 std::string TrackingParameters::asString() const
 {
-  std::string str = std::format("NZb:{} NPhB:{} PerVtx:{} DropFail:{} ClSh:{} TrklMinPt:{:.2f} MinCl:{}", ZBins, PhiBins, PerPrimaryVertexProcessing, DropTFUponFailure, ClusterSharing, TrackletMinPt, MinTrackLength);
+  std::string str = std::format("NZb:{} NPhB:{} PerVtx:{} DropFail:{} ClSh:{} TrklMinPt:{:.2f} MinCl:{} NSigma:{}", ZBins, PhiBins, PerPrimaryVertexProcessing, DropTFUponFailure, ClusterSharing, TrackletMinPt, MinTrackLength, NSigmaCut);
   bool first = true;
   for (int il = NLayers; il >= MinTrackLength; il--) {
     int slot = NLayers - il;

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

@@ -18,6 +18,7 @@
 #include <oneapi/tbb/parallel_sort.h>
 
 #include <algorithm>
+#include <cstddef>
 #include <iterator>
 #include <limits>
 #include <ranges>
@@ -40,6 +41,7 @@
 #include "ITStracking/BoundedAllocator.h"
 #include "ITStracking/IndexTableUtils.h"
 #include "ITStracking/Tracklet.h"
+#include "ITStracking/Utils.h"
 #include "ReconstructionDataFormats/Track.h"
 
 /// optimization output
@@ -58,7 +60,7 @@ namespace o2::its
 {
 namespace
 {
-utils::TreeStreamRedirector* sDBGOut{nullptr};
+o2::utils::TreeStreamRedirector* sDBGOut{nullptr};
 }
 
 struct PassMode {
@@ -570,8 +572,9 @@ void TrackerTraits<NLayers>::findCellSeeds(const int iteration)
         int rof = mTimeFrame->getClusterROF(i, cls[i]);
         int rofStartBC = mTimeFrame->getROFOverlapTableView().getLayer(i).getROFStartInBC(rof);
         int rofEndBC = mTimeFrame->getROFOverlapTableView().getLayer(i).getROFEndInBC(rof);
-        startBC = o2::gpu::CAMath::Min(startBC, rofStartBC);
-        endBC = o2::gpu::CAMath::Max(endBC, rofEndBC);
+        // if the start/end of the cluster is after/before the current bracket need to enlarge
+        startBC = (rofStartBC <= startBC) ? rofStartBC : o2::gpu::CAMath::Max(startBC, rofStartBC);
+        endBC = (rofEndBC >= endBC) ? rofEndBC : o2::gpu::CAMath::Min(endBC, rofEndBC);
       }
       if (endBC - startBC < 0) { // this should not happen
         ltracks[iCell].markDead();
@@ -946,10 +949,10 @@ void TrackerTraits<NLayers>::processNeighbours(int iteration, int iLayer, int iL
         if (neighbourCell.getSecondTrackletIndex() != currentCell.getFirstTrackletIndex()) {
           continue;
         }
-        if (mTimeFrame->isClusterUsed(iLayer - 1, neighbourCell.getFirstClusterIndex())) {
+        if (currentCell.getLevel() - 1 != neighbourCell.getLevel()) {
           continue;
         }
-        if (currentCell.getLevel() - 1 != neighbourCell.getLevel()) {
+        if (mTimeFrame->isClusterUsed(iLayer - 1, neighbourCell.getFirstClusterIndex())) {
           continue;
         }
 
@@ -1143,15 +1146,17 @@ void TrackerTraits<NLayers>::findRoads(const int iteration)
       }
 
       deepVectorClear(trackSeeds);
+      // sort tracks in quality (accounting for 1. length; 2. chi2)
+      // needed since then tracks with shared clusters can be marked/discarded
       tbb::parallel_sort(tracks.begin(), tracks.end(), [](const auto& a, const auto& b) {
-        return a.getChi2() < b.getChi2();
+        return a.isBetter(b);
       });
     });
 
     for (auto& track : tracks) {
       int nShared = 0;
       bool isFirstShared{false};
-      int firstLayer{-1}, firstCluster{-1};
+      int firstLayer{constants::UnusedIndex}, firstCluster{constants::UnusedIndex};
       for (int iLayer{0}; iLayer < mRecoParams[iteration].params.NLayers; ++iLayer) {
         if (track.getClusterIndex(iLayer) == constants::UnusedIndex) {
           continue;
@@ -1172,20 +1177,32 @@ void TrackerTraits<NLayers>::findRoads(const int iteration)
 
       // here we can do the calculation of the time bracket simply
       // by checkig in which rofs the clusters are
-      int bcStart{0}, bcEnd{std::numeric_limits<int>::max()};
+      std::array<utils::Bracket, NLayers> brackets;
+      brackets.fill(utils::InvalidBracket);
       for (int iLayer{0}; iLayer < mRecoParams[iteration].params.NLayers; ++iLayer) {
         if (track.getClusterIndex(iLayer) == constants::UnusedIndex) {
           continue;
         }
-        mTimeFrame->markUsedCluster(iLayer, track.getClusterIndex(iLayer));
         int currentROF = mTimeFrame->getClusterROF(iLayer, track.getClusterIndex(iLayer));
-        int bcClsSta = mTimeFrame->getROFOverlapTableView().getLayer(iLayer).getROFStartInBC(currentROF);
-        int bcClsEnd = mTimeFrame->getROFOverlapTableView().getLayer(iLayer).getROFEndInBC(currentROF);
-        bcStart = std::max(bcStart, bcClsSta);
-        bcEnd = std::min(bcEnd, bcClsEnd);
+        // need to account for the imposed delay
+        int bcClsSta = mTimeFrame->getROFOverlapTableView().getLayer(iLayer).getROFStartInBC(currentROF, true);
+        int bcClsEnd = mTimeFrame->getROFOverlapTableView().getLayer(iLayer).getROFEndInBC(currentROF, true);
+        brackets[iLayer] = utils::Bracket{bcClsSta, bcClsEnd};
+        LOGP(debug, "\tlay:{} bcClsSta={} bcClsEnd={}", iLayer, bcClsSta, bcClsEnd);
+      }
+      const auto best = utils::computeSmallestBracket(brackets);
+      if (best == utils::InvalidBracket) {
+        continue; // track has an impossible span, discard
+      }
+      // mark used clusters for the next iteration
+      for (int iLayer{0}; iLayer < mRecoParams[iteration].params.NLayers; ++iLayer) {
+        if (track.getClusterIndex(iLayer) == constants::UnusedIndex) {
+          continue;
+        }
+        mTimeFrame->markUsedCluster(iLayer, track.getClusterIndex(iLayer));
       }
-      track.getTimeStamp().setTimeStamp(bcStart);
-      track.getTimeStamp().setTimeStampError(bcEnd - bcStart + 1);
+      track.getTimeStamp().setTimeStamp(best.first);
+      track.getTimeStamp().setTimeStampError(best.second - best.first + 1);
       track.setUserField(0);
       track.getParamOut().setUserField(0);
       mTimeFrame->getTracks().emplace_back(track);

Detectors/ITSMFT/ITS/tracking/test/CMakeLists.txt

@@ -20,3 +20,9 @@ o2_add_test(roflookuptables
             COMPONENT_NAME its-tracking
             LABELS "its;tracking"
             PUBLIC_LINK_LIBRARIES O2::ITStracking)
+
+o2_add_test(utils
+            SOURCES testUtils.cxx
+            COMPONENT_NAME its-tracking
+            LABELS "its;tracking"
+            PUBLIC_LINK_LIBRARIES O2::ITStracking)

Detectors/ITSMFT/ITS/tracking/test/testUtils.cxx

@@ -0,0 +1,91 @@
+// 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 <boost/test/tools/old/interface.hpp>
+#define BOOST_TEST_MODULE ITS Utils
+#define BOOST_TEST_MAIN
+#define BOOST_TEST_DYN_LINK
+
+#include <boost/test/unit_test.hpp>
+
+#include "ITStracking/Utils.h"
+using namespace o2::its::utils;
+
+BOOST_AUTO_TEST_CASE(test_1)
+{
+  const std::array<Bracket, 2> brackets{{{594, 1188},
+                                         {297, 891}}};
+  const auto result = computeSmallestBracket(brackets);
+  BOOST_CHECK_EQUAL(result.first, 594);
+  BOOST_CHECK_EQUAL(result.second, 891);
+}
+
+BOOST_AUTO_TEST_CASE(test_2)
+{
+  const std::array<Bracket, 7> brackets{{{594, 1188},
+                                         {594, 1188},
+                                         {594, 1188},
+                                         {297, 891},
+                                         {297, 891},
+                                         {297, 891},
+                                         {891, 1485}}};
+  const auto result = computeSmallestBracket(brackets);
+  BOOST_CHECK_EQUAL(result.first, 594);
+  BOOST_CHECK_EQUAL(result.second, 1188);
+}
+
+BOOST_AUTO_TEST_CASE(test_3)
+{
+  std::array<Bracket, 7> brackets;
+  brackets.fill(InvalidBracket);
+  brackets[0] = {500, 600};
+  brackets[3] = {500, 600};
+  brackets[6] = {500, 600};
+  const auto result = computeSmallestBracket(brackets);
+  BOOST_CHECK_EQUAL(result.first, 500);
+  BOOST_CHECK_EQUAL(result.second, 600);
+}
+
+BOOST_AUTO_TEST_CASE(test_4)
+{
+  std::array<Bracket, 3> brackets;
+  brackets.fill(InvalidBracket);
+  brackets[0] = {100, 200};
+  brackets[2] = {500, 600};
+  const auto result = computeSmallestBracket(brackets);
+  BOOST_CHECK(result == InvalidBracket);
+}
+
+BOOST_AUTO_TEST_CASE(test_5)
+{
+  std::array<Bracket, 7> brackets{{{5346, 5940},
+                                   {5346, 5940},
+                                   {5346, 5940},
+                                   {5049, 5643},
+                                   {4455, 5049},
+                                   {5049, 5643},
+                                   {5049, 5643}}};
+  const auto result = computeSmallestBracket(brackets);
+  BOOST_CHECK(result != InvalidBracket);
+}
+
+BOOST_AUTO_TEST_CASE(test_6)
+{
+  std::array<Bracket, 7> brackets{{{46926, 47520},
+                                   {46926, 47520},
+                                   {46926, 47520},
+                                   {47817, 48411},
+                                   {47817, 48411},
+                                   {47817, 48411},
+                                   {47817, 48411}}};
+  const auto result = computeSmallestBracket(brackets);
+  BOOST_CHECK(result == InvalidBracket);
+}