Skip to content

Commit

Permalink
DCAFitterN: log-throttling for err.messages + user policy for bad CovMat
Browse files Browse the repository at this point in the history 
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).
  • Loading branch information
@shahor02
shahor02 committed Dec 10, 2024
1 parent d06c2cf commit c5cbdc4
Show file tree
Hide file tree
Showing 2 changed files with 114 additions and 22 deletions.
19 changes: 17 additions & 2 deletions Common/DCAFitter/README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.

Expand Down Expand Up @@ -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).

117 changes: 97 additions & 20 deletions Common/DCAFitter/include/DCAFitter/DCAFitterN.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
}
};

Expand All @@ -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
{
Expand All @@ -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;
Expand Down Expand Up @@ -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;
Expand All @@ -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
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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
Expand All @@ -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;
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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;
}
Expand All @@ -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;
}
Expand Down

0 comments on commit c5cbdc4

Please sign in to comment.