Implemented comments of David and Ruben

Author: tubagundem

Detectors/TPC/calibration/include/TPCCalibration/CalculatedEdx.h

@@ -213,7 +213,7 @@ class CalculatedEdx
   /// \return returns the number of rows where refit/propagation failed (row.propagationFailed) since the last resetDebugCounters(); with setRefit(), this only counts rows where the propagation fallback was also unable to recover the row
   long getNPropagationFailed() const { return mNPropagationFailed; }
 
-  /// \return returns the number of rows where setRefit()'s RefitTrackAsGPU() could not reach the row and the row's track state instead came from the propagation fallback since the last resetDebugCounters(); always 0 outside setRefit() mode
+  /// \return returns the number of rows of tracks where setRefit()'s RefitTrackAsGPU() failed and the row's track state instead came from propagating the original track since the last resetDebugCounters(); always 0 outside setRefit() mode
   long getNRefitFallback() const { return mNRefitFallback; }
 
   /// \return returns the number of rows gathered by gatherRowClusterData() (processed for refit/propagation) since the last resetDebugCounters()
@@ -400,7 +400,7 @@ class CalculatedEdx
     int stackNumber;
     StackID stackID;
     bool propagationFailed;               ///< true if refit/propagation to this row failed, or the resulting track param is NaN, and no fallback recovered it either
-    bool refitFellBack;                   ///< true if setRefit() mode's RefitTrackAsGPU() could not reach this row
+    bool refitFellBack;                   ///< true if setRefit() mode's refit of this track failed and the row's track state came from propagating the original track
     int missingClusters;                  ///< number of skipped rows since the previous entry in rowData (i.e. rowIndex - previous rowIndex - 1); same for every settings entry since rowOrder does not depend on the settings
     bool sameSectorAsPrevRow;             ///< true if this row's sector equals the previous entry in rowData's sector
     bool missingClusterGapDeadOrEdge;     ///< true if any of the missingClusters skipped row(s) would land on a dead channel or off the padrow edge
@@ -433,11 +433,13 @@ class CalculatedEdx
   /// \param sectorIndexOld sectorIndex of the previous entry appended to rowData (255 if this is the first row)
   /// \param occupancyROC per-region occupancy accumulator, updated in place
   /// \param rowData output per-row data; the new row is appended, and rowData.back() (if non-empty) is read as the previous row for the missing-cluster-gap check
-  /// \param refitAbandoned setRefit() mode only: false as long as RefitTrackAsGPU() keeps succeeding; the first
-  ///        time it fails for this track, set to true and stays true for the rest of the track. On that first
-  ///        failure, the propagation fallback resumes from the track's state just before the failed attempt
-  ///        (i.e. its state after the last successfully refit row), not from the track's pristine pre-loop state.
-  void gatherRowClusterDataForRow(o2::tpc::TrackTPC& track, const std::vector<o2::tpc::ClusterNative>& fragmentClusters, const std::vector<bool>& fragmentIsShared, unsigned char sectorIndex, unsigned char rowIndex, bool mergeable, unsigned char rowIndexOld, unsigned char sectorIndexOld, std::array<std::vector<unsigned int>, 4>& occupancyROC, std::vector<RowClusterData>& rowData, bool& refitAbandoned);
+  /// \param refitFailed setRefit() mode only: true if the single per-track refit in refitTrack() failed, in which case the rows are propagated from the track's original state
+  void gatherRowClusterDataForRow(o2::tpc::TrackTPC& track, const std::vector<o2::tpc::ClusterNative>& fragmentClusters, const std::vector<bool>& fragmentIsShared, unsigned char sectorIndex, unsigned char rowIndex, bool mergeable, unsigned char rowIndexOld, unsigned char sectorIndexOld, std::array<std::vector<unsigned int>, 4>& occupancyROC, std::vector<RowClusterData>& rowData, bool refitFailed);
+
+  /// setRefit() mode only: refit the whole track once, outward, so that it ends at the outermost cluster; the track is then propagated row by row in gatherRowClusterDataForRow()
+  /// \param track input track, mutated in place by the refit; left unchanged if the refit fails
+  /// \return true if the refit failed
+  bool refitTrack(o2::tpc::TrackTPC& track);
 
   /// geometrically propagate track (rotating into the row's sector frame first) to xPosition
   /// \return true if any of the three attempts succeeded (track left at xPosition); false if all failed (track left unchanged, at its state on entry)
@@ -479,7 +481,7 @@ class CalculatedEdx
   std::unordered_map<std::string, std::unique_ptr<o2::utils::TreeStreamRedirector>> mStreamers; ///< debug streamers, keyed by output file name so each debugRootFile gets its own tree
   long mDebugTrackIndex{-1};                                                                    ///< running index of the track being processed, written to the debug trees so per-cluster rows can be grouped back into tracks
   long mNPropagationFailed{0};                                                                  ///< number of rows where refit/propagation failed (and the fallback below, if applicable, also failed) since the last resetDebugCounters()
-  long mNRefitFallback{0};                                                                      ///< number of rows where setRefit()'s RefitTrackAsGPU() failed but the propagation fallback recovered the row, since the last resetDebugCounters()
+  long mNRefitFallback{0};                                                                      ///< number of rows of tracks where setRefit()'s RefitTrackAsGPU() failed but propagating the original track recovered the row, since the last resetDebugCounters()
   long mNRowsProcessed{0};                                                                      ///< number of rows gathered by gatherRowClusterData() since the last resetDebugCounters()
   std::vector<long> mNSubThresholdFilledPerSettings;                                            ///< number of row gaps filled as subthreshold clusters, per dEdxSettings list index, since the last resetDebugCounters()
 

Detectors/TPC/calibration/macro/calculatedEdx.C

@@ -326,7 +326,7 @@ void calculatedEdx(const std::string dir = ".",
     // per-event summary: refit/propagation failures and how many row gaps were filled as subthreshold clusters per settingsList entry
     const long nPropagationFailed = calcdEdx.getNPropagationFailed();
     const long nRowsProcessed = calcdEdx.getNRowsProcessed();
-    const auto& nSubThresholdFilledPerSettings = calcdEdx.getNSubThresholdFilledPerSettings();
+    const auto nSubThresholdFilledPerSettings = calcdEdx.getNSubThresholdFilledPerSettings(); // copy: resetDebugCounters() below clears the member
     calcdEdx.resetDebugCounters();
     std::string subThresholdBreakdown;
     for (size_t i = 0; i < nSubThresholdFilledPerSettings.size(); i++) {

Detectors/TPC/calibration/src/CalculatedEdx.cxx

@@ -261,7 +261,7 @@ void CalculatedEdx::gatherRowClusterData(o2::tpc::TrackTPC& track, std::vector<R
 {
   rowData.clear();
 
-  bool refitAbandoned = false;
+  const bool refitFailed = refitTrack(track);
 
   // handle same (sector, row) clusters
   std::vector<std::pair<unsigned char, unsigned char>> rowOrder;
@@ -298,7 +298,7 @@ void CalculatedEdx::gatherRowClusterData(o2::tpc::TrackTPC& track, std::vector<R
     }
 
     const bool mergeable = mergeableRows.count(rowKey) > 0;
-    gatherRowClusterDataForRow(track, fragmentClusters, fragmentIsShared, sectorIndex, rowIndex, mergeable, rowIndexOld, sectorIndexOld, occupancyROC, rowData, refitAbandoned);
+    gatherRowClusterDataForRow(track, fragmentClusters, fragmentIsShared, sectorIndex, rowIndex, mergeable, rowIndexOld, sectorIndexOld, occupancyROC, rowData, refitFailed);
     rowIndexOld = rowIndex;
     sectorIndexOld = sectorIndex;
   }
@@ -320,6 +320,21 @@ void CalculatedEdx::gatherRowClusterData(o2::tpc::TrackTPC& track, std::vector<R
   }
 }
 
+bool CalculatedEdx::refitTrack(o2::tpc::TrackTPC& track)
+{
+  if (!mRefit) {
+    return false;
+  }
+  // RefitTrackAsGPU() refits the whole track from its clusters (not just up to the current row) and may write into the track even if it fails, so keep the original state to fall back from
+  const o2::track::TrackParCov trackBeforeRefit = track;
+  // outward refit: the track ends at the outermost cluster, i.e. at the start of the cluster reference list that is traversed row by row
+  const bool failed = (mRefit->RefitTrackAsGPU(track, true, true) < 0) || std::isnan(track.getParam(1));
+  if (failed) {
+    static_cast<o2::track::TrackParCov&>(track) = trackBeforeRefit;
+  }
+  return failed;
+}
+
 bool CalculatedEdx::propagateTrackToX(o2::track::TrackParCov& track, float xPosition, unsigned char sectorIndex) const
 {
   const o2::track::TrackParCov trackBackup = track;
@@ -347,7 +362,7 @@ bool CalculatedEdx::propagateTrackToX(o2::track::TrackParCov& track, float xPosi
   return check;
 }
 
-void CalculatedEdx::gatherRowClusterDataForRow(o2::tpc::TrackTPC& track, const std::vector<o2::tpc::ClusterNative>& fragmentClusters, const std::vector<bool>& fragmentIsShared, unsigned char sectorIndex, unsigned char rowIndex, bool mergeable, unsigned char rowIndexOld, unsigned char sectorIndexOld, std::array<std::vector<unsigned int>, 4>& occupancyROC, std::vector<RowClusterData>& rowData, bool& refitAbandoned)
+void CalculatedEdx::gatherRowClusterDataForRow(o2::tpc::TrackTPC& track, const std::vector<o2::tpc::ClusterNative>& fragmentClusters, const std::vector<bool>& fragmentIsShared, unsigned char sectorIndex, unsigned char rowIndex, bool mergeable, unsigned char rowIndexOld, unsigned char sectorIndexOld, std::array<std::vector<unsigned int>, 4>& occupancyROC, std::vector<RowClusterData>& rowData, bool refitFailed)
 {
   RowClusterData row;
   row.sectorIndex = sectorIndex;
@@ -418,32 +433,9 @@ void CalculatedEdx::gatherRowClusterDataForRow(o2::tpc::TrackTPC& track, const s
   // get the x position of the track
   const float xPosition = Mapper::instance().getPadCentre(PadPos(rowIndex, 0)).X();
   bool check = true;
-  bool refitFellBack = false;
-  if (mRefit) {
-    if (!refitAbandoned) {
-      // snapshot the track's state as it stood before this row's refit attempt (i.e. after the previous row's
-      // successful refit) -- RefitTrackAsGPU() writes back into `track` even on failure, so on failure this is
-      // the most recent known-good state to fall back from, not the track's pristine pre-loop state
-      const o2::track::TrackParCov trackBeforeRefit = track;
-      // refit this track
-      mRefit->setTrackReferenceX(xPosition);
-      // RefitTrackAsGPU() returns < 0 when it fails; reachedReference is false when the fit succeeded but the final move-to-reference step could not reach xPosition, so both trigger the fallback.
-      bool reachedReference = true;
-      check = (mRefit->RefitTrackAsGPU(track, false, true, &reachedReference) < 0) ? false : reachedReference;
-      if (!check || std::isnan(track.getParam(1))) {
-        refitAbandoned = true;
-        static_cast<o2::track::TrackParCov&>(track) = trackBeforeRefit;
-        check = propagateTrackToX(track, xPosition, sectorIndex);
-        refitFellBack = check;
-      }
-    } else {
-      // already abandoned refit for this track (the previous row's RefitTrackAsGPU() failed)
-      // keep propagating incrementally from `track`'s current state
-      check = propagateTrackToX(track, xPosition, sectorIndex);
-      refitFellBack = check;
-    }
-  } else if (mPropagateTrack) {
-    // propagate this track to the plane X=xk (cm) in the field "b" (kG)
+  if (mRefit || mPropagateTrack) {
+    // mRefit: the track was refit once in gatherRowClusterData() and is only propagated from row to row here, if that refit failed, the propagation starts from the track's original state
+    // mPropagateTrack: propagate this track to the plane X=xk (cm) in the field "b" (kG)
     check = propagateTrackToX(track, xPosition, sectorIndex);
   } else if (mPropagateParams) {
     // propagate the params of the track instead of full propagation; same rollback rationale as mPropagateTrack above
@@ -458,7 +450,7 @@ void CalculatedEdx::gatherRowClusterDataForRow(o2::tpc::TrackTPC& track, const s
   }
 
   row.propagationFailed = (!check || std::isnan(track.getParam(1)));
-  row.refitFellBack = refitFellBack && !row.propagationFailed;
+  row.refitFellBack = mRefit && refitFailed && !row.propagationFailed;
   ++mNRowsProcessed;
   if (row.propagationFailed) {
     ++mNPropagationFailed;
@@ -509,7 +501,7 @@ void CalculatedEdx::gatherRowClusterData(o2::tpc::TrackTPC& track, const std::ve
 {
   rowData.clear();
 
-  bool refitAbandoned = false;
+  const bool refitFailed = refitTrack(track);
 
   // handle same (sector, row) clusters
   std::vector<std::pair<unsigned char, unsigned char>> rowOrder;
@@ -544,7 +536,7 @@ void CalculatedEdx::gatherRowClusterData(o2::tpc::TrackTPC& track, const std::ve
     }
 
     const bool mergeable = mergeableRows.count(rowKey) > 0;
-    gatherRowClusterDataForRow(track, fragmentClusters, fragmentIsShared, sectorIndex, rowIndex, mergeable, rowIndexOld, sectorIndexOld, occupancyROC, rowData, refitAbandoned);
+    gatherRowClusterDataForRow(track, fragmentClusters, fragmentIsShared, sectorIndex, rowIndex, mergeable, rowIndexOld, sectorIndexOld, occupancyROC, rowData, refitFailed);
     rowData.back().inputClusterIndices = clusterIndices; // positions in the externally supplied clusters vector grouped into this row
     rowIndexOld = rowIndex;
     sectorIndexOld = sectorIndex;

GPU/GPUTracking/Interface/GPUO2InterfaceRefit.cxx

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

GPU/GPUTracking/Interface/GPUO2InterfaceRefit.h

@@ -64,10 +64,10 @@ class GPUO2InterfaceRefit
   GPUO2InterfaceRefit(const o2::tpc::ClusterNativeAccess* cl, const o2::gpu::TPCFastTransformPOD* trans, float bzNominalGPU, const o2::tpc::TPCClRefElem* trackRef, uint32_t nHbfPerTf = 0, const uint8_t* sharedmap = nullptr, const uint32_t* occupancymap = nullptr, int32_t occupancyMapSize = -1, const std::vector<o2::tpc::TrackTPC>* trks = nullptr, o2::base::Propagator* p = nullptr);
   ~GPUO2InterfaceRefit();
 
-  int32_t RefitTrackAsGPU(o2::tpc::TrackTPC& trk, bool outward = false, bool resetCov = false, bool* reachedReferenceOut = nullptr);
-  int32_t RefitTrackAsTrackParCov(o2::tpc::TrackTPC& trk, bool outward = false, bool resetCov = false, bool* reachedReferenceOut = nullptr);
-  int32_t RefitTrackAsGPU(o2::track::TrackParCov& trk, const o2::tpc::TrackTPCClusRef& clusRef, float time0, float* chi2 = nullptr, bool outward = false, bool resetCov = false, bool* reachedReferenceOut = nullptr);
-  int32_t RefitTrackAsTrackParCov(o2::track::TrackParCov& trk, const o2::tpc::TrackTPCClusRef& clusRef, float time0, float* chi2 = nullptr, bool outward = false, bool resetCov = false, bool* reachedReferenceOut = nullptr);
+  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);
   void setTrackReferenceX(float v);
   void setIgnoreErrorsAtTrackEnds(bool v);
   void updateCalib(const o2::gpu::TPCFastTransformPOD* trans, float bzNominalGPU);

GPU/GPUTracking/Merger/GPUTPCGMTrackParam.cxx

@@ -304,11 +304,11 @@ GPUd() bool GPUTPCGMTrackParam::Fit(GPUTPCGMMerger* GPUrestrict() merger, int32_
               relTime /= clusterCount;
               relTime = relTime - CAMath::Round(relTime);
               if (acc) {
-                dEdx.fillCluster(qtot, qmax, cluster.row, cluster.sector, mP[2], mP[3], merger->GetConstantMem()->calibObjects, zz, pad, relTime, iTrk, clusterState);
+                dEdx.fillCluster(qtot, qmax, cluster.row, cluster.sector, mP[2], mP[3], merger->GetConstantMem()->calibObjects, zz, pad, relTime GPUCA_DEBUG_STREAMER_CHECK(, iTrk, clusterState));
               }
               if GPUCA_RTC_CONSTEXPR (GPUCA_GET_CONSTEXPR(param.rec.tpc, dEdxClusterRejectionFlagMask) != GPUCA_GET_CONSTEXPR(param.rec.tpc, dEdxClusterRejectionFlagMaskAlt)) {
                 if (accAlt) {
-                  dEdxAlt.fillCluster(qtot, qmax, cluster.row, cluster.sector, mP[2], mP[3], merger->GetConstantMem()->calibObjects, zz, pad, relTime, iTrk, clusterState);
+                  dEdxAlt.fillCluster(qtot, qmax, cluster.row, cluster.sector, mP[2], mP[3], merger->GetConstantMem()->calibObjects, zz, pad, relTime GPUCA_DEBUG_STREAMER_CHECK(, iTrk, clusterState));
                 }
               }
             }
@@ -364,16 +364,14 @@ GPUd() bool GPUTPCGMTrackParam::Fit(GPUTPCGMMerger* GPUrestrict() merger, int32_
   return true;
 }
 
-GPUdni() bool GPUTPCGMTrackParam::MoveToReference(GPUTPCGMPropagator& prop, const GPUParam& param, float& Alpha)
+GPUdni() void GPUTPCGMTrackParam::MoveToReference(GPUTPCGMPropagator& prop, const GPUParam& param, float& Alpha)
 {
   static constexpr float kDeg2Rad = M_PI / 180.f;
   static constexpr float kSectAngle = 2 * M_PI / 18.f;
 
-  bool reachedReference = true;
   if (param.rec.tpc.trackReferenceX <= 500) {
     GPUTPCGMTrackParam save = *this;
     float saveAlpha = Alpha;
-    reachedReference = false;
     for (int32_t attempt = 0; attempt < 3; attempt++) {
       float dAngle = CAMath::Round(CAMath::ATan2(mP[0], mX) / kDeg2Rad / 20.f) * kSectAngle;
       Alpha += dAngle;
@@ -382,7 +380,7 @@ GPUdni() bool GPUTPCGMTrackParam::MoveToReference(GPUTPCGMPropagator& prop, cons
       }
       ConstrainSinPhi();
       if (CAMath::Abs(mP[0]) <= mX * CAMath::Tan(kSectAngle / 2.f)) {
-        return true;
+        return;
       }
     }
     *this = save;
@@ -394,7 +392,6 @@ GPUdni() bool GPUTPCGMTrackParam::MoveToReference(GPUTPCGMPropagator& prop, cons
     ConstrainSinPhi();
     Alpha += dAngle;
   }
-  return reachedReference;
 }
 
 GPUd() void GPUTPCGMTrackParam::MirrorTo(GPUTPCGMPropagator& GPUrestrict() prop, float toY, float toZ, bool inFlyDirection, const GPUParam& param, uint8_t row, uint8_t clusterState, bool mirrorParameters, int8_t sector)