Skip to content

Please consider the following formatting changes to #14132 #38

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
alibuild wants to merge 2 commits into from
Closed

Please consider the following formatting changes to #14132 #38

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
ba169e6
Common: DCAFitter add fit status code
f3sch Mar 31, 2025
14db0d4
Please consider the following formatting changes
alibuild Apr 1, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
80 changes: 67 additions & 13 deletions Common/DCAFitter/include/DCAFitter/DCAFitterN.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -118,10 +118,26 @@ 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)
enum BadCovPolicy : uint8_t { // 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)
};

enum FitStatus : uint8_t { // fit status of crossing hypothesis
None = 0, // no status set
Converged = 1, // fit converged
NoCrossing = 2, // no reasaonable crossing was found
RejRadius = 3, // radius of crossing was not acceptable
RejTrackX = 4, // one candidate track x was below the mimimum required radius
RejTrackRoughZ = 5, // rejected by rough cut on tracks Z difference
RejChi2Max = 6, // rejected by maximum chi2 cut
FailProp = 7, // propagation of at least prong to PCA failed
FailInvCov = 8, // inversion of cov.-matrix failed
FailInvWeight = 9, // inversion of Ti weight matrix failed
FailInv2ndDeriv = 10, // inversion of 2nd derivatives failed
FailCorrTracks = 11, // correction of tracks to updated x failed
FailCloserAlt = 12, // alternative PCA is closer
};

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

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

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

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

GPUdi() FitStatus getFitStatus(int cand = 0) const noexcept { return mFitStatus[mOrder[cand]]; }

///< return number of iterations during minimization (no check for its validity)
GPUdi() int getNIterations(int cand = 0) const { return mNIters[mOrder[cand]]; }
GPUdi() void setPropagateToPCA(bool v = true) { mPropagateToPCA = v; }
Expand Down Expand Up @@ -315,6 +333,8 @@ class DCAFitterN
{
mCurHyp = 0;
mAllowAltPreference = true;
mOrder.fill(0);
mFitStatus.fill(FitStatus::None);
}

GPUdi() static void setTrackPos(Vec3D& pnt, const Track& tr)
Expand Down Expand Up @@ -362,12 +382,13 @@ class DCAFitterN
LogLogThrottler mLoggerBadCov{};
LogLogThrottler mLoggerBadInv{};
LogLogThrottler mLoggerBadProp{};
MatSym3D mWeightInv; // inverse weight of single track, [sum{M^T E M}]^-1 in EQ.T
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
o2::gpu::gpustd::array<FitStatus, MAXHYP> mFitStatus{}; // fit status of each hypothesis fit
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 All @@ -390,7 +411,7 @@ class DCAFitterN
float mMaxStep = 2.0; // Max step for propagation with Propagator
int mFitterID = 0; // locat fitter ID (mostly for debugging)
size_t mCallID = 0;
ClassDefNV(DCAFitterN, 2);
ClassDefNV(DCAFitterN, 3);
};

///_________________________________________________________________________
Expand All @@ -407,7 +428,8 @@ GPUd() int DCAFitterN<N, Args...>::process(const Tr&... args)
mTrAux[i].set(*mOrigTrPtr[i], mBz);
}
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
return 0; // no crossing
mFitStatus[mCurHyp] = FitStatus::NoCrossing;
return 0;
}
for (int ih = 0; ih < MAXHYP; ih++) {
mPropFailed[ih] = false;
Expand All @@ -428,6 +450,7 @@ GPUd() int DCAFitterN<N, Args...>::process(const Tr&... args)
for (int ic = 0; ic < mCrossings.nDCA; ic++) {
// check if radius is acceptable
if (mCrossings.xDCA[ic] * mCrossings.xDCA[ic] + mCrossings.yDCA[ic] * mCrossings.yDCA[ic] > mMaxR2) {
mFitStatus[mCurHyp] = FitStatus::RejRadius;
continue;
}
mCrossIDCur = ic;
Expand Down Expand Up @@ -468,6 +491,7 @@ GPUd() bool DCAFitterN<N, Args...>::calcPCACoefs()
{
//< calculate Ti matrices for global vertex decomposition to V = sum_{0<i<N} Ti pi, see EQ.T in the ref
if (!calcInverseWeight()) {
mFitStatus[mCurHyp] = FitStatus::FailInvWeight;
return false;
}
for (int i = N; i--;) { // build Mi*Ei matrix
Expand Down Expand Up @@ -727,6 +751,7 @@ GPUd() bool DCAFitterN<N, Args...>::recalculatePCAWithErrors(int cand)
#endif
mLoggerBadCov.evCountPrev = mLoggerBadCov.evCount;
}
mFitStatus[mCurHyp] = FitStatus::FailInvCov;
if (mBadCovPolicy == Discard) {
return false;
} else if (mBadCovPolicy == OverrideAndFlag) {
Expand Down Expand Up @@ -935,10 +960,14 @@ GPUd() bool DCAFitterN<N, Args...>::minimizeChi2()
for (int i = N; i--;) {
mCandTr[mCurHyp][i] = *mOrigTrPtr[i];
auto x = mTrAux[i].c * mPCA[mCurHyp][0] + mTrAux[i].s * mPCA[mCurHyp][1]; // X of PCA in the track frame
if (x < mMinXSeed || !propagateToX(mCandTr[mCurHyp][i], x)) {
if (x < mMinXSeed) {
mFitStatus[mCurHyp] = FitStatus::RejTrackX;
return false;
}
setTrackPos(mTrPos[mCurHyp][i], mCandTr[mCurHyp][i]); // prepare positions
if (!propagateToX(mCandTr[mCurHyp][i], x)) {
return false;
}
setTrackPos(mTrPos[mCurHyp][i], mCandTr[mCurHyp][i]); // prepare positions
if (!mTrcEInv[mCurHyp][i].set(mCandTr[mCurHyp][i], XerrFactor)) { // prepare inverse cov.matrices at starting point
if (mLoggerBadCov.needToLog()) {
#ifndef GPUCA_GPUCODE
Expand All @@ -950,6 +979,7 @@ GPUd() bool DCAFitterN<N, Args...>::minimizeChi2()
#endif
mLoggerBadCov.evCountPrev = mLoggerBadCov.evCount;
}
mFitStatus[mCurHyp] = FitStatus::FailInvCov;
if (mBadCovPolicy == Discard) {
return false;
} else if (mBadCovPolicy == OverrideAndFlag) {
Expand All @@ -959,6 +989,7 @@ GPUd() bool DCAFitterN<N, Args...>::minimizeChi2()
}

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

Expand All @@ -979,28 +1010,36 @@ GPUd() bool DCAFitterN<N, Args...>::minimizeChi2()
printf("fitter %d: error (%ld muted): Inversion failed\n", mFitterID, mLoggerBadCov.getNMuted());
mLoggerBadInv.evCountPrev = mLoggerBadInv.evCount;
}
mFitStatus[mCurHyp] = FitStatus::FailInv2ndDeriv;
return false;
}
VecND dx = mD2Chi2Dx2 * mDChi2Dx;
if (!correctTracks(dx)) {
mFitStatus[mCurHyp] = FitStatus::FailCorrTracks;
return false;
}
calcPCA(); // updated PCA
if (mCrossIDAlt >= 0 && closerToAlternative()) {
mFitStatus[mCurHyp] = FitStatus::FailCloserAlt;
mAllowAltPreference = false;
return false;
}
calcTrackResiduals(); // updated residuals
chi2Upd = calcChi2(); // updated chi2
if (getAbsMax(dx) < mMinParamChange || chi2Upd > chi2 * mMinRelChi2Change) {
chi2 = chi2Upd;
mFitStatus[mCurHyp] = FitStatus::Converged;
break; // converged
}
chi2 = chi2Upd;
} while (++mNIters[mCurHyp] < mMaxIter);
//
mChi2[mCurHyp] = chi2 * NInv;
return mChi2[mCurHyp] < mMaxChi2;
if (mChi2[mCurHyp] >= mMaxChi2) {
mFitStatus[mCurHyp] = FitStatus::RejChi2Max;
return false;
}
return true;
}

//___________________________________________________________________
Expand All @@ -1012,12 +1051,17 @@ GPUd() bool DCAFitterN<N, Args...>::minimizeChi2NoErr()
for (int i = N; i--;) {
mCandTr[mCurHyp][i] = *mOrigTrPtr[i];
auto x = mTrAux[i].c * mPCA[mCurHyp][0] + mTrAux[i].s * mPCA[mCurHyp][1]; // X of PCA in the track frame
if (x < mMinXSeed || !propagateParamToX(mCandTr[mCurHyp][i], x)) {
if (x < mMinXSeed) {
mFitStatus[mCurHyp] = FitStatus::RejTrackX;
return false;
}
if (!propagateParamToX(mCandTr[mCurHyp][i], x)) {
return false;
}
setTrackPos(mTrPos[mCurHyp][i], mCandTr[mCurHyp][i]); // prepare positions
}
if (mMaxDZIni > 0 && !roughDZCut()) { // apply rough cut on tracks Z difference
mFitStatus[mCurHyp] = FitStatus::RejTrackX;
return false;
}

Expand All @@ -1035,28 +1079,36 @@ GPUd() bool DCAFitterN<N, Args...>::minimizeChi2NoErr()
printf("itter %d: error (%ld muted): Inversion failed\n", mFitterID, mLoggerBadCov.getNMuted());
mLoggerBadInv.evCountPrev = mLoggerBadInv.evCount;
}
mFitStatus[mCurHyp] = FitStatus::FailInv2ndDeriv;
return false;
}
VecND dx = mD2Chi2Dx2 * mDChi2Dx;
if (!correctTracks(dx)) {
mFitStatus[mCurHyp] = FitStatus::FailCorrTracks;
return false;
}
calcPCANoErr(); // updated PCA
if (mCrossIDAlt >= 0 && closerToAlternative()) {
mFitStatus[mCurHyp] = FitStatus::FailCloserAlt;
mAllowAltPreference = false;
return false;
}
calcTrackResiduals(); // updated residuals
chi2Upd = calcChi2NoErr(); // updated chi2
if (getAbsMax(dx) < mMinParamChange || chi2Upd > chi2 * mMinRelChi2Change) {
chi2 = chi2Upd;
mFitStatus[mCurHyp] = FitStatus::Converged;
break; // converged
}
chi2 = chi2Upd;
} while (++mNIters[mCurHyp] < mMaxIter);
//
mChi2[mCurHyp] = chi2 * NInv;
return mChi2[mCurHyp] < mMaxChi2;
if (mChi2[mCurHyp] >= mMaxChi2) {
mFitStatus[mCurHyp] = FitStatus::RejChi2Max;
return false;
}
return true;
}

//___________________________________________________________________
Expand Down Expand Up @@ -1182,6 +1234,7 @@ GPUdi() bool DCAFitterN<N, Args...>::propagateParamToX(o2::track::TrackPar& t, f
res = t.propagateParamTo(x, mBz);
}
if (!res) {
mFitStatus[mCurHyp] = FitStatus::FailProp;
mPropFailed[mCurHyp] = true;
if (mLoggerBadProp.needToLog()) {
#ifndef GPUCA_GPUCODE
Expand All @@ -1208,6 +1261,7 @@ GPUdi() bool DCAFitterN<N, Args...>::propagateToX(o2::track::TrackParCov& t, flo
res = t.propagateTo(x, mBz);
}
if (!res) {
mFitStatus[mCurHyp] = FitStatus::FailProp;
mPropFailed[mCurHyp] = true;
if (mLoggerBadProp.needToLog()) {
#ifndef GPUCA_GPUCODE
Expand Down