Skip to content

Commit bbb9289

Browse files
f3schf3sch
committed
Common: DCAFitter add fit status code
1 parent 5a7fba3 commit bbb9289

File tree

1 file changed

+67
-13
lines changed

1 file changed

+67
-13
lines changed

‎Common/DCAFitter/include/DCAFitter/DCAFitterN.h‎

Lines changed: 67 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -118,10 +118,26 @@ class DCAFitterN
118118
using ArrTrPos = o2::gpu::gpustd::array<Vec3D, N>; // container of Track positions
119119

120120
public:
121-
enum BadCovPolicy { // if encountering non-positive defined cov. matrix, the choice is:
122-
Discard = 0, // stop evaluation
123-
Override = 1, // override correlation coef. to have cov.matrix pos.def and continue
124-
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)
121+
enum BadCovPolicy : uint8_t { // if encountering non-positive defined cov. matrix, the choice is:
122+
Discard = 0, // stop evaluation
123+
Override = 1, // override correlation coef. to have cov.matrix pos.def and continue
124+
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)
125+
};
126+
127+
enum FitStatus : uint8_t { // fit status of crossing hypothesis
128+
None = 0, // no status set
129+
Converged = 1, // fit converged
130+
NoCrossing = 2, // no reasaonable crossing was found
131+
RejRadius = 3, // radius of crossing was not acceptable
132+
RejTrackX = 4, // one candidate track x was below the mimimum required radius
133+
RejTrackRoughZ = 5, // rejected by rough cut on tracks Z difference
134+
RejChi2Max = 6, // rejected by maximum chi2 cut
135+
FailProp = 7, // propagation of at least prong to PCA failed
136+
FailInvCov = 8, // inversion of cov.-matrix failed
137+
FailInvWeight = 9, // inversion of Ti weight matrix failed
138+
FailInv2ndDeriv = 10, // inversion of 2nd derivatives failed
139+
FailCorrTracks = 11, // correction of tracks to updated x failed
140+
FailCloserAlt = 12, // alternative PCA is closer
125141
};
126142

127143
static constexpr int getNProngs() { return N; }
@@ -154,7 +170,7 @@ class DCAFitterN
154170
///< check if propagation of tracks to candidate vertex was done
155171
GPUd() bool isPropagateTracksToVertexDone(int cand = 0) const { return mTrPropDone[mOrder[cand]]; }
156172

157-
///< check if propagation of tracks to candidate vertex was done
173+
///< check if propagation of tracks to candidate vertex failed
158174
bool isPropagationFailure(int cand = 0) const { return mPropFailed[mOrder[cand]]; }
159175

160176
///< track param propagated to V0 candidate (no check for the candidate validity)
@@ -201,6 +217,8 @@ class DCAFitterN
201217

202218
const Track* getOrigTrackPtr(int i) const { return mOrigTrPtr[i]; }
203219

220+
GPUdi() FitStatus getFitStatus(int cand = 0) const noexcept { return mFitStatus[mOrder[cand]]; }
221+
204222
///< return number of iterations during minimization (no check for its validity)
205223
GPUdi() int getNIterations(int cand = 0) const { return mNIters[mOrder[cand]]; }
206224
GPUdi() void setPropagateToPCA(bool v = true) { mPropagateToPCA = v; }
@@ -315,6 +333,8 @@ class DCAFitterN
315333
{
316334
mCurHyp = 0;
317335
mAllowAltPreference = true;
336+
mOrder.fill(0);
337+
mFitStatus.fill(FitStatus::None);
318338
}
319339

320340
GPUdi() static void setTrackPos(Vec3D& pnt, const Track& tr)
@@ -362,12 +382,13 @@ class DCAFitterN
362382
LogLogThrottler mLoggerBadCov{};
363383
LogLogThrottler mLoggerBadInv{};
364384
LogLogThrottler mLoggerBadProp{};
365-
MatSym3D mWeightInv; // inverse weight of single track, [sum{M^T E M}]^-1 in EQ.T
385+
MatSym3D mWeightInv; // inverse weight of single track, [sum{M^T E M}]^-1 in EQ.T
366386
o2::gpu::gpustd::array<int, MAXHYP> mOrder{0};
367387
int mCurHyp = 0;
368388
int mCrossIDCur = 0;
369389
int mCrossIDAlt = -1;
370390
BadCovPolicy mBadCovPolicy{BadCovPolicy::Discard}; // what to do in case of non-pos-def. cov. matrix, see BadCovPolicy enum
391+
o2::gpu::gpustd::array<FitStatus, MAXHYP> mFitStatus{}; // fit status of each hypothesis fit
371392
bool mAllowAltPreference = true; // if the fit converges to alternative PCA seed, abandon the current one
372393
bool mUseAbsDCA = false; // use abs. distance minimization rather than chi2
373394
bool mWeightedFinalPCA = false; // recalculate PCA as a cov-matrix weighted mean, even if absDCA method was used
@@ -390,7 +411,7 @@ class DCAFitterN
390411
float mMaxStep = 2.0; // Max step for propagation with Propagator
391412
int mFitterID = 0; // locat fitter ID (mostly for debugging)
392413
size_t mCallID = 0;
393-
ClassDefNV(DCAFitterN, 2);
414+
ClassDefNV(DCAFitterN, 3);
394415
};
395416

396417
///_________________________________________________________________________
@@ -407,7 +428,8 @@ GPUd() int DCAFitterN<N, Args...>::process(const Tr&... args)
407428
mTrAux[i].set(*mOrigTrPtr[i], mBz);
408429
}
409430
if (!mCrossings.set(mTrAux[0], *mOrigTrPtr[0], mTrAux[1], *mOrigTrPtr[1], mMaxDXYIni, mIsCollinear)) { // even for N>2 it should be enough to test just 1 loop
410-
return 0; // no crossing
431+
mFitStatus[mCurHyp] = FitStatus::NoCrossing;
432+
return 0;
411433
}
412434
for (int ih = 0; ih < MAXHYP; ih++) {
413435
mPropFailed[ih] = false;
@@ -428,6 +450,7 @@ GPUd() int DCAFitterN<N, Args...>::process(const Tr&... args)
428450
for (int ic = 0; ic < mCrossings.nDCA; ic++) {
429451
// check if radius is acceptable
430452
if (mCrossings.xDCA[ic] * mCrossings.xDCA[ic] + mCrossings.yDCA[ic] * mCrossings.yDCA[ic] > mMaxR2) {
453+
mFitStatus[mCurHyp] = FitStatus::RejRadius;
431454
continue;
432455
}
433456
mCrossIDCur = ic;
@@ -468,6 +491,7 @@ GPUd() bool DCAFitterN<N, Args...>::calcPCACoefs()
468491
{
469492
//< calculate Ti matrices for global vertex decomposition to V = sum_{0<i<N} Ti pi, see EQ.T in the ref
470493
if (!calcInverseWeight()) {
494+
mFitStatus[mCurHyp] = FitStatus::FailInvWeight;
471495
return false;
472496
}
473497
for (int i = N; i--;) { // build Mi*Ei matrix
@@ -727,6 +751,7 @@ GPUd() bool DCAFitterN<N, Args...>::recalculatePCAWithErrors(int cand)
727751
#endif
728752
mLoggerBadCov.evCountPrev = mLoggerBadCov.evCount;
729753
}
754+
mFitStatus[mCurHyp] = FitStatus::FailInvCov;
730755
if (mBadCovPolicy == Discard) {
731756
return false;
732757
} else if (mBadCovPolicy == OverrideAndFlag) {
@@ -935,10 +960,14 @@ GPUd() bool DCAFitterN<N, Args...>::minimizeChi2()
935960
for (int i = N; i--;) {
936961
mCandTr[mCurHyp][i] = *mOrigTrPtr[i];
937962
auto x = mTrAux[i].c * mPCA[mCurHyp][0] + mTrAux[i].s * mPCA[mCurHyp][1]; // X of PCA in the track frame
938-
if (x < mMinXSeed || !propagateToX(mCandTr[mCurHyp][i], x)) {
963+
if (x < mMinXSeed) {
964+
mFitStatus[mCurHyp] = FitStatus::RejTrackX;
939965
return false;
940966
}
941-
setTrackPos(mTrPos[mCurHyp][i], mCandTr[mCurHyp][i]); // prepare positions
967+
if (!propagateToX(mCandTr[mCurHyp][i], x)) {
968+
return false;
969+
}
970+
setTrackPos(mTrPos[mCurHyp][i], mCandTr[mCurHyp][i]); // prepare positions
942971
if (!mTrcEInv[mCurHyp][i].set(mCandTr[mCurHyp][i], XerrFactor)) { // prepare inverse cov.matrices at starting point
943972
if (mLoggerBadCov.needToLog()) {
944973
#ifndef GPUCA_GPUCODE
@@ -950,6 +979,7 @@ GPUd() bool DCAFitterN<N, Args...>::minimizeChi2()
950979
#endif
951980
mLoggerBadCov.evCountPrev = mLoggerBadCov.evCount;
952981
}
982+
mFitStatus[mCurHyp] = FitStatus::FailInvCov;
953983
if (mBadCovPolicy == Discard) {
954984
return false;
955985
} else if (mBadCovPolicy == OverrideAndFlag) {
@@ -959,6 +989,7 @@ GPUd() bool DCAFitterN<N, Args...>::minimizeChi2()
959989
}
960990

961991
if (mMaxDZIni > 0 && !roughDZCut()) { // apply rough cut on tracks Z difference
992+
mFitStatus[mCurHyp] = FitStatus::RejTrackX;
962993
return false;
963994
}
964995

@@ -979,28 +1010,36 @@ GPUd() bool DCAFitterN<N, Args...>::minimizeChi2()
9791010
printf("fitter %d: error (%ld muted): Inversion failed\n", mFitterID, mLoggerBadCov.getNMuted());
9801011
mLoggerBadInv.evCountPrev = mLoggerBadInv.evCount;
9811012
}
1013+
mFitStatus[mCurHyp] = FitStatus::FailInv2ndDeriv;
9821014
return false;
9831015
}
9841016
VecND dx = mD2Chi2Dx2 * mDChi2Dx;
9851017
if (!correctTracks(dx)) {
1018+
mFitStatus[mCurHyp] = FitStatus::FailCorrTracks;
9861019
return false;
9871020
}
9881021
calcPCA(); // updated PCA
9891022
if (mCrossIDAlt >= 0 && closerToAlternative()) {
1023+
mFitStatus[mCurHyp] = FitStatus::FailCloserAlt;
9901024
mAllowAltPreference = false;
9911025
return false;
9921026
}
9931027
calcTrackResiduals(); // updated residuals
9941028
chi2Upd = calcChi2(); // updated chi2
9951029
if (getAbsMax(dx) < mMinParamChange || chi2Upd > chi2 * mMinRelChi2Change) {
9961030
chi2 = chi2Upd;
1031+
mFitStatus[mCurHyp] = FitStatus::Converged;
9971032
break; // converged
9981033
}
9991034
chi2 = chi2Upd;
10001035
} while (++mNIters[mCurHyp] < mMaxIter);
10011036
//
10021037
mChi2[mCurHyp] = chi2 * NInv;
1003-
return mChi2[mCurHyp] < mMaxChi2;
1038+
if (mChi2[mCurHyp] >= mMaxChi2) {
1039+
mFitStatus[mCurHyp] = FitStatus::RejChi2Max;
1040+
return false;
1041+
}
1042+
return true;
10041043
}
10051044

10061045
//___________________________________________________________________
@@ -1012,12 +1051,17 @@ GPUd() bool DCAFitterN<N, Args...>::minimizeChi2NoErr()
10121051
for (int i = N; i--;) {
10131052
mCandTr[mCurHyp][i] = *mOrigTrPtr[i];
10141053
auto x = mTrAux[i].c * mPCA[mCurHyp][0] + mTrAux[i].s * mPCA[mCurHyp][1]; // X of PCA in the track frame
1015-
if (x < mMinXSeed || !propagateParamToX(mCandTr[mCurHyp][i], x)) {
1054+
if (x < mMinXSeed) {
1055+
mFitStatus[mCurHyp] = FitStatus::RejTrackX;
1056+
return false;
1057+
}
1058+
if (!propagateParamToX(mCandTr[mCurHyp][i], x)) {
10161059
return false;
10171060
}
10181061
setTrackPos(mTrPos[mCurHyp][i], mCandTr[mCurHyp][i]); // prepare positions
10191062
}
10201063
if (mMaxDZIni > 0 && !roughDZCut()) { // apply rough cut on tracks Z difference
1064+
mFitStatus[mCurHyp] = FitStatus::RejTrackX;
10211065
return false;
10221066
}
10231067

@@ -1035,28 +1079,36 @@ GPUd() bool DCAFitterN<N, Args...>::minimizeChi2NoErr()
10351079
printf("itter %d: error (%ld muted): Inversion failed\n", mFitterID, mLoggerBadCov.getNMuted());
10361080
mLoggerBadInv.evCountPrev = mLoggerBadInv.evCount;
10371081
}
1082+
mFitStatus[mCurHyp] = FitStatus::FailInv2ndDeriv;
10381083
return false;
10391084
}
10401085
VecND dx = mD2Chi2Dx2 * mDChi2Dx;
10411086
if (!correctTracks(dx)) {
1087+
mFitStatus[mCurHyp] = FitStatus::FailCorrTracks;
10421088
return false;
10431089
}
10441090
calcPCANoErr(); // updated PCA
10451091
if (mCrossIDAlt >= 0 && closerToAlternative()) {
1092+
mFitStatus[mCurHyp] = FitStatus::FailCloserAlt;
10461093
mAllowAltPreference = false;
10471094
return false;
10481095
}
10491096
calcTrackResiduals(); // updated residuals
10501097
chi2Upd = calcChi2NoErr(); // updated chi2
10511098
if (getAbsMax(dx) < mMinParamChange || chi2Upd > chi2 * mMinRelChi2Change) {
10521099
chi2 = chi2Upd;
1100+
mFitStatus[mCurHyp] = FitStatus::Converged;
10531101
break; // converged
10541102
}
10551103
chi2 = chi2Upd;
10561104
} while (++mNIters[mCurHyp] < mMaxIter);
10571105
//
10581106
mChi2[mCurHyp] = chi2 * NInv;
1059-
return mChi2[mCurHyp] < mMaxChi2;
1107+
if (mChi2[mCurHyp] >= mMaxChi2) {
1108+
mFitStatus[mCurHyp] = FitStatus::RejChi2Max;
1109+
return false;
1110+
}
1111+
return true;
10601112
}
10611113

10621114
//___________________________________________________________________
@@ -1182,6 +1234,7 @@ GPUdi() bool DCAFitterN<N, Args...>::propagateParamToX(o2::track::TrackPar& t, f
11821234
res = t.propagateParamTo(x, mBz);
11831235
}
11841236
if (!res) {
1237+
mFitStatus[mCurHyp] = FitStatus::FailProp;
11851238
mPropFailed[mCurHyp] = true;
11861239
if (mLoggerBadProp.needToLog()) {
11871240
#ifndef GPUCA_GPUCODE
@@ -1208,6 +1261,7 @@ GPUdi() bool DCAFitterN<N, Args...>::propagateToX(o2::track::TrackParCov& t, flo
12081261
res = t.propagateTo(x, mBz);
12091262
}
12101263
if (!res) {
1264+
mFitStatus[mCurHyp] = FitStatus::FailProp;
12111265
mPropFailed[mCurHyp] = true;
12121266
if (mLoggerBadProp.needToLog()) {
12131267
#ifndef GPUCA_GPUCODE

0 commit comments

Comments
 (0)