TPC TimeSeries: fix silent track loss from binning overflow (#15658)

* TPC TimeSeries: fix silent track loss from binning overflow

Bin indices for tgl, phi, qPt, and multiplicity were used as implicit
track selection cuts: tracks outside histogram range were silently
dropped (return). Replace with std::clamp — edge bins become overflow
bins (standard ROOT convention). No change for tracks within range.

Bug: changing --max-qPt or --mult-max removed tracks from ALL outputs
(DCA, dEdx, etc.), not just the binned histograms.

* TPC TimeSeries: clamp fix + ITS cluster sizes + TRD tracklets + TRD matching

Phase 0.2 — binning overflow fix:
  - Replace bounds-check-and-return with std::clamp on all 4 bin indices
  - Edge bins act as saturated overflow; no tracks silently dropped

Phase 0.3 D1 — ITS cluster sizes (per-track, unbinned):
  - itsClusterSizes: packed 4-bit per layer (bit 28 kSharedClusters masked)
  - itsHasSharedClusters, itsPattern: 7-bit layer hit pattern

Phase 0.3 D2 — TRD tracklet objects (per-track, unbinned):
  - Native Tracklet64[6] and CalibratedTracklet[6] per layer
  - trdPattern (6-bit validity mask), nTRDTracklets
  - requestTRDTracklets added to DataRequest

Phase 0.3 D3 — TRD matching fraction (per-TF):
  - nITSTPCBasedPVContributors, nITSTPCWithTRDPVContributors, fracTRD
  - NaN for zero denominator. ClassDefNV 7 -> 8.

* TPC TimeSeries: store TRD tracklets as native objects (std::vector)

Replace flat primitive arrays with std::vector<Tracklet64> and
std::vector<CalibratedTracklet>. std::array failed ROOT serialization
(missing ShowMember); std::vector with ROOT dictionary works.

* Clang

---------

Co-authored-by: miranov25 <marian.ivanov@cern.cg>

Author: miranov25

Detectors/Calibration/include/DetectorsCalibration/IntegratedClusterCalibrator.h

@@ -365,6 +365,10 @@ struct TimeSeriesITSTPC {
   std::vector<float> vertexY_ITSTPC_RMS;                ///< vertex y RMS with ITS-TPC cut (nContributorsITS + nContributorsITSTPC)<0.95
   std::vector<float> vertexZ_ITSTPC_RMS;                ///< vertex z RMS with ITS-TPC cut (nContributorsITS + nContributorsITSTPC)<0.95
 
+  std::vector<float> nITSTPCBasedPVContributors;   ///< number of ITS-TPC-based PV contributors (denominator for TRD matching fraction)
+  std::vector<float> nITSTPCWithTRDPVContributors; ///< number of ITS-TPC-TRD PV contributors (numerator for TRD matching fraction)
+  std::vector<float> fracTRD;                      ///< fraction of ITS-TPC PV contributors with TRD match (NaN if denominator=0)
+
   int quantileValues = 23;                          ///<! number of values in quantiles + truncated mean (hardcoded for the moment)
   std::vector<float> nVertexContributors_Quantiles; ///< number of primary vertices for quantiles 0.1, 0.2, ... 0.9 and truncated mean values 0.05->0.95, 0.1->0.9, 0.2->0.8
 
@@ -498,12 +502,15 @@ struct TimeSeriesITSTPC {
     vertexX_ITSTPC_RMS.resize(nTotalVtx);
     vertexY_ITSTPC_RMS.resize(nTotalVtx);
     vertexZ_ITSTPC_RMS.resize(nTotalVtx);
+    nITSTPCBasedPVContributors.resize(nTotalVtx);
+    nITSTPCWithTRDPVContributors.resize(nTotalVtx);
+    fracTRD.resize(nTotalVtx);
 
     const int nTotalQ = quantileValues * nTotal / mTSTPC.getNBins();
     nVertexContributors_Quantiles.resize(nTotalQ);
   }
 
-  ClassDefNV(TimeSeriesITSTPC, 7);
+  ClassDefNV(TimeSeriesITSTPC, 8);
 };
 
 } // end namespace tpc

Detectors/TPC/workflow/src/TPCTimeSeriesSpec.cxx

@@ -41,6 +41,9 @@
 #include <chrono>
 #include "DataFormatsTPC/PIDResponse.h"
 #include "DataFormatsITS/TrackITS.h"
+#include "DataFormatsTRD/TrackTRD.h"
+#include "DataFormatsTRD/Tracklet64.h"
+#include "DataFormatsTRD/CalibratedTracklet.h"
 #include "TROOT.h"
 #include "ReconstructionDataFormats/MatchInfoTOF.h"
 #include "DataFormatsTOF/Cluster.h"
@@ -63,6 +66,13 @@ namespace tpc
 class TPCTimeSeries : public Task
 {
  public:
+  /// D2: per-track TRD tracklet lookup data
+  struct TRDTrackletData {
+    uint8_t trdPattern = 0;
+    uint8_t nTRDTracklets = 0;
+    int trackletIndices[6] = {-1, -1, -1, -1, -1, -1};
+  };
+
   /// \constructor
   TPCTimeSeries(std::shared_ptr<o2::base::GRPGeomRequest> req, const bool disableWriter, const o2::base::Propagator::MatCorrType matType, const bool enableUnbinnedWriter, const bool tpcOnly, std::shared_ptr<o2::globaltracking::DataRequest> dr) : mCCDBRequest(req), mDisableWriter(disableWriter), mMatType(matType), mUnbinnedWriter(enableUnbinnedWriter), mTPCOnly(tpcOnly), mDataRequest(dr) {};
 
@@ -303,6 +313,38 @@ class TPCTimeSeries : public Task
     // find nearest vertex of tracks which have no vertex assigned
     findNearesVertex(tracksTPC, vertices);
 
+    // D2: build TPC track index → TRD tracklet data map (for unbinned output)
+    // For each TPC track that has a TRD match, store the TrackTRD tracklet indices
+    std::unordered_map<unsigned int, TRDTrackletData> tpcToTRDMap;
+    auto trdTracklets = mTPCOnly ? gsl::span<const o2::trd::Tracklet64>() : recoData.getTRDTracklets();
+    auto trdCalibTracklets = mTPCOnly ? gsl::span<const o2::trd::CalibratedTracklet>() : recoData.getTRDCalibratedTracklets();
+    if (mUnbinnedWriter && !mTPCOnly) {
+      // scan ITS-TPC-TRD tracks
+      auto itstpctrdTracks = recoData.getITSTPCTRDTracks<o2::trd::TrackTRD>();
+      for (unsigned int ig = 0; ig < itstpctrdTracks.size(); ++ig) {
+        auto gid = GTrackID(ig, GTrackID::ITSTPCTRD);
+        auto refTPC = recoData.getTPCContributorGID(gid);
+        if (!refTPC.isIndexSet()) {
+          continue;
+        }
+        auto refTRD = recoData.getSingleDetectorRefs(gid)[GTrackID::TRD];
+        if (!refTRD.isIndexSet()) {
+          continue;
+        }
+        const auto& trdTrack = recoData.getTrack<o2::trd::TrackTRD>(refTRD);
+        TRDTrackletData trdData;
+        for (int iLay = 0; iLay < 6; ++iLay) {
+          auto trkltId = trdTrack.getTrackletIndex(iLay);
+          if (trkltId >= 0) {
+            trdData.trdPattern |= (1 << iLay);
+            trdData.nTRDTracklets++;
+            trdData.trackletIndices[iLay] = trkltId;
+          }
+        }
+        tpcToTRDMap[refTPC] = trdData;
+      }
+    }
+
     // getting cluster references for cluster bitmask
     if (mUnbinnedWriter) {
       mTPCTrackClIdx = pc.inputs().get<gsl::span<o2::tpc::TPCClRefElem>>("trackTPCClRefs");
@@ -472,7 +514,7 @@ class TPCTimeSeries : public Task
       auto myThread = [&](int iThread) {
         for (size_t i = iThread; i < loopEnd; i += mNThreads) {
           if (acceptTrack(tracksTPC[i])) {
-            fillDCA(tracksTPC, tracksITSTPC, vertices, i, iThread, indicesITSTPC, tracksITS, idxTPCTrackToTOFCluster, tofClusters);
+            fillDCA(tracksTPC, tracksITSTPC, vertices, i, iThread, indicesITSTPC, tracksITS, idxTPCTrackToTOFCluster, tofClusters, tpcToTRDMap, trdTracklets, trdCalibTracklets);
           }
         }
       };
@@ -489,7 +531,7 @@ class TPCTimeSeries : public Task
       auto myThread = [&](int iThread) {
         for (size_t i = iThread; i < loopEnd; i += mNThreads) {
           if (acceptTrack(tracksTPC[i])) {
-            fillDCA(tracksTPC, tracksITSTPC, vertices, i, iThread, indicesITSTPC, tracksITS, idxTPCTrackToTOFCluster, tofClusters);
+            fillDCA(tracksTPC, tracksITSTPC, vertices, i, iThread, indicesITSTPC, tracksITS, idxTPCTrackToTOFCluster, tofClusters, tpcToTRDMap, trdTracklets, trdCalibTracklets);
           }
         }
       };
@@ -1133,7 +1175,7 @@ class TPCTimeSeries : public Task
     return isGoodTrack;
   }
 
-  void fillDCA(const gsl::span<const TrackTPC> tracksTPC, const gsl::span<const o2::dataformats::TrackTPCITS> tracksITSTPC, const gsl::span<const o2::dataformats::PrimaryVertex> vertices, const int iTrk, const int iThread, const std::unordered_map<unsigned int, std::array<int, 2>>& indicesITSTPC, const gsl::span<const o2::its::TrackITS> tracksITS, const std::vector<std::tuple<int, float, float, o2::track::TrackLTIntegral, double, float, unsigned int, unsigned short>>& idxTPCTrackToTOFCluster, const gsl::span<const o2::tof::Cluster> tofClusters)
+  void fillDCA(const gsl::span<const TrackTPC> tracksTPC, const gsl::span<const o2::dataformats::TrackTPCITS> tracksITSTPC, const gsl::span<const o2::dataformats::PrimaryVertex> vertices, const int iTrk, const int iThread, const std::unordered_map<unsigned int, std::array<int, 2>>& indicesITSTPC, const gsl::span<const o2::its::TrackITS> tracksITS, const std::vector<std::tuple<int, float, float, o2::track::TrackLTIntegral, double, float, unsigned int, unsigned short>>& idxTPCTrackToTOFCluster, const gsl::span<const o2::tof::Cluster> tofClusters, const std::unordered_map<unsigned int, TRDTrackletData>& tpcToTRDMap, const gsl::span<const o2::trd::Tracklet64> trdTracklets, const gsl::span<const o2::trd::CalibratedTracklet> trdCalibTracklets)
   {
     const auto& trackFull = tracksTPC[iTrk];
     const bool isGoodTrack = checkTrack(trackFull);
@@ -1179,21 +1221,22 @@ class TPCTimeSeries : public Task
       return;
     }
 
-    const int tglBin = mTglBins * std::abs(trackTmp.getTgl()) / mMaxTgl + mPhiBins;
-    const int phiBin = mPhiBins * trackTmp.getPhi() / o2::constants::math::TwoPI;
+    // Saturate bin indices — edge bins act as overflow (Phase 0.2 fix)
+    const int tglBin = std::clamp(static_cast<int>(mTglBins * std::abs(trackTmp.getTgl()) / mMaxTgl) + mPhiBins,
+                                  mPhiBins, mPhiBins + mTglBins - 1);
+    const int phiBin = std::clamp(static_cast<int>(mPhiBins * trackTmp.getPhi() / o2::constants::math::TwoPI),
+                                  0, mPhiBins - 1);
 
     const int offsQPtBin = mPhiBins + mTglBins;
-    const int qPtBin = offsQPtBin + mQPtBins * (trackTmp.getQ2Pt() + mMaxQPt) / (2 * mMaxQPt);
+    const int qPtBin = std::clamp(offsQPtBin + static_cast<int>(mQPtBins * (trackTmp.getQ2Pt() + mMaxQPt) / (2 * mMaxQPt)),
+                                  offsQPtBin, offsQPtBin + mQPtBins - 1);
     const int localMult = mNTracksWindow[iTrk];
 
     const int offsMult = offsQPtBin + mQPtBins;
-    const int multBin = offsMult + mMultBins * localMult / mMultMax;
+    const int multBin = std::clamp(offsMult + static_cast<int>(mMultBins * localMult / mMultMax),
+                                   offsMult, offsMult + mMultBins - 1);
     const int nBins = getNBins();
 
-    if ((phiBin < 0) || (phiBin > mPhiBins) || (tglBin < mPhiBins) || (tglBin > offsQPtBin) || (qPtBin < offsQPtBin) || (qPtBin > offsMult) || (multBin < offsMult) || (multBin > offsMult + mMultBins)) {
-      return;
-    }
-
     float sigmaY2 = 0;
     float sigmaZ2 = 0;
     const int sector = o2::math_utils::angle2Sector(trackTmp.getPhiPos());
@@ -1354,6 +1397,30 @@ class TPCTimeSeries : public Task
         const float chi2match_ITSTPC = hasITSTPC ? tracksITSTPC[idxITSTPC.front()].getChi2Match() : -1;
         const int nClITS = idxITSCheck ? tracksITS[idxITSTrack].getNClusters() : -1;
         const int chi2ITS = idxITSCheck ? tracksITS[idxITSTrack].getChi2() : -1;
+        // D1: ITS cluster sizes (4-bit per layer, mask bit 28 = kSharedClusters)
+        const uint32_t itsClusterSizes = idxITSCheck ? (static_cast<uint32_t>(tracksITS[idxITSTrack].getClusterSizes()) & 0x0FFFFFFFu) : 0u;
+        const bool itsHasSharedClusters = idxITSCheck ? tracksITS[idxITSTrack].hasSharedClusters() : false;
+        const uint32_t itsPattern = idxITSCheck ? (tracksITS[idxITSTrack].getPattern() & 0x7Fu) : 0u;
+
+        // D2: TRD tracklet data — native objects per layer
+        uint8_t trdPattern = 0;
+        uint8_t nTRDTracklets = 0;
+        std::vector<o2::trd::Tracklet64> trdTrackletVec(6);
+        std::vector<o2::trd::CalibratedTracklet> trdCalibVec(6);
+        auto itTRD = tpcToTRDMap.find(iTrk);
+        if (itTRD != tpcToTRDMap.end()) {
+          const auto& trdData = itTRD->second;
+          trdPattern = trdData.trdPattern;
+          nTRDTracklets = trdData.nTRDTracklets;
+          for (int iLay = 0; iLay < 6; ++iLay) {
+            if (trdData.trackletIndices[iLay] >= 0) {
+              trdTrackletVec[iLay] = trdTracklets[trdData.trackletIndices[iLay]];
+              if (trdData.trackletIndices[iLay] < static_cast<int>(trdCalibTracklets.size())) {
+                trdCalibVec[iLay] = trdCalibTracklets[trdData.trackletIndices[iLay]];
+              }
+            }
+          }
+        }
         int typeSide = 2; // A- and C-Side cluster
         if (trackFull.hasASideClustersOnly()) {
           typeSide = 0;
@@ -1488,6 +1555,14 @@ class TPCTimeSeries : public Task
                             << "mX_ITS=" << mx_ITS
                             << "nClITS=" << nClITS
                             << "chi2ITS=" << chi2ITS
+                            << "itsClusterSizes=" << itsClusterSizes
+                            << "itsHasSharedClusters=" << itsHasSharedClusters
+                            << "itsPattern=" << itsPattern
+                            // D2: TRD tracklet data
+                            << "trdPattern=" << trdPattern
+                            << "nTRDTracklets=" << nTRDTracklets
+                            << "trdTracklets=" << trdTrackletVec
+                            << "trdCalibTracklets=" << trdCalibVec
                             << "chi2match_ITSTPC=" << chi2match_ITSTPC
                             << "PID=" << trkOrig.getPID().getID()
                             // TPC cov at vertex (without vertex constrained)
@@ -1680,6 +1755,7 @@ class TPCTimeSeries : public Task
 
     std::unordered_map<int, int> nContributors_ITS;    // ITS: vertex ID -> n contributors
     std::unordered_map<int, int> nContributors_ITSTPC; // ITS-TPC (and ITS-TPC-TRD, ITS-TPC-TOF, ITS-TPC-TRD-TOF): vertex ID -> n contributors
+    std::unordered_map<int, int> nContributors_TRD;    // ITS-TPC-TRD (and ITS-TPC-TRD-TOF): vertex ID -> n TRD-matched PV contributors
 
     // loop over collisions
     if (!vertices.empty()) {
@@ -1700,6 +1776,10 @@ class TPCTimeSeries : public Task
               if (refITSTPC.isIndexSet()) {
                 indicesITSTPC_vtx[refITSTPC] = vID;
                 ++nContributors_ITSTPC[vID];
+                // count TRD-matched PV contributors
+                if (source == TrkSrc::ITSTPCTRD || source == TrkSrc::ITSTPCTRDTOF) {
+                  ++nContributors_TRD[vID];
+                }
               } else {
                 ++nContributors_ITS[vID];
               }
@@ -1761,6 +1841,17 @@ class TPCTimeSeries : public Task
     mBufferDCA.vertexY_ITSTPC_RMS.front() = avgVtxITSTPC[1].getStdDev();
     mBufferDCA.vertexZ_ITSTPC_RMS.front() = avgVtxITSTPC[2].getStdDev();
 
+    // TRD matching fraction (summed over all vertices in this TF)
+    int sumITSTPCBased = 0;
+    int sumWithTRD = 0;
+    for (int ivtx = 0; ivtx < vertices.size(); ++ivtx) {
+      sumITSTPCBased += nContributors_ITSTPC[ivtx];
+      sumWithTRD += nContributors_TRD[ivtx];
+    }
+    mBufferDCA.nITSTPCBasedPVContributors.front() = sumITSTPCBased;
+    mBufferDCA.nITSTPCWithTRDPVContributors.front() = sumWithTRD;
+    mBufferDCA.fracTRD.front() = (sumITSTPCBased > 0) ? static_cast<float>(sumWithTRD) / sumITSTPCBased : std::nanf("");
+
     // quantiles and truncated mean
     RobustAverage avg(vertices.size(), false);
     for (const auto& vtx : vertices) {
@@ -1850,6 +1941,10 @@ o2::framework::DataProcessorSpec getTPCTimeSeriesSpec(const bool disableWriter,
   if (src[GTrackID::TPC]) {
     dataRequest->requestClusters(GTrackID::getSourcesMask("TPC"), useMC);
   }
+  // D2: request TRD tracklets for tracks with TRD contribution
+  if (srcTracks[GTrackID::ITSTPCTRD] || srcTracks[GTrackID::ITSTPCTRDTOF]) {
+    dataRequest->requestTRDTracklets(useMC);
+  }
 
   bool tpcOnly = srcTracks == GTrackID::getSourcesMask("TPC");
   if (srcTracks.any() && !tpcOnly) {