diff --git a/Common/DCAFitter/DCAFitterN_derivation.md b/Common/DCAFitter/DCAFitterN_derivation.md new file mode 100644 index 0000000000000..eb1c7e36c1fe2 --- /dev/null +++ b/Common/DCAFitter/DCAFitterN_derivation.md @@ -0,0 +1,552 @@ +# 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. 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. +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..d02efa68526db 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,19 @@ 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; + // 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; } }; @@ -98,7 +110,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; @@ -213,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 @@ -313,17 +325,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) @@ -389,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; @@ -502,13 +504,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 +534,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 +672,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 +688,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 +709,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 +755,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", @@ -797,30 +808,106 @@ GPUd() void DCAFitterN::calcPCANoErr() //___________________________________________________________________ template -GPUd() o2::math_utils::SMatrix> DCAFitterN::calcPCACovMatrix(int cand) const +GPUd() double DCAFitterN::calcCollinearInflation(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++; + 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; } } - if (nAdded && covm.Invert()) { - return covm; + sin2Mean /= npairs; + + constexpr double Sin2Ref = 1.e-5; + constexpr double MaxInflation = 1.e4; + if (sin2Mean <= 0.) { + return MaxInflation; } - // correct way has failed, use simple sum - MatSym3D covmSum; + return sin2Mean < Sin2Ref ? std::min(MaxInflation, Sin2Ref / sin2Mean) : 1.; +} + +//___________________________________________________________________ +template +GPUd() o2::math_utils::SMatrix> DCAFitterN::calcPCACovMatrix(int cand) const +{ + // 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; } - return covmSum; + const double maxDiag = o2::gpu::GPUCommonMath::Max(o2::gpu::GPUCommonMath::Max(info(0, 0), info(1, 1)), info(2, 2)); + const double det2 = info(0, 0) * info(1, 1) - info(1, 0) * info(1, 0); + const double det3 = info(0, 0) * (info(1, 1) * info(2, 2) - info(2, 1) * info(2, 1)) - + info(1, 0) * (info(1, 0) * info(2, 2) - info(2, 1) * info(2, 0)) + + info(2, 0) * (info(1, 0) * info(2, 1) - info(1, 1) * info(2, 0)); + constexpr double MinRelDet = 1.e-12; + constexpr double InflateRelDet = 1.e-6; + constexpr double MaxInflation = 1.e4; + const bool isWellConditionedInfo = maxDiag > 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 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) = 4.; + info(1, 1) = 4.; + info(2, 2) = 4.; + return info; } //___________________________________________________________________ @@ -856,7 +943,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 +966,26 @@ 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; + /* + // 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); + if (!propagated) { + return false; + } + setTrackPos(mTrPos[mCurHyp][i], trc); } return true; } @@ -972,7 +1073,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 +1277,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]; } @@ -1245,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 } } @@ -1271,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 } } 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);