Merge 9be6fb1 into sapling-pr-archive-ktf

Author: ktf

Common/Constants/include/CommonConstants/PhysicsConstants.h

@@ -92,7 +92,9 @@ enum Pdg {
   kHyperHelium4 = 1010020040,
   kHyperHelium5 = 1010020050,
   kHyperHelium4Sigma = 1110020040,
-  kLambda1520_Py = 102134
+  kLambda1520_Py = 102134,
+  kK1_1270_0 = 10313,
+  kK1_1270Plus = 10323
 };
 
 /// \brief Declarations of masses for additional particles
@@ -158,6 +160,8 @@ constexpr double MassHyperHelium4 = 3.921728;
 constexpr double MassHyperHelium5 = 4.839961;
 constexpr double MassHyperHelium4Sigma = 3.995;
 constexpr double MassLambda1520_Py = 1.5195;
+constexpr double MassK1_1270_0 = 1.253;
+constexpr double MassK1_1270Plus = 1.272;
 
 /// \brief Declarations of masses for particles in ROOT PDG_t
 constexpr double MassDown = 0.00467;

Common/Constants/include/CommonConstants/make_pdg_header.py

@@ -151,6 +151,8 @@ class Pdg(Enum):
     kHyperHelium5 = 1010020050
     kHyperHelium4Sigma = 1110020040
     kLambda1520_Py = 102134  # PYTHIA code different from PDG
+    kK1_1270_0 = 10313
+    kK1_1270Plus = 10323
 
 
 dbPdg = o2.O2DatabasePDG

Detectors/GlobalTracking/src/MatchTPCITS.cxx

@@ -880,9 +880,13 @@ void MatchTPCITS::doMatching(int sec)
   int idxMinTPC = timeStartTPC[minROFITS];                             // index of 1st cached TPC track within cached ITS ROFrames
   auto t2nbs = tpcTimeBin2MUS(mZ2TPCBin * mParams->tpcTimeICMatchingNSigma);
   bool checkInteractionCandidates = mUseFT0 && mParams->validateMatchByFIT != MatchTPCITSParams::Disable;
+  LOGP(alarm, "doMatching sec={} nTPC={} nITS={} idxMinTPC={} mTPCWork.size={} mITSWork.size={}", sec, nTracksTPC, nTracksITS, idxMinTPC, mTPCWork.size(), mITSWork.size());
 
   int itsROBin = 0;
   for (int itpc = idxMinTPC; itpc < nTracksTPC; itpc++) {
+    if (cacheTPC[itpc] >= (int)mTPCWork.size()) {
+      LOGP(alarm, "OOB TPC: itpc={} cacheTPC[itpc]={} mTPCWork.size={}", itpc, cacheTPC[itpc], mTPCWork.size());
+    }
     auto& trefTPC = mTPCWork[cacheTPC[itpc]];
     // estimate ITS 1st ROframe bin this track may match to: TPC track are sorted according to their
     // timeMax, hence the timeMax - MaxmNTPCBinsFullDrift are non-decreasing
@@ -895,6 +899,9 @@ void MatchTPCITS::doMatching(int sec)
     int iits0 = timeStartITS[itsROBin];
     nCheckTPCControl++;
     for (auto iits = iits0; iits < nTracksITS; iits++) {
+      if (cacheITS[iits] >= (int)mITSWork.size()) {
+        LOGP(alarm, "OOB ITS: iits={} cacheITS[iits]={} mITSWork.size={}", iits, cacheITS[iits], mITSWork.size());
+      }
       auto& trefITS = mITSWork[cacheITS[iits]];
       // compare if the ITS and TPC tracks may overlap in time
       LOG(debug) << "TPC bracket: " << trefTPC.tBracket.asString() << " ITS bracket: " << trefITS.tBracket.asString() << " TPCtgl: " << trefTPC.getTgl() << " ITStgl: " << trefITS.getTgl();
@@ -934,6 +941,9 @@ void MatchTPCITS::doMatching(int sec)
       }
 
       nCheckITSControl++;
+      if (nCheckITSControl <= 3) {
+        LOGP(alarm, "doMatching sec={} itpc={} iits={} nCheckITS={}", sec, itpc, iits, nCheckITSControl);
+      }
       float chi2 = -1;
       int rejFlag = compareTPCITSTracks(trefITS, trefTPC, chi2);
 
@@ -1006,6 +1016,7 @@ void MatchTPCITS::doMatching(int sec)
               << " N TPC tracks checked: " << nCheckTPCControl << " (starting from " << idxMinTPC
               << "), checks: " << nCheckITSControl << ", matches:" << nMatchesControl;
   }
+  LOGP(alarm, "doMatching sec={} done: nCheckTPC={} nCheckITS={} nMatches={}", sec, nCheckTPCControl, nCheckITSControl, nMatchesControl);
   mNMatchesControl += nMatchesControl;
 }
 

Detectors/ITSMFT/ITS/macros/test/CheckStaggering.C

@@ -80,9 +80,9 @@ void CheckStaggering(int runNumber, int max = -1, const std::string& dir = "")
   auto& ccdbmgr = o2::ccdb::BasicCCDBManager::instance();
   ccdbmgr.setURL("https://alice-ccdb.cern.ch");
   auto runDuration = ccdbmgr.getRunDuration(runNumber);
-  auto tRun = runDuration.first + (runDuration.second - runDuration.first) / 2; // time stamp for the middle of the run duration
+  auto tRun = runDuration.first + ((runDuration.second - runDuration.first) / 2); // time stamp for the middle of the run duration
   ccdbmgr.setTimestamp(tRun);
-  printf("Run %d has TS %lld", runNumber, tRun);
+  printf("Run %d has TS %ld", runNumber, tRun);
   auto geoAligned = ccdbmgr.get<TGeoManager>("GLO/Config/GeometryAligned");
   auto magField = ccdbmgr.get<o2::parameters::GRPMagField>("GLO/Config/GRPMagField");
   auto grpLHC = ccdbmgr.get<o2::parameters::GRPLHCIFData>("GLO/Config/GRPLHCIF");
@@ -134,7 +134,7 @@ void CheckStaggering(int runNumber, int max = -1, const std::string& dir = "")
   auto hVtxZ = new TH1F("hVtxZ", "seeding vertices Z", 200, -16, 16);
   auto hVtxNCont = new TH1F("hVtxNCont", "seeding vertices contributors", 100, 0, 100);
   auto hVtxZNCont = new TProfile("hVtxZNCont", "seeding vertices z-contributors", 200, -16, 16);
-  auto hVtxCls = new TProfile("hVtxCls", ";Cls/TF;Cls/Vtx", 400, 20000, 60000);
+  auto hVtxCls = new TProfile("hVtxCls", ";Cls/TF;Cls/Vtx", 2000, 600000, 900000);
   auto hVtxTS = new TH1D("hVtxTS", "vtx time t0;t0 (BC)", o2::constants::lhc::LHCMaxBunches, 0, o2::constants::lhc::LHCMaxBunches);
 
   const float minVtxWeight{5};
@@ -183,8 +183,12 @@ void CheckStaggering(int runNumber, int max = -1, const std::string& dir = "")
 
     tTrks->SetBranchAddress("ITSTrack", &trkArrPtr);
     tTrks->SetBranchAddress("Vertices", &vtxArrPtr);
-    for (int i{0}; i < 7; ++i) {
-      tCls->SetBranchAddress(Form("ITSClusterComp_%d", i), &clsArr[i]);
+    if (tCls->GetBranchStatus("ITSClusterComp")) {
+      tCls->SetBranchAddress("ITSClusterComp", &clsArr[0]);
+    } else {
+      for (int i{0}; i < 7; ++i) {
+        tCls->SetBranchAddress(Form("ITSClusterComp_%d", i), &clsArr[i]);
+      }
     }
 
     for (int iTF{0}; tTrks->LoadTree(iTF) >= 0; ++iTF) {
@@ -193,7 +197,9 @@ void CheckStaggering(int runNumber, int max = -1, const std::string& dir = "")
 
       size_t ncls = 0;
       for (int i{0}; i < 7; ++i) {
-        ncls += clsArr[i]->size();
+        if (clsArr[i]) {
+          ncls += clsArr[i]->size();
+        }
       }
 
       // for each TF built pool of positive and negaitve tracks

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

@@ -47,20 +47,6 @@ namespace o2::its
 namespace gpu
 {
 
-struct sort_tracklets {
-  GPUhd() bool operator()(const Tracklet& a, const Tracklet& b)
-  {
-    if (a.firstClusterIndex != b.firstClusterIndex) {
-      return a.firstClusterIndex < b.firstClusterIndex;
-    }
-    return a.secondClusterIndex < b.secondClusterIndex;
-  }
-};
-
-struct equal_tracklets {
-  GPUhd() bool operator()(const Tracklet& a, const Tracklet& b) { return a.firstClusterIndex == b.firstClusterIndex && a.secondClusterIndex == b.secondClusterIndex; }
-};
-
 template <typename T1, typename T2>
 struct sort_by_second {
   GPUhd() bool operator()(const gpuPair<T1, T2>& a, const gpuPair<T1, T2>& b) const { return a.second < b.second; }
@@ -408,7 +394,7 @@ GPUg() void __launch_bounds__(256, 1) computeLayerTrackletsMultiROFKernel(
         const float sqInverseDeltaZ0{1.f / (math_utils::Sq(currentCluster.zCoordinate - primaryVertex.getZ()) + constants::Tolerance)}; /// protecting from overflows adding the detector resolution
         const float sigmaZ{o2::gpu::CAMath::Sqrt(math_utils::Sq(resolution) * math_utils::Sq(tanLambda) * ((math_utils::Sq(inverseR0) + sqInverseDeltaZ0) * math_utils::Sq(meanDeltaR) + 1.f) + math_utils::Sq(meanDeltaR * MSAngle))};
         const int4 selectedBinsRect{o2::its::getBinsRect(currentCluster, layerIndex + 1, zAtRmin, zAtRmax, sigmaZ * NSigmaCut, phiCut, *utils)};
-        if (selectedBinsRect.x == 0 && selectedBinsRect.y == 0 && selectedBinsRect.z == 0 && selectedBinsRect.w == 0) {
+        if (selectedBinsRect.x < 0) {
           continue;
         }
         int phiBinsNum{selectedBinsRect.w - selectedBinsRect.y + 1};
@@ -687,8 +673,8 @@ void computeTrackletsInROFsHandler(const IndexTableUtils<NLayers>* utils,
     mulScatAng[layer]);
   thrust::device_ptr<Tracklet> tracklets_ptr(spanTracklets[layer]);
   auto nosync_policy = THRUST_NAMESPACE::par_nosync(gpu::TypedAllocator<char>(alloc)).on(streams[layer].get());
-  thrust::sort(nosync_policy, tracklets_ptr, tracklets_ptr + nTracklets[layer], gpu::sort_tracklets());
-  auto unique_end = thrust::unique(nosync_policy, tracklets_ptr, tracklets_ptr + nTracklets[layer], gpu::equal_tracklets());
+  thrust::sort(nosync_policy, tracklets_ptr, tracklets_ptr + nTracklets[layer]);
+  auto unique_end = thrust::unique(nosync_policy, tracklets_ptr, tracklets_ptr + nTracklets[layer]);
   nTracklets[layer] = unique_end - tracklets_ptr;
   if (layer) {
     GPUChkErrS(cudaMemsetAsync(trackletsLUTsHost[layer], 0, (nClusters[layer] + 1) * sizeof(int), streams[layer].get()));

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

@@ -56,18 +56,17 @@ class SeedBase : public o2::track::TrackParCovF
   GPUhd() const auto& clustersRaw() const { return mClusters; }
 
  private:
-  float mChi2 = -999.f;
-  int mLevel = constants::UnusedIndex;
+  float mChi2{constants::UnsetValue};
+  int mLevel{constants::UnusedIndex};
   std::array<int, 2> mTracklets = constants::helpers::initArray<int, 2, constants::UnusedIndex>();
   std::array<int, NClusters> mClusters = constants::helpers::initArray<int, NClusters, constants::UnusedIndex>();
   TimeEstBC mTime;
 };
 
 /// CellSeed: connections of three clusters
-class CellSeed final : public SeedBase<3>
+class CellSeed final : public SeedBase<constants::ClustersPerCell>
 {
-  static constexpr int NStoredClusters = 3;
-  using Base = SeedBase<NStoredClusters>;
+  using Base = SeedBase<constants::ClustersPerCell>;
 
  public:
   GPUhdDefault() CellSeed() = default;
@@ -98,7 +97,7 @@ class CellSeed final : public SeedBase<3>
   GPUhd() int getCluster(int layer) const
   {
     const int rel = layer - getInnerLayer();
-    return (rel >= 0 && rel < NStoredClusters) ? this->clustersRaw()[rel] : constants::UnusedIndex;
+    return (rel >= 0 && rel < constants::ClustersPerCell) ? this->clustersRaw()[rel] : constants::UnusedIndex;
   }
 };
 

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

@@ -19,9 +19,6 @@
 #include <array>
 #include <utility>
 
-#include "GPUCommonDef.h"
-#include "GPUCommonDefAPI.h"
-
 namespace o2::its::constants
 {
 
@@ -30,12 +27,14 @@ constexpr float MB = KB * KB;
 constexpr float GB = MB * KB;
 constexpr bool DoTimeBenchmarks = true;
 constexpr bool SaveTimeBenchmarks = false;
-constexpr float Tolerance = 1e-12; // numerical tolerance
-constexpr int ClustersPerCell = 3;
-constexpr int UnusedIndex = -1;
-constexpr float Radl = 9.36f; // Radiation length of Si [cm]
-constexpr float Rho = 2.33f;  // Density of Si [g/cm^3]
-constexpr int MaxIter = 4;    // Max. supported iterations
+constexpr float Tolerance = 1e-12;                  // numerical tolerance
+constexpr int ClustersPerCell = 3;                  // number of clusters for a cell
+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]
+constexpr float Rho = 2.33f;                        // Density of Si [g/cm^3]
+constexpr int MaxIter = 4;                          // Max. supported iterations
+constexpr int MaxSelectedTrackletsPerCluster = 100; // vertexer: max lines per cluster
 
 namespace helpers
 {

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

@@ -50,6 +50,19 @@ struct LogLogThrottler {
     return false;
   }
 };
+
+struct TimingStats {
+  std::uint64_t calls = 0;
+  double totalTimeMs = 0.;
+
+  void add(double timeMs)
+  {
+    ++calls;
+    totalTimeMs += timeMs;
+  }
+  double averageTimeMs() const { return calls ? totalTimeMs / static_cast<double>(calls) : 0.; }
+};
+
 } // namespace o2::its
 
-#endif
+#endif

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

@@ -124,7 +124,7 @@ GPUhdi() int4 getBinsRect(const Cluster& currentCluster, const int layerIndex,
 
   if (zRangeMax < -utils.getLayerZ(layerIndex) ||
       zRangeMin > utils.getLayerZ(layerIndex) || zRangeMin > zRangeMax) {
-    return int4{0, 0, 0, 0};
+    return int4{-1, -1, -1, -1};
   }
 
   return int4{o2::gpu::GPUCommonMath::Max(0, utils.getZBinIndex(layerIndex, zRangeMin)),

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

@@ -89,6 +89,12 @@ GPUhdi() float smallestAngleDifference(float a, float b)
   return o2::gpu::CAMath::Remainderf(b - a, o2::constants::math::TwoPI);
 }
 
+GPUhdi() bool isPhiDifferenceBelow(const float phiA, const float phiB, const float phiCut)
+{
+  const float deltaPhi = o2::gpu::CAMath::Abs(phiA - phiB);
+  return deltaPhi < phiCut || deltaPhi > o2::constants::math::TwoPI - phiCut;
+}
+
 GPUhdi() constexpr float Sq(float v)
 {
   return v * v;