From 4758e53daed59396c206bca0173e81ef7fa91f90 Mon Sep 17 00:00:00 2001 From: shahoian Date: Thu, 27 Nov 2025 02:01:31 +0100 Subject: [PATCH] Methods for Kalman filter linearized wrt reference track --- .../TrackParametrization.h | 1 + .../TrackParametrizationWithError.h | 8 +- .../src/TrackParametrization.cxx | 32 ++ .../src/TrackParametrizationWithError.cxx | 373 +++++++++++++++++- .../Base/include/DetectorsBase/Propagator.h | 4 + Detectors/Base/src/Propagator.cxx | 69 +++- 6 files changed, 461 insertions(+), 26 deletions(-) diff --git a/DataFormats/Reconstruction/include/ReconstructionDataFormats/TrackParametrization.h b/DataFormats/Reconstruction/include/ReconstructionDataFormats/TrackParametrization.h index 8cb22efd39e38..6d117efd37a44 100644 --- a/DataFormats/Reconstruction/include/ReconstructionDataFormats/TrackParametrization.h +++ b/DataFormats/Reconstruction/include/ReconstructionDataFormats/TrackParametrization.h @@ -229,6 +229,7 @@ class TrackParametrization // parameters manipulation GPUd() bool correctForELoss(value_t xrho, bool anglecorr = false); GPUd() bool rotateParam(value_t alpha); + GPUd() bool rotateParamCS(value_t& alpha, value_t& ca, value_t& sa); GPUd() bool propagateParamTo(value_t xk, value_t b); GPUd() bool propagateParamTo(value_t xk, const dim3_t& b); GPUd() void invertParam(); diff --git a/DataFormats/Reconstruction/include/ReconstructionDataFormats/TrackParametrizationWithError.h b/DataFormats/Reconstruction/include/ReconstructionDataFormats/TrackParametrizationWithError.h index cd9d1517a81b1..acfc217cb1dfe 100644 --- a/DataFormats/Reconstruction/include/ReconstructionDataFormats/TrackParametrizationWithError.h +++ b/DataFormats/Reconstruction/include/ReconstructionDataFormats/TrackParametrizationWithError.h @@ -89,9 +89,11 @@ class TrackParametrizationWithError : public TrackParametrization // parameters + covmat manipulation GPUd() bool testRotate(value_t alpha) const; GPUd() bool rotate(value_t alpha); - GPUd() bool propagateTo(value_t xk, value_t b); + GPUd() bool rotateWithRef(value_t alpha, TrackParametrization& linRef, value_t bz); + GPUd() bool propagateTo(value_t xk, value_t bz); + GPUd() bool propagateToWithRef(value_t xk, TrackParametrization& linRef0, value_t bz); GPUd() bool propagateTo(value_t xk, const dim3_t& b); - GPUd() bool propagateToDCA(const o2::dataformats::VertexBase& vtx, value_t b, o2::dataformats::DCA* dca = nullptr, value_t maxD = 999.f); + GPUd() bool propagateToDCA(const o2::dataformats::VertexBase& vtx, value_t bz, o2::dataformats::DCA* dca = nullptr, value_t maxD = 999.f); GPUd() void invert(); GPUd() value_t getPredictedChi2(const dim2_t& p, const dim3_t& cov) const; GPUd() value_t getPredictedChi2Quiet(const dim2_t& p, const dim3_t& cov) const; @@ -118,7 +120,7 @@ class TrackParametrizationWithError : public TrackParametrization GPUd() bool update(const BaseCluster& p); GPUd() bool correctForMaterial(value_t x2x0, value_t xrho, bool anglecorr = false); - + GPUd() bool correctForMaterial(TrackParametrization& linRef, value_t x2x0, value_t xrho, bool anglecorr = false); GPUd() void resetCovariance(value_t s2 = 0); GPUd() void checkCovariance(); GPUd() void checkCorrelations(); diff --git a/DataFormats/Reconstruction/src/TrackParametrization.cxx b/DataFormats/Reconstruction/src/TrackParametrization.cxx index 7086e4d93cec8..74ff7a476c5a4 100644 --- a/DataFormats/Reconstruction/src/TrackParametrization.cxx +++ b/DataFormats/Reconstruction/src/TrackParametrization.cxx @@ -188,6 +188,38 @@ GPUd() bool TrackParametrization::rotateParam(value_t alpha) return true; } +//______________________________________________________________ +template +GPUd() bool TrackParametrization::rotateParamCS(value_t& alpha, value_t& ca, value_t& sa) +{ + // rotate to alpha frame + if (gpu::CAMath::Abs(getSnp()) > constants::math::Almost1) { + LOGP(debug, "Precondition is not satisfied: |sin(phi)|>1 ! {:f}", getSnp()); + return false; + } + // + math_utils::detail::bringToPMPi(alpha); + math_utils::detail::sincos(alpha - getAlpha(), sa, ca); + value_t snp = getSnp(), csp = gpu::CAMath::Sqrt((1.f - snp) * (1.f + snp)); // Improve precision + // RS: check if rotation does no invalidate track model (cos(local_phi)>=0, i.e. particle direction in local frame is along the X axis + if ((csp * ca + snp * sa) < 0) { + // LOGF(warning,"Rotation failed: local cos(phi) would become {:.2f}", csp * ca + snp * sa); + return false; + } + // + value_t tmp = snp * ca - csp * sa; + if (gpu::CAMath::Abs(tmp) > constants::math::Almost1) { + LOGP(debug, "Rotation failed: new snp {:.2f}", tmp); + return false; + } + value_t xold = getX(), yold = getY(); + mAlpha = alpha; + mX = xold * ca + yold * sa; + mP[kY] = -xold * sa + yold * ca; + mP[kSnp] = tmp; + return true; +} + //____________________________________________________________ template GPUd() bool TrackParametrization::propagateParamTo(value_t xk, const dim3_t& b) diff --git a/DataFormats/Reconstruction/src/TrackParametrizationWithError.cxx b/DataFormats/Reconstruction/src/TrackParametrizationWithError.cxx index 01849bd0c9e8f..a0b9487daa993 100644 --- a/DataFormats/Reconstruction/src/TrackParametrizationWithError.cxx +++ b/DataFormats/Reconstruction/src/TrackParametrizationWithError.cxx @@ -43,7 +43,7 @@ GPUd() void TrackParametrizationWithError::invert() //______________________________________________________________ template -GPUd() bool TrackParametrizationWithError::propagateTo(value_t xk, value_t b) +GPUd() bool TrackParametrizationWithError::propagateTo(value_t xk, value_t bz) { //---------------------------------------------------------------- // propagate this track to the plane X=xk (cm) in the field "b" (kG) @@ -52,7 +52,7 @@ GPUd() bool TrackParametrizationWithError::propagateTo(value_t xk, valu if (gpu::CAMath::Abs(dx) < constants::math::Almost0) { return true; } - value_t crv = this->getCurvature(b); + value_t crv = this->getCurvature(bz); value_t x2r = crv * dx; value_t f1 = this->getSnp(), f2 = f1 + x2r; if ((gpu::CAMath::Abs(f1) > constants::math::Almost1) || (gpu::CAMath::Abs(f2) > constants::math::Almost1)) { @@ -66,7 +66,8 @@ GPUd() bool TrackParametrizationWithError::propagateTo(value_t xk, valu if (gpu::CAMath::Abs(r2) < constants::math::Almost0) { return false; } - double dy2dx = (f1 + f2) / (r1 + r2); + double r1pr2Inv = 1. / (r1 + r2); + double dy2dx = (f1 + f2) * r1pr2Inv; bool arcz = gpu::CAMath::Abs(x2r) > 0.05f; params_t dP{0.f}; if (arcz) { @@ -106,14 +107,110 @@ GPUd() bool TrackParametrizationWithError::propagateTo(value_t xk, valu &c44 = mC[kSigQ2Pt2]; // evaluate matrix in double prec. - double rinv = 1. / r1; - double r3inv = rinv * rinv * rinv; - double f24 = dx * b * constants::math::B2C; // x2r/mP[kQ2Pt]; - double f02 = dx * r3inv; - double f04 = 0.5 * f24 * f02; - double f12 = f02 * this->getTgl() * f1; - double f14 = 0.5 * f24 * f12; // 0.5*f24*f02*getTgl()*f1; - double f13 = dx * rinv; + value_t kb = bz * constants::math::B2C; + double r2inv = 1. / r2, r1inv = 1. / r1; + double dx2r1pr2 = dx * r1pr2Inv; + + double hh = dx2r1pr2 * r2inv * (1. + r1 * r2 + f1 * f2), jj = dx * (dy2dx - f2 * r2inv); + double f02 = hh * r1inv; + double f04 = hh * dx2r1pr2 * kb; + double f24 = dx * kb; // x2r/mP[kQ2Pt]; + double f12 = this->getTgl() * (f02 * f2 + jj); + double f13 = dx * (r2 + f2 * dy2dx); + double f14 = this->getTgl() * (f04 * f2 + jj * f24); + + // b = C*ft + double b00 = f02 * c20 + f04 * c40, b01 = f12 * c20 + f14 * c40 + f13 * c30; + double b02 = f24 * c40; + double b10 = f02 * c21 + f04 * c41, b11 = f12 * c21 + f14 * c41 + f13 * c31; + double b12 = f24 * c41; + double b20 = f02 * c22 + f04 * c42, b21 = f12 * c22 + f14 * c42 + f13 * c32; + double b22 = f24 * c42; + double b40 = f02 * c42 + f04 * c44, b41 = f12 * c42 + f14 * c44 + f13 * c43; + double b42 = f24 * c44; + double b30 = f02 * c32 + f04 * c43, b31 = f12 * c32 + f14 * c43 + f13 * c33; + double b32 = f24 * c43; + + // a = f*b = f*C*ft + double a00 = f02 * b20 + f04 * b40, a01 = f02 * b21 + f04 * b41, a02 = f02 * b22 + f04 * b42; + double a11 = f12 * b21 + f14 * b41 + f13 * b31, a12 = f12 * b22 + f14 * b42 + f13 * b32; + double a22 = f24 * b42; + + // F*C*Ft = C + (b + bt + a) + c00 += b00 + b00 + a00; + c10 += b10 + b01 + a01; + c20 += b20 + b02 + a02; + c30 += b30; + c40 += b40; + c11 += b11 + b11 + a11; + c21 += b21 + b12 + a12; + c31 += b31; + c41 += b41; + c22 += b22 + b22 + a22; + c32 += b32; + c42 += b42; + + checkCovariance(); + + return true; +} + +//______________________________________________________________ +template +GPUd() bool TrackParametrizationWithError::propagateToWithRef(value_t xk, TrackParametrization& linRef0, value_t bz) +{ + //---------------------------------------------------------------- + // propagate this track to the plane X=xk (cm) in the field "b" (kG), using linRef as linearization point + //---------------------------------------------------------------- + if (this->getAbsCharge() == 0) { + bz = 0; + } + value_t dx = xk - this->getX(); + if (gpu::CAMath::Abs(dx) < constants::math::Almost0) { + this->setX(xk); + linRef0.setX(xk); + return true; + } + // propagate reference track + TrackParametrization linRef1 = linRef0; + if (!linRef1.propagateTo(xk, bz)) { + return false; + } + value_t kb = bz * constants::math::B2C; + // evaluate in double prec. + double snpRef0 = linRef0.getSnp(), cspRef0 = gpu::CAMath::Sqrt((1 - snpRef0) * (1 + snpRef0)); + double snpRef1 = linRef1.getSnp(), cspRef1 = gpu::CAMath::Sqrt((1 - snpRef1) * (1 + snpRef1)); + double cspRef0Inv = 1 / cspRef0, cspRef1Inv = 1 / cspRef1, cc = cspRef0 + cspRef1, ccInv = 1 / cc, dy2dx = (snpRef0 + snpRef1) * ccInv; + double dxccInv = dx * ccInv, hh = dxccInv * cspRef1Inv * (1 + cspRef0 * cspRef1 + snpRef0 * snpRef1), jj = dx * (dy2dx - snpRef1 * cspRef1Inv); + + double f02 = hh * cspRef0Inv; + double f04 = hh * dxccInv * kb; + double f24 = dx * kb; + double f12 = linRef0.getTgl() * (f02 * snpRef1 + jj); + double f13 = dx * (cspRef1 + snpRef1 * dy2dx); // dS + double f14 = linRef0.getTgl() * (f04 * snpRef1 + jj * f24); + + // difference between the current and reference state + value_t diff[5]; + for (int i = 0; i < 5; i++) { + diff[i] = this->getParam(i) - linRef0.getParam(i); + } + value_t snpUpd = snpRef1 + diff[kSnp] + f24 * diff[kQ2Pt]; + if (gpu::CAMath::Abs(snpUpd) > constants::math::Almost1) { + return false; + } + linRef0 = linRef1; // update reference track + this->setX(xk); + this->setY(linRef1.getY() + diff[kY] + f02 * diff[kSnp] + f04 * diff[kQ2Pt]); + this->setZ(linRef1.getZ() + diff[kZ] + f13 * diff[kTgl] + f14 * diff[kQ2Pt]); + this->setSnp(snpUpd); + this->setTgl(linRef1.getTgl() + diff[kTgl]); + this->setQ2Pt(linRef1.getQ2Pt() + diff[kQ2Pt]); + + value_t &c00 = mC[kSigY2], &c10 = mC[kSigZY], &c11 = mC[kSigZ2], &c20 = mC[kSigSnpY], &c21 = mC[kSigSnpZ], + &c22 = mC[kSigSnp2], &c30 = mC[kSigTglY], &c31 = mC[kSigTglZ], &c32 = mC[kSigTglSnp], &c33 = mC[kSigTgl2], + &c40 = mC[kSigQ2PtY], &c41 = mC[kSigQ2PtZ], &c42 = mC[kSigQ2PtSnp], &c43 = mC[kSigQ2PtTgl], + &c44 = mC[kSigQ2Pt2]; // b = C*ft double b00 = f02 * c20 + f04 * c40, b01 = f12 * c20 + f14 * c40 + f13 * c30; @@ -158,6 +255,7 @@ GPUd() bool TrackParametrizationWithError::testRotate(value_t) const // no ops return true; } + //______________________________________________________________ template GPUd() bool TrackParametrizationWithError::rotate(value_t alpha) @@ -213,6 +311,101 @@ GPUd() bool TrackParametrizationWithError::rotate(value_t alpha) return true; } +//______________________________________________________________ +template +GPUd() bool TrackParametrizationWithError::rotateWithRef(value_t alpha, TrackParametrization& linRef0, value_t bz) +{ + // RS: similar to int32_t GPUTPCGMPropagator::RotateToAlpha(float newAlpha), i.e. rotate the track to new frame alpha, using linRef as linearization point + // rotate to alpha frame the reference (linearization point) trackParam, then align the current track to it + if (gpu::CAMath::Abs(this->getSnp()) > constants::math::Almost1) { + LOGP(debug, "Precondition is not satisfied: |sin(phi)|>1 ! {:f}", this->getSnp()); + return false; + } + // + math_utils::detail::bringToPMPi(alpha); + // + value_t ca = 0, sa = 0; + TrackParametrization linRef1 = linRef0; + // rotate the reference, adjusting alpha to +-pi, return precalculated cos and sin of alpha - alphaOld + if (!linRef1.rotateParamCS(alpha, ca, sa)) { + return false; + } + + value_t trackX = this->getX() * ca + this->getY() * sa; // X of the rotated current track + if (!linRef1.propagateParamTo(trackX, bz)) { + return false; + } + + // now rotate the current track + value_t snp = this->getSnp(), csp = gpu::CAMath::Sqrt((1.f - snp) * (1.f + snp)), updSnp = snp * ca - csp * sa; + if ((csp * ca + snp * sa) < 0 || gpu::CAMath::Abs(updSnp) > constants::math::Almost1) { + // LOGP(warning,"Rotation failed: local cos(phi) would become {:.2f}", csp * ca + snp * sa); + return false; + } + this->setY(-sa * this->getX() + ca * this->getY()); + this->setX(trackX); + this->setSnp(updSnp); + this->setAlpha(alpha); + + // rotate covariance, accounting for the extra error from the rotated X + value_t snpRef0 = linRef0.getSnp(), cspRef0 = gpu::CAMath::Sqrt((value_t(1) - snpRef0) * (value_t(1) + snpRef0)); // original reference + value_t snpRef1 = linRef1.getSnp(), cspRef1 = ca * cspRef0 + sa * snpRef0; // rotated reference + value_t rr = cspRef1 / cspRef0; // cos1_ref / cos0_ref + + // "extra row" of the lower triangle of cov. matrix + value_t cXSigY = mC[kSigY2] * ca * sa; + value_t cXSigZ = mC[kSigZY] * sa; + value_t cXSigSnp = mC[kSigSnpY] * rr * sa; + value_t cXSigTgl = mC[kSigTglY] * sa; + value_t cXSigQ2Pt = mC[kSigQ2PtY] * sa; + value_t cSigX2 = mC[kSigY2] * sa * sa; + + // plane rotation of existing cov matrix + mC[kSigY2] *= ca * ca; + mC[kSigZY] *= ca; + mC[kSigSnpY] *= ca * rr; + mC[kSigSnpZ] *= rr; + mC[kSigSnp2] *= rr * rr; + mC[kSigTglY] *= ca; + mC[kSigTglSnp] *= rr; + mC[kSigQ2PtY] *= ca; + mC[kSigQ2PtSnp] *= rr; + + // transport covariance from pseudo 6x6 matrix to usual 5x5, Jacobian (trust to Sergey): + auto cspRef1Inv = value_t(1) / cspRef1; + auto j3 = -snpRef1 * cspRef1Inv; // -pYmod/pXmod = -tg_pho = -sin_phi_mod / cos_phi_mod + auto j4 = -linRef1.getTgl() * cspRef1Inv; // -pZmod/pXmod = -tgl_mod / cos_phi_mod + auto j5 = linRef1.getCurvature(bz); + // Y Z Sin DzDs q/p X + // { { 1, 0, 0, 0, 0, j3 }, // Y + // { 0, 1, 0, 0, 0, j4 }, // Z + // { 0, 0, 1, 0, 0, j5 }, // snp + // { 0, 0, 0, 1, 0, 0 }, // tgl + // { 0, 0, 0, 0, 1, 0 } }; // q/pt + auto hXSigY = cXSigY + cSigX2 * j3; + auto hXSigZ = cXSigZ + cSigX2 * j4; + auto hXSigSnp = cXSigSnp + cSigX2 * j5; + + mC[kSigY2] += j3 * (cXSigY + hXSigY); + mC[kSigZ2] += j4 * (cXSigZ + hXSigZ); + mC[kSigSnpY] += cXSigSnp * j3 + hXSigY * j5; + mC[kSigSnp2] += j5 * (cXSigSnp + hXSigSnp); + mC[kSigTglZ] += cXSigTgl * j4; + mC[kSigQ2PtY] += cXSigQ2Pt * j3; + mC[kSigQ2PtSnp] += cXSigQ2Pt * j5; + + mC[kSigZY] += cXSigZ * j3 + hXSigY * j4; + mC[kSigSnpZ] += cXSigSnp * j4 + hXSigZ * j5; + mC[kSigTglY] += cXSigTgl * j3; + mC[kSigTglSnp] += cXSigTgl * j5; + mC[kSigQ2PtZ] += cXSigQ2Pt * j4; + + checkCovariance(); + linRef0 = linRef1; + + return true; +} + //_______________________________________________________________________ template GPUd() bool TrackParametrizationWithError::propagateToDCA(const o2::dataformats::VertexBase& vtx, value_t b, o2::dataformats::DCA* dca, value_t maxD) @@ -476,8 +669,8 @@ GPUd() bool TrackParametrizationWithError::propagateTo(value_t xk, cons if (gpu::CAMath::Abs(r2) < constants::math::Almost0) { return false; } - - value_t dy2dx = (f1 + f2) / (r1 + r2); + double r1pr2Inv = 1. / (r1 + r2), r2inv = 1. / r2, r1inv = 1. / r1; + double dy2dx = (f1 + f2) * r1pr2Inv, dx2r1pr2 = dx * r1pr2Inv; value_t step = (gpu::CAMath::Abs(x2r) < 0.05f) ? dx * gpu::CAMath::Abs(r2 + f2 * dy2dx) // chord : 2.f * gpu::CAMath::ASin(0.5f * dx * gpu::CAMath::Sqrt(1.f + dy2dx * dy2dx) * crv) / crv; // arc step *= gpu::CAMath::Sqrt(1.f + this->getTgl() * this->getTgl()); @@ -493,15 +686,16 @@ GPUd() bool TrackParametrizationWithError::propagateTo(value_t xk, cons &c22 = mC[kSigSnp2], &c30 = mC[kSigTglY], &c31 = mC[kSigTglZ], &c32 = mC[kSigTglSnp], &c33 = mC[kSigTgl2], &c40 = mC[kSigQ2PtY], &c41 = mC[kSigQ2PtZ], &c42 = mC[kSigQ2PtSnp], &c43 = mC[kSigQ2PtTgl], &c44 = mC[kSigQ2Pt2]; + // evaluate matrix in double prec. - double rinv = 1. / r1; - double r3inv = rinv * rinv * rinv; - double f24 = dx * b[2] * constants::math::B2C; // x2r/track[kQ2Pt]; - double f02 = dx * r3inv; - double f04 = 0.5 * f24 * f02; - double f12 = f02 * this->getTgl() * f1; - double f14 = 0.5 * f24 * f12; // 0.5*f24*f02*getTgl()*f1; - double f13 = dx * rinv; + value_t kb = b[2] * constants::math::B2C; + double hh = dx2r1pr2 * r2inv * (1. + r1 * r2 + f1 * f2), jj = dx * (dy2dx - f2 * r2inv); + double f02 = hh * r1inv; + double f04 = hh * dx2r1pr2 * kb; + double f24 = dx * kb; // x2r/mP[kQ2Pt]; + double f12 = this->getTgl() * (f02 * f2 + jj); + double f13 = dx * (r2 + f2 * dy2dx); + double f14 = this->getTgl() * (f04 * f2 + jj * f24); // b = C*ft double b00 = f02 * c20 + f04 * c40, b01 = f12 * c20 + f14 * c40 + f13 * c30; @@ -1121,6 +1315,143 @@ GPUd() bool TrackParametrizationWithError::correctForMaterial(value_t x return true; } +//______________________________________________ +template +GPUd() bool TrackParametrizationWithError::correctForMaterial(TrackParametrization& linRef, value_t x2x0, value_t xrho, bool anglecorr) +{ + //------------------------------------------------------------------ + // This function corrects the reference and current track parameters for the crossed material + // "x2x0" - X/X0, the thickness in units of the radiation length. + // "xrho" - is the product length*density (g/cm^2). + // It should be passed as negative when propagating tracks + // from the intreaction point to the outside of the central barrel. + // "dedx" - mean enery loss (GeV/(g/cm^2), if <=kCalcdEdxAuto : calculate on the fly + // "anglecorr" - switch for the angular correction + //------------------------------------------------------------------ + constexpr value_t kMSConst2 = 0.0136f * 0.0136f; + constexpr value_t kMinP = 0.01f; // kill below this momentum + + value_t csp2 = (1.f - linRef.getSnp()) * (1.f + linRef.getSnp()); // cos(phi)^2 + value_t cst2I = (1.f + linRef.getTgl() * linRef.getTgl()); // 1/cos(lambda)^2 + if (anglecorr) { // Apply angle correction, if requested + value_t angle = gpu::CAMath::Sqrt(cst2I / csp2); + x2x0 *= angle; + xrho *= angle; + } + auto pid = linRef.getPID(); + auto m = pid.getMass(); + int charge2 = linRef.getAbsCharge() * linRef.getAbsCharge(); + value_t p = linRef.getP(), p0 = p, p02 = p * p, e2 = p02 + pid.getMass2(), massInv = 1. / m, bg = p * massInv, dETot = 0.; + value_t e = gpu::CAMath::Sqrt(e2), e0 = e; + if (m > 0 && xrho != 0.f) { + value_t ekin = e - m, dedx = this->getdEdxBBOpt(bg); +#ifdef _BB_NONCONST_CORR_ + value_t dedxDer = 0., dedx1 = dedx; +#endif + if (charge2 != 1) { + dedx *= charge2; + } + value_t dE = dedx * xrho; + int na = 1 + int(gpu::CAMath::Abs(dE) / ekin * ELoss2EKinThreshInv); + if (na > MaxELossIter) { + na = MaxELossIter; + } + if (na > 1) { + dE /= na; + xrho /= na; +#ifdef _BB_NONCONST_CORR_ + dedxDer = this->getBetheBlochSolidDerivativeApprox(dedx1, bg); // require correction for non-constantness of dedx vs betagamma + if (charge2 != 1) { + dedxDer *= charge2; + } +#endif + } + while (na--) { +#ifdef _BB_NONCONST_CORR_ + if (dedxDer != 0.) { // correction for non-constantness of dedx vs beta*gamma (in linear approximation): for a single step dE -> dE * [(exp(dedxDer) - 1)/dedxDer] + if (xrho < 0) { + dedxDer = -dedxDer; // E.loss ( -> positive derivative) + } + auto corrC = (gpu::CAMath::Exp(dedxDer) - 1.) / dedxDer; + dE *= corrC; + } +#endif + e += dE; + if (e > m) { // stopped + p = gpu::CAMath::Sqrt(e * e - pid.getMass2()); + } else { + return false; + } + if (na) { + bg = p * massInv; + dedx = this->getdEdxBBOpt(bg); +#ifdef _BB_NONCONST_CORR_ + dedxDer = this->getBetheBlochSolidDerivativeApprox(dedx, bg); +#endif + if (charge2 != 1) { + dedx *= charge2; +#ifdef _BB_NONCONST_CORR_ + dedxDer *= charge2; +#endif + } + dE = dedx * xrho; + } + } + + if (p < kMinP) { + return false; + } + dETot = e - e0; + } // end of e.loss correction + + // Calculating the multiple scattering corrections****************** + value_t& fC22 = mC[kSigSnp2]; + value_t& fC33 = mC[kSigTgl2]; + value_t& fC43 = mC[kSigQ2PtTgl]; + value_t& fC44 = mC[kSigQ2Pt2]; + // + value_t cC22(0.f), cC33(0.f), cC43(0.f), cC44(0.f); + if (x2x0 != 0.f) { + value_t beta2 = p02 / e2, theta2 = kMSConst2 / (beta2 * p02) * gpu::CAMath::Abs(x2x0); + value_t fp34 = linRef.getTgl(); + if (charge2 != 1) { + theta2 *= charge2; + fp34 *= linRef.getCharge2Pt(); + } + if (theta2 > constants::math::PI * constants::math::PI) { + return false; + } + value_t t2c2I = theta2 * cst2I; + cC22 = t2c2I * csp2; + cC33 = t2c2I * cst2I; + cC43 = t2c2I * fp34; + cC44 = theta2 * fp34 * fp34; + // optimize this + // cC22 = theta2*((1.-getSnp())*(1.+getSnp()))*(1. + this->getTgl()*getTgl()); + // cC33 = theta2*(1. + this->getTgl()*getTgl())*(1. + this->getTgl()*getTgl()); + // cC43 = theta2*getTgl()*this->getQ2Pt()*(1. + this->getTgl()*getTgl()); + // cC44 = theta2*getTgl()*this->getQ2Pt()*getTgl()*this->getQ2Pt(); + } + + // the energy loss correction contribution to cov.matrix: approximate energy loss fluctuation (M.Ivanov) + constexpr value_t knst = 0.0007f; // To be tuned. + value_t sigmadE = knst * gpu::CAMath::Sqrt(gpu::CAMath::Abs(dETot)) * e0 / p02 * linRef.getCharge2Pt(); + cC44 += sigmadE * sigmadE; + + // Applying the corrections***************************** + fC22 += cC22; + fC33 += cC33; + fC43 += cC43; + fC44 += cC44; + auto pscale = p0 / p; + linRef.setQ2Pt(linRef.getQ2Pt() * pscale); + this->setQ2Pt(this->getQ2Pt() * pscale); + + checkCovariance(); + + return true; +} + //______________________________________________________________ template GPUd() bool TrackParametrizationWithError::getCovXYZPxPyPzGlo(std::array& cv) const diff --git a/Detectors/Base/include/DetectorsBase/Propagator.h b/Detectors/Base/include/DetectorsBase/Propagator.h index d9b1522f4295b..b41fe21e3dd53 100644 --- a/Detectors/Base/include/DetectorsBase/Propagator.h +++ b/Detectors/Base/include/DetectorsBase/Propagator.h @@ -84,6 +84,10 @@ class PropagatorImpl value_type maxSnp = MAX_SIN_PHI, value_type maxStep = MAX_STEP, MatCorrType matCorr = MatCorrType::USEMatCorrLUT, track::TrackLTIntegral* tofInfo = nullptr, int signCorr = 0) const; + GPUd() bool propagateToX(TrackParCov_t& track, TrackPar_t& linRef, value_type x, value_type bZ, + value_type maxSnp = MAX_SIN_PHI, value_type maxStep = MAX_STEP, MatCorrType matCorr = MatCorrType::USEMatCorrLUT, + track::TrackLTIntegral* tofInfo = nullptr, int signCorr = 0) const; + GPUd() bool propagateToX(TrackPar_t& track, value_type x, value_type bZ, value_type maxSnp = MAX_SIN_PHI, value_type maxStep = MAX_STEP, MatCorrType matCorr = MatCorrType::USEMatCorrLUT, track::TrackLTIntegral* tofInfo = nullptr, int signCorr = 0) const; diff --git a/Detectors/Base/src/Propagator.cxx b/Detectors/Base/src/Propagator.cxx index 0763eb48ff474..494019d5b6a00 100644 --- a/Detectors/Base/src/Propagator.cxx +++ b/Detectors/Base/src/Propagator.cxx @@ -295,8 +295,7 @@ GPUd() bool PropagatorImpl::propagateToX(TrackParCov_t& track, value_ty //---------------------------------------------------------------- // // Propagates the track to the plane X=xk (cm) - // taking into account all the three components of the magnetic field - // and correcting for the crossed material. + // Use bz only and correct for the crossed material. // // maxStep - maximal step for propagation // tofInfo - optional container for track length and PID-dependent TOF integration @@ -352,6 +351,72 @@ GPUd() bool PropagatorImpl::propagateToX(TrackParCov_t& track, value_ty return true; } +//_______________________________________________________________________ +template +GPUd() bool PropagatorImpl::propagateToX(TrackParCov_t& track, TrackPar_t& linRef, value_type xToGo, value_type bZ, value_type maxSnp, value_type maxStep, + PropagatorImpl::MatCorrType matCorr, track::TrackLTIntegral* tofInfo, int signCorr) const +{ + //---------------------------------------------------------------- + // + // Propagates the track to the plane X=xk (cm), using linRef as a Kalman linearisation point. + // Use bz only and correct for the crossed material if requested. + // + // maxStep - maximal step for propagation + // tofInfo - optional container for track length and PID-dependent TOF integration + // + // matCorr - material correction type, it is up to the user to make sure the pointer is attached (if LUT is requested) + //---------------------------------------------------------------- + auto dx = xToGo - track.getX(); + int dir = dx > 0.f ? 1 : -1; + if (!signCorr) { + signCorr = -dir; // sign of eloss correction is not imposed + } + + while (math_utils::detail::abs(dx) > Epsilon) { + auto step = math_utils::detail::min(math_utils::detail::abs(dx), maxStep); + if (dir < 0) { + step = -step; + } + auto x = track.getX() + step; + auto xyz0 = linRef.getXYZGlo(); + + auto correct = [&track, &linRef, &xyz0, tofInfo, matCorr, signCorr, this]() { + bool res = true; + if (matCorr != MatCorrType::USEMatCorrNONE) { + auto xyz1 = linRef.getXYZGlo(); + auto mb = this->getMatBudget(matCorr, xyz0, xyz1); + if (!track.correctForMaterial(linRef, mb.meanX2X0, mb.getXRho(signCorr))) { + res = false; + } + if (tofInfo) { + tofInfo->addStep(mb.length, linRef.getQ2P2()); // fill L,ToF info using already calculated step length + tofInfo->addX2X0(mb.meanX2X0); + tofInfo->addXRho(mb.getXRho(signCorr)); + } + } else if (tofInfo) { // if tofInfo filling was requested w/o material correction, we need to calculate the step lenght + auto xyz1 = linRef.getXYZGlo(); + math_utils::Vector3D stepV(xyz1.X() - xyz0.X(), xyz1.Y() - xyz0.Y(), xyz1.Z() - xyz0.Z()); + tofInfo->addStep(stepV.R(), linRef.getQ2P2()); + } + return res; + }; + + if (!track.propagateToWithRef(x, linRef, bZ)) { // linRef also updated + return false; + } + if (maxSnp > 0 && math_utils::detail::abs(track.getSnp()) >= maxSnp) { + correct(); + return false; + } + if (!correct()) { + return false; + } + dx = xToGo - track.getX(); + } + track.setX(xToGo); + return true; +} + //_______________________________________________________________________ template GPUd() bool PropagatorImpl::propagateToX(TrackPar_t& track, value_type xToGo, value_type bZ, value_type maxSnp, value_type maxStep,