Add residuals for external detectors

shahor02 · shahor02 · commit a23c62ae4c9a · 2025-08-18T00:08:45.000+02:00
diff --git a/Detectors/TPC/calibration/SpacePoints/include/SpacePoints/TrackInterpolation.h b/Detectors/TPC/calibration/SpacePoints/include/SpacePoints/TrackInterpolation.h
@@ -121,12 +121,14 @@ struct TrackData {
   float chi2TPC{};                           ///< chi2 of TPC track
   float chi2ITS{};                           ///< chi2 of ITS track
   float chi2TRD{};                           ///< chi2 of TRD track
+
   unsigned short nClsTPC{};                  ///< number of attached TPC clusters
   unsigned short nClsITS{};                  ///< number of attached ITS clusters
   unsigned short nTrkltsTRD{};               ///< number of attached TRD tracklets
   unsigned short clAvailTOF{};               ///< whether or not track seed has a matched TOF cluster
-  o2::dataformats::RangeReference<> clIdx{}; ///< index of first cluster residual and total number of cluster residuals of this track
-  ClassDefNV(TrackData, 6);
+  uint8_t nExtDetResid = 0;                  ///< number of external detectors (to TPC) residuals stored, on top of clIdx.getEntries
+  o2::dataformats::RangeReference<> clIdx{}; ///< index of first cluster residual and total number of TPC cluster residuals of this track
+  ClassDefNV(TrackData, 7);
 };
 
 /// \class TrackInterpolation
@@ -265,6 +267,8 @@ class TrackInterpolation
   /// Set the centre of mass energy required for pT downsampling Tsalis function
   void setSqrtS(float s) { mSqrtS = s; }
 
+  int processTRDLayer(const o2::trd::TrackTRD& trkTRD, int iLayer, o2::track::TracParCov& trkWork, std::array<float, 2>* trkltTRDYZ = nullptr, std::array<float, 3>* trkltTRDCov = nullptr);
+
   // --------------------------------- output ---------------------------------------------
   std::vector<UnbinnedResid>& getClusterResiduals() { return mClRes; }
   std::vector<TrackDataCompact>& getTrackDataCompact() { return mTrackDataCompact; }
@@ -297,6 +301,7 @@ class TrackInterpolation
   std::vector<o2::globaltracking::RecoContainer::GlobalIDSet> mGIDtables{}; ///< GIDs of contributors from single detectors for each seed
   std::vector<float> mTrackTimes{};                                         ///< time estimates for all input tracks in micro seconds
   std::vector<o2::track::TrackParCov> mSeeds{};                             ///< seeding track parameters (ITS tracks)
+  std::vector<int> mParentID{};                                             ///< entry of more global parent track for skimmed seeds (-1: no parent)
   std::map<int, int> mTrackTypes;                                           ///< mapping of track source to array index in mTrackIndices
   std::array<std::vector<uint32_t>, 4> mTrackIndices;                       ///< keep GIDs of input tracks separately for each track type
   gsl::span<const TPCClRefElem> mTPCTracksClusIdx;                          ///< input TPC cluster indices from span
diff --git a/Detectors/TPC/calibration/SpacePoints/src/TrackInterpolation.cxx b/Detectors/TPC/calibration/SpacePoints/src/TrackInterpolation.cxx
@@ -18,6 +18,7 @@
 #include "SpacePoints/TrackInterpolation.h"
 #include "SpacePoints/TrackResiduals.h"
 #include "ITStracking/IOUtils.h"
+#include "ITSBase/GeometryTGeo.h"
 #include "TPCBase/ParameterElectronics.h"
 #include "DataFormatsTPC/TrackTPC.h"
 #include "DataFormatsTPC/Defs.h"
@@ -65,6 +66,9 @@ void TrackInterpolation::init(o2::dataformats::GlobalTrackID::mask_t src, o2::da
   mTrackTypes.insert({GTrackID::ITSTPCTOF, 2});
   mTrackTypes.insert({GTrackID::ITSTPCTRDTOF, 3});
 
+  o2::its::GeometryTGeo* geom = GeometryTGeo::Instance();
+  geom->fillMatrixCache(o2::math_utils::bit2Mask(o2::math_utils::TransformType::T2L, o2::math_utils::TransformType::L2G));
+
   mInitDone = true;
   LOGP(info, "Done initializing TrackInterpolation. Configured track input: {}. Track input specifically for map: {}",
        GTrackID::getSourcesNames(mSourcesConfigured), mSingleSourcesConfigured ? "identical" : GTrackID::getSourcesNames(mSourcesConfiguredMap));
@@ -273,6 +277,9 @@ void TrackInterpolation::process()
   trackIndices.insert(trackIndices.end(), mTrackIndices[mTrackTypes[GTrackID::ITSTPC]].begin(), mTrackIndices[mTrackTypes[GTrackID::ITSTPC]].end());
 
   int nSeeds = mSeeds.size(), lastChecked = 0;
+  mParentID.clear();
+  mParentID.resize(nSeeds, -1);
+
   int maxOutputTracks = (mMaxTracksPerTF >= 0) ? mMaxTracksPerTF + mAddTracksForMapPerTF : nSeeds;
   mTrackData.reserve(maxOutputTracks);
   mClRes.reserve(maxOutputTracks * param::NPadRows);
@@ -292,6 +299,7 @@ void TrackInterpolation::process()
       this->mGIDtables.push_back(this->mRecoCont->getSingleDetectorRefs(this->mGIDs.back()));
       this->mTrackTimes.push_back(this->mTrackTimes[seedIndex]);
       this->mSeeds.push_back(this->mSeeds[seedIndex]);
+      this->mParentID.push_back(seedIndex); // store parent seed id
     };
 
     GTrackID::mask_t partsAdded;
@@ -457,41 +465,15 @@ void TrackInterpolation::interpolateTrack(int iSeed)
       (*trackDataExtended).trkTRD = trkTRD;
     }
     for (int iLayer = o2::trd::constants::NLAYER - 1; iLayer >= 0; --iLayer) {
-      int trkltIdx = trkTRD.getTrackletIndex(iLayer);
-      if (trkltIdx < 0) {
-        // no TRD tracklet in this layer
+      std::array<float, 2> trkltTRDYZ{};
+      std::array<float, 3> trkltTRDCov{};
+      int res = processTRDLayer(trkTRD, iLayer, trkWork, &trkltTRDYZ, &trkltTRDCov);
+      if (res == 0) { // no TRD tracklet in this layer
         continue;
       }
-      const auto& trdSP = mRecoCont->getTRDCalibratedTracklets()[trkltIdx];
-      const auto& trdTrklt = mRecoCont->getTRDTracklets()[trkltIdx];
-      if (mDumpTrackPoints) {
-        (*trackDataExtended).trkltTRD.push_back(trdTrklt);
-        (*trackDataExtended).clsTRD.push_back(trdSP);
-      }
-      auto trkltDet = trdTrklt.getDetector();
-      auto trkltSec = trkltDet / (o2::trd::constants::NLAYER * o2::trd::constants::NSTACK);
-      if (trkltSec != o2::math_utils::angle2Sector(trkWork.getAlpha())) {
-        if (!trkWork.rotate(o2::math_utils::sector2Angle(trkltSec))) {
-          LOG(debug) << "Track could not be rotated in TRD tracklet coordinate system in layer " << iLayer;
-          return;
-        }
-      }
-      if (!propagator->PropagateToXBxByBz(trkWork, trdSP.getX(), mParams->maxSnp, mParams->maxStep, mMatCorr)) {
-        LOG(debug) << "Failed propagation to TRD layer " << iLayer;
+      if (res < 0) { // failed to reach this layer
         return;
       }
-
-      const auto* pad = mGeoTRD->getPadPlane(trkltDet);
-      float tilt = tan(TMath::DegToRad() * pad->getTiltingAngle()); // tilt is signed! and returned in degrees
-      float tiltCorrUp = tilt * (trdSP.getZ() - trkWork.getZ());
-      float zPosCorrUp = trdSP.getZ() + mRecoParam.getZCorrCoeffNRC() * trkWork.getTgl(); // maybe Z can be corrected on avarage already by the tracklet transformer?
-      float padLength = pad->getRowSize(trdTrklt.getPadRow());
-      if (!((trkWork.getSigmaZ2() < (padLength * padLength / 12.f)) && (std::fabs(trdSP.getZ() - trkWork.getZ()) < padLength))) {
-        tiltCorrUp = 0.f;
-      }
-      std::array<float, 2> trkltTRDYZ{trdSP.getY() - tiltCorrUp, zPosCorrUp};
-      std::array<float, 3> trkltTRDCov;
-      mRecoParam.recalcTrkltCov(tilt, trkWork.getSnp(), pad->getRowSize(trdTrklt.getPadRow()), trkltTRDCov);
       if (!trkWork.update(trkltTRDYZ, trkltTRDCov)) {
         LOG(debug) << "Failed to update track at TRD layer " << iLayer;
         return;
@@ -611,6 +593,46 @@ void TrackInterpolation::interpolateTrack(int iSeed)
   }
 }
 
+int TrackInterpolation::processTRDLayer(const o2::trd::TrackTRD& trkTRD, int iLayer, o2::track::TracParCov& trkWork,
+                                        std::array<float, 2>* trkltTRDYZ, std::array<float, 3>* trkltTRDCov)
+{
+  // return 1 in case of successful processing, 0 if there is no TRD tracklet at given layer, -1 if processing failed
+  int trkltIdx = trkTRD.getTrackletIndex(iLayer);
+  if (trkltIdx < 0) {
+    return 0; // no TRD tracklet in this layer
+  }
+  const auto& trdSP = mRecoCont->getTRDCalibratedTracklets()[trkltIdx];
+  const auto& trdTrklt = mRecoCont->getTRDTracklets()[trkltIdx];
+  auto trkltDet = trdTrklt.getDetector();
+  auto trkltSec = trkltDet / (o2::trd::constants::NLAYER * o2::trd::constants::NSTACK);
+  if (trkltSec != o2::math_utils::angle2Sector(trkWork.getAlpha())) {
+    if (!trkWork.rotate(o2::math_utils::sector2Angle(trkltSec))) {
+      LOG(debug) << "Track could not be rotated in TRD tracklet coordinate system in layer " << iLayer;
+      return -1;
+    }
+  }
+  if (!propagator->PropagateToXBxByBz(trkWork, trdSP.getX(), mParams->maxSnp, mParams->maxStep, mMatCorr)) {
+    LOG(debug) << "Failed propagation to TRD layer " << iLayer;
+    return -1;
+  }
+  if (trkltTRDYZ) {
+    const auto* pad = mGeoTRD->getPadPlane(trkltDet);
+    float tilt = tan(TMath::DegToRad() * pad->getTiltingAngle()); // tilt is signed! and returned in degrees
+    float tiltCorrUp = tilt * (trdSP.getZ() - trkWork.getZ());
+    float zPosCorrUp = trdSP.getZ() + mRecoParam.getZCorrCoeffNRC() * trkWork.getTgl(); // maybe Z can be corrected on avarage already by the tracklet transformer?
+    float padLength = pad->getRowSize(trdTrklt.getPadRow());
+    if (!((trkWork.getSigmaZ2() < (padLength * padLength / 12.f)) && (std::fabs(trdSP.getZ() - trkWork.getZ()) < padLength))) {
+      tiltCorrUp = 0.f;
+    }
+    (*trkltTRDYZ)[0] = trdSP.getY() - tiltCorrUp;
+    (*trkltTRDYZ)[1] = zPosCorrUp;
+    if (trkltTRDCov) {
+      mRecoParam.recalcTrkltCov(tilt, trkWork.getSnp(), pad->getRowSize(trdTrklt.getPadRow()), *trkltTRDCov);
+    }
+  }
+  return 1;
+}
+
 void TrackInterpolation::extrapolateTrack(int iSeed)
 {
   // extrapolate ITS-only track through TPC and store residuals to TPC clusters in the output vectors
@@ -638,7 +660,7 @@ void TrackInterpolation::extrapolateTrack(int iSeed)
   trackData.gid = mGIDs[iSeed];
   trackData.par = mSeeds[iSeed];
 
-  auto& trkWork = mSeeds[iSeed];
+  auto trkWork = mSeeds[iSeed];
   float clusterTimeBinOffset = mTrackTimes[iSeed] / mTPCTimeBinMUS;
   auto propagator = o2::base::Propagator::Instance();
   unsigned short rowPrev = 0; // used to calculate dRow of two consecutive cluster residuals
@@ -681,6 +703,13 @@ void TrackInterpolation::extrapolateTrack(int iSeed)
     rowPrev = row;
     ++nMeasurements;
   }
+
+  TrackParams params; // for refitted track parameters and flagging rejected clusters
+  if (clusterResiduals.size() > constants::MAXGLOBALPADROW) {
+    LOGP(warn, "Extrapolated ITS-TPC track and found more reesiduals than possible ({})", clusterResiduals.size());
+    return;
+  }
+
   trackData.chi2TPC = trkTPC.getChi2();
   trackData.chi2ITS = trkITS.getChi2();
   trackData.nClsTPC = trkTPC.getNClusterReferences();
@@ -691,19 +720,14 @@ void TrackInterpolation::extrapolateTrack(int iSeed)
     (*trackDataExtended).trkOuter = trkWork;
   }
 
-  TrackParams params; // for refitted track parameters and flagging rejected clusters
-  if (clusterResiduals.size() > constants::MAXGLOBALPADROW) {
-    LOGP(warn, "Extrapolated ITS-TPC track and found more reesiduals than possible ({})", clusterResiduals.size());
-    return;
-  }
   if (mParams->skipOutlierFiltering || validateTrack(trackData, params, clusterResiduals)) {
-    // track is good
-    int nClValidated = 0;
-    int iRow = 0;
-    for (unsigned int iCl = 0; iCl < clusterResiduals.size(); ++iCl) {
+    // track is good, store TPC part
+
+    int nClValidated = 0, iRow = 0;
+    unsigned int iCl = 0;
+    for (iCl = 0; iCl < clusterResiduals.size(); ++iCl) {
       iRow += clusterResiduals[iCl].dRow;
-      if (params.flagRej[iCl]) {
-        // skip masked cluster residual
+      if (iRow < param::NPadRows && params.flagRej[iCl]) { // skip masked cluster residual
         continue;
       }
       ++nClValidated;
@@ -720,9 +744,94 @@ void TrackInterpolation::extrapolateTrack(int iSeed)
       }
     }
     trackData.clIdx.setEntries(nClValidated);
+
+    // do we have TRD residuals to add?
+    int iSeedFull = mParentID[iSeed] == -1 ? iSeed : mParentID[iSeed];
+    auto gidFull = mGIDs[iSeedFull];
+    const auto& gidTableFull = mGIDtables[iSeedFull];
+    bool stopPropagation = false;
+    if (gidTableFull[GTrackID::TRD].isIndexSet()) {
+      const auto& trkTRD = mRecoCont->getITSTPCTRDTrack<o2::trd::TrackTRD>(gidTableFull[GTrackID::ITSTPCTRD]);
+      for (int iLayer = 0; iLayer < o2::trd::constants::NLAYER; iLayer++) {
+        std::array<float, 2> trkltTRDYZ{};
+        int res = processTRDLayer(trkTRD, iLayer, trkWork, &trkltTRDYZ);
+        if (res == 0) { // no traklet on this layer
+          continue;
+        }
+        if (res < 0) { // failed to reach this layer
+          stopPropagation = true;
+          break;
+        }
+
+        float tgPhi = trkWork.getSnp() / std::sqrt((1.f - trkWork.getSnp()) * (1.f + trkWork.getSnp()));
+        auto dy = trkltTRDYZ[0] - trkWork.getY();
+        auto dz = trkltTRDYZ[1] - trkWork.getZ();
+        const auto sec = clusterResiduals[iCl].sec;
+        if ((std::abs(dy) < param::MaxResid) && (std::abs(dz) < param::MaxResid) && (std::abs(trkWork.getY()) < param::MaxY) && (std::abs(trkWork.getZ()) < param::MaxZ) && (std::abs(tgPhi) < param::MaxTgSlp)) {
+          mClRes.emplace_back(dy, dz, tgPhi, trkWork.getY(), trkWork.getZ(), 160 + iLayer, o2::math_utils::angle2Sector(trkWork.getAlpha()));
+          trackData.nTrkltsTRD++;
+          trackData.nExtDetResid++;
+          // NEED TO STORE: TRD X
+        }
+      }
+    }
+
+    // do we have TOF residual to add?
+    while (gidTableFull[GTrackID::TOF].isIndexSet() && !stopPropagation) {
+      const auto& clTOF = mRecoCont->getTOFClusters()[gidTableFull[GTrackID::TOF]];
+      const float clTOFAlpha = o2::math_utils::sector2Angle(clTOF.getCount());
+      if (trkWork.getAlpha() != clTOFAlpha && !trkWork.rotate(clTOFAlpha)) {
+        LOG(debug) << "Failed to rotate into TOF cluster sector frame";
+        stopPropagation = true;
+        break;
+      }
+      if (!propagator->PropagateToXBxByBz(trkWork, clTOF.getX(), mParams->maxSnp, mParams->maxStep, mMatCorr)) {
+        LOG(debug) << "Failed final propagation to TOF radius";
+        break;
+      }
+
+      float tgPhi = trkWork.getSnp() / std::sqrt((1.f - trkWork.getSnp()) * (1.f + trkWork.getSnp()));
+      auto dy = clTOF.getY() - trkWork.getY();
+      auto dz = clTOF.getZ() - trkWork.getZ();
+      if ((std::abs(dy) < param::MaxResid) && (std::abs(dz) < param::MaxResid) && (std::abs(trkWork.getY()) < param::MaxY) && (std::abs(trkWork.getZ()) < param::MaxZ) && (std::abs(tgPhi) < param::MaxTgSlp)) {
+        mClRes.emplace_back(dy, dz, tgPhi, trkWork.getY(), trkWork.getZ(), 170, clTOF.getCount());
+        trackData.clAvailTOF = 1;
+        trackData.nExtDetResid++;
+        // NEED TO STORE: TOF X
+      }
+      break;
+    }
+
+    // add ITS residuals
+    while (!stopPropagation) {
+      auto trkWorkITS = trackData.par; // this is ITS outer param
+      auto nClITS = trkITS.getNumberOfClusters();
+      auto clEntryITS = trkITS.getFirstClusterEntry();
+      o2::its::GeometryTGeo* geom = GeometryTGeo::Instance();
+      for (int iCl = 0; iCl < nCl; iCl++) { // clusters are stored from outer to inner layers
+        const auto& clsITS = mITSClustersArray[mITSTrackClusIdx[clEntry + iCl]];
+        int chip = clsITS.getSensorID();
+        float chipX, chipAlpha;
+        geom->getSensorXAlphaRefPlane(clsITS.getSensorID(), chipX, chipAlpha);
+        if (!trkWorkITS.rotate(chipAlpha) || !propagator->PropagateToXBxByBz(trkWorkITS, chipX, mParams->maxSnp, mParams->maxStep, mMatCorr)) {
+          LOGP(debug, "Failed final propagation to ITS X={} alpha={}", chipX, chipAlpha);
+          stopPropagation = true;
+          break;
+        }
+        float tgPhi = trkWorkITS.getSnp() / std::sqrt((1.f - trkWorkITS.getSnp()) * (1.f + trkWorkITS.getSnp()));
+        auto dy = clsITS.getY() - trkWorkITS.getY();
+        auto dz = clsITS.getZ() - trkWorkITS.getZ();
+        if ((std::abs(dy) < param::MaxResid) && (std::abs(dz) < param::MaxResid) && (std::abs(trkWorkITS.getY()) < param::MaxY) && (std::abs(trkWorkITS.getZ()) < param::MaxZ) && (std::abs(tgPhi) < param::MaxTgSlp)) {
+          mClRes.emplace_back(dy, dz, tgPhi, trkWorkITS.getY(), trkWorkITS.getZ(), 180 + geom->getLayer(clsITS.getSensorID()), 0);
+          trackData.nExtDetResid++;
+          // NEED TO STORE: ITS X, alpha
+        }
+      }
+      break;
+    }
     mTrackData.push_back(std::move(trackData));
     mGIDsSuccess.push_back(mGIDs[iSeed]);
-    mTrackDataCompact.emplace_back(mClRes.size() - nClValidated, nClValidated, mGIDs[iSeed].getSource());
+    mTrackDataCompact.emplace_back(mClRes.size() - nClValidated, nClValidated, mGIDs[iSeed].getSource()); // RS Check
     if (mDumpTrackPoints) {
       (*trackDataExtended).clIdx.setEntries(nClValidated);
       mTrackDataExtended.push_back(std::move(*trackDataExtended));
