diff --git a/RecoPixelVertexing/PixelTrackFitting/interface/BrokenLine.h b/RecoPixelVertexing/PixelTrackFitting/interface/BrokenLine.h index 86fe6a278777c..7bf9d5fa45e9e 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; } /*! diff --git a/RecoPixelVertexing/PixelTrackFitting/interface/RiemannFit.h b/RecoPixelVertexing/PixelTrackFitting/interface/RiemannFit.h index 52cf4b637fb37..63f139c70463c 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); } @@ -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) { @@ -254,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; @@ -514,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 @@ -908,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; 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));