MatchCosmics: energy-loss sign along the muon's flight; optional energy-loss sign in the TPC refit

The propagator applies the energy-loss correction along the direction of the track parameters (loss forward, gain backward), and
refitWinners fits a cosmic from the bottom leg to the top one: after the inward refit of the bottom leg its parameters are inverted
and carried up through the top leg, i.e. they describe a particle moving upward. A cosmic muon flies downward, from the top leg to
the bottom leg, so only the inward refit of the bottom leg (backward along its parameters = backward along the flight) had the
physical sign; the propagation of the inverted bottom leg to the top leg, the outward refit of the top leg, the inward refit of the
top leg, its propagation to the bottom parameter and the propagation of the top seeds to the DCA had the opposite one. The sign is
now imposed along the flight in all of them, with the propagator's signCorr argument.

GPUO2InterfaceRefit / GPUTrackingRefit::RefitTrackAsTrackParCov (TrackParCov + cluster references) get the analogous optional last
argument eLossSign, passed to their propagations (0 = along the parameters as before, +1 gain, -1 loss; flipped on the optional
way to the reference X when that runs against the refit direction); the default keeps the behaviour for all other users.

MC truth (cosmic MC, 1073 muons crossing the inner region; true parameters at the bottom TPC inner radius propagated back through
the closest approach to the top one, compared to the true momentum there, q16 / q50 / q84 in MeV): true loss between the two
14.5 / 16.1 / 19.9 MeV (median for |z| at the closest approach 0-25 / 25-50 / 50-75 cm); predicted - true with the sign of the
upward fit -14.5 / -15.2 / -18.6 MeV (as without material: the gain inward and the loss outward cancel), with the sign along the
flight +1.2 / +2.1 / +1.0 MeV (spread q16-q84 about -2 ... +9 MeV). The matched cosmics hardly change because the two halves'
errors cancel in their average: cosmic MC 7003 -> 7006 correct pairs (0 wrong), chi2Match median unchanged; cosmics run 562658
358 cosmics before and after; refitted q/pt within +-0.05 %. The halves themselves (e.g. the momentum at the TPC entrance) change
by the full loss.

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

Author: matthias-kleiner

Detectors/GlobalTracking/src/MatchCosmics.cxx

@@ -45,6 +45,15 @@ using namespace o2::globaltracking;
 using GTrackID = o2d::GlobalTrackID;
 using MatCorrType = o2::base::Propagator::MatCorrType;
 
+namespace
+{
+// energy-loss sign of the propagations: the propagator applies the loss along the direction of the track parameters, and the refit runs
+// from the bottom leg up through the top leg (parameters of an upward-moving particle); the cosmic muon flies from the top leg to the
+// bottom leg, so the sign is imposed along its flight
+constexpr int ELossGain = 1;
+constexpr int ELossLoss = -1;
+} // namespace
+
 //________________________________________________________
 void MatchCosmics::process(const o2::globaltracking::RecoContainer& data)
 {
@@ -64,7 +73,12 @@ void MatchCosmics::process(const o2::globaltracking::RecoContainer& data)
   for (int i = 0; i < ntr; i++) {
     auto& trc = mSeeds[i];
     if (trc.matchID != Reject) {
-      if (!prop->propagateToDCABxByBz(v, trc, mMatchParams->maxStep, mMatchParams->matCorr)) {
+      // a cosmic muon flies downward: along a leg whose outward direction points up (top leg) the inward propagation follows the flight,
+      // so the energy is lost; along the bottom leg it goes back in the flight, so the energy is gained
+      std::array<float, 3> momentum{};
+      trc.getPxPyPzGlo(momentum);
+      const int eLossSign = momentum[1] > 0.f ? ELossLoss : ELossGain;
+      if (!prop->propagateToDCABxByBz(v, trc, mMatchParams->maxStep, mMatchParams->matCorr, nullptr, nullptr, eLossSign)) {
         trc.matchID = Reject; // reject track
         continue;
       }
@@ -154,7 +168,7 @@ void MatchCosmics::refitWinners(const o2::globaltracking::RecoContainer& data)
     }
   }
 
-  auto refitITSTrack = [this, &data, &itsTracksROF, &itsClusters](o2::track::TrackParCov& trFit, GTrackID gidx, float& chi2, bool inward = false) {
+  auto refitITSTrack = [this, &data, &itsTracksROF, &itsClusters](o2::track::TrackParCov& trFit, GTrackID gidx, float& chi2, bool inward, int eLossSign) {
     const auto& itsTrOrig = data.getITSTrack(gidx);
     int nclRefit = 0, ncl = itsTrOrig.getNumberOfClusters(), rof = itsTracksROF[gidx.getIndex()];
     const auto& itsTrackClusRefs = data.getITSTracksClusterRefs();
@@ -170,7 +184,7 @@ void MatchCosmics::refitWinners(const o2::globaltracking::RecoContainer& data)
     for (int icl = from; icl != to; icl += step) { // ITS clusters are referred in layer decreasing order
       const auto& clus = itsClusters[itsTrackClusRefs[clEntry + icl]];
       float alpha = geomITS->getSensorRefAlpha(clus.getSensorID()), x = clus.getX();
-      if (!trFit.rotate(alpha) || !propagator->propagateToX(trFit, x, propagator->getNominalBz(), this->mMatchParams->maxSnp, this->mMatchParams->maxStep, this->mMatchParams->matCorr)) {
+      if (!trFit.rotate(alpha) || !propagator->propagateToX(trFit, x, propagator->getNominalBz(), this->mMatchParams->maxSnp, this->mMatchParams->maxStep, this->mMatchParams->matCorr, nullptr, eLossSign)) {
         break;
       }
       chi2 += trFit.getPredictedChi2(clus);
@@ -221,7 +235,7 @@ void MatchCosmics::refitWinners(const o2::globaltracking::RecoContainer& data)
       if (!mFieldON) {
         trCosm.setQ2Pt(-o2::track::kMostProbablePt);
       }
-      int retVal = tpcRefitter->RefitTrackAsTrackParCov(trCosm, tpcTrOrig.getClusterRef(), t0 * tpcTBinMUSInv, &chi2, false, false); // inward refit, reset
+      int retVal = tpcRefitter->RefitTrackAsTrackParCov(trCosm, tpcTrOrig.getClusterRef(), t0 * tpcTBinMUSInv, &chi2, false, false, ELossGain); // inward refit, reset
       if (retVal < 0) {                                                                                                             // refit failed
         LOG(debug) << "Inward refit of btm TPC track failed.";
         continue;
@@ -247,7 +261,7 @@ void MatchCosmics::refitWinners(const o2::globaltracking::RecoContainer& data)
     }
     trCosm.invert();
     if (!trCosm.rotate(mSeeds[poolEntryID[top]].getAlpha()) ||
-        !o2::base::Propagator::Instance()->PropagateToXBxByBz(trCosm, mSeeds[poolEntryID[top]].getX(), mMatchParams->maxSnp, mMatchParams->maxStep, mMatchParams->matCorr)) {
+        !o2::base::Propagator::Instance()->PropagateToXBxByBz(trCosm, mSeeds[poolEntryID[top]].getX(), mMatchParams->maxSnp, mMatchParams->maxStep, mMatchParams->matCorr, nullptr, ELossGain)) {
       LOG(debug) << "Rotation/propagation of btm-track to top-track frame failed.";
       continue;
     }
@@ -260,7 +274,7 @@ void MatchCosmics::refitWinners(const o2::globaltracking::RecoContainer& data)
 
     // is there ITS sub-track?
     if (gidxListTop[GTrackID::ITS].isIndexSet()) {
-      auto nclfit = refitITSTrack(trCosm, gidxListTop[GTrackID::ITS], chi2, false);
+      auto nclfit = refitITSTrack(trCosm, gidxListTop[GTrackID::ITS], chi2, false, ELossGain);
       if (nclfit < 0) {
         continue;
       }
@@ -273,13 +287,13 @@ void MatchCosmics::refitWinners(const o2::globaltracking::RecoContainer& data)
       if (trCosm.getX() * trCosm.getX() + trCosm.getY() * trCosm.getY() <= o2::constants::geom::XTPCInnerRef * o2::constants::geom::XTPCInnerRef) {
         float xtogo = 0;
         if (!trCosm.getXatLabR(o2::constants::geom::XTPCInnerRef, xtogo, mBz, o2::track::DirOutward) ||
-            !o2::base::Propagator::Instance()->PropagateToXBxByBz(trCosm, xtogo, mMatchParams->maxSnp, mMatchParams->maxStep, mMatchParams->matCorr)) {
+            !o2::base::Propagator::Instance()->PropagateToXBxByBz(trCosm, xtogo, mMatchParams->maxSnp, mMatchParams->maxStep, mMatchParams->matCorr, nullptr, ELossGain)) {
           LOG(debug) << "Propagation to inner TPC boundary X=" << xtogo << " failed";
           continue;
         }
       }
       const auto& tpcTrOrig = data.getTPCTrack(gidxListTop[GTrackID::TPC]);
-      int retVal = tpcRefitter->RefitTrackAsTrackParCov(trCosm, tpcTrOrig.getClusterRef(), t0 * tpcTBinMUSInv, &chi2, true, false); // outward refit, no reset
+      int retVal = tpcRefitter->RefitTrackAsTrackParCov(trCosm, tpcTrOrig.getClusterRef(), t0 * tpcTBinMUSInv, &chi2, true, false, ELossGain); // outward refit, no reset
       if (retVal < 0) {                                                                                                             // refit failed
         LOG(debug) << "Outward refit of top TPC track failed.";
         continue;
@@ -293,23 +307,23 @@ void MatchCosmics::refitWinners(const o2::globaltracking::RecoContainer& data)
     auto trCosmTop = outerLegs[top];
     if (gidxListTop[GTrackID::TPC].isIndexSet()) { // inward refit in TPC
       const auto& tpcTrOrig = data.getTPCTrack(gidxListTop[GTrackID::TPC]);
-      int retVal = tpcRefitter->RefitTrackAsTrackParCov(trCosmTop, tpcTrOrig.getClusterRef(), t0 * tpcTBinMUSInv, &chi2Dummy, false, true); // inward refit, reset
+      int retVal = tpcRefitter->RefitTrackAsTrackParCov(trCosmTop, tpcTrOrig.getClusterRef(), t0 * tpcTBinMUSInv, &chi2Dummy, false, true, ELossLoss); // inward refit, reset
       if (retVal < 0) {                                                                                                                     // refit failed
-        LOG(debug) << "Outward refit of top TPC track failed.";
+        LOG(debug) << "Inward refit of top TPC track failed.";
         continue;
       } // inward refit in TPC
     }
     // is there ITS sub-track ?
     if (gidxListTop[GTrackID::ITS].isIndexSet()) {
-      auto nclfit = refitITSTrack(trCosmTop, gidxListTop[GTrackID::ITS], chi2Dummy, true);
+      auto nclfit = refitITSTrack(trCosmTop, gidxListTop[GTrackID::ITS], chi2Dummy, true, ELossLoss);
       if (nclfit < 0) {
         continue;
       }
       nclTot += nclfit;
     } // ITS refit
     // propagate to bottom param
     if (!trCosmTop.rotate(trCosmBtm.getAlpha()) ||
-        !o2::base::Propagator::Instance()->PropagateToXBxByBz(trCosmTop, trCosmBtm.getX(), mMatchParams->maxSnp, mMatchParams->maxStep, mMatchParams->matCorr)) {
+        !o2::base::Propagator::Instance()->PropagateToXBxByBz(trCosmTop, trCosmBtm.getX(), mMatchParams->maxSnp, mMatchParams->maxStep, mMatchParams->matCorr, nullptr, ELossLoss)) {
       LOG(debug) << "Rotation/propagation of top-track to bottom-track frame failed.";
       continue;
     }

GPU/GPUTracking/Interface/GPUO2InterfaceRefit.cxx

@@ -136,7 +136,7 @@ void GPUO2InterfaceRefit::updateCalib(const TPCFastTransformPOD* trans, float bz
 int32_t GPUO2InterfaceRefit::RefitTrackAsGPU(o2::tpc::TrackTPC& trk, bool outward, bool resetCov) { return mRefit->RefitTrackAsGPU(trk, outward, resetCov); }
 int32_t GPUO2InterfaceRefit::RefitTrackAsTrackParCov(o2::tpc::TrackTPC& trk, bool outward, bool resetCov) { return mRefit->RefitTrackAsTrackParCov(trk, outward, resetCov); }
 int32_t GPUO2InterfaceRefit::RefitTrackAsGPU(o2::track::TrackParCov& trk, const o2::tpc::TrackTPCClusRef& clusRef, float time0, float* chi2, bool outward, bool resetCov) { return mRefit->RefitTrackAsGPU(trk, clusRef, time0, chi2, outward, resetCov); }
-int32_t GPUO2InterfaceRefit::RefitTrackAsTrackParCov(o2::track::TrackParCov& trk, const o2::tpc::TrackTPCClusRef& clusRef, float time0, float* chi2, bool outward, bool resetCov) { return mRefit->RefitTrackAsTrackParCov(trk, clusRef, time0, chi2, outward, resetCov); }
+int32_t GPUO2InterfaceRefit::RefitTrackAsTrackParCov(o2::track::TrackParCov& trk, const o2::tpc::TrackTPCClusRef& clusRef, float time0, float* chi2, bool outward, bool resetCov, int32_t eLossSign) { return mRefit->RefitTrackAsTrackParCov(trk, clusRef, time0, chi2, outward, resetCov, eLossSign); }
 void GPUO2InterfaceRefit::setIgnoreErrorsAtTrackEnds(bool v) { mRefit->mIgnoreErrorsOnTrackEnds = v; }
 void GPUO2InterfaceRefit::setTrackReferenceX(float v) { mParam->rec.tpc.trackReferenceX = v; }
 

GPU/GPUTracking/Interface/GPUO2InterfaceRefit.h

@@ -67,7 +67,8 @@ class GPUO2InterfaceRefit
   int32_t RefitTrackAsGPU(o2::tpc::TrackTPC& trk, bool outward = false, bool resetCov = false);
   int32_t RefitTrackAsTrackParCov(o2::tpc::TrackTPC& trk, bool outward = false, bool resetCov = false);
   int32_t RefitTrackAsGPU(o2::track::TrackParCov& trk, const o2::tpc::TrackTPCClusRef& clusRef, float time0, float* chi2 = nullptr, bool outward = false, bool resetCov = false);
-  int32_t RefitTrackAsTrackParCov(o2::track::TrackParCov& trk, const o2::tpc::TrackTPCClusRef& clusRef, float time0, float* chi2 = nullptr, bool outward = false, bool resetCov = false);
+  // eLossSign: energy-loss sign of the propagations, 0 along the track parameters (loss forward, gain backward), +1 gain, -1 loss
+  int32_t RefitTrackAsTrackParCov(o2::track::TrackParCov& trk, const o2::tpc::TrackTPCClusRef& clusRef, float time0, float* chi2 = nullptr, bool outward = false, bool resetCov = false, int32_t eLossSign = 0);
   void setTrackReferenceX(float v);
   void setIgnoreErrorsAtTrackEnds(bool v);
   void updateCalib(const o2::gpu::TPCFastTransformPOD* trans, float bzNominalGPU);

GPU/GPUTracking/Refit/GPUTrackingRefit.cxx

@@ -221,6 +221,7 @@ GPUd() int32_t GPUTrackingRefit::RefitTrack(T& trkX, bool outward, bool resetCov
   convertTrack<S, T, typename internal::refitTrackTypes<S>::propagator>(trk, trkX, prop, &TrackParCovChi2);
   int32_t begin = 0, count;
   float tOffset;
+  [[maybe_unused]] int32_t eLossSign = 0; // energy-loss sign of the TrackParCov propagations (0: from the direction, outward = loss)
   if constexpr (std::is_same_v<T, GPUTPCGMMergedTrack>) {
     count = trkX.NClusters();
     tOffset = trkX.GetParam().GetTOffset();
@@ -230,6 +231,7 @@ GPUd() int32_t GPUTrackingRefit::RefitTrack(T& trkX, bool outward, bool resetCov
   } else if constexpr (std::is_same_v<T, TrackParCovWithArgs>) {
     count = trkX.clusRef.getEntries();
     tOffset = trkX.time0;
+    eLossSign = trkX.eLossSign;
   } else {
     static_assert("Invalid template");
   }
@@ -357,7 +359,7 @@ GPUd() int32_t GPUTrackingRefit::RefitTrack(T& trkX, bool outward, bool resetCov
         IgnoreErrors(trk.getSnp());
         return -1;
       }
-      if (!prop->PropagateToXBxByBz(trk, x, constants::MAX_SIN_PHI_LOW)) {
+      if (!prop->PropagateToXBxByBz(trk, x, constants::MAX_SIN_PHI_LOW, Propagator::MAX_STEP, Propagator::MatCorrType::USEMatCorrLUT, nullptr, eLossSign)) {
         IgnoreErrors(trk.getSnp());
         return -2;
       }
@@ -401,11 +403,13 @@ GPUd() int32_t GPUTrackingRefit::RefitTrack(T& trkX, bool outward, bool resetCov
     constexpr float kDeg2Rad = M_PI / 180.f;
     constexpr float kSectAngle = 2 * M_PI / 18.f;
     if (mPparam->rec.tpc.trackReferenceX <= 500) {
-      if (prop->PropagateToXBxByBz(trk, mPparam->rec.tpc.trackReferenceX)) {
+      // a forced energy-loss sign holds along the refit direction: the way to the reference X can run against it
+      const int32_t eLossSignRef = (mPparam->rec.tpc.trackReferenceX > trk.getX()) == outward ? eLossSign : -eLossSign;
+      if (prop->PropagateToXBxByBz(trk, mPparam->rec.tpc.trackReferenceX, Propagator::MAX_SIN_PHI, Propagator::MAX_STEP, Propagator::MatCorrType::USEMatCorrLUT, nullptr, eLossSignRef)) {
         if (CAMath::Abs(trk.getY()) > trk.getX() * CAMath::Tan(kSectAngle / 2.f)) {
           float newAlpha = trk.getAlpha() + CAMath::Round(CAMath::ATan2(trk.getY(), trk.getX()) / kDeg2Rad / 20.f) * kSectAngle;
           GPUTPCGMTrackParam::NormalizeAlpha(newAlpha);
-          trk.rotate(newAlpha) && prop->PropagateToXBxByBz(trk, mPparam->rec.tpc.trackReferenceX);
+          trk.rotate(newAlpha) && prop->PropagateToXBxByBz(trk, mPparam->rec.tpc.trackReferenceX, Propagator::MAX_SIN_PHI, Propagator::MAX_STEP, Propagator::MatCorrType::USEMatCorrLUT, nullptr, eLossSignRef);
         }
       }
     }

GPU/GPUTracking/Refit/GPUTrackingRefit.h

@@ -73,15 +73,16 @@ class GPUTrackingRefit
     const o2::tpc::TrackTPCClusRef& clusRef;
     float time0;
     float* chi2;
+    int32_t eLossSign = 0; // energy-loss sign of the TrackParCov propagations: 0 along the track parameters (loss forward), +1 gain, -1 loss
   };
   GPUd() int32_t RefitTrackAsGPU(o2::track::TrackParCov& trk, const o2::tpc::TrackTPCClusRef& clusRef, float time0, float* chi2 = nullptr, bool outward = false, bool resetCov = false)
   {
     TrackParCovWithArgs x{trk, clusRef, time0, chi2};
     return RefitTrack<TrackParCovWithArgs, GPUTPCGMTrackParam>(x, outward, resetCov);
   }
-  GPUd() int32_t RefitTrackAsTrackParCov(o2::track::TrackParCov& trk, const o2::tpc::TrackTPCClusRef& clusRef, float time0, float* chi2 = nullptr, bool outward = false, bool resetCov = false)
+  GPUd() int32_t RefitTrackAsTrackParCov(o2::track::TrackParCov& trk, const o2::tpc::TrackTPCClusRef& clusRef, float time0, float* chi2 = nullptr, bool outward = false, bool resetCov = false, int32_t eLossSign = 0)
   {
-    TrackParCovWithArgs x{trk, clusRef, time0, chi2};
+    TrackParCovWithArgs x{trk, clusRef, time0, chi2, eLossSign};
     return RefitTrack<TrackParCovWithArgs, o2::track::TrackParCov>(x, outward, resetCov);
   }