From 4ce08135db46a3335702c1f70950f9269747aa89 Mon Sep 17 00:00:00 2001 From: shahoian Date: Sat, 18 Jul 2026 21:36:49 +0200 Subject: [PATCH 1/3] Fix x-axis related error treatments Inspired by the refining (@galocco) of similar errors in the na60+ adaptation of DCAFitter (@fprino). List of changes: 1. TrackCovI now stores a full symmetric 3D information matrix. The old dummy X variance is removed. The X information is derived from the track slopes dY/dX and dZ/dX via H^T C_YZ^{-1} H. 2. XerrFactor was removed. The fit no longer creates an artificial longitudinal uncertainty from sigma_Y. 3. getTrackCovMatrix() was removed. PCA covariance is now computed from summed track information matrices, not by rotating an artificial 3D covariance with dummy X error. 4. calcPCACoefs() now multiplies the full local information matrix. Terms involving local XY and XZ correlations are included in the Ti coefficient matrices. 5. calcInverseWeight() now rotates and sums the full information matrix. All six local information elements contribute to the global 3x3 weight. 6. calcChi2Derivatives() now contracts residual derivatives with the full information matrix. The Hessian curvature term is added only to diagonal Hessian elements, matching the fact that each trajectory has second derivative only with respect to its own running parameter. 7. calcChi2DerivativesNoErr() was corrected analogously. The Gauss-Newton term contributes to all Hessian elements, while the residual-times-second-derivative term contributes only to diagonal elements. 8. calcPCACovMatrix() now returns the inverse of summed track information. The longitudinal vertex uncertainty comes from track slopes rather than from a dummy X covariance. 9. calcChi2() now uses all six information-matrix elements. Residual correlations involving X are included in the chi2 value. 10. correctTracks() now propagates mCandTr analytically to the corrected X. This keeps the candidate track state, slopes, covariance, and mTrPos synchronized between Newton iterations while avoiding full-field Propagator calls for the small Newton corrections. 11. createParentTrackParCov() now builds the parent momentum covariance from the native O2 momentum parameters (snp,tgl,q/pt). The Jacobian includes the track-frame alpha rotation, and daughter momentum covariances are summed in lab px,py,pz before the existing O2 parent TrackParCov constructor converts to the parent frame. --- Common/DCAFitter/DCAFitterN_derivation.md | 447 ++++++++++++++++++ .../DCAFitter/include/DCAFitter/DCAFitterN.h | 214 +++++---- 2 files changed, 577 insertions(+), 84 deletions(-) create mode 100644 Common/DCAFitter/DCAFitterN_derivation.md diff --git a/Common/DCAFitter/DCAFitterN_derivation.md b/Common/DCAFitter/DCAFitterN_derivation.md new file mode 100644 index 0000000000000..83ba1646dfa20 --- /dev/null +++ b/Common/DCAFitter/DCAFitterN_derivation.md @@ -0,0 +1,447 @@ +# DCAFitterN derivation notes + +This file combines the extracted content of +`~/Downloads/DCAFitterN_derivation_for_codex.md` with the O2-specific covariance +updates made in `include/DCAFitter/DCAFitterN.h`. + +## 1. Convention check from the extracted scan + +The extracted scan starts with a matrix named `M_i` and states that it maps local +track coordinates to global coordinates: + +```text +p_i^g = M_i p_i +``` + +but the printed 2D block + +```text +[ cos(alpha_i) sin(alpha_i) ] +[ -sin(alpha_i) cos(alpha_i) ] +``` + +is the inverse of the O2 local-to-global rotation. O2 uses + +```text +R_i = +[ cos(alpha_i) -sin(alpha_i) 0 ] +[ sin(alpha_i) cos(alpha_i) 0 ] +[ 0 0 1 ] +``` + +with + +```text +p_i^g = R_i p_i, p_i = R_i^T p_i^g . +``` + +The later extracted formulas note this possible sign reversal. The code in +`DCAFitterN.h` matches the O2 convention above. In the rest of this document +`R_i` is always the O2 local-to-global rotation. + +## 2. Weighted N-prong vertex fit + +For track `i`, let the current point in its local frame be + +```text +p_i = (x_i, y_i, z_i)^T . +``` + +Let `B_i` be the inverse covariance, or information matrix, for this local +point. The fitted vertex `V` is in the global frame. Its representation in +track `i`'s local frame is `R_i^T V`, so the residual is + +```text +Delta_i = p_i - R_i^T V . +``` + +The weighted objective is + +```text +chi2 = 1/2 sum_i Delta_i^T B_i Delta_i . +``` + +For fixed track points, differentiating with respect to `V` gives + +```text +A V = sum_i R_i B_i p_i +``` + +where + +```text +A = sum_i R_i B_i R_i^T . +``` + +Therefore + +```text +V = A^{-1} sum_i R_i B_i p_i . +``` + +Define + +```text +T_i = A^{-1} R_i B_i . +``` + +Then + +```text +V = sum_i T_i p_i . +``` + +This is the `calcInverseWeight`, `calcPCACoefs`, and `calcPCA` structure in the +code. The useful identity is + +```text +sum_i T_i R_i^T = I . +``` + +## 3. Residuals after eliminating the vertex + +Substituting `V = sum_k T_k p_k` into the residual gives + +```text +Delta_i = sum_k D_ik p_k +``` + +with + +```text +D_ik = delta_ik I - R_i^T T_k . +``` + +The fit parameters are the local running coordinates `x_k`. During one Newton +linearization, `R_i`, `B_i`, `A`, and `T_i` are treated as fixed. Each track +point depends only on its own parameter: + +```text +d p_i / d x_k = delta_ik p_i' +``` + +so + +```text +d Delta_i / d x_k = D_ik p_k' +``` + +and + +```text +d^2 Delta_i / d x_k d x_l = delta_kl D_ik p_k'' . +``` + +## 4. Track derivatives in the O2 local frame + +O2 central-barrel parameters are + +```text +(Y, Z, snp, tgl, q/pt) +``` + +where + +```text +snp = sin(phi_local), csp = sqrt(1 - snp^2), +tgl = tan(lambda), kappa = curvature = (q/pt) Bz B2C . +``` + +For the fast constant-`Bz` helix model: + +```text +d snp / dX = kappa +dY / dX = snp / csp +dZ / dX = tgl / csp +``` + +and + +```text +d2Y / dX2 = kappa / csp^3 +d2Z / dX2 = kappa tgl snp / csp^3 . +``` + +Thus + +```text +p_i' = (1, dY/dX, dZ/dX)^T +p_i'' = (0, d2Y/dX2, d2Z/dX2)^T . +``` + +This matches `TrackDeriv::set`. + +## 5. Gradient and Hessian + +With symmetric `B_i`, + +```text +g_k = d chi2 / d x_k + = sum_i Delta_i^T B_i D_ik p_k' . +``` + +The exact Hessian is + +```text +H_kl = + sum_i (D_il p_l')^T B_i (D_ik p_k') + + delta_kl sum_i Delta_i^T B_i D_ik p_k'' . +``` + +The first term is the Gauss-Newton term. It contributes to diagonal and mixed +Hessian elements. The second term is the residual-curvature term. Since each +trajectory has an intrinsic second derivative only with respect to its own +running coordinate, this term contributes only when `k == l`. This is the +reason for the code condition + +```cpp +if (i == j) { + ... +} +``` + +when computing the Hessian element `H_ij`. + +The implementation solves + +```text +H dX = g +``` + +and then applies + +```text +X_new = X_old - dX . +``` + +This is equivalent to the more usual Newton notation `deltaX = -H^{-1} g`. + +## 6. No-error special case + +For the absolute-distance fit, take all information matrices as identity: + +```text +B_i = I . +``` + +Then + +```text +A = N I, A^{-1} = (1/N) I, +T_i = (1/N) R_i . +``` + +The vertex is the average of global track points: + +```text +V = (1/N) sum_i R_i p_i . +``` + +The residual is + +```text +Delta_i = p_i - (1/N) sum_j R_i^T R_j p_j . +``` + +Define + +```text +R_ij = (1/N) R_i^T R_j . +``` + +With the O2 convention, + +```text +R_i^T R_j = +[ cos(ai-aj) sin(ai-aj) 0 ] +[ -sin(ai-aj) cos(ai-aj) 0 ] +[ 0 0 1 ] . +``` + +This matches the code definitions + +```cpp +mCosDif[i][j] = (ci*cj + si*sj) / N; +mSinDif[i][j] = (si*cj - ci*sj) / N; +``` + +and the residual derivative components in `calcResidDerivativesNoErr`. + +## 7. Track information matrix involving the local X axis + +The old DCAFitterN code assigned a dummy uncertainty to local `X`, derived from +the local `Y` uncertainty. The corrected treatment derives longitudinal vertex +information from the track geometry. + +At fixed local `X`, the track measures `(Y,Z)` with covariance + +```text +C_YZ = +[ C_YY C_YZ ] +[ C_YZ C_ZZ ] . +``` + +Let + +```text +D = C_YY C_ZZ - C_YZ^2 +``` + +and + +```text +W = C_YZ^{-1} + = 1/D [ C_ZZ -C_YZ ] + [ -C_YZ C_YY ] . +``` + +Writing + +```text +wYY = C_ZZ / D +wYZ = -C_YZ / D +wZZ = C_YY / D +``` + +and + +```text +y' = dY/dX = snp/csp +z' = dZ/dX = tgl/csp +``` + +the vertex measurement matrix in local `(X,Y,Z)` coordinates is + +```text +H = +[ -y' 1 0 ] +[ -z' 0 1 ] . +``` + +The local 3D information matrix is + +```text +I = H^T W H . +``` + +Its independent elements are + +```text +I_YY = wYY +I_YZ = wYZ +I_ZZ = wZZ +I_XY = -(wYY y' + wYZ z') +I_XZ = -(wYZ y' + wZZ z') +I_XX = y'^2 wYY + 2 y' z' wYZ + z'^2 wZZ . +``` + +These are the six members of `TrackCovI`: + +```text +sxx, sxy, sxz, syy, syz, szz . +``` + +To combine tracks in the global frame, rotate the local information: + +```text +I_i^g = R_i I_i R_i^T . +``` + +The vertex covariance is then + +```text +C_V = (sum_i I_i^g)^{-1} . +``` + +## 8. Parent momentum covariance + +The parent momentum is the sum of independent daughter momenta: + +```text +P = sum_i p_i, C_P = sum_i C_{p_i} . +``` + +For one O2 daughter track, + +```text +px = pt (csp cos(alpha) - snp sin(alpha)) +py = pt (snp cos(alpha) + csp sin(alpha)) +pz = pt tgl +``` + +with `pt = |q|/|q/pt|`. The charge is treated as exact and discrete. The +derivatives with respect to native momentum parameters `(snp, tgl, q/pt)` are + +```text +dpx/dsnp = -pt (snp cos(alpha)/csp + sin(alpha)) +dpy/dsnp = pt (cos(alpha) - snp sin(alpha)/csp) +dpz/dtgl = pt + +dpx/d(q/pt) = -px / (q/pt) +dpy/d(q/pt) = -py / (q/pt) +dpz/d(q/pt) = -pz / (q/pt) +``` + +and the omitted mixed derivatives in this Jacobian are zero: + +```text +dpx/dtgl = 0, dpy/dtgl = 0, +dpz/dsnp = 0 . +``` + +Let `A` be this 3 by 3 Jacobian and let `C_a` be the daughter covariance +submatrix in the native momentum-parameter order `(snp,tgl,q/pt)`. Then + +```text +C_p = A C_a A^T . +``` + +The six independent elements of `C_p` are accumulated into the O2 lab covariance +slots + +```text +cov[9], cov[13], cov[14], cov[18], cov[19], cov[20] +``` + +which correspond to + +```text +Cov(px,px), Cov(py,px), Cov(py,py), +Cov(pz,px), Cov(pz,py), Cov(pz,pz). +``` + +## 9. Parent track covariance + +The final lab covariance passed to the O2 parent constructor contains the fitted +vertex covariance and the summed parent momentum covariance: + +```text +C_lab = +[ C_V 0 ] +[ 0 C_P ] . +``` + +Position-momentum cross-covariances are not included in this approximation. +The existing O2 constructor + +```cpp +TrackParCov(xyz, pxpypz, cov, charge, sectorAlpha) +``` + +then transforms this lab covariance to the selected parent track frame, +including the parent `alpha` convention. + +## 10. Code changes summarized + +1. `TrackCovI` now stores a full symmetric local 3D information matrix. +2. The dummy `XerrFactor` approximation was removed. +3. PCA weights and chi2 derivatives now use `sxx,sxy,sxz,syy,syz,szz`. +4. PCA covariance is the inverse of the summed global information matrix. +5. The Hessian residual-curvature term is restricted to diagonal Hessian + elements. +6. `correctTracks()` now propagates the actual candidate track state with O2's + analytic constant-`Bz` transport, keeping `mCandTr` and `mTrPos` + synchronized. +7. `createParentTrackParCov()` now propagates daughter momentum covariance from + native O2 `(snp,tgl,q/pt)` parameters to lab `(px,py,pz)` before constructing + the parent `TrackParCov`. diff --git a/Common/DCAFitter/include/DCAFitter/DCAFitterN.h b/Common/DCAFitter/include/DCAFitter/DCAFitterN.h index 2641dec84aed9..a78bfd7cdca30 100644 --- a/Common/DCAFitter/include/DCAFitter/DCAFitterN.h +++ b/Common/DCAFitter/include/DCAFitter/DCAFitterN.h @@ -12,7 +12,8 @@ /// \file DCAFitterN.h /// \brief Defintions for N-prongs secondary vertex fit /// \author ruben.shahoyan@cern.ch -/// For the formulae derivation see /afs/cern.ch/user/s/shahoian/public/O2/DCAFitter/DCAFitterN.pdf +/// For the original derivation see /afs/cern.ch/user/s/shahoian/public/O2/DCAFitter/DCAFitterN.pdf +/// The AI-assisted readme is in DCAFitterN_derivation.md #ifndef _ALICEO2_DCA_FITTERN_ #define _ALICEO2_DCA_FITTERN_ @@ -28,17 +29,19 @@ namespace vertexing { ///__________________________________________________________________________________ -///< Inverse cov matrix (augmented by a dummy X error) of the point defined by the track +///< Inverse covariance matrix of the point defined by the track struct TrackCovI { - float sxx, syy, syz, szz; + // Independent elements of the symmetric 3D information matrix + // H^T Cyz^{-1} H. A track constrains Y and Z at a given X through + // H = {{-dY/dX, 1, 0}, {-dZ/dX, 0, 1}}. + float sxx, sxy, sxz, syy, syz, szz; GPUdDefault() TrackCovI() = default; - GPUd() bool set(const o2::track::TrackParCov& trc, float xerrFactor = 1.f) + GPUd() bool set(const o2::track::TrackParCov& trc) { - // 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; + // Invert the 2D covariance of the measured track position (Y,Z). + float cyy = trc.getSigmaY2(), czz = trc.getSigmaZ2(), cyz = trc.getSigmaZY(); float detYZ = cyy * czz - cyz * cyz; bool res = true; if (detYZ <= 0.) { @@ -47,10 +50,15 @@ struct TrackCovI { res = false; } auto detYZI = 1. / detYZ; - sxx = 1. / cxx; syy = czz * detYZI; syz = -cyz * detYZI; szz = cyy * detYZI; + const float cspI = 1.f / trc.getCsp(); + const float dydx = trc.getSnp() * cspI; + const float dzdx = trc.getTgl() * cspI; + sxy = -(syy * dydx + syz * dzdx); + sxz = -(syz * dydx + szz * dzdx); + sxx = dydx * dydx * syy + 2.f * dydx * dzdx * syz + dzdx * dzdx * szz; return res; } }; @@ -98,7 +106,6 @@ class DCAFitterN static constexpr double NMax = 4; static constexpr double NInv = 1. / N; static constexpr int MAXHYP = 2; - static constexpr float XerrFactor = 5.; // factor for conversion of track covYY to dummy covXX using Track = o2::track::TrackParCov; using TrackAuxPar = o2::track::TrackAuxPar; using CrossInfo = o2::track::CrossInfo; @@ -313,17 +320,6 @@ class DCAFitterN return mat; } - GPUd() MatSym3D getTrackCovMatrix(int i, int cand = 0) const // generate covariance matrix of track position, adding fake X error - { - const auto& trc = mCandTr[mOrder[cand]][i]; - MatSym3D mat; - mat(0, 0) = trc.getSigmaY2() * XerrFactor; - mat(1, 1) = trc.getSigmaY2(); - mat(2, 2) = trc.getSigmaZ2(); - mat(2, 1) = trc.getSigmaZY(); - return mat; - } - GPUd() void assign(int) {} template GPUd() void assign(int i, const T& t, const Tr&... args) @@ -502,13 +498,13 @@ GPUd() bool DCAFitterN::calcPCACoefs() const auto& taux = mTrAux[i]; const auto& tcov = mTrcEInv[mCurHyp][i]; MatStd3D miei; - miei[0][0] = taux.c * tcov.sxx; - miei[0][1] = -taux.s * tcov.syy; - miei[0][2] = -taux.s * tcov.syz; - miei[1][0] = taux.s * tcov.sxx; - miei[1][1] = taux.c * tcov.syy; - miei[1][2] = taux.c * tcov.syz; - miei[2][0] = 0; + miei[0][0] = taux.c * tcov.sxx - taux.s * tcov.sxy; + miei[0][1] = taux.c * tcov.sxy - taux.s * tcov.syy; + miei[0][2] = taux.c * tcov.sxz - taux.s * tcov.syz; + miei[1][0] = taux.s * tcov.sxx + taux.c * tcov.sxy; + miei[1][1] = taux.s * tcov.sxy + taux.c * tcov.syy; + miei[1][2] = taux.s * tcov.sxz + taux.c * tcov.syz; + miei[2][0] = tcov.sxz; miei[2][1] = tcov.syz; miei[2][2] = tcov.szz; mTrCFVT[mCurHyp][i] = mWeightInv * miei; @@ -532,11 +528,11 @@ GPUd() bool DCAFitterN::calcInverseWeight() for (int i = N; i--;) { const auto& taux = mTrAux[i]; const auto& tcov = mTrcEInv[mCurHyp][i]; - arrmat[XX] += taux.cc * tcov.sxx + taux.ss * tcov.syy; - arrmat[XY] += taux.cs * (tcov.sxx - tcov.syy); - arrmat[XZ] += -taux.s * tcov.syz; - arrmat[YY] += taux.cc * tcov.syy + taux.ss * tcov.sxx; - arrmat[YZ] += taux.c * tcov.syz; + arrmat[XX] += taux.cc * tcov.sxx - 2. * taux.cs * tcov.sxy + taux.ss * tcov.syy; + arrmat[XY] += taux.cs * (tcov.sxx - tcov.syy) + (taux.cc - taux.ss) * tcov.sxy; + arrmat[XZ] += taux.c * tcov.sxz - taux.s * tcov.syz; + arrmat[YY] += taux.ss * tcov.sxx + 2. * taux.cs * tcov.sxy + taux.cc * tcov.syy; + arrmat[YZ] += taux.s * tcov.sxz + taux.c * tcov.syz; arrmat[ZZ] += tcov.szz; } // invert 3x3 symmetrix matrix @@ -670,9 +666,9 @@ GPUd() void DCAFitterN::calcChi2Derivatives() const auto& covI = mTrcEInv[mCurHyp][j]; // inverse cov matrix of track j const auto& dr1 = mDResidDx[j][i]; // vector of j-th residuals 1st derivative over X param of track i auto& cidr = covIDrDx[i][j]; // vector covI_j * dres_j/dx_i, save for 2nd derivative calculation - cidr[0] = covI.sxx * dr1[0]; - cidr[1] = covI.syy * dr1[1] + covI.syz * dr1[2]; - cidr[2] = covI.syz * dr1[1] + covI.szz * dr1[2]; + cidr[0] = covI.sxx * dr1[0] + covI.sxy * dr1[1] + covI.sxz * dr1[2]; + cidr[1] = covI.sxy * dr1[0] + covI.syy * dr1[1] + covI.syz * dr1[2]; + cidr[2] = covI.sxz * dr1[0] + covI.syz * dr1[1] + covI.szz * dr1[2]; // calculate res_i * covI_j * dres_j/dx_i dchi1 += o2::math_utils::Dot(res, cidr); } @@ -686,11 +682,13 @@ GPUd() void DCAFitterN::calcChi2Derivatives() const auto& dr1j = mDResidDx[k][j]; // vector of k-th residuals 1st derivative over X param of track j const auto& cidrkj = covIDrDx[i][k]; // vector covI_k * dres_k/dx_i dchi2 += o2::math_utils::Dot(dr1j, cidrkj); - if (k == j) { + if (i == j) { const auto& res = mTrRes[mCurHyp][k]; // vector of residuals of track k const auto& covI = mTrcEInv[mCurHyp][k]; // inverse cov matrix of track k - const auto& dr2ij = mD2ResidDx2[k][j]; // vector of k-th residuals 2nd derivative over X params of track j - dchi2 += res[0] * covI.sxx * dr2ij[0] + res[1] * (covI.syy * dr2ij[1] + covI.syz * dr2ij[2]) + res[2] * (covI.syz * dr2ij[1] + covI.szz * dr2ij[2]); + const auto& dr2ij = mD2ResidDx2[k][i]; // vector of k-th residuals 2nd derivative over X param i + dchi2 += res[0] * (covI.sxx * dr2ij[0] + covI.sxy * dr2ij[1] + covI.sxz * dr2ij[2]) + + res[1] * (covI.sxy * dr2ij[0] + covI.syy * dr2ij[1] + covI.syz * dr2ij[2]) + + res[2] * (covI.sxz * dr2ij[0] + covI.syz * dr2ij[1] + covI.szz * dr2ij[2]); } } } @@ -705,16 +703,23 @@ GPUd() void DCAFitterN::calcChi2DerivativesNoErr() for (int i = N; i--;) { auto& dchi1 = mDChi2Dx[i]; // DChi2/Dx_i = sum_j { res_j * Dres_j/Dx_i } dchi1 = 0; // chi2 1st derivative - for (int j = N; j--;) { - const auto& res = mTrRes[mCurHyp][j]; // vector of residuals of track j - const auto& dr1 = mDResidDx[j][i]; // vector of j-th residuals 1st derivative over X param of track i + for (int k = N; k--;) { + const auto& res = mTrRes[mCurHyp][k]; // vector of residuals of track k + const auto& dr1 = mDResidDx[k][i]; // vector of k-th residuals 1st derivative over X param of track i dchi1 += o2::math_utils::Dot(res, dr1); - if (i >= j) { // symmetrix matrix - // chi2 2nd derivative - auto& dchi2 = mD2Chi2Dx2[i][j]; // D2Chi2/Dx_i/Dx_j = sum_k { Dres_k/Dx_j * covI_k * Dres_k/Dx_i + res_k * covI_k * D2res_k/Dx_i/Dx_j } - dchi2 = o2::math_utils::Dot(mTrRes[mCurHyp][i], mD2ResidDx2[i][j]); - for (int k = N; k--;) { - dchi2 += o2::math_utils::Dot(mDResidDx[k][i], mDResidDx[k][j]); + } + } + for (int i = N; i--;) { + for (int j = i + 1; j--;) { + auto& dchi2 = mD2Chi2Dx2[i][j]; + dchi2 = 0.; + for (int k = N; k--;) { + // Gauss-Newton term, present for diagonal and mixed elements. + dchi2 += o2::math_utils::Dot(mDResidDx[k][i], mDResidDx[k][j]); + // A trajectory has a second derivative only with respect to its own + // X parameter, hence the curvature term contributes only to H_ii. + if (i == j) { + dchi2 += o2::math_utils::Dot(mTrRes[mCurHyp][k], mD2ResidDx2[k][i]); } } } @@ -744,7 +749,7 @@ GPUd() bool DCAFitterN::recalculatePCAWithErrors(int cand) mCurHyp = mOrder[cand]; if (mUseAbsDCA) { for (int i = N; i--;) { - if (!mTrcEInv[mCurHyp][i].set(mCandTr[mCurHyp][i], XerrFactor)) { // prepare inverse cov.matrices at starting point + if (!mTrcEInv[mCurHyp][i].set(mCandTr[mCurHyp][i])) { // 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", @@ -799,28 +804,40 @@ GPUd() void DCAFitterN::calcPCANoErr() template GPUd() o2::math_utils::SMatrix> DCAFitterN::calcPCACovMatrix(int cand) const { - // calculate covariance matrix for the point of closest approach - MatSym3D covm; - int nAdded = 0; - for (int i = N; i--;) { // calculate sum of inverses - // MatSym3D covTr = o2::math_utils::Similarity(mUseAbsDCA ? getTrackRotMatrix(i) : mTrCFVT[mOrder[cand]][i], getTrackCovMatrix(i, cand)); - // RS by using Similarity(mTrCFVT[mOrder[cand]][i], getTrackCovMatrix(i, cand)) we underestimate the error, use simple rotation - MatSym3D covTr = o2::math_utils::Similarity(getTrackRotMatrix(i), getTrackCovMatrix(i, cand)); - if (covTr.Invert()) { - covm += covTr; - nAdded++; - } - } - if (nAdded && covm.Invert()) { - return covm; - } - // correct way has failed, use simple sum - MatSym3D covmSum; + // Each track measures Y and Z at the vertex X. With the local slopes + // sy = dY/dX and sz = dZ/dX, its vertex measurement matrix is + // H = {{-sy, 1, 0}, {-sz, 0, 1}}. The longitudinal information must come + // from the track geometry, not from a dummy X variance. + MatSym3D info; + auto* arrmat = info.Array(); + memset(arrmat, 0, sizeof(info)); + enum { XX, + XY, + YY, + XZ, + YZ, + ZZ }; + const int ord = mOrder[cand]; for (int i = N; i--;) { - MatSym3D covTr = o2::math_utils::Similarity(getTrackRotMatrix(i), getTrackCovMatrix(i, cand)); - covmSum += covTr; + const auto& taux = mTrAux[i]; + TrackCovI tcov; + tcov.set(mCandTr[ord][i]); + arrmat[XX] += taux.cc * tcov.sxx - 2. * taux.cs * tcov.sxy + taux.ss * tcov.syy; + arrmat[XY] += taux.cs * (tcov.sxx - tcov.syy) + (taux.cc - taux.ss) * tcov.sxy; + arrmat[XZ] += taux.c * tcov.sxz - taux.s * tcov.syz; + arrmat[YY] += taux.ss * tcov.sxx + 2. * taux.cs * tcov.sxy + taux.cc * tcov.syy; + arrmat[YZ] += taux.s * tcov.sxz + taux.c * tcov.syz; + arrmat[ZZ] += tcov.szz; + } + if (info.Invert()) { + return info; } - return covmSum; + // Fall back on identity diagonal matrix with 1 cm error. + memset(arrmat, 0, sizeof(info)); + info(0, 0) = 1.; + info(1, 1) = 1.; + info(2, 2) = 1.; + return info; } //___________________________________________________________________ @@ -856,7 +873,8 @@ GPUdi() double DCAFitterN::calcChi2() const for (int i = N; i--;) { const auto& res = mTrRes[mCurHyp][i]; const auto& covI = mTrcEInv[mCurHyp][i]; - chi2 += res[0] * res[0] * covI.sxx + res[1] * res[1] * covI.syy + res[2] * res[2] * covI.szz + 2. * res[1] * res[2] * covI.syz; + chi2 += res[0] * res[0] * covI.sxx + res[1] * res[1] * covI.syy + res[2] * res[2] * covI.szz + + 2. * (res[0] * res[1] * covI.sxy + res[0] * res[2] * covI.sxz + res[1] * res[2] * covI.syz); } return chi2; } @@ -878,13 +896,17 @@ GPUdi() double DCAFitterN::calcChi2NoErr() const template GPUd() bool DCAFitterN::correctTracks(const VecND& corrX) { - // propagate tracks to updated X + // Propagate the actual candidate tracks to the updated X. Use the analytic + // constant-Bz transport here: Newton corrections are small, but the track + // state must stay synchronized with mTrPos for the next derivative update. for (int i = N; i--;) { - const auto& trDer = mTrDer[mCurHyp][i]; - auto dx2h = 0.5 * corrX[i] * corrX[i]; - mTrPos[mCurHyp][i][0] -= corrX[i]; - mTrPos[mCurHyp][i][1] -= trDer.dydx * corrX[i] - dx2h * trDer.d2ydx2; - mTrPos[mCurHyp][i][2] -= trDer.dzdx * corrX[i] - dx2h * trDer.d2zdx2; + auto& trc = mCandTr[mCurHyp][i]; + const float x = static_cast(mTrPos[mCurHyp][i][0] - corrX[i]); + const bool propagated = mUseAbsDCA ? trc.propagateParamTo(x, mBz) : trc.propagateTo(x, mBz); + if (!propagated) { + return false; + } + setTrackPos(mTrPos[mCurHyp][i], trc); } return true; } @@ -972,7 +994,7 @@ GPUd() bool DCAFitterN::minimizeChi2() 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 (!mTrcEInv[mCurHyp][i].set(mCandTr[mCurHyp][i])) { // 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", @@ -1176,21 +1198,45 @@ GPUd() void DCAFitterN::print() const template GPUd() o2::track::TrackParCov DCAFitterN::createParentTrackParCov(int cand, bool sectorAlpha) const { - const auto& trP = getTrack(0, cand); - const auto& trN = getTrack(1, cand); std::array covV = {0.}; std::array pvecV = {0.}; int q = 0; for (int it = 0; it < N; it++) { const auto& trc = getTrack(it, cand); std::array pvecT = {0.}; - std::array covT = {0.}; - trc.getPxPyPzGlo(pvecT); - trc.getCovXYZPxPyPzGlo(covT); - constexpr int MomInd[6] = {9, 13, 14, 18, 19, 20}; // cov matrix elements for momentum component - for (int i = 0; i < 6; i++) { - covV[MomInd[i]] += covT[MomInd[i]]; + const bool hasMomentum = trc.getPxPyPzGlo(pvecT); + + // Propagate only the native momentum-parameter covariance + // (snp,tgl,q/pt) to the lab momentum covariance. The final constructor + // below rotates the summed lab covariance to the parent alpha frame. + if (hasMomentum) { + const double snp = trc.getSnp(); + const double csp = trc.getCsp(); + const double pt = trc.getPt(); + const double alpha = trc.getAlpha(); + double sna = 0., csa = 0.; + o2::math_utils::detail::sincos(alpha, sna, csa); + const double dPxdSnp = -pt * (snp * csa / csp + sna); + const double dPydSnp = pt * (csa - snp * sna / csp); + const double dPzdTgl = pt; + const double q2ptI = 1. / trc.getQ2Pt(); + const double dPxdQ = -pvecT[0] * q2ptI; + const double dPydQ = -pvecT[1] * q2ptI; + const double dPzdQ = -pvecT[2] * q2ptI; + const double cSnpSnp = trc.getSigmaSnp2(); + const double cTglSnp = trc.getSigmaTglSnp(); + const double cTglTgl = trc.getSigmaTgl2(); + const double cQSnp = trc.getSigma1PtSnp(); + const double cQTgl = trc.getSigma1PtTgl(); + const double cQQ = trc.getSigma1Pt2(); + covV[9] += dPxdSnp * dPxdSnp * cSnpSnp + 2. * dPxdSnp * dPxdQ * cQSnp + dPxdQ * dPxdQ * cQQ; + covV[13] += dPydSnp * (dPxdSnp * cSnpSnp + dPxdQ * cQSnp) + dPydQ * (dPxdSnp * cQSnp + dPxdQ * cQQ); + covV[14] += dPydSnp * dPydSnp * cSnpSnp + 2. * dPydSnp * dPydQ * cQSnp + dPydQ * dPydQ * cQQ; + covV[18] += dPzdTgl * (dPxdSnp * cTglSnp + dPxdQ * cQTgl) + dPzdQ * (dPxdSnp * cQSnp + dPxdQ * cQQ); + covV[19] += dPzdTgl * (dPydSnp * cTglSnp + dPydQ * cQTgl) + dPzdQ * (dPydSnp * cQSnp + dPydQ * cQQ); + covV[20] += dPzdTgl * dPzdTgl * cTglTgl + 2. * dPzdTgl * dPzdQ * cQTgl + dPzdQ * dPzdQ * cQQ; } + for (int i = 0; i < 3; i++) { pvecV[i] += pvecT[i]; } From be983611f3972e31a853968a011824e72aa0f47d Mon Sep 17 00:00:00 2001 From: shahoian Date: Sun, 19 Jul 2026 21:15:07 +0200 Subject: [PATCH 2/3] Add regularization for degenerate matrices --- Common/DCAFitter/DCAFitterN_derivation.md | 107 +++++++++++++++++- .../DCAFitter/include/DCAFitter/DCAFitterN.h | 105 ++++++++++++++--- 2 files changed, 198 insertions(+), 14 deletions(-) diff --git a/Common/DCAFitter/DCAFitterN_derivation.md b/Common/DCAFitter/DCAFitterN_derivation.md index 83ba1646dfa20..eb1c7e36c1fe2 100644 --- a/Common/DCAFitter/DCAFitterN_derivation.md +++ b/Common/DCAFitter/DCAFitterN_derivation.md @@ -431,7 +431,112 @@ TrackParCov(xyz, pxpypz, cov, charge, sectorAlpha) then transforms this lab covariance to the selected parent track frame, including the parent `alpha` convention. -## 10. Code changes summarized +## 10. Regularization of singular or weakly constrained covariance matrices + +### Single-track `TrackCovI` + +The matrix + +```text +I = H^T W H +``` + +has rank at most two for one track, because one track measures only the two +coordinates transverse to its own trajectory. The null direction is the local +track tangent + +```text +t = (1, y_prime, z_prime)^T . +``` + +Indeed + +```text +H t = 0, I t = 0 . +``` + +Therefore a single-track `TrackCovI` is positive semidefinite, not positive +definite. Inverting this 3D matrix by itself is not meaningful; the apparent +negative diagonal elements seen in such an inverse are numerical symptoms of +trying to invert a rank-deficient information matrix. The physically meaningful +inverse at single-track level is only the original 2 by 2 `(Y,Z)` covariance. + +For the multi-track vertex fit, different track directions usually make + +```text +A = sum_i R_i I_i R_i^T +``` + +full rank. However, nearly parallel or otherwise weak geometries can still make +the longitudinal eigenvalue very small. To avoid an exactly singular +single-track contribution while preserving the measured `(Y,Z)` block, the code +may add a weak positive longitudinal prior only to `I_XX`: + +```text +I_XX -> I_XX + 1/sigma_X,prior^2 . +``` + +In code this is implemented as + +```cpp +constexpr float XRegErrFactor = 10.f; +const float sigmaX2 = C_YY * XRegErrFactor; +sxx += 1.f / sigmaX2; +``` + +This is different from multiplying `I_XX` by a number below one. A reduction of +`I_XX` can make the matrix indefinite, while adding a positive diagonal term +keeps the information matrix positive definite in the tangent direction and does +not alter `I_YY`, `I_YZ`, or `I_ZZ`. + +The regularization should remain weak. It is a numerical stabilizer for badly +conditioned geometries, not an additional detector measurement of the local +track `X` coordinate. + +### `calcPCACovMatrix()` + +The vertex covariance is + +```text +C_V = A^-1, A = sum_i R_i I_i R_i^T . +``` + +If `A` is singular or ill-conditioned, returning a small fallback covariance +would incorrectly shrink the uncertainty along the weakly constrained direction. +The conservative behavior is: + +1. Check that `A` is compatible with a positive-definite information matrix. +2. Reject cases where the determinant is too small compared with the diagonal + scale. +3. Invert only when these checks pass. +4. If inversion fails, or if the inverted covariance has non-positive diagonal + elements, return a deliberately loose fallback covariance. + +The code uses the leading-minor checks + +```text +A_XX > 0, +det(A_XY block) > 0, +det(A) > epsilon max(diag(A))^3 +``` + +with + +```text +epsilon = 1e-12 . +``` + +On failure, it returns a diagonal covariance with + +```text +sigma^2 = 1e6 cm^2 . +``` + +This corresponds to a 10 m uncertainty in each coordinate. It is intentionally +loose: the fallback marks the parent vertex as poorly constrained rather than +creating an artificially precise parent covariance. + +## 11. Code changes summarized 1. `TrackCovI` now stores a full symmetric local 3D information matrix. 2. The dummy `XerrFactor` approximation was removed. diff --git a/Common/DCAFitter/include/DCAFitter/DCAFitterN.h b/Common/DCAFitter/include/DCAFitter/DCAFitterN.h index a78bfd7cdca30..d02efa68526db 100644 --- a/Common/DCAFitter/include/DCAFitter/DCAFitterN.h +++ b/Common/DCAFitter/include/DCAFitter/DCAFitterN.h @@ -59,6 +59,10 @@ struct TrackCovI { sxy = -(syy * dydx + syz * dzdx); sxz = -(syz * dydx + szz * dzdx); sxx = dydx * dydx * syy + 2.f * dydx * dzdx * syz + dzdx * dzdx * szz; + // The matrix is degenerate by construction, regularize sxx term to preserve the original YZ block exactly + constexpr float XRegErrFactor = 10.f; + const float sigmaX2 = cyy * XRegErrFactor; + sxx += 1.f / sigmaX2; return res; } }; @@ -220,6 +224,7 @@ class DCAFitterN ///< recalculate PCA as a cov-matrix weighted mean, even if absDCA method was used GPUd() bool recalculatePCAWithErrors(int cand = 0); + GPUd() double calcCollinearInflation(int cand) const; GPUd() MatSym3D calcPCACovMatrix(int cand = 0) const; std::array calcPCACovMatrixFlat(int cand = 0) const @@ -385,9 +390,10 @@ class DCAFitterN std::array mNIters; // number of iterations for each seed std::array mTrPropDone{}; // Flag that the tracks are fully propagated to PCA std::array mPropFailed{}; // Flag that some propagation failed for this PCA candidate - LogLogThrottler mLoggerBadCov{}; - LogLogThrottler mLoggerBadInv{}; - LogLogThrottler mLoggerBadProp{}; + mutable LogLogThrottler mLoggerBadCov{}; + mutable LogLogThrottler mLoggerBadInv{}; + mutable LogLogThrottler mLoggerBadProp{}; + mutable LogLogThrottler mLoggerBadPCACov{}; MatSym3D mWeightInv; // inverse weight of single track, [sum{M^T E M}]^-1 in EQ.T std::array mOrder{0}; int mCurHyp = 0; @@ -800,6 +806,50 @@ GPUd() void DCAFitterN::calcPCANoErr() pca[2] *= NInv; } +//___________________________________________________________________ +template +GPUd() double DCAFitterN::calcCollinearInflation(int cand) const +{ + std::array, N> u{}; + int nu = 0; + + for (int i = 0; i < N; ++i) { + std::array p{}; + if (!getTrack(i, cand).getPxPyPzGlo(p)) { + continue; + } + const double p2 = p[0] * p[0] + p[1] * p[1] + p[2] * p[2]; + if (p2 <= 0.) { + continue; + } + const double pI = 1. / std::sqrt(p2); + u[nu++] = {p[0] * pI, p[1] * pI, p[2] * pI}; + } + + if (nu < 2) { + return 1.; + } + + double sin2Mean = 0.; + int npairs = 0; + for (int i = 0; i < nu; ++i) { + for (int j = i + 1; j < nu; ++j) { + double cij = u[i][0] * u[j][0] + u[i][1] * u[j][1] + u[i][2] * u[j][2]; + cij = std::clamp(cij, -1., 1.); + sin2Mean += std::max(0., 1. - cij * cij); + ++npairs; + } + } + sin2Mean /= npairs; + + constexpr double Sin2Ref = 1.e-5; + constexpr double MaxInflation = 1.e4; + if (sin2Mean <= 0.) { + return MaxInflation; + } + return sin2Mean < Sin2Ref ? std::min(MaxInflation, Sin2Ref / sin2Mean) : 1.; +} + //___________________________________________________________________ template GPUd() o2::math_utils::SMatrix> DCAFitterN::calcPCACovMatrix(int cand) const @@ -829,14 +879,34 @@ GPUd() o2::math_utils::SMatrix 0. && info(0, 0) > 0. && det2 > 0. && det3 > MinRelDet * maxDiag * maxDiag * maxDiag; + if (isWellConditionedInfo) { + auto cov = info; + if (cov.Invert() && cov(0, 0) > 0. && cov(1, 1) > 0. && cov(2, 2) > 0.) { + // if (mIsCollinear) { + // cov *= calcCollinearInflation(cand); + // } + return cov; + } + } + if (mLoggerBadPCACov.needToLog()) { + printf("fitter %d: error (%ld muted): override ill-conditioned PCACovMatrix by dummy matrix", mFitterID, mLoggerBadPCACov.evCount); } - // Fall back on identity diagonal matrix with 1 cm error. + // Fall back on a deliberately loose vertex covariance. Returning a tight + // identity covariance for a singular or ill-conditioned information matrix + // would shrink the uncertainty in the weakly constrained direction. memset(arrmat, 0, sizeof(info)); - info(0, 0) = 1.; - info(1, 1) = 1.; - info(2, 2) = 1.; + info(0, 0) = 4.; + info(1, 1) = 4.; + info(2, 2) = 4.; return info; } @@ -900,6 +970,15 @@ GPUd() bool DCAFitterN::correctTracks(const VecND& corrX) // constant-Bz transport here: Newton corrections are small, but the track // state must stay synchronized with mTrPos for the next derivative update. for (int i = N; i--;) { + /* + // Updating only mTrPos by Taylor expansion leaves mCandTr at the previous X, + // leaving calcTrackDerivatives() insensitive to the update. Use full fast propagation instead. + const auto& trDer = mTrDer[mCurHyp][i]; + auto dx2h = 0.5 * corrX[i] * corrX[i]; + mTrPos[mCurHyp][i][0] -= corrX[i]; + mTrPos[mCurHyp][i][1] -= trDer.dydx * corrX[i] - dx2h * trDer.d2ydx2; + mTrPos[mCurHyp][i][2] -= trDer.dzdx * corrX[i] - dx2h * trDer.d2zdx2; + */ auto& trc = mCandTr[mCurHyp][i]; const float x = static_cast(mTrPos[mCurHyp][i][0] - corrX[i]); const bool propagated = mUseAbsDCA ? trc.propagateParamTo(x, mBz) : trc.propagateTo(x, mBz); @@ -1291,9 +1370,9 @@ 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.evCount, t.asString().c_str()); + printf("fitter %d: error (%ld muted): propagation to %.4f failed for %s\n", mFitterID, mLoggerBadProp.evCount, x, t.asString().c_str()); #else - printf("fitter %d: error (%ld muted): propagation failed\n", mFitterID, mLoggerBadProp.evCount); + printf("fitter %d: error (%ld muted): propagation to %.4f failed\n", mFitterID, mLoggerBadProp.evCount, x); #endif } } @@ -1317,9 +1396,9 @@ 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.evCount, t.asString().c_str()); + printf("fitter %d: error (%ld muted): propagation to %.4f failed for %s\n", mFitterID, mLoggerBadProp.evCount, x, t.asString().c_str()); #else - printf("fitter %d: error (%ld muted): propagation failed\n", mFitterID, mLoggerBadProp.evCount); + printf("fitter %d: error (%ld muted): propagation to %.4f failed\n", mFitterID, mLoggerBadProp.evCount, x); #endif } } From 95ea7b2b3aa8871ca1e087e8786883ff6be2b3ec Mon Sep 17 00:00:00 2001 From: shahoian Date: Sun, 19 Jul 2026 21:15:44 +0200 Subject: [PATCH 3/3] Add pulls and gen.parent TrackParCov to test output --- Common/DCAFitter/test/testDCAFitterN.cxx | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/Common/DCAFitter/test/testDCAFitterN.cxx b/Common/DCAFitter/test/testDCAFitterN.cxx index bd00b5bed841e..8ab2c7f4ca323 100644 --- a/Common/DCAFitter/test/testDCAFitterN.cxx +++ b/Common/DCAFitter/test/testDCAFitterN.cxx @@ -56,11 +56,21 @@ float checkResults(o2::utils::TreeStreamRedirector& outs, std::string& treeName, double dst = TMath::Sqrt(df[0] * df[0] + df[1] * df[1] + df[2] * df[2]); distMin = dst < distMin ? dst : distMin; auto parentTrack = fitter.createParentTrackParCov(ic); - // float genX + const std::array genPos{static_cast(vgen[0]), static_cast(vgen[1]), static_cast(vgen[2])}; + const std::array genMom{static_cast(genPar.Px()), static_cast(genPar.Py()), static_cast(genPar.Pz())}; + o2::track::TrackPar genParentTrack(genPos, genMom, parentTrack.getCharge(), false); + genParentTrack.rotateParam(parentTrack.getAlpha()); + std::array parentCovGlo{}; + const bool hasParentCovGlo = parentTrack.getCovXYZPxPyPzGlo(parentCovGlo); + const double pullX = hasParentCovGlo && parentCovGlo[0] > 0.f ? df[0] / TMath::Sqrt(parentCovGlo[0]) : 0.; + const double pullY = hasParentCovGlo && parentCovGlo[2] > 0.f ? df[1] / TMath::Sqrt(parentCovGlo[2]) : 0.; + const double pullZ = hasParentCovGlo && parentCovGlo[5] > 0.f ? df[2] / TMath::Sqrt(parentCovGlo[5]) : 0.; outs << treeName.c_str() << "cand=" << ic << "ncand=" << nCand << "nIter=" << nIter << "chi2=" << chi2 << "genPart=" << genPar << "recPart=" << moth << "genX=" << vgen[0] << "genY=" << vgen[1] << "genZ=" << vgen[2] << "dx=" << df[0] << "dy=" << df[1] << "dz=" << df[2] << "dst=" << dst + << "pullX=" << pullX << "pullY=" << pullY << "pullZ=" << pullZ + << "genParentTrack=" << genParentTrack << "useAbsDCA=" << absDCA << "useWghDCA=" << useWghDCA << "parent=" << parentTrack; for (int i = 0; i < fitter.getNProngs(); i++) { outs << treeName.c_str() << fmt::format("prong{}=", i).c_str() << fitter.getTrack(i, ic);