TPC SCD: MC truth at the TRD (entrance and ideal tracklet positions)

For tracks with TRD residuals (TRD seeds and ITS-TPC map seeds whose most complete track has TRD), TrackDataMC
now also holds the truth at the TRD entrance (first TRD track reference in time, sector frame) and, per layer, the
true y and z at the x of the tracklet in the tracklet's sector frame, i.e. the frame of the stored TRD residual
(nearest TRD track reference propagated with the workflow's material correction), for the true particle of the
ITS-TPC part (ideal tracklets). TrackInterpolation records which ITS-TPC-TRD track gave the TRD residuals of each
stored track (getTRDGIDsSuccess, aligned with the track data).

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

Author: matthias-kleiner

Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/TPCInterpolationSpec.cxx

@@ -177,15 +177,16 @@ void TPCInterpolationDPL::fillMCTruth(const RecoContainer& recoData)
 {
   // MC truth for every stored TrackData: labels of the ITS-TPC part of the seed and of its ITS and TPC parts, the truth at
   // the ITS outer parameters (the ITS track reference nearest to TrackData::par, propagated to its x with the material
-  // correction of the workflow and the mass of the true particle) and the truth at the TPC entrance (the first TPC track
-  // reference in time, in the sector frame)
+  // correction of the workflow and the mass of the true particle), the truth at the TPC entrance (the first TPC track
+  // reference in time, in the sector frame) and, for TRD-matched seeds, the truth at the TRD entrance and the true
+  // positions at the x of the TRD tracklets of the track (ideal tracklets)
   const auto& trkData = mInterpolation.getReferenceTracks();
   mTrackDataMC.clear();
   mTrackDataMC.resize(trkData.size());
   struct Lookup {
-    o2::MCCompLabel lbl;
-    uint32_t idx;
-    bool its; // ITS outer (true) or TPC entrance (false)
+    o2::MCCompLabel lbl{};
+    uint32_t idx{0};
+    uint8_t kind{0}; // 0: ITS outer, 1: TPC entrance, 2: TRD
   };
   std::vector<Lookup> lookups;
   lookups.reserve(2 * trkData.size());
@@ -208,10 +209,13 @@ void TPCInterpolationDPL::fillMCTruth(const RecoContainer& recoData)
     const auto& lblITS = mc.labelITS.isValid() ? mc.labelITS : mc.label;
     const auto& lblTPC = mc.labelTPC.isValid() ? mc.labelTPC : mc.label;
     if (lblITS.isValid()) {
-      lookups.push_back({lblITS, uint32_t(i), true});
+      lookups.push_back({lblITS, uint32_t(i), 0});
     }
     if (lblTPC.isValid()) {
-      lookups.push_back({lblTPC, uint32_t(i), false});
+      lookups.push_back({lblTPC, uint32_t(i), 1});
+    }
+    if (mInterpolation.getTRDGIDsSuccess()[i].isIndexSet() && lblTPC.isValid()) { // track with TRD residuals: the true particle of the ITS-TPC part at the TRD
+      lookups.push_back({mc.label.isValid() ? mc.label : lblTPC, uint32_t(i), 2});
     }
   }
   // the reader loads the kinematics of a whole event (can be >100 MB): process event by event and release it right after
@@ -265,7 +269,7 @@ void TPCInterpolationDPL::fillMCTruth(const RecoContainer& recoData)
     const auto* pPDG = mcTrk ? O2DatabasePDG::Instance()->GetParticle(pdg) : nullptr;
     int charge = pPDG ? int(std::lround(pPDG->Charge() / 3.)) : 0; // TParticlePDG charge is in units of |e|/3
     auto refs = mMCReader->getTrackRefs(lbl.getSourceID(), lbl.getEventID(), lbl.getTrackID());
-    if (lk.its) { // ITS outer: track reference of the ITS part nearest to TrackData::par
+    if (lk.kind == 0) { // ITS outer: track reference of the ITS part nearest to TrackData::par
       mc.pdg = pdg;
       const o2::TrackReference* best = nullptr;
       float bestD2 = 1e30f;
@@ -291,7 +295,7 @@ void TPCInterpolationDPL::fillMCTruth(const RecoContainer& recoData)
           mc.flags |= TrackDataMC::HasITSOut;
         }
       }
-    } else { // TPC entrance: first TPC track reference in time of the TPC part
+    } else if (lk.kind == 1) { // TPC entrance: first TPC track reference in time of the TPC part
       if (!mc.labelITS.isValid() && !mc.label.isValid()) {
         mc.pdg = pdg; // no ITS lookup for this track
       }
@@ -321,6 +325,56 @@ void TPCInterpolationDPL::fillMCTruth(const RecoContainer& recoData)
           }
         }
       }
+    } else if (charge) { // TRD: entrance (first TRD track reference in time) and the true positions at the tracklet x of each layer
+      const o2::TrackReference* first = nullptr;
+      for (const auto& ref : refs) {
+        if (ref.getDetectorId() == DetID::TRD && (!first || ref.getTime() < first->getTime())) {
+          first = &ref;
+        }
+      }
+      if (!first) {
+        continue;
+      }
+      mc.parTRDIn = refToPar(*first, charge, pdg, true);
+      mc.flags |= TrackDataMC::HasTRDIn;
+      const auto& trkTRD = recoData.getITSTPCTRDTrack<o2::trd::TrackTRD>(mInterpolation.getTRDGIDsSuccess()[lk.idx]); // the TRD track of the stored TRD residuals
+      const auto trkltsCalib = recoData.getTRDCalibratedTracklets();
+      const auto tracklets = recoData.getTRDTracklets();
+      for (int iLayer = 0; iLayer < o2::trd::constants::NLAYER; iLayer++) {
+        int trkltIdx = trkTRD.getTrackletIndex(iLayer);
+        if (trkltIdx < 0) {
+          continue;
+        }
+        const auto& sp = trkltsCalib[trkltIdx]; // tracklet x, y, z in its sector frame, as used for the TRD residual
+        int sec = tracklets[trkltIdx].getDetector() / (o2::trd::constants::NLAYER * o2::trd::constants::NSTACK);
+        float alpha = o2::math_utils::sector2Angle(sec);
+        float cs = std::cos(alpha);
+        float sn = std::sin(alpha);
+        float gx = sp.getX() * cs - sp.getY() * sn;
+        float gy = sp.getX() * sn + sp.getY() * cs;
+        float gz = sp.getZ();
+        const o2::TrackReference* best = nullptr;
+        float bestD2 = 1e30f;
+        for (const auto& ref : refs) {
+          if (ref.getDetectorId() != DetID::TRD) {
+            continue;
+          }
+          float dx = ref.X() - gx;
+          float dy = ref.Y() - gy;
+          float dz = ref.Z() - gz;
+          float d2 = dx * dx + dy * dy + dz * dz;
+          if (d2 < bestD2) {
+            bestD2 = d2;
+            best = &ref;
+          }
+        }
+        auto par = refToPar(*best, charge, pdg, false);
+        if (par.rotateParam(alpha) && prop->PropagateToXBxByBz(par, sp.getX(), 0.999f, o2::base::Propagator::MAX_STEP, matCorr)) {
+          mc.yTRD[iLayer] = par.getY();
+          mc.zTRD[iLayer] = par.getZ();
+          mc.trdLayerMask |= uint8_t(1) << iLayer;
+        }
+      }
     }
   }
   if (curSrc >= 0) {

Detectors/TPC/calibration/SpacePoints/include/SpacePoints/TrackInterpolation.h

@@ -249,20 +249,26 @@ struct TrackData {
 struct TrackDataMC {
   enum Flags : uint8_t { HasITSOut = 0x1,    ///< parITSOut is filled
                          HasTPCIn = 0x2,     ///< parTPCIn is filled
-                         FakeITSTPC = 0x4 }; ///< ITS and TPC parts of the track have different MC labels
+                         FakeITSTPC = 0x4,   ///< ITS and TPC parts of the track have different MC labels
+                         HasTRDIn = 0x8 };   ///< parTRDIn is filled
   o2::MCCompLabel label{};                   ///< MC label of the ITS-TPC part of the seeding track
   o2::MCCompLabel labelITS{};                ///< MC label of its ITS part
   o2::MCCompLabel labelTPC{};                ///< MC label of its TPC part
   o2::track::TrackPar parITSOut{};           ///< truth at x and alpha of TrackData::par, from the nearest ITS track reference (propagated with the material correction, true mass)
   o2::track::TrackPar parTPCIn{};            ///< truth at the first TPC track reference (sector frame)
   float distITSRef{-1.f};                    ///< 3D distance between the ITS track reference used and TrackData::par in cm
   float distTPCRef{-1.f};                    ///< distance (y,z) between parTPCIn propagated to the innermost TPC cluster of the track and that cluster in cm (large: wrong leg, looper, fake)
+  o2::track::TrackPar parTRDIn{};            ///< truth at the first TRD track reference (sector frame), TRD-matched seeds only
+  float yTRD[6] = {};                        ///< truth y at the x of the TRD tracklet of each layer (tracklet sector frame), see trdLayerMask
+  float zTRD[6] = {};                        ///< truth z at the x of the TRD tracklet of each layer (tracklet sector frame), see trdLayerMask
+  uint8_t trdLayerMask{0};                   ///< bit i set: yTRD[i], zTRD[i] filled
   int pdg{0};                                ///< PDG code of the MC particle of the ITS part (TPC part if no ITS label)
   uint8_t flags{0};
   bool hasITSOut() const { return flags & HasITSOut; }
   bool hasTPCIn() const { return flags & HasTPCIn; }
   bool isFakeITSTPC() const { return flags & FakeITSTPC; }
-  ClassDefNV(TrackDataMC, 1);
+  bool hasTRDIn() const { return flags & HasTRDIn; }
+  ClassDefNV(TrackDataMC, 2);
 };
 
 /// \class TrackInterpolation
@@ -445,6 +451,8 @@ class TrackInterpolation
   std::vector<TrackDataCompact>& getTrackDataCompact() { return mTrackDataCompact; }
   std::vector<TrackDataExtended>& getTrackDataExtended() { return mTrackDataExtended; }
   std::vector<TrackData>& getReferenceTracks() { return mTrackData; }
+  /// ITS-TPC-TRD track whose tracklets gave the TRD residuals of each stored track (not set if none), aligned with getReferenceTracks()
+  const std::vector<o2::dataformats::GlobalTrackID>& getTRDGIDsSuccess() const { return mTRDGIDsSuccess; }
 
   void setLane(int lID, int nL)
   {
@@ -512,6 +520,7 @@ class TrackInterpolation
   // cache
   std::array<CacheStruct, constants::MAXGLOBALPADROW> mCache{{}}; ///< caching positions, covariances and angles for track extrapolations and interpolation
   std::vector<o2::dataformats::GlobalTrackID> mGIDsSuccess;       ///< keep track of the GIDs which could be processed successfully
+  std::vector<o2::dataformats::GlobalTrackID> mTRDGIDsSuccess;    ///< ITS-TPC-TRD track used for the TRD residuals of each stored track (not set if none)
 
   TrackValidationData mTrackValidation;
 

Detectors/TPC/calibration/SpacePoints/src/TrackInterpolation.cxx

@@ -807,10 +807,12 @@ void TrackInterpolation::interpolateTrack(int iSeed)
     }
 
     bool stopPropagation = !mExtDetResid;
+    GTrackID gidTRDUsed{};
     if (!stopPropagation) {
       // do we have TRD residuals to add?
       trkWork = trkOuter;
       if (!allLost && gidTable[GTrackID::TRD].isIndexSet()) { // allLost: trkOuter is not a valid outer param
+        gidTRDUsed = gidTable[GTrackID::ITSTPCTRD];
         const auto& trkTRD = mRecoCont->getITSTPCTRDTrack<o2::trd::TrackTRD>(gidTable[GTrackID::ITSTPCTRD]);
         for (int iLayer = 0; iLayer < o2::trd::constants::NLAYER; iLayer++) {
           std::array<float, 2> trkltTRDYZ{};
@@ -926,6 +928,7 @@ void TrackInterpolation::interpolateTrack(int iSeed)
     }
 
     mGIDsSuccess.push_back(mGIDs[iSeed]);
+    mTRDGIDsSuccess.push_back(gidTRDUsed);
     mTrackDataCompact.emplace_back(trackData.clIdx.getFirstEntry(), trackData.multStack, nClValidated, mGIDs[iSeed].getSource(), trackData.nExtDetResid, trackData.filterFlag);
     mTrackData.push_back(std::move(trackData));
     stored = true;
@@ -1187,12 +1190,14 @@ void TrackInterpolation::extrapolateTrack(int iSeed)
     }
 
     bool stopPropagation = !mExtDetResid;
+    GTrackID gidTRDUsed{};
     if (!stopPropagation) {
       // do we have TRD residuals to add?
       int iSeedFull = mParentID[iSeed] == -1 ? iSeed : mParentID[iSeed];
       auto gidFull = mGIDs[iSeedFull];
       const auto& gidTableFull = mGIDtables[iSeedFull];
       if (!refLost && gidTableFull[GTrackID::TRD].isIndexSet()) { // refLost: trkWork did not reach the TPC outer end
+        gidTRDUsed = gidTableFull[GTrackID::ITSTPCTRD];
         const auto& trkTRD = mRecoCont->getITSTPCTRDTrack<o2::trd::TrackTRD>(gidTableFull[GTrackID::ITSTPCTRD]);
         trackData.nTrkltsTRD = trkTRD.getNtracklets();
         trackData.chi2TRD = trkTRD.getChi2();
@@ -1319,6 +1324,7 @@ void TrackInterpolation::extrapolateTrack(int iSeed)
     mTrackData.push_back(std::move(trackData));
     stored = true;
     mGIDsSuccess.push_back(mGIDs[iSeed]);
+    mTRDGIDsSuccess.push_back(gidTRDUsed);
     mTrackDataCompact.emplace_back(trackData.clIdx.getFirstEntry(), trackData.multStack, nClValidated, mGIDs[iSeed].getSource(), trackData.nExtDetResid, trackData.filterFlag);
     if (mDumpTrackPoints) {
       (*trackDataExtended).clIdx.setEntries(nClValidated);
@@ -1694,6 +1700,7 @@ void TrackInterpolation::reset()
   mClRes.clear();
   mDetInfoRes.clear();
   mGIDsSuccess.clear();
+  mTRDGIDsSuccess.clear();
   for (auto& vec : mTrackIndices) {
     vec.clear();
   }