Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Code performance improvements in follow-up to #31722 #33864

Closed
wants to merge 7 commits into from
Closed

Code performance improvements in follow-up to #31722 #33864

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
2cfb0d6
Reduced duplicate computations with temporary variables. (#614)
ericcano Apr 8, 2021
55be17f
Minor simplification (#614)
ericcano Apr 8, 2021
55b7e36
Reduced computations with a local variable.
ericcano Apr 9, 2021
1141a92
Simplified computations.
ericcano Apr 9, 2021
f4a03f3
Update RecoPixelVertexing/PixelTrackFitting/interface/RiemannFit.h
ericcano Apr 9, 2021
03c09bd
Removed unneeded range checking.
ericcano Apr 13, 2021
aee2dfb
Fixed code format
ericcano May 27, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
41 changes: 26 additions & 15 deletions RecoPixelVertexing/PixelTrackFitting/interface/BrokenLine.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -208,31 +208,42 @@ namespace brokenline {
const riemannFit::VectorNd<n>& varBeta) {
riemannFit::MatrixNd<n> c_uMat = riemannFit::MatrixNd<n>::Zero();
for (u_int i = 0; i < n; i++) {
// Pre compute differences used multiple time. Notation is t<j><k> = 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;
}

/*!
Expand Down
32 changes: 18 additions & 14 deletions RecoPixelVertexing/PixelTrackFitting/interface/RiemannFit.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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);
}
Expand Down Expand Up @@ -125,10 +125,13 @@ namespace riemannFit {
VectorNd<N> 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<N> s_values;
VectorNd<N> rad_lengths;
const Vector2d oVec(fast_fit(0), fast_fit(1));
Expand All @@ -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<N> scatter_cov_rad = MatrixNd<N>::Zero();
VectorNd<N> 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) {
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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<double, 2, 2> 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<N> res = p2D_rot.row(1).transpose() - aMat.transpose() * sol;
double chi2 = res.transpose() * vyInvMat * res;
Expand Down
4 changes: 2 additions & 2 deletions RecoPixelVertexing/PixelTrackFitting/plugins/storeTracks.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand All @@ -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));
Expand Down