In MkFinder::selectHitIndicesV2() use track errors for extending bin …

…search and for doing hit preselection.

RecoTracker/MkFitCore/src/Matriplex/Matriplex.h

@@ -115,6 +115,24 @@ namespace Matriplex {
       return *this;
     }
 
+    Matriplex operator-() {
+      Matriplex t;
+      for (idx_t i = 0; i < kTotSize; ++i)
+        t.fArray[i] = -fArray[i];
+      return t;
+    }
+
+    Matriplex& abs(const Matriplex& a) {
+      for (idx_t i = 0; i < kTotSize; ++i)
+        fArray[i] = std::abs(a.fArray[i]);
+      return *this;
+    }
+    Matriplex& abs() {
+      for (idx_t i = 0; i < kTotSize; ++i)
+        fArray[i] = std::abs(fArray[i]);
+      return *this;
+    }
+
     Matriplex& sqrt(const Matriplex& a) {
       for (idx_t i = 0; i < kTotSize; ++i)
         fArray[i] = std::sqrt(a.fArray[i]);
@@ -401,6 +419,12 @@ namespace Matriplex {
     return t;
   }
 
+  template <typename T, idx_t D1, idx_t D2, idx_t N>
+  MPlex<T, D1, D2, N> abs(const MPlex<T, D1, D2, N>& a) {
+    MPlex<T, D1, D2, N> t;
+    return t.abs(a);
+  }
+
   template <typename T, idx_t D1, idx_t D2, idx_t N>
   MPlex<T, D1, D2, N> sqrt(const MPlex<T, D1, D2, N>& a) {
     MPlex<T, D1, D2, N> t;
@@ -410,7 +434,7 @@ namespace Matriplex {
   template <typename T, idx_t D1, idx_t D2, idx_t N>
   MPlex<T, D1, D2, N> sqr(const MPlex<T, D1, D2, N>& a) {
     MPlex<T, D1, D2, N> t;
-    return t.sqrt(a);
+    return t.sqr(a);
   }
 
   template <typename T, idx_t D1, idx_t D2, idx_t N>

RecoTracker/MkFitCore/src/Matriplex/MatriplexSym.h

@@ -263,6 +263,14 @@ namespace Matriplex {
         a[5 * N + n] = s * c22;
       }
     }
+
+    Matriplex<T, 1, 1, N> ReduceFixedIJ(idx_t i, idx_t j) const {
+      Matriplex<T, 1, 1, N> t;
+      for (idx_t n = 0; n < N; ++n) {
+        t[n] = constAt(n, i, j);
+      }
+      return t;
+    }
   };
 
   template <typename T, idx_t D, idx_t N>

RecoTracker/MkFitCore/src/MkFinder.cc

@@ -789,6 +789,8 @@ namespace mkfit {
       mp::StatePlex sp1, sp2;
       int n_proc;
 
+      MPlexQF dphi_track, dq_track;  // 3 sigma track errors at initial state
+
       // debug & ntuple dump -- to be local in functions
       MPlexQF phi_c, dphi;
       MPlexQF q_c, qmin, qmax;
@@ -808,10 +810,10 @@ namespace mkfit {
         }
       }
 
-      void find_bin_ranges(const LayerInfo &li, const LayerOfHits &loh) {
+      void find_bin_ranges(const LayerInfo &li, const LayerOfHits &loh, const MPlexLS &err) {
         // Below made members for debugging
         // MPlexQF phi_c, dphi_min, dphi_max;
-        phi_c = mp::fast_atan2(isp.y, isp.x);
+        // phi_c = mp::fast_atan2(isp.y, isp.x);  // calculated below as difference
 
         // Matriplex::min_max(sp1.dphi, sp2.dphi, dphi_min, dphi_max);
         // the above is wrong: dalpha is not dphi --> renamed variable in State
@@ -833,33 +835,47 @@ namespace mkfit {
           //       phi_c[ii], xp1[ii], xp2[ii], pmin[ii], pmax[ii], dphi[ii]);
         }
 
+        const auto calc_err_xy = [&](const MPlexQF &x, const MPlexQF &y) {
+          return x * x * err.ReduceFixedIJ(0, 0) + y * y * err.ReduceFixedIJ(1, 1) +
+                 2.0f * x * y * err.ReduceFixedIJ(0, 1);
+        };
+
+        // Calculate dphi_track, dq_track differs for barrel/endcap
+        MPlexQF r2_c = isp.x * isp.x + isp.y * isp.y;
+        MPlexQF r2inv_c = 1.0f / r2_c;
+        MPlexQF dphidx_c = -isp.y * r2inv_c;
+        MPlexQF dphidy_c = isp.x * r2inv_c;
+        dphi_track = 3.0f * calc_err_xy(dphidx_c, dphidy_c).abs().sqrt();
+
         // MPlexQF qmin, qmax;
         if (li.is_barrel()) {
           Matriplex::min_max(sp1.z, sp2.z, qmin, qmax);
           q_c = isp.z;
+          dq_track = 3.0f * err.ReduceFixedIJ(2, 2).abs().sqrt();
         } else {
           Matriplex::min_max(Matriplex::hypot(sp1.x, sp1.y), Matriplex::hypot(sp2.x, sp2.y), qmin, qmax);
-          q_c = Matriplex::hypot(isp.x, isp.y);
+          q_c = Matriplex::sqrt(r2_c);
+          dq_track = 3.0f * (r2inv_c * calc_err_xy(isp.x, isp.y).abs()).sqrt();
         }
 
         for (int i = 0; i < p1.kTotSize; ++i) {
           // Clamp crazy sizes. This actually only happens when prop-fail flag is set.
           // const float dphi_clamp = 0.1;
           // if (dphi_min[i] > 0.0f || dphi_min[i] < -dphi_clamp) dphi_min[i] = -dphi_clamp;
           // if (dphi_max[i] < 0.0f || dphi_max[i] > dphi_clampf) dphi_max[i] = dphi_clamp;
-          p1[i] = loh.phiBinChecked(pmin[i] - PHI_BIN_EXTRA_FAC * 0.0123f);
-          p2[i] = loh.phiBinChecked(pmax[i] + PHI_BIN_EXTRA_FAC * 0.0123f);
+          p1[i] = loh.phiBinChecked(pmin[i] - dphi_track[i] - PHI_BIN_EXTRA_FAC * 0.0123f);
+          p2[i] = loh.phiBinChecked(pmax[i] + dphi_track[i] + PHI_BIN_EXTRA_FAC * 0.0123f);
 
           q0[i] = loh.qBinChecked(q_c[i]);
-          q1[i] = loh.qBinChecked(qmin[i] - Q_BIN_EXTRA_FAC * 0.5f * li.q_bin());
-          q2[i] = loh.qBinChecked(qmax[i] + Q_BIN_EXTRA_FAC * 0.5f * li.q_bin()) + 1;
+          q1[i] = loh.qBinChecked(qmin[i] - dq_track[i] - Q_BIN_EXTRA_FAC * 0.5f * li.q_bin());
+          q2[i] = loh.qBinChecked(qmax[i] + dq_track[i] + Q_BIN_EXTRA_FAC * 0.5f * li.q_bin()) + 1;
         }
       }
     };
 
     Bins B(m_Par[iI], m_Chg, N_proc);
     B.prop_to_limits(LI);
-    B.find_bin_ranges(LI, L);
+    B.find_bin_ranges(LI, L, m_Err[iI]);
 
     for (int i = 0; i < N_proc; ++i) {
       m_XHitSize[i] = 0;
@@ -982,8 +998,8 @@ namespace mkfit {
             new_ddphi = cdist(std::abs(new_phi - L.hit_phi(hi)));
             new_ddq = std::abs(new_q - L.hit_q(hi));
 
-            bool dqdphi_presel =
-                new_ddq < DDQ_PRESEL_FAC * L.hit_q_half_length(hi) && new_ddphi < DDPHI_PRESEL_FAC * 0.0123f;
+            bool dqdphi_presel = new_ddq < B.dq_track[itrack] + DDQ_PRESEL_FAC * L.hit_q_half_length(hi) &&
+                                 new_ddphi < B.dphi_track[itrack] + DDPHI_PRESEL_FAC * 0.0123f;
 
             // clang-format off
             dprintf("     SHI %3u %4u %5u  %6.3f %6.3f %6.4f %7.5f  PROP-%s  %s\n",