From bbb92897fda664881db853dd547eb1870ad7c905 Mon Sep 17 00:00:00 2001 From: Felix Schlepper Date: Mon, 31 Mar 2025 13:20:01 +0200 Subject: [PATCH 1/5] Common: DCAFitter add fit status code --- .../DCAFitter/include/DCAFitter/DCAFitterN.h | 80 ++++++++++++++++--- 1 file changed, 67 insertions(+), 13 deletions(-) diff --git a/Common/DCAFitter/include/DCAFitter/DCAFitterN.h b/Common/DCAFitter/include/DCAFitter/DCAFitterN.h index 97ea6d206247b..a8bf1bb63bfb8 100644 --- a/Common/DCAFitter/include/DCAFitter/DCAFitterN.h +++ b/Common/DCAFitter/include/DCAFitter/DCAFitterN.h @@ -118,10 +118,26 @@ class DCAFitterN using ArrTrPos = o2::gpu::gpustd::array; // 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; } @@ -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) @@ -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; } @@ -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) @@ -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 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 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 @@ -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); }; ///_________________________________________________________________________ @@ -407,7 +428,8 @@ GPUd() int DCAFitterN::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; @@ -428,6 +450,7 @@ GPUd() int DCAFitterN::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; @@ -468,6 +491,7 @@ GPUd() bool DCAFitterN::calcPCACoefs() { //< calculate Ti matrices for global vertex decomposition to V = sum_{0::recalculatePCAWithErrors(int cand) #endif mLoggerBadCov.evCountPrev = mLoggerBadCov.evCount; } + mFitStatus[mCurHyp] = FitStatus::FailInvCov; if (mBadCovPolicy == Discard) { return false; } else if (mBadCovPolicy == OverrideAndFlag) { @@ -935,10 +960,14 @@ GPUd() bool DCAFitterN::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 @@ -950,6 +979,7 @@ GPUd() bool DCAFitterN::minimizeChi2() #endif mLoggerBadCov.evCountPrev = mLoggerBadCov.evCount; } + mFitStatus[mCurHyp] = FitStatus::FailInvCov; if (mBadCovPolicy == Discard) { return false; } else if (mBadCovPolicy == OverrideAndFlag) { @@ -959,6 +989,7 @@ GPUd() bool DCAFitterN::minimizeChi2() } if (mMaxDZIni > 0 && !roughDZCut()) { // apply rough cut on tracks Z difference + mFitStatus[mCurHyp] = FitStatus::RejTrackX; return false; } @@ -979,14 +1010,17 @@ GPUd() bool DCAFitterN::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; } @@ -994,13 +1028,18 @@ GPUd() bool DCAFitterN::minimizeChi2() 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; } //___________________________________________________________________ @@ -1012,12 +1051,17 @@ GPUd() bool DCAFitterN::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; } @@ -1035,14 +1079,17 @@ GPUd() bool DCAFitterN::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; } @@ -1050,13 +1097,18 @@ GPUd() bool DCAFitterN::minimizeChi2NoErr() 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; } //___________________________________________________________________ @@ -1182,6 +1234,7 @@ GPUdi() bool DCAFitterN::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 @@ -1208,6 +1261,7 @@ GPUdi() bool DCAFitterN::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 From 2a86f4ee9dfa64a7fc06bc4b0f28369351ada1c8 Mon Sep 17 00:00:00 2001 From: Felix Schlepper Date: Tue, 1 Apr 2025 21:08:42 +0200 Subject: [PATCH 2/5] Common: DCAFitter add test/summary of stats Signed-off-by: Felix Schlepper --- Common/DCAFitter/test/testDCAFitterN.cxx | 49 +++++++++++++++++++++++- 1 file changed, 48 insertions(+), 1 deletion(-) diff --git a/Common/DCAFitter/test/testDCAFitterN.cxx b/Common/DCAFitter/test/testDCAFitterN.cxx index a102a0a4253e3..bd00b5bed841e 100644 --- a/Common/DCAFitter/test/testDCAFitterN.cxx +++ b/Common/DCAFitter/test/testDCAFitterN.cxx @@ -134,7 +134,7 @@ TLorentzVector generate(Vec3D& vtx, std::vector& vctr, f float rad = forceQ[i] == 0 ? 600. : TMath::Abs(1. / trc.getCurvature(bz)); if (!trc.propagateTo(trc.getX() + (gRandom->Rndm() - 0.5) * rad * 0.05, bz) || !trc.rotate(trc.getAlpha() + (gRandom->Rndm() - 0.5) * 0.2)) { - printf("Failed to randomize "); + LOGP(error, "Failed to randomize "); trc.print(); } } @@ -143,6 +143,22 @@ TLorentzVector generate(Vec3D& vtx, std::vector& vctr, f return parent; } +static constexpr int NFitStatus{14}; +using FitStatusArray = std::array, NFitStatus>; +static constexpr const char* FitStatusNames[NFitStatus] = { + "None", "Converged", "MaxIter", "NoCrossing", "RejRadius", "RejTrackX", "RejTrackRoughZ", "RejChi2Max", + "FailProp", "FailInvConv", "FailInvWeight", "FailInv2ndDeriv", "FailCorrTracks", "FailCloserAlt"}; +inline void printStat(const FitStatusArray& a) +{ + LOGP(info, "FitStatus summary : ....A / ..AWD / ...WD (A=abs.dist;AWD=abs.wghPCA.dist;WD=wgh.dist)"); + for (int i{0}; i < NFitStatus; ++i) { + LOGP(info, "{:2d}={:20s}: {:5d} / {:5d} / {:5d}", i, FitStatusNames[i], a[i][0], a[i][1], a[i][2]); + } + BOOST_CHECK(a[0][0] == 0); // ensure coverage of all possible states + BOOST_CHECK(a[0][1] == 0); + BOOST_CHECK(a[0][2] == 0); +} + BOOST_AUTO_TEST_CASE(DCAFitterNProngs) { constexpr int NTest = 10000; @@ -159,6 +175,7 @@ BOOST_AUTO_TEST_CASE(DCAFitterNProngs) std::vector k0dec = {pion, pion}; std::vector dchdec = {pion, kch, pion}; std::vector vctracks; + FitStatusArray fitstat; Vec3D vtxGen; double bz = 5.0; @@ -166,6 +183,7 @@ BOOST_AUTO_TEST_CASE(DCAFitterNProngs) { LOG(info) << "\n\nProcessing 2-prong Helix - Helix case"; std::vector forceQ{1, 1}; + std::memset(fitstat.data(), 0, sizeof(fitstat)); o2::vertexing::DCAFitterN<2> ft; // 2 prong fitter ft.setBz(bz); @@ -196,6 +214,7 @@ BOOST_AUTO_TEST_CASE(DCAFitterNProngs) meanDA += minD; nfoundA++; } + ++fitstat[ft.getFitStatus()][0]; ft.setUseAbsDCA(true); ft.setWeightedFinalPCA(true); @@ -208,6 +227,7 @@ BOOST_AUTO_TEST_CASE(DCAFitterNProngs) meanDAW += minD; nfoundAW++; } + ++fitstat[ft.getFitStatus()][1]; ft.setUseAbsDCA(false); ft.setWeightedFinalPCA(false); @@ -220,6 +240,7 @@ BOOST_AUTO_TEST_CASE(DCAFitterNProngs) meanDW += minD; nfoundW++; } + ++fitstat[ft.getFitStatus()][2]; } // ft.print(); meanDA /= nfoundA ? nfoundA : 1; @@ -232,6 +253,7 @@ BOOST_AUTO_TEST_CASE(DCAFitterNProngs) << " mean.dist to truth: " << meanDAW << " CPU time: " << swAW.CpuTime() * 1000 << " ms"; LOG(info) << "2-prongs with wgh.dist minization: eff= " << float(nfoundW) / NTest << " mean.dist to truth: " << meanDW << " CPU time: " << swW.CpuTime() * 1000 << " ms"; + printStat(fitstat); BOOST_CHECK(nfoundA > 0.99 * NTest); BOOST_CHECK(nfoundAW > 0.99 * NTest); BOOST_CHECK(nfoundW > 0.99 * NTest); @@ -245,6 +267,7 @@ BOOST_AUTO_TEST_CASE(DCAFitterNProngs) { LOG(info) << "\n\nProcessing 2-prong Helix - Helix case gamma conversion"; std::vector forceQ{1, 1}; + std::memset(fitstat.data(), 0, sizeof(fitstat)); o2::vertexing::DCAFitterN<2> ft; // 2 prong fitter ft.setBz(bz); @@ -254,6 +277,8 @@ BOOST_AUTO_TEST_CASE(DCAFitterNProngs) ft.setMaxDXYIni(4); // do not consider V0 seeds with tracks XY-distance exceeding this. This is default anyway ft.setMinParamChange(1e-3); // stop iterations if max correction is below this value. This is default anyway ft.setMinRelChi2Change(0.9); // stop iterations if chi2 improves by less that this factor + ft.setMaxChi2(); + ft.setCollinear(true); std::string treeName2A = "gpr2a", treeName2AW = "gpr2aw", treeName2W = "gpr2w"; TStopwatch swA, swAW, swW; @@ -275,6 +300,7 @@ BOOST_AUTO_TEST_CASE(DCAFitterNProngs) meanDA += minD; nfoundA++; } + ++fitstat[ft.getFitStatus()][0]; ft.setUseAbsDCA(true); ft.setWeightedFinalPCA(true); @@ -287,6 +313,7 @@ BOOST_AUTO_TEST_CASE(DCAFitterNProngs) meanDAW += minD; nfoundAW++; } + ++fitstat[ft.getFitStatus()][1]; ft.setUseAbsDCA(false); ft.setWeightedFinalPCA(false); @@ -299,6 +326,7 @@ BOOST_AUTO_TEST_CASE(DCAFitterNProngs) meanDW += minD; nfoundW++; } + ++fitstat[ft.getFitStatus()][2]; } // ft.print(); meanDA /= nfoundA ? nfoundA : 1; @@ -311,6 +339,7 @@ BOOST_AUTO_TEST_CASE(DCAFitterNProngs) << " mean.dist to truth: " << meanDAW << " CPU time: " << swAW.CpuTime() * 1000 << " ms"; LOG(info) << "2-prongs with wgh.dist minization: eff= " << float(nfoundW) / NTest << " mean.dist to truth: " << meanDW << " CPU time: " << swW.CpuTime() * 1000 << " ms"; + printStat(fitstat); BOOST_CHECK(nfoundA > 0.99 * NTest); BOOST_CHECK(nfoundAW > 0.99 * NTest); BOOST_CHECK(nfoundW > 0.99 * NTest); @@ -324,6 +353,8 @@ BOOST_AUTO_TEST_CASE(DCAFitterNProngs) { std::vector forceQ{1, 1}; LOG(info) << "\n\nProcessing 2-prong Helix - Line case"; + std::memset(fitstat.data(), 0, sizeof(fitstat)); + o2::vertexing::DCAFitterN<2> ft; // 2 prong fitter ft.setBz(bz); ft.setPropagateToPCA(true); // After finding the vertex, propagate tracks to the DCA. This is default anyway @@ -354,6 +385,7 @@ BOOST_AUTO_TEST_CASE(DCAFitterNProngs) meanDA += minD; nfoundA++; } + ++fitstat[ft.getFitStatus()][0]; ft.setUseAbsDCA(true); ft.setWeightedFinalPCA(true); @@ -366,6 +398,7 @@ BOOST_AUTO_TEST_CASE(DCAFitterNProngs) meanDAW += minD; nfoundAW++; } + ++fitstat[ft.getFitStatus()][1]; ft.setUseAbsDCA(false); ft.setWeightedFinalPCA(false); @@ -378,6 +411,7 @@ BOOST_AUTO_TEST_CASE(DCAFitterNProngs) meanDW += minD; nfoundW++; } + ++fitstat[ft.getFitStatus()][2]; } // ft.print(); meanDA /= nfoundA ? nfoundA : 1; @@ -390,6 +424,7 @@ BOOST_AUTO_TEST_CASE(DCAFitterNProngs) << " mean.dist to truth: " << meanDAW << " CPU time: " << swAW.CpuTime() * 1000 << " ms"; LOG(info) << "2-prongs with wgh.dist minization: eff= " << float(nfoundW) / NTest << " mean.dist to truth: " << meanDW << " CPU time: " << swW.CpuTime() * 1000 << " ms"; + printStat(fitstat); BOOST_CHECK(nfoundA > 0.99 * NTest); BOOST_CHECK(nfoundAW > 0.99 * NTest); BOOST_CHECK(nfoundW > 0.99 * NTest); @@ -403,6 +438,8 @@ BOOST_AUTO_TEST_CASE(DCAFitterNProngs) { std::vector forceQ{0, 0}; LOG(info) << "\n\nProcessing 2-prong Line - Line case"; + std::memset(fitstat.data(), 0, sizeof(fitstat)); + o2::vertexing::DCAFitterN<2> ft; // 2 prong fitter ft.setBz(bz); ft.setPropagateToPCA(true); // After finding the vertex, propagate tracks to the DCA. This is default anyway @@ -432,6 +469,7 @@ BOOST_AUTO_TEST_CASE(DCAFitterNProngs) meanDA += minD; nfoundA++; } + ++fitstat[ft.getFitStatus()][0]; ft.setUseAbsDCA(true); ft.setWeightedFinalPCA(true); @@ -444,6 +482,7 @@ BOOST_AUTO_TEST_CASE(DCAFitterNProngs) meanDAW += minD; nfoundAW++; } + ++fitstat[ft.getFitStatus()][1]; ft.setUseAbsDCA(false); ft.setWeightedFinalPCA(false); @@ -456,6 +495,7 @@ BOOST_AUTO_TEST_CASE(DCAFitterNProngs) meanDW += minD; nfoundW++; } + ++fitstat[ft.getFitStatus()][2]; } // ft.print(); meanDA /= nfoundA ? nfoundA : 1; @@ -468,6 +508,7 @@ BOOST_AUTO_TEST_CASE(DCAFitterNProngs) << " mean.dist to truth: " << meanDAW << " CPU time: " << swAW.CpuTime() * 1000 << " ms"; LOG(info) << "2-prongs with wgh.dist minization: eff= " << float(nfoundW) / NTest << " mean.dist to truth: " << meanDW << " CPU time: " << swW.CpuTime() * 1000 << " ms"; + printStat(fitstat); BOOST_CHECK(nfoundA > 0.99 * NTest); BOOST_CHECK(nfoundAW > 0.99 * NTest); BOOST_CHECK(nfoundW > 0.99 * NTest); @@ -481,6 +522,8 @@ BOOST_AUTO_TEST_CASE(DCAFitterNProngs) { LOG(info) << "\n\nProcessing 3-prong vertices"; std::vector forceQ{1, 1, 1}; + std::memset(fitstat.data(), 0, sizeof(fitstat)); + o2::vertexing::DCAFitterN<3> ft; // 3 prong fitter ft.setBz(bz); ft.setPropagateToPCA(true); // After finding the vertex, propagate tracks to the DCA. This is default anyway @@ -509,6 +552,7 @@ BOOST_AUTO_TEST_CASE(DCAFitterNProngs) meanDA += minD; nfoundA++; } + ++fitstat[ft.getFitStatus()][0]; ft.setUseAbsDCA(true); ft.setWeightedFinalPCA(true); @@ -521,6 +565,7 @@ BOOST_AUTO_TEST_CASE(DCAFitterNProngs) meanDAW += minD; nfoundAW++; } + ++fitstat[ft.getFitStatus()][1]; ft.setUseAbsDCA(false); ft.setWeightedFinalPCA(false); @@ -533,6 +578,7 @@ BOOST_AUTO_TEST_CASE(DCAFitterNProngs) meanDW += minD; nfoundW++; } + ++fitstat[ft.getFitStatus()][2]; } // ft.print(); meanDA /= nfoundA ? nfoundA : 1; @@ -545,6 +591,7 @@ BOOST_AUTO_TEST_CASE(DCAFitterNProngs) << " mean.dist to truth: " << meanDAW << " CPU time: " << swAW.CpuTime() * 1000 << " ms"; LOG(info) << "3-prongs with wgh.dist minization: eff= " << float(nfoundW) / NTest << " mean.dist to truth: " << meanDW << " CPU time: " << swW.CpuTime() * 1000 << " ms"; + printStat(fitstat); BOOST_CHECK(nfoundA > 0.99 * NTest); BOOST_CHECK(nfoundAW > 0.99 * NTest); BOOST_CHECK(nfoundW > 0.99 * NTest); From 87fdd2ff8bddaeb948f01631a8577eac4b1ab319 Mon Sep 17 00:00:00 2001 From: Felix Schlepper Date: Wed, 2 Apr 2025 09:57:45 +0200 Subject: [PATCH 3/5] Common: DCAFitter catch maxIter reached Signed-off-by: Felix Schlepper --- .../DCAFitter/include/DCAFitter/DCAFitterN.h | 47 +++++++++++-------- 1 file changed, 28 insertions(+), 19 deletions(-) diff --git a/Common/DCAFitter/include/DCAFitter/DCAFitterN.h b/Common/DCAFitter/include/DCAFitter/DCAFitterN.h index a8bf1bb63bfb8..a967b8d7e3ca3 100644 --- a/Common/DCAFitter/include/DCAFitter/DCAFitterN.h +++ b/Common/DCAFitter/include/DCAFitter/DCAFitterN.h @@ -125,19 +125,24 @@ class DCAFitterN }; 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 + None, // no status set (should not be possible!) + + /* Good Conditions */ + Converged, // fit converged + MaxIter, // max iterations reached before fit convergence + + /* Error Conditions */ + NoCrossing, // no reasaonable crossing was found + RejRadius, // radius of crossing was not acceptable + RejTrackX, // one candidate track x was below the mimimum required radius + RejTrackRoughZ, // rejected by rough cut on tracks Z difference + RejChi2Max, // rejected by maximum chi2 cut + FailProp, // propagation of at least prong to PCA failed + FailInvCov, // inversion of cov.-matrix failed + FailInvWeight, // inversion of Ti weight matrix failed + FailInv2ndDeriv, // inversion of 2nd derivatives failed + FailCorrTracks, // correction of tracks to updated x failed + FailCloserAlt, // alternative PCA is closer }; static constexpr int getNProngs() { return N; } @@ -334,6 +339,10 @@ class DCAFitterN mCurHyp = 0; mAllowAltPreference = true; mOrder.fill(0); + mPropFailed.fill(false); + mTrPropDone.fill(false); + mNIters.fill(0); + mChi2.fill(-1); mFitStatus.fill(FitStatus::None); } @@ -431,9 +440,6 @@ GPUd() int DCAFitterN::process(const Tr&... args) mFitStatus[mCurHyp] = FitStatus::NoCrossing; return 0; } - for (int ih = 0; ih < MAXHYP; ih++) { - mPropFailed[ih] = false; - } if (mUseAbsDCA) { calcRMatrices(); // needed for fast residuals derivatives calculation in case of abs. distance minimization } @@ -455,9 +461,6 @@ GPUd() int DCAFitterN::process(const Tr&... args) } mCrossIDCur = ic; mCrossIDAlt = (mCrossings.nDCA == 2 && mAllowAltPreference) ? 1 - ic : -1; // works for max 2 crossings - mNIters[mCurHyp] = 0; - mTrPropDone[mCurHyp] = false; - mChi2[mCurHyp] = -1.; mPCA[mCurHyp][0] = mCrossings.xDCA[ic]; mPCA[mCurHyp][1] = mCrossings.yDCA[ic]; @@ -1033,6 +1036,9 @@ GPUd() bool DCAFitterN::minimizeChi2() } chi2 = chi2Upd; } while (++mNIters[mCurHyp] < mMaxIter); + if (mNIters[mCurHyp] == mMaxIter) { + mFitStatus[mCurHyp] = FitStatus::MaxIter; + } // mChi2[mCurHyp] = chi2 * NInv; if (mChi2[mCurHyp] >= mMaxChi2) { @@ -1102,6 +1108,9 @@ GPUd() bool DCAFitterN::minimizeChi2NoErr() } chi2 = chi2Upd; } while (++mNIters[mCurHyp] < mMaxIter); + if (mNIters[mCurHyp] == mMaxIter) { + mFitStatus[mCurHyp] = FitStatus::MaxIter; + } // mChi2[mCurHyp] = chi2 * NInv; if (mChi2[mCurHyp] >= mMaxChi2) { From 0b2ae984d1465c781236e19f3e0b39a5c0f9527d Mon Sep 17 00:00:00 2001 From: Felix Schlepper Date: Wed, 2 Apr 2025 10:40:58 +0200 Subject: [PATCH 4/5] Common: DCAFitter fix spell Signed-off-by: Felix Schlepper --- Common/DCAFitter/include/DCAFitter/DCAFitterN.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Common/DCAFitter/include/DCAFitter/DCAFitterN.h b/Common/DCAFitter/include/DCAFitter/DCAFitterN.h index a967b8d7e3ca3..a3f930de8508a 100644 --- a/Common/DCAFitter/include/DCAFitter/DCAFitterN.h +++ b/Common/DCAFitter/include/DCAFitter/DCAFitterN.h @@ -1082,7 +1082,7 @@ GPUd() bool DCAFitterN::minimizeChi2NoErr() // do Newton-Rapson iteration with corrections = - dchi2/d{x0..xN} * [ d^2chi2/d{x0..xN}^2 ]^-1 if (!mD2Chi2Dx2.Invert()) { if (mLoggerBadInv.needToLog()) { - printf("itter %d: error (%ld muted): Inversion failed\n", mFitterID, mLoggerBadCov.getNMuted()); + printf("fitter %d: error (%ld muted): Inversion failed\n", mFitterID, mLoggerBadCov.getNMuted()); mLoggerBadInv.evCountPrev = mLoggerBadInv.evCount; } mFitStatus[mCurHyp] = FitStatus::FailInv2ndDeriv; From 4774ec24a0f678f42aa26f5aef2d1d42031f20ae Mon Sep 17 00:00:00 2001 From: Felix Schlepper Date: Wed, 2 Apr 2025 11:00:17 +0200 Subject: [PATCH 5/5] Common: DCAFitter simplify exp.-backoff of logging Signed-off-by: Felix Schlepper --- .../DCAFitter/include/DCAFitter/DCAFitterN.h | 43 ++++++++----------- 1 file changed, 18 insertions(+), 25 deletions(-) diff --git a/Common/DCAFitter/include/DCAFitter/DCAFitterN.h b/Common/DCAFitter/include/DCAFitter/DCAFitterN.h index a3f930de8508a..569b3ea49e515 100644 --- a/Common/DCAFitter/include/DCAFitter/DCAFitterN.h +++ b/Common/DCAFitter/include/DCAFitter/DCAFitterN.h @@ -75,21 +75,20 @@ struct TrackDeriv { ///< Log log-throttling helper struct LogLogThrottler { size_t evCount{0}; - size_t evCountPrev{0}; - size_t logCount{0}; - + size_t nextLog{1}; GPUdi() bool needToLog() { - if (size_t(o2::gpu::GPUCommonMath::Log(++evCount)) + 1 > logCount) { - logCount++; + if (++evCount > nextLog) { + nextLog *= 2; return true; } return false; } - - GPUdi() size_t getNMuted() const { return evCount - evCountPrev - 1; } - - GPUdi() void clear() { evCount = evCountPrev = logCount = 0; } + GPUdi() void clear() + { + evCount = 0; + nextLog = 1; + } }; template @@ -747,12 +746,11 @@ GPUd() bool DCAFitterN::recalculatePCAWithErrors(int cand) 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()); + mFitterID, mLoggerBadCov.evCount, 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()); + mFitterID, mLoggerBadCov.evCount, mCandTr[mCurHyp][i].getSigmaY2(), mCandTr[mCurHyp][i].getSigmaZ2(), mCandTr[mCurHyp][i].getSigmaZY()); #endif - mLoggerBadCov.evCountPrev = mLoggerBadCov.evCount; } mFitStatus[mCurHyp] = FitStatus::FailInvCov; if (mBadCovPolicy == Discard) { @@ -975,12 +973,11 @@ GPUd() bool DCAFitterN::minimizeChi2() 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()); + mFitterID, mLoggerBadCov.evCount, 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()); + mFitterID, mLoggerBadCov.evCount, mCandTr[mCurHyp][i].getSigmaY2(), mCandTr[mCurHyp][i].getSigmaZ2(), mCandTr[mCurHyp][i].getSigmaZY()); #endif - mLoggerBadCov.evCountPrev = mLoggerBadCov.evCount; } mFitStatus[mCurHyp] = FitStatus::FailInvCov; if (mBadCovPolicy == Discard) { @@ -1010,8 +1007,7 @@ GPUd() bool DCAFitterN::minimizeChi2() // do Newton-Rapson iteration with corrections = - dchi2/d{x0..xN} * [ d^2chi2/d{x0..xN}^2 ]^-1 if (!mD2Chi2Dx2.Invert()) { if (mLoggerBadInv.needToLog()) { - printf("fitter %d: error (%ld muted): Inversion failed\n", mFitterID, mLoggerBadCov.getNMuted()); - mLoggerBadInv.evCountPrev = mLoggerBadInv.evCount; + printf("fitter %d: error (%ld muted): Inversion failed\n", mFitterID, mLoggerBadCov.evCount); } mFitStatus[mCurHyp] = FitStatus::FailInv2ndDeriv; return false; @@ -1082,8 +1078,7 @@ GPUd() bool DCAFitterN::minimizeChi2NoErr() // do Newton-Rapson iteration with corrections = - dchi2/d{x0..xN} * [ d^2chi2/d{x0..xN}^2 ]^-1 if (!mD2Chi2Dx2.Invert()) { if (mLoggerBadInv.needToLog()) { - printf("fitter %d: error (%ld muted): Inversion failed\n", mFitterID, mLoggerBadCov.getNMuted()); - mLoggerBadInv.evCountPrev = mLoggerBadInv.evCount; + printf("fitter %d: error (%ld muted): Inversion failed\n", mFitterID, mLoggerBadCov.evCount); } mFitStatus[mCurHyp] = FitStatus::FailInv2ndDeriv; return false; @@ -1247,11 +1242,10 @@ GPUdi() bool DCAFitterN::propagateParamToX(o2::track::TrackPar& t, f 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()); + printf("fitter %d: error (%ld muted): propagation failed for %s\n", mFitterID, mLoggerBadProp.evCount, t.asString().c_str()); #else - printf("fitter %d: error (%ld muted): propagation failed\n", mFitterID, mLoggerBadProp.getNMuted()); + printf("fitter %d: error (%ld muted): propagation failed\n", mFitterID, mLoggerBadProp.evCount); #endif - mLoggerBadProp.evCountPrev = mLoggerBadProp.evCount; } } return res; @@ -1274,11 +1268,10 @@ GPUdi() bool DCAFitterN::propagateToX(o2::track::TrackParCov& t, flo 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()); + printf("fitter %d: error (%ld muted): propagation failed for %s\n", mFitterID, mLoggerBadProp.evCount, t.asString().c_str()); #else - printf("fitter %d: error (%ld muted): propagation failed\n", mFitterID, mLoggerBadProp.getNMuted()); + printf("fitter %d: error (%ld muted): propagation failed\n", mFitterID, mLoggerBadProp.evCount); #endif - mLoggerBadProp.evCountPrev = mLoggerBadProp.evCount; } } return res;