Unbinned residuals with ITS refit and PV point added

With scdcalib.refitITS=true (default) the ITS track outer param will be refitted from scratch using
stable Kalman filter staring from the inner param (and using it as as linearization point) and imposing
PID of the global track (if any).
Note that the reconstruction is done with extra syst.errors on the ITS clusters (20 microns): they can be
added also for this refit with usual ITSCATrackerParam.sysErrY2... settings (though this is not necessary).

The track residuals wrt the PV are added to the unbinned residuals with identified of padrow 190.
Note that PV residuals interpretation differs from the rest: they are provided at the PCA of the track to PV,
with PV rotated to the frame of the track. The X of the vertex in this frame is stored in the channel slot
mapped from [-0.5:0.5] to short; the alpha of the track frame is stored as the angle in [-pi : pi] mapped to short:

auto dy = yv - trkAtPCA.getY(), auto dz = zv - trkAtPCA.getZ();
short compXV = static_cast<short>(xv * 0x7fff / param::MaxVtxX); // MaxVtxX = 0.5 cm
mClRes.emplace_back(dy, dz, trkAtPCA.getAlpha() / TMath::Pi(), trkWorkITS.getY(), trkWorkITS.getZ(), 190, -1, compXV);

Author: shahor02

Detectors/TPC/calibration/SpacePoints/CMakeLists.txt

@@ -8,6 +8,7 @@
 # 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.
+# add_compile_options(-O0 -g -fPIC)
 
 o2_add_library(SpacePoints
                SOURCES src/SpacePointsCalibParam.cxx

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

@@ -41,6 +41,7 @@ struct SpacePointsCalibConfParam : public o2::conf::ConfigurableParamHelper<Spac
   int minTOFTRDPVContributors = 2;      ///< min contributors from TRD or TOF (fast detectors) to consider tracks of this PV
   bool ignoreNonPVContrib = true;       ///< flag if tracks which did not contribute to the PV should be ignored or not
   bool enableTrackDownsampling = false; ///< flag if track sampling shall be enabled or not
+  bool refitITS = true;                 ///< refit ITS tracks with PID attached to the seed
   float tsalisThreshold = 1.f;          ///< in case the sampling functions returns a value smaller than this the track is discarded (1. means no downsampling)
 
   // other settings for track interpolation

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

@@ -58,6 +58,7 @@ static constexpr float MaxY = 50.f;      ///< max value for y position (sector c
 static constexpr float MaxZ = 300.f;     ///< max value for z position
 static constexpr float MaxTgSlp = 1.f;   ///< max value for phi (from snp, converted to tangens)
 static constexpr float MaxTRDSlope = 5.; ///< max value for the TRD tracklet getDy
+static constexpr float MaxVtxX = 0.5;    ///< max X value of the PV rotated to the frame of track at its DCA to PV
 // miscellaneous
 static constexpr float sEps = 1e-6f; ///< small number for float comparisons
 

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

@@ -239,6 +239,8 @@ class TrackInterpolation
   /// Reset cache and output vectors
   void reset();
 
+  // refit ITS track taking PID (unless already refitted) from the seed and reassign to the seed
+  bool refITSTrack(o2::dataformats::GlobalTrackID, int iSeed);
   // -------------------------------------- outlier rejection --------------------------------------------------
 
   /// Validates the given input track and its residuals
@@ -328,6 +330,7 @@ class TrackInterpolation
   std::vector<o2::dataformats::GlobalTrackID> mGIDs{};                      ///< GIDs of input tracks
   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<int> mTrackPVID{};                                            ///< track vertex index (if any)
   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
@@ -337,6 +340,7 @@ class TrackInterpolation
   // ITS specific input only needed for debugging
   gsl::span<const int> mITSTrackClusIdx;                    ///< input ITS track cluster indices span
   std::vector<o2::BaseCluster<float>> mITSClustersArray;    ///< ITS clusters created in run() method from compact clusters
+  std::vector<int> mITSRefitSeedID;                         ///< seed ID first using refitted ITS track
   const o2::itsmft::TopologyDictionary* mITSDict = nullptr; ///< cluster patterns dictionary
 
   // output

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

@@ -221,6 +221,9 @@ void TrackInterpolation::prepareInputTrackSample(const o2::globaltracking::RecoC
   int nv = vtxRefs.size() - 1;
   GTrackID::mask_t allowedSources = GTrackID::getSourcesMask("ITS-TPC,ITS-TPC-TRD,ITS-TPC-TOF,ITS-TPC-TRD-TOF") & mSourcesConfigured;
   constexpr std::array<int, 3> SrcFast = {int(GTrackID::ITSTPCTRD), int(GTrackID::ITSTPCTOF), int(GTrackID::ITSTPCTRDTOF)};
+  if (mParams->refitITS) {
+    mITSRefitSeedID.resize(mRecoCont->getITSTracks().size(), -1);
+  }
 
   for (int iv = 0; iv < nv; iv++) {
     LOGP(debug, "processing PV {} of {}", iv, nv);
@@ -281,6 +284,7 @@ void TrackInterpolation::prepareInputTrackSample(const o2::globaltracking::RecoC
         mGIDtables.push_back(gidTable);
         mTrackTimes.push_back(pv.getTimeStamp().getTimeStamp());
         mTrackIndices[mTrackTypes[vid.getSource()]].push_back(nTrackSeeds++);
+        mTrackPVID.push_back(iv);
       }
     }
   }
@@ -360,13 +364,13 @@ void TrackInterpolation::process()
     if (mParams->enableTrackDownsampling && !isTrackSelected(mSeeds[seedIndex])) {
       continue;
     }
-
     auto addPart = [this, seedIndex](GTrackID::Source src) {
       this->mGIDs.push_back(this->mGIDtables[seedIndex][src]);
       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
+      this->mTrackPVID.push_back(this->mTrackPVID[seedIndex]);
     };
 
     GTrackID::mask_t partsAdded;
@@ -450,9 +454,12 @@ void TrackInterpolation::interpolateTrack(int iSeed)
       (*trackDataExtended).clsITS.push_back(clsITS);
     }
   }
+  if (mParams->refitITS && !refITSTrack(gidTable[GTrackID::ITS], iSeed)) {
+    return;
+  }
   trackData.gid = mGIDs[iSeed];
   trackData.par = mSeeds[iSeed];
-  auto& trkWork = mSeeds[iSeed];
+  auto trkWork = mSeeds[iSeed];
   o2::track::TrackPar trkInner{trkWork};
   // reset the cache array (sufficient to set cluster available to zero)
   for (auto& elem : mCache) {
@@ -734,6 +741,27 @@ void TrackInterpolation::interpolateTrack(int iSeed)
             trackData.nExtDetResid++;
           }
         }
+        if (!stopPropagation) { // add residual to PV
+          const auto& pv = mRecoCont->getPrimaryVertices()[mTrackPVID[iSeed]];
+          o2::math_utils::Point3D<float> vtx{pv.getX(), pv.getY(), pv.getZ()};
+          if (!propagator->propagateToDCA(vtx, trkWorkITS, mBz, mParams->maxStep, mMatCorr)) {
+            LOGP(debug, "Failed propagation to DCA to PV ({} {} {}), {}", pv.getX(), pv.getY(), pv.getZ(), trkWorkITS.asString());
+            stopPropagation = true;
+            break;
+          }
+          // rotate PV to the track frame
+          float sn, cs, alpha = trkWorkITS.getAlpha();
+          math_utils::detail::bringToPMPi(alpha);
+          math_utils::detail::sincos<float>(alpha, sn, cs);
+          float xv = vtx.X() * cs + vtx.Y() * sn, yv = -vtx.X() * sn + vtx.Y() * cs, zv = vtx.Z();
+          auto dy = yv - trkWorkITS.getY();
+          auto dz = zv - 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) && abs(xv) < param::MaxVtxX) {
+            short compXV = static_cast<short>(xv * 0x7fff / param::MaxVtxX);
+            mClRes.emplace_back(dy, dz, alpha / TMath::Pi(), trkWorkITS.getY(), trkWorkITS.getZ(), 190, -1, compXV);
+            trackData.nExtDetResid++;
+          }
+        }
         break;
       }
     }
@@ -826,6 +854,9 @@ void TrackInterpolation::extrapolateTrack(int iSeed)
       (*trackDataExtended).clsITS.push_back(clsITS);
     }
   }
+  if (mParams->refitITS && !refITSTrack(gidTable[GTrackID::ITS], iSeed)) {
+    return;
+  }
   trackData.gid = mGIDs[iSeed];
   trackData.par = mSeeds[iSeed];
 
@@ -987,7 +1018,7 @@ void TrackInterpolation::extrapolateTrack(int iSeed)
           int chip = cls.getSensorID();
           float chipX, chipAlpha;
           geom->getSensorXAlphaRefPlane(cls.getSensorID(), chipX, chipAlpha);
-          if (!trkWorkITS.rotate(chipAlpha) || !propagator->PropagateToXBxByBz(trkWorkITS, chipX, mParams->maxSnp, mParams->maxStep, mMatCorr)) {
+          if (!trkWorkITS.rotate(chipAlpha) || !propagator->propagateToX(trkWorkITS, chipX, mBz, mParams->maxSnp, mParams->maxStep, mMatCorr)) {
             LOGP(debug, "Failed final propagation to ITS X={} alpha={}", chipX, chipAlpha);
             stopPropagation = true;
             break;
@@ -1000,6 +1031,27 @@ void TrackInterpolation::extrapolateTrack(int iSeed)
             trackData.nExtDetResid++;
           }
         }
+        if (!stopPropagation) { // add residual to PV
+          const auto& pv = mRecoCont->getPrimaryVertices()[mTrackPVID[iSeed]];
+          o2::math_utils::Point3D<float> vtx{pv.getX(), pv.getY(), pv.getZ()};
+          if (!propagator->propagateToDCA(vtx, trkWorkITS, mBz, mParams->maxStep, mMatCorr)) {
+            LOGP(debug, "Failed propagation to DCA to PV ({} {} {}), {}", pv.getX(), pv.getY(), pv.getZ(), trkWorkITS.asString());
+            stopPropagation = true;
+            break;
+          }
+          // rotate PV to the track frame
+          float sn, cs, alpha = trkWorkITS.getAlpha();
+          math_utils::detail::bringToPMPi(alpha);
+          math_utils::detail::sincos<float>(alpha, sn, cs);
+          float xv = vtx.X() * cs + vtx.Y() * sn, yv = -vtx.X() * sn + vtx.Y() * cs, zv = vtx.Z();
+          auto dy = yv - trkWorkITS.getY();
+          auto dz = zv - 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) && abs(xv) < param::MaxVtxX) {
+            short compXV = static_cast<short>(xv * 0x7fff / param::MaxVtxX);
+            mClRes.emplace_back(dy, dz, alpha / TMath::Pi(), trkWorkITS.getY(), trkWorkITS.getZ(), 190, -1, compXV);
+            trackData.nExtDetResid++;
+          }
+        }
         break;
       }
     }
@@ -1403,6 +1455,8 @@ void TrackInterpolation::reset()
   mGIDtables.clear();
   mTrackTimes.clear();
   mSeeds.clear();
+  mITSRefitSeedID.clear();
+  mTrackPVID.clear();
 }
 
 //______________________________________________
@@ -1416,3 +1470,50 @@ void TrackInterpolation::setTPCVDrift(const o2::tpc::VDriftCorrFact& v)
     o2::tpc::TPCFastTransformHelperO2::instance()->updateCalibration(*mFastTransform, 0, 1.0, mTPCVDriftRef, mTPCDriftTimeOffsetRef);
   }
 }
+
+//______________________________________________
+bool TrackInterpolation::refITSTrack(o2::dataformats::GlobalTrackID gid, int seedID)
+{
+  // refit ITS track outwards taking PID (unless already refitted) from the seed and reassign to the seed
+  auto& seed = mSeeds[seedID];
+  int refitID = mITSRefitSeedID[gid.getIndex()];
+  if (refitID >= 0) { // track was already refitted
+    if (mSeeds[refitID].getPID() == seed.getPID()) {
+      seed = mSeeds[refitID];
+    }
+    return true;
+  }
+  const auto& trkITS = mRecoCont->getITSTrack(gid);
+  // fetch clusters
+  auto nCl = trkITS.getNumberOfClusters();
+  auto clEntry = trkITS.getFirstClusterEntry();
+  o2::track::TrackParCov track(trkITS); // start from the inner param
+  track.setPID(seed.getPID());
+  o2::track::TrackPar refLin(track); // and use it also as linearization reference
+  auto geom = o2::its::GeometryTGeo::Instance();
+  auto prop = o2::base::Propagator::Instance();
+  for (int iCl = nCl - 1; iCl >= 0; iCl--) { // clusters are stored from outer to inner layers
+    const auto& cls = mITSClustersArray[mITSTrackClusIdx[clEntry + iCl]];
+    int chip = cls.getSensorID();
+    float chipX, chipAlpha;
+    geom->getSensorXAlphaRefPlane(cls.getSensorID(), chipX, chipAlpha);
+    if (!track.rotate(chipAlpha, refLin, mBz)) {
+      LOGP(debug, "failed to rotate ITS tracks to alpha={} for the refit: {}", chipAlpha, track.asString());
+      return false;
+    }
+    if (!prop->propagateToX(track, refLin, cls.getX(), mBz, o2::base::PropagatorImpl<float>::MAX_SIN_PHI, o2::base::PropagatorImpl<float>::MAX_STEP, o2::base::PropagatorF::MatCorrType::USEMatCorrLUT)) {
+      LOGP(debug, "failed to propagate ITS tracks to X={}: {}", cls.getX(), track.asString());
+      return false;
+    }
+    std::array<float, 2> posTF{cls.getY(), cls.getZ()};
+    std::array<float, 3> covTF{cls.getSigmaY2(), cls.getSigmaYZ(), cls.getSigmaZ2()};
+    if (!track.update(posTF, covTF)) {
+      LOGP(debug, "failed to update ITS tracks by cluster ({},{})/({},{},{})", track.asString(), cls.getY(), cls.getZ(), cls.getSigmaY2(), cls.getSigmaYZ(), cls.getSigmaZ2());
+      return false;
+    }
+  }
+  seed = track;
+  // memorize that this ITS track was already refitted
+  mITSRefitSeedID[gid.getIndex()] = seedID;
+  return true;
+}