DCAFitterN: log-throttling for err.messages + user policy for bad CovMat

Due to the linearization errors the covariance matrix of the track propagated to some point may become non-positive defined.
In this case an error will be logged (logarithmically throttled), the relevant correlation coefficient of the cov.matrix is
redefined to cure the position part of the cov.matrix and further program flow depends on the user settings for
DCAFitterN::setBadCovPolicy(v):

DCAFitterN::setBadCovPolicy(DCAFitterN::Discard) : abandon fit (default)

DCAFitterN::setBadCovPolicy(DCAFitterN::Override) : continue fit with overridden cov.matrix

DCAFitterN::setBadCovPolicy(DCAFitterN::OverrideAnFlag): continue fit with overridden cov.matrix but set the propagation failure flag (can be checked using the
isPropagationFailure(int cand = 0) method).

Author: shahor02

Common/DCAFitter/README.md

@@ -2,7 +2,7 @@
 \page refDetectorsVertexing DCAFitter
 /doxy -->
 
-## DCAFitterN
+# DCAFitterN
 
 Templated class to fit the Point of Closest Approach (PCA) of secondary vertex with N prongs. Allows minimization of either absolute or weighted Distances of Closest Approach (DCA) of N tracks to their common PCA.
 
@@ -74,7 +74,22 @@ Extra method `setWeightedFinalPCA(bool)` is provided for the "mixed" mode: if `s
 but the final V0 position will be calculated using weighted average. One can also recalculate the V0 position by the weighted average method by calling explicitly
 `ft.recalculatePCAWithErrors(int icand=0)`, w/o prior call of `setWeightedFinalPCA(true)`: this will update the position returned by the `getPCACandidate(int cand = 0)`.
 
-The covariance matrix of the V0 position is calculated as an inversed sum of tracks inversed covariances at respective `X_dca` points.
+The covariance matrix of the V0 position is calculated as an inverted sum of tracks inversed covariances at respective `X_dca` points.
 
 See ``O2/Common/DCAFitter/test/testDCAFitterN.cxx`` for more extended example.
 Currently only 2 and 3 prongs permitted, thought this can be changed by modifying ``DCAFitterN::NMax`` constant.
+
+## Error handling
+
+It may happen that the track propagation to the the proximity of the PCA fails at the various stage of the fit. In this case the fit is abandoned and the failure flag is set, it can be checked using
+isPropagationFailure(int cand = 0)` method.
+
+Also, due to the linearization errors the covariance matrix of the track propagated to some point may become non-positive defined.
+In this case the relevant correlation coefficient of the cov.matrix is redefined to cure the position part of the cov.matrix and further program flow depends on the user settings for `DCAFitterN::setBadCovPolicy(v)`:
+
+`DCAFitterN::setBadCovPolicy(DCAFitterN::Discard);` : abandon fit (default)
+
+`DCAFitterN::setBadCovPolicy(DCAFitterN::Override);` : continue fit with overridden cov.matrix
+
+`DCAFitterN::setBadCovPolicy(DCAFitterN::OverrideAnFlag);` continue fit with overridden cov.matrix but set the propagation failure flag (can be checked using the same `isPropagationFailure(int cand = 0)` method).
+

Common/DCAFitter/include/DCAFitter/DCAFitterN.h

@@ -26,35 +26,32 @@ namespace o2
 {
 namespace vertexing
 {
+
 ///__________________________________________________________________________________
 ///< Inverse cov matrix (augmented by a dummy X error) of the point defined by the track
 struct TrackCovI {
   float sxx, syy, syz, szz;
 
-  GPUd() TrackCovI(const o2::track::TrackParCov& trc, float xerrFactor = 1.) { set(trc, xerrFactor); }
-
   GPUdDefault() TrackCovI() = default;
 
-  GPUd() void set(const o2::track::TrackParCov& trc, float xerrFactor = 1)
+  GPUd() bool set(const o2::track::TrackParCov& trc, float xerrFactor = 1.f)
   {
     // we assign Y error to X for DCA calculation
     // (otherwise for quazi-collinear tracks the X will not be constrained)
     float cyy = trc.getSigmaY2(), czz = trc.getSigmaZ2(), cyz = trc.getSigmaZY(), cxx = cyy * xerrFactor;
     float detYZ = cyy * czz - cyz * cyz;
+    bool res = true;
     if (detYZ <= 0.) {
-#ifndef GPUCA_GPUCODE
-      printf("overriding invalid track covariance from %s\n", trc.asString().c_str());
-#else
-      printf("overriding invalid track covariance cyy:%e czz:%e cyz:%e\n", cyy, czz, cyz);
-#endif
       cyz = o2::gpu::GPUCommonMath::Sqrt(cyy * czz) * (cyz > 0 ? 0.98f : -0.98f);
       detYZ = cyy * czz - cyz * cyz;
+      res = false;
     }
     auto detYZI = 1. / detYZ;
     sxx = 1. / cxx;
     syy = czz * detYZI;
     syz = -cyz * detYZI;
     szz = cyy * detYZI;
+    return res;
   }
 };
 
@@ -74,6 +71,27 @@ struct TrackDeriv {
   }
 };
 
+///__________________________________________________________________________
+///< Log log-throttling helper
+struct LogLogThrottler {
+  size_t evCount{0};
+  size_t evCountPrev{0};
+  size_t logCount{0};
+
+  GPUdi() bool needToLog()
+  {
+    if (size_t(o2::gpu::GPUCommonMath::Log(++evCount)) + 1 > logCount) {
+      logCount++;
+      return true;
+    }
+    return false;
+  }
+
+  GPUdi() size_t getNMuted() const { return evCount - evCountPrev - 1; }
+
+  GPUdi() void clear() { evCount = evCountPrev = logCount = 0; }
+};
+
 template <int N, typename... Args>
 class DCAFitterN
 {
@@ -100,6 +118,12 @@ class DCAFitterN
   using ArrTrPos = o2::gpu::gpustd::array<Vec3D, N>;         // container of Track positions
 
  public:
+  enum BadCovPolicy {   // if encountering non-positive defined cov. matrix, the choice is:
+    Discard = 0,        // stop evaluation
+    Override = 1,       // override correlation coef. to have cov.matrix pos.def and continue
+    OverrideAndFlag = 2 // override correlation coef. to have cov.matrix pos.def, set mPropFailed flag of corresponding candidate to true and continue (up to the user to check the flag)
+  };
+
   static constexpr int getNProngs() { return N; }
 
   DCAFitterN() = default;
@@ -300,6 +324,9 @@ class DCAFitterN
     pnt[2] = tr.getZ();
   }
 
+  void setBadCovPolicy(BadCovPolicy v) { mBadCovPolicy = v; }
+  BadCovPolicy getBadCovPolicy() const { return mBadCovPolicy; }
+
  private:
   // vectors of 1st derivatives of track local residuals over X parameters
   o2::gpu::gpustd::array<o2::gpu::gpustd::array<Vec3D, N>, N> mDResidDx;
@@ -325,11 +352,15 @@ class DCAFitterN
   o2::gpu::gpustd::array<int, MAXHYP> mNIters;           // number of iterations for each seed
   o2::gpu::gpustd::array<bool, MAXHYP> mTrPropDone{};    // Flag that the tracks are fully propagated to PCA
   o2::gpu::gpustd::array<bool, MAXHYP> mPropFailed{};    // Flag that some propagation failed for this PCA candidate
+  LogLogThrottler mLoggerBadCov{};
+  LogLogThrottler mLoggerBadInv{};
+  LogLogThrottler mLoggerBadProp{};
   MatSym3D mWeightInv;                                   // inverse weight of single track, [sum{M^T E M}]^-1 in EQ.T
   o2::gpu::gpustd::array<int, MAXHYP> mOrder{0};
   int mCurHyp = 0;
   int mCrossIDCur = 0;
   int mCrossIDAlt = -1;
+  BadCovPolicy mBadCovPolicy{BadCovPolicy::Discard};                                              // what to do in case of non-pos-def. cov. matrix, see BadCovPolicy enum
   bool mAllowAltPreference = true;                                                                // if the fit converges to alternative PCA seed, abandon the current one
   bool mUseAbsDCA = false;                                                                        // use abs. distance minimization rather than chi2
   bool mWeightedFinalPCA = false;                                                                 // recalculate PCA as a cov-matrix weighted mean, even if absDCA method was used
@@ -678,7 +709,23 @@ GPUd() bool DCAFitterN<N, Args...>::recalculatePCAWithErrors(int cand)
   mCurHyp = mOrder[cand];
   if (mUseAbsDCA) {
     for (int i = N; i--;) {
-      mTrcEInv[mCurHyp][i].set(mCandTr[mCurHyp][i], XerrFactor); // prepare inverse cov.matrices at starting point
+      if (!mTrcEInv[mCurHyp][i].set(mCandTr[mCurHyp][i], XerrFactor)) { // prepare inverse cov.matrices at starting point
+        if (mLoggerBadCov.needToLog()) {
+#ifndef GPUCA_GPUCODE
+          printf("fitter %d: error (%ld muted): overrode invalid track covariance from %s\n",
+                 mFitterID, mLoggerBadCov.getNMuted(), mCandTr[mCurHyp][i].asString().c_str());
+#else
+          printf("fitter %d: error (%ld muted): overrode invalid track covariance cyy:%e czz:%e cyz:%e\n",
+                 mFitterID, mLoggerBadCov.getNMuted(), mCandTr[mCurHyp][i].getSigmaY2(), mCandTr[mCurHyp][i].getSigmaZ2(), mCandTr[mCurHyp][i].getSigmaZY());
+#endif
+          mLoggerBadCov.evCountPrev = mLoggerBadCov.evCount;
+        }
+        if (mBadCovPolicy == Discard) {
+          return false;
+        } else if (mBadCovPolicy == OverrideAndFlag) {
+          mPropFailed[mCurHyp] = true;
+        } // otherwise, just use overridden errors w/o flagging
+      }
     }
     if (!calcPCACoefs()) {
       mCurHyp = saveCurHyp;
@@ -885,7 +932,23 @@ GPUd() bool DCAFitterN<N, Args...>::minimizeChi2()
       return false;
     }
     setTrackPos(mTrPos[mCurHyp][i], mCandTr[mCurHyp][i]);      // prepare positions
-    mTrcEInv[mCurHyp][i].set(mCandTr[mCurHyp][i], XerrFactor); // prepare inverse cov.matrices at starting point
+    if (!mTrcEInv[mCurHyp][i].set(mCandTr[mCurHyp][i], XerrFactor)) { // prepare inverse cov.matrices at starting point
+      if (mLoggerBadCov.needToLog()) {
+#ifndef GPUCA_GPUCODE
+        printf("fitter %d: error (%ld muted): overrode invalid track covariance from %s\n",
+               mFitterID, mLoggerBadCov.getNMuted(), mCandTr[mCurHyp][i].asString().c_str());
+#else
+        printf("fitter %d: error (%ld muted): overrode invalid track covariance cyy:%e czz:%e cyz:%e\n",
+               mFitterID, mLoggerBadCov.getNMuted(), mCandTr[mCurHyp][i].getSigmaY2(), mCandTr[mCurHyp][i].getSigmaZ2(), mCandTr[mCurHyp][i].getSigmaZY());
+#endif
+        mLoggerBadCov.evCountPrev = mLoggerBadCov.evCount;
+      }
+      if (mBadCovPolicy == Discard) {
+        return false;
+      } else if (mBadCovPolicy == OverrideAndFlag) {
+        mPropFailed[mCurHyp] = true;
+      } // otherwise, just use overridden errors w/o flagging
+    }
   }
 
   if (mMaxDZIni > 0 && !roughDZCut()) { // apply rough cut on tracks Z difference
@@ -905,11 +968,10 @@ GPUd() bool DCAFitterN<N, Args...>::minimizeChi2()
 
     // do Newton-Rapson iteration with corrections = - dchi2/d{x0..xN} * [ d^2chi2/d{x0..xN}^2 ]^-1
     if (!mD2Chi2Dx2.Invert()) {
-#ifndef GPUCA_GPUCODE_DEVICE
-      LOG(error) << "InversionFailed";
-#else
-      printf("InversionFailed\n");
-#endif
+      if (mLoggerBadInv.needToLog()) {
+        printf("fitter %d: error (%ld muted): Inversion failed\n", mFitterID, mLoggerBadCov.getNMuted());
+        mLoggerBadInv.evCountPrev = mLoggerBadInv.evCount;
+      }
       return false;
     }
     VecND dx = mD2Chi2Dx2 * mDChi2Dx;
@@ -962,11 +1024,10 @@ GPUd() bool DCAFitterN<N, Args...>::minimizeChi2NoErr()
 
     // do Newton-Rapson iteration with corrections = - dchi2/d{x0..xN} * [ d^2chi2/d{x0..xN}^2 ]^-1
     if (!mD2Chi2Dx2.Invert()) {
-#ifndef GPUCA_GPUCODE_DEVICE
-      LOG(error) << "InversionFailed";
-#else
-      printf("InversionFailed\n");
-#endif
+      if (mLoggerBadInv.needToLog()) {
+        printf("itter %d: error (%ld muted): Inversion failed\n", mFitterID, mLoggerBadCov.getNMuted());
+        mLoggerBadInv.evCountPrev = mLoggerBadInv.evCount;
+      }
       return false;
     }
     VecND dx = mD2Chi2Dx2 * mDChi2Dx;
@@ -1109,6 +1170,14 @@ GPUdi() bool DCAFitterN<N, Args...>::propagateParamToX(o2::track::TrackPar& t, f
   }
   if (!res) {
     mPropFailed[mCurHyp] = true;
+    if (mLoggerBadProp.needToLog()) {
+#ifndef GPUCA_GPUCODE
+      printf("fitter %d: error (%ld muted): propagation failed for %s\n", mFitterID, mLoggerBadProp.getNMuted(), t.asString().c_str());
+#else
+      printf("fitter %d: error (%ld muted): propagation failed\n", mFitterID, mLoggerBadProp.getNMuted());
+#endif
+      mLoggerBadProp.evCountPrev = mLoggerBadProp.evCount;
+    }
   }
   return res;
 }
@@ -1127,6 +1196,14 @@ GPUdi() bool DCAFitterN<N, Args...>::propagateToX(o2::track::TrackParCov& t, flo
   }
   if (!res) {
     mPropFailed[mCurHyp] = true;
+    if (mLoggerBadProp.needToLog()) {
+#ifndef GPUCA_GPUCODE
+      printf("fitter %d: error (%ld muted): propagation failed for %s\n", mFitterID, mLoggerBadProp.getNMuted(), t.asString().c_str());
+#else
+      printf("fitter %d: error (%ld muted): propagation failed\n", mFitterID, mLoggerBadProp.getNMuted());
+#endif
+      mLoggerBadProp.evCountPrev = mLoggerBadProp.evCount;
+    }
   }
   return res;
 }