From b46cacfdaa908ef80a5bc93aa6c1f43f6180af9b Mon Sep 17 00:00:00 2001
From: Eric Cano <eric.cano@cern.ch>
Date: Fri, 9 Apr 2021 17:21:32 +0200
Subject: [PATCH] Reduced computations with a local variable.
 (cms-patatrack/cmssw#614)

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<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));
@@ -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) {