From 2cfb0d694ce1bffc05aa2af77f0a0963ab77f5fc Mon Sep 17 00:00:00 2001 From: Eric Cano Date: Thu, 8 Apr 2021 15:17:53 +0200 Subject: [PATCH 1/7] Reduced duplicate computations with temporary variables. (#614) Also removed addition of transposed matrix by directly copying the needed elements. https://github.com/cms-sw/cmssw/pull/31722#discussion_r603395679 https://github.com/cms-sw/cmssw/pull/31722#discussion_r603396371 --- .../PixelTrackFitting/interface/BrokenLine.h | 41 ++++++++++++------- 1 file changed, 26 insertions(+), 15 deletions(-) diff --git a/RecoPixelVertexing/PixelTrackFitting/interface/BrokenLine.h b/RecoPixelVertexing/PixelTrackFitting/interface/BrokenLine.h index 86fe6a278777c..76b4bd19fc2cd 100644 --- a/RecoPixelVertexing/PixelTrackFitting/interface/BrokenLine.h +++ b/RecoPixelVertexing/PixelTrackFitting/interface/BrokenLine.h @@ -208,31 +208,42 @@ namespace brokenline { const riemannFit::VectorNd& varBeta) { riemannFit::MatrixNd c_uMat = riemannFit::MatrixNd::Zero(); for (u_int i = 0; i < n; i++) { + // Pre compute differences used multiple time. Notation is t = sTotal(i + k) - sTotal(i + k - j). + // We pre-compute only when possible/needed + // The compiler with unroll this loop with a known boundary, and no tests should + // be actually executed. + const auto t10 = (i > 0) ? sTotal(i) - sTotal(i - 1) : 0; + const auto t21 = ((i > 0) && (i < (n - 1))) ? sTotal(i + 1) - sTotal(i - 1) : 0; + const auto t11 = (i < (n - 1)) ? sTotal(i + 1) - sTotal(i) : 0; + const auto t12 = (i < (n - 2)) ? sTotal(i + 2) - sTotal(i + 1) : 0; + const auto t22 = (i < (n - 2)) ? sTotal(i + 2) - sTotal(i) : 0; + const auto t11xt10 = t11 * t10; + const auto t12xt11 = t12 * t11; c_uMat(i, i) = weights(i); if (i > 1) - c_uMat(i, i) += 1. / (varBeta(i - 1) * riemannFit::sqr(sTotal(i) - sTotal(i - 1))); + c_uMat(i, i) += 1. / (varBeta(i - 1) * riemannFit::sqr(t10)); + if (i > 0 && i < n - 1) - c_uMat(i, i) += - (1. / varBeta(i)) * riemannFit::sqr((sTotal(i + 1) - sTotal(i - 1)) / - ((sTotal(i + 1) - sTotal(i)) * (sTotal(i) - sTotal(i - 1)))); + c_uMat(i, i) += riemannFit::sqr(t21 / t11xt10) / varBeta(i); + if (i < n - 2) - c_uMat(i, i) += 1. / (varBeta(i + 1) * riemannFit::sqr(sTotal(i + 1) - sTotal(i))); + c_uMat(i, i) += 1. / (varBeta(i + 1) * riemannFit::sqr(t11)); if (i > 0 && i < n - 1) - c_uMat(i, i + 1) = - 1. / (varBeta(i) * (sTotal(i + 1) - sTotal(i))) * - (-(sTotal(i + 1) - sTotal(i - 1)) / ((sTotal(i + 1) - sTotal(i)) * (sTotal(i) - sTotal(i - 1)))); - if (i < n - 2) - c_uMat(i, i + 1) += - 1. / (varBeta(i + 1) * (sTotal(i + 1) - sTotal(i))) * - (-(sTotal(i + 2) - sTotal(i)) / ((sTotal(i + 2) - sTotal(i + 1)) * (sTotal(i + 1) - sTotal(i)))); + c_uMat(i, i + 1) = -t21 / (varBeta(i) * t11 * t11xt10); if (i < n - 2) - c_uMat(i, i + 2) = 1. / (varBeta(i + 1) * (sTotal(i + 2) - sTotal(i + 1)) * (sTotal(i + 1) - sTotal(i))); + c_uMat(i, i + 1) += -t22 / (varBeta(i + 1) * t11 * t12xt11); - c_uMat(i, i) *= 0.5; + if (i < n - 1) // for 2 previous steps, non diagonal element copy for this symmetric matrix + c_uMat(i + 1, i) = c_uMat(i, i + 1); + + if (i < n - 2) { + c_uMat(i, i + 2) = 1. / (varBeta(i + 1) * t12xt11); + c_uMat(i + 2, i) = c_uMat(i, i + 2); // non-diagonal element copy. + } } - return c_uMat + c_uMat.transpose(); + return c_uMat; } /*! From 55be17f76f0cbf3f2bc5bdee8a10cacede2ec6f5 Mon Sep 17 00:00:00 2001 From: Eric Cano Date: Thu, 8 Apr 2021 16:38:29 +0200 Subject: [PATCH 2/7] Minor simplification (#614) https://github.com/cms-sw/cmssw/pull/31722#discussion_r603646172 --- RecoPixelVertexing/PixelTrackFitting/interface/RiemannFit.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/RecoPixelVertexing/PixelTrackFitting/interface/RiemannFit.h b/RecoPixelVertexing/PixelTrackFitting/interface/RiemannFit.h index 52cf4b637fb37..3405840be2d94 100644 --- a/RecoPixelVertexing/PixelTrackFitting/interface/RiemannFit.h +++ b/RecoPixelVertexing/PixelTrackFitting/interface/RiemannFit.h @@ -93,7 +93,7 @@ namespace riemannFit { for (uint k = 0; k < n; ++k) { for (uint l = k; l < n; ++l) { for (uint i = 0; i < std::min(k, l); ++i) { - tmp(k + n, l + n) += std::abs(s_values(k) - s_values(i)) * std::abs(s_values(l) - s_values(i)) * sig2_S(i); + tmp(k + n, l + n) += std::abs((s_values(k) - s_values(i)) * (s_values(l) - s_values(i))) * sig2_S(i); } tmp(l + n, k + n) = tmp(k + n, l + n); } From 55b7e36a25488a3b2bf1d0b4aa0396623708da2b Mon Sep 17 00:00:00 2001 From: Eric Cano Date: Fri, 9 Apr 2021 17:21:32 +0200 Subject: [PATCH 3/7] Reduced computations with a local variable. Also introduced a computation simplification suggested in: https://github.com/cms-sw/cmssw/pull/31722#discussion_r603648859 --- .../PixelTrackFitting/interface/RiemannFit.h | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/RecoPixelVertexing/PixelTrackFitting/interface/RiemannFit.h b/RecoPixelVertexing/PixelTrackFitting/interface/RiemannFit.h index 3405840be2d94..7e35295586569 100644 --- a/RecoPixelVertexing/PixelTrackFitting/interface/RiemannFit.h +++ b/RecoPixelVertexing/PixelTrackFitting/interface/RiemannFit.h @@ -125,10 +125,13 @@ namespace riemannFit { VectorNd const& rad, double B) { constexpr uint n = N; - double p_t = std::min(20., fast_fit(2) * B); // limit pt to avoid too small error!!! - double p_2 = p_t * p_t * (1. + 1. / sqr(fast_fit(3))); - double theta = atan(fast_fit(3)); - theta = theta < 0. ? theta + M_PI : theta; + const double p_t = std::min(20., fast_fit(2) * B); // limit pt to avoid too small error!!! + // fast_fit(3) = tan(theta) => + // 1 / sqr(sin(theta)) = (sqr(sin(theta) + sqr(cos(theta))) / sqr(sin(theta)) + // = 1 + 1 / sqr(tan(theta)) + // = 1 + 1 / sqr(fast_fit(3)) + const double invSqrSinTheta = 1. + 1. / sqr(fast_fit(3)); + const double p_2 = sqr(p_t) * invSqrSinTheta; VectorNd s_values; VectorNd rad_lengths; const Vector2d oVec(fast_fit(0), fast_fit(1)); @@ -141,10 +144,10 @@ namespace riemannFit { const double tempAtan2 = atan2(cross, dot); s_values(i) = std::abs(tempAtan2 * fast_fit(2)); } - computeRadLenUniformMaterial(s_values * sqrt(1. + 1. / sqr(fast_fit(3))), rad_lengths); + computeRadLenUniformMaterial(s_values * sqrt(invSqrSinTheta), rad_lengths); MatrixNd scatter_cov_rad = MatrixNd::Zero(); VectorNd sig2 = (1. + 0.038 * rad_lengths.array().log()).abs2() * rad_lengths.array(); - sig2 *= 0.000225 / (p_2 * sqr(sin(theta))); + sig2 *= 0.000225 / p_2 * invSqrSinTheta; for (uint k = 0; k < n; ++k) { for (uint l = k; l < n; ++l) { for (uint i = 0; i < std::min(k, l); ++i) { From 1141a92bbd59c7cb28881779d9d435a77e850549 Mon Sep 17 00:00:00 2001 From: Eric Cano Date: Fri, 9 Apr 2021 19:46:50 +0200 Subject: [PATCH 4/7] Simplified computations. https://github.com/cms-sw/cmssw/pull/31722#discussion_r603650560 https://github.com/cms-sw/cmssw/pull/31722#discussion_r603653076 https://github.com/cms-sw/cmssw/pull/31722#discussion_r603654860 --- .../PixelTrackFitting/interface/RiemannFit.h | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/RecoPixelVertexing/PixelTrackFitting/interface/RiemannFit.h b/RecoPixelVertexing/PixelTrackFitting/interface/RiemannFit.h index 7e35295586569..d435a157d8e11 100644 --- a/RecoPixelVertexing/PixelTrackFitting/interface/RiemannFit.h +++ b/RecoPixelVertexing/PixelTrackFitting/interface/RiemannFit.h @@ -257,7 +257,7 @@ namespace riemannFit { const double tan_c = -y2 / x2; const double tan_c2 = sqr(tan_c); cov_rad(i) = - 1. / (1. + tan_c2) * (cov_cart(i, i) + cov_cart(i + n, i + n) * tan_c2 + 2 * cov_cart(i, i + n) * tan_c); + (cov_cart(i, i) + cov_cart(i + n, i + n) * tan_c2 + 2 * cov_cart(i, i + n) * tan_c) / (1. + tan_c2) ; } } return cov_rad; @@ -517,7 +517,7 @@ namespace riemannFit { // scale const double tempQ = mc.squaredNorm(); - const double tempS = sqrt(n * 1. / tempQ); // scaling factor + const double tempS = sqrt(n / tempQ); // scaling factor p3D *= tempS; // project on paraboloid @@ -911,13 +911,14 @@ namespace riemannFit { // We need now to transfer back the results in the original s-z plane const auto sinTheta = sin(theta); const auto cosTheta = cos(theta); - auto common_factor = 1. / (sinTheta - sol(1, 0) * cosTheta); + const auto common_factor = 1. / (sinTheta - sol(1, 0) * cosTheta); + const auto sqr_common_factor = sqr(common_factor); Eigen::Matrix jMat; - jMat << 0., common_factor * common_factor, common_factor, sol(0, 0) * cosTheta * common_factor * common_factor; + jMat << 0., sqr_common_factor, common_factor, sol(0, 0) * cosTheta * sqr_common_factor; - double tempM = common_factor * (sol(1, 0) * sinTheta + cosTheta); - double tempQ = common_factor * sol(0, 0); - auto cov_mq = jMat * covParamsMat * jMat.transpose(); + const double tempM = common_factor * (sol(1, 0) * sinTheta + cosTheta); + const double tempQ = common_factor * sol(0, 0); + const auto cov_mq = jMat * covParamsMat * jMat.transpose(); VectorNd res = p2D_rot.row(1).transpose() - aMat.transpose() * sol; double chi2 = res.transpose() * vyInvMat * res; From f4a03f3d81719dcdad5efa54a974e8ce92c32e8b Mon Sep 17 00:00:00 2001 From: Eric Cano <37585813+ericcano@users.noreply.github.com> Date: Fri, 9 Apr 2021 17:34:10 +0200 Subject: [PATCH 5/7] Update RecoPixelVertexing/PixelTrackFitting/interface/RiemannFit.h Co-authored-by: Andrea Bocci --- RecoPixelVertexing/PixelTrackFitting/interface/RiemannFit.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/RecoPixelVertexing/PixelTrackFitting/interface/RiemannFit.h b/RecoPixelVertexing/PixelTrackFitting/interface/RiemannFit.h index d435a157d8e11..8968e646fa5fa 100644 --- a/RecoPixelVertexing/PixelTrackFitting/interface/RiemannFit.h +++ b/RecoPixelVertexing/PixelTrackFitting/interface/RiemannFit.h @@ -127,7 +127,7 @@ namespace riemannFit { constexpr uint n = N; const double p_t = std::min(20., fast_fit(2) * B); // limit pt to avoid too small error!!! // fast_fit(3) = tan(theta) => - // 1 / sqr(sin(theta)) = (sqr(sin(theta) + sqr(cos(theta))) / sqr(sin(theta)) + // 1 / sqr(sin(theta)) = (sqr(sin(theta)) + sqr(cos(theta))) / sqr(sin(theta)) // = 1 + 1 / sqr(tan(theta)) // = 1 + 1 / sqr(fast_fit(3)) const double invSqrSinTheta = 1. + 1. / sqr(fast_fit(3)); From 03c09bd511f005bab45526d35d7376c245eff503 Mon Sep 17 00:00:00 2001 From: Eric Cano Date: Tue, 13 Apr 2021 12:03:53 +0200 Subject: [PATCH 6/7] Removed unneeded range checking. https://github.com/cms-sw/cmssw/pull/31722#discussion_r603666335 https://github.com/cms-sw/cmssw/pull/31722#discussion_r603666574 --- RecoPixelVertexing/PixelTrackFitting/plugins/storeTracks.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/RecoPixelVertexing/PixelTrackFitting/plugins/storeTracks.h b/RecoPixelVertexing/PixelTrackFitting/plugins/storeTracks.h index 59101b6ba5214..09cf7eefd0197 100644 --- a/RecoPixelVertexing/PixelTrackFitting/plugins/storeTracks.h +++ b/RecoPixelVertexing/PixelTrackFitting/plugins/storeTracks.h @@ -47,7 +47,7 @@ void storeTracks(Ev& ev, const TWH& tracksWithHits, const TrackerTopology& ttopo reco::TrackExtra theTrackExtra{}; //fill the TrackExtra with TrackingRecHitRef - unsigned int nHits = tracks->at(k).numberOfValidHits(); + unsigned int nHits = (*tracks)[k].numberOfValidHits(); theTrackExtra.setHits(hitCollProd, cc, nHits); cc += nHits; AlgebraicVector5 v = AlgebraicVector5(0, 0, 0, 0, 0); @@ -63,7 +63,7 @@ void storeTracks(Ev& ev, const TWH& tracksWithHits, const TrackerTopology& ttopo for (int k = 0; k < nTracks; k++) { const reco::TrackExtraRef theTrackExtraRef(ohTE, k); - (tracks->at(k)).setExtra(theTrackExtraRef); + (*tracks)[k].setExtra(theTrackExtraRef); } ev.put(std::move(tracks)); From aee2dfb9276f04839813266811ef8cf52c4c5b6e Mon Sep 17 00:00:00 2001 From: Eric Cano Date: Thu, 27 May 2021 15:46:10 +0200 Subject: [PATCH 7/7] Fixed code format --- .../PixelTrackFitting/interface/BrokenLine.h | 14 +++++++------- .../PixelTrackFitting/interface/RiemannFit.h | 2 +- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/RecoPixelVertexing/PixelTrackFitting/interface/BrokenLine.h b/RecoPixelVertexing/PixelTrackFitting/interface/BrokenLine.h index 76b4bd19fc2cd..7bf9d5fa45e9e 100644 --- a/RecoPixelVertexing/PixelTrackFitting/interface/BrokenLine.h +++ b/RecoPixelVertexing/PixelTrackFitting/interface/BrokenLine.h @@ -212,11 +212,11 @@ namespace brokenline { // We pre-compute only when possible/needed // The compiler with unroll this loop with a known boundary, and no tests should // be actually executed. - const auto t10 = (i > 0) ? sTotal(i) - sTotal(i - 1) : 0; + const auto t10 = (i > 0) ? sTotal(i) - sTotal(i - 1) : 0; const auto t21 = ((i > 0) && (i < (n - 1))) ? sTotal(i + 1) - sTotal(i - 1) : 0; - const auto t11 = (i < (n - 1)) ? sTotal(i + 1) - sTotal(i) : 0; - const auto t12 = (i < (n - 2)) ? sTotal(i + 2) - sTotal(i + 1) : 0; - const auto t22 = (i < (n - 2)) ? sTotal(i + 2) - sTotal(i) : 0; + const auto t11 = (i < (n - 1)) ? sTotal(i + 1) - sTotal(i) : 0; + const auto t12 = (i < (n - 2)) ? sTotal(i + 2) - sTotal(i + 1) : 0; + const auto t22 = (i < (n - 2)) ? sTotal(i + 2) - sTotal(i) : 0; const auto t11xt10 = t11 * t10; const auto t12xt11 = t12 * t11; c_uMat(i, i) = weights(i); @@ -230,17 +230,17 @@ namespace brokenline { c_uMat(i, i) += 1. / (varBeta(i + 1) * riemannFit::sqr(t11)); if (i > 0 && i < n - 1) - c_uMat(i, i + 1) = -t21 / (varBeta(i) * t11 * t11xt10); + c_uMat(i, i + 1) = -t21 / (varBeta(i) * t11 * t11xt10); if (i < n - 2) c_uMat(i, i + 1) += -t22 / (varBeta(i + 1) * t11 * t12xt11); - if (i < n - 1) // for 2 previous steps, non diagonal element copy for this symmetric matrix + if (i < n - 1) // for 2 previous steps, non diagonal element copy for this symmetric matrix c_uMat(i + 1, i) = c_uMat(i, i + 1); if (i < n - 2) { c_uMat(i, i + 2) = 1. / (varBeta(i + 1) * t12xt11); - c_uMat(i + 2, i) = c_uMat(i, i + 2); // non-diagonal element copy. + c_uMat(i + 2, i) = c_uMat(i, i + 2); // non-diagonal element copy. } } return c_uMat; diff --git a/RecoPixelVertexing/PixelTrackFitting/interface/RiemannFit.h b/RecoPixelVertexing/PixelTrackFitting/interface/RiemannFit.h index 8968e646fa5fa..63f139c70463c 100644 --- a/RecoPixelVertexing/PixelTrackFitting/interface/RiemannFit.h +++ b/RecoPixelVertexing/PixelTrackFitting/interface/RiemannFit.h @@ -257,7 +257,7 @@ namespace riemannFit { const double tan_c = -y2 / x2; const double tan_c2 = sqr(tan_c); cov_rad(i) = - (cov_cart(i, i) + cov_cart(i + n, i + n) * tan_c2 + 2 * cov_cart(i, i + n) * tan_c) / (1. + tan_c2) ; + (cov_cart(i, i) + cov_cart(i + n, i + n) * tan_c2 + 2 * cov_cart(i, i + n) * tan_c) / (1. + tan_c2); } } return cov_rad;