mkFit: overlap area candidate re-check and start using the new hit selection V2 in pixelLess iteration

RecoTracker/MkFitCMS/standalone/Shell.cc

@@ -325,7 +325,9 @@ namespace mkfit {
     reco tracks labels are seed indices.
     seed labels are sim track indices
     --
-    mkfit labels are seed indices in given iteration after cleaning (at seed load-time)
+    mkfit labels are seed indices in given iteration after cleaning (at seed load-time).
+          This is no longer true -- was done like that in branch where this code originated from.
+          It seems the label is the same as seed label.
   */
 
   int Shell::LabelFromHits(Track &t, bool replace, float good_frac) {

RecoTracker/MkFitCore/interface/IterationConfig.h

@@ -90,6 +90,8 @@ namespace mkfit {
     float chi2Cut_min = 15.0;
     float chi2CutOverlap = 3.5;
     float pTCutOverlap = 0.0;
+    bool recheckOverlap = false;
+    bool useHitSelectionV2 = false;
 
     //quality filter params
     int minHitsQF = 4;

RecoTracker/MkFitCore/interface/TrackStructures.h

@@ -224,6 +224,7 @@ namespace mkfit {
     }
 
     void addHitIdx(int hitIdx, int hitLyr, float chi2);
+    bool popOverlap();
 
     HoTNode& refLastHoTNode();              // for filling up overlap info
     const HoTNode& refLastHoTNode() const;  // for dump traversal
@@ -564,6 +565,25 @@ namespace mkfit {
     }
   }
 
+  inline bool TrackCand::popOverlap() {
+    auto popHitIdx = getLastHitIdx();
+    auto popHitLyr = getLastHitLyr();
+    auto popPrev = refLastHoTNode().m_prev_idx;
+    auto popChi2 = refLastHoTNode().m_chi2;
+    // sanity checks first, then just shift lastHitIdx_ to popPrev
+    if (lastHitIdx_ == 0 || popHitIdx < 0)
+      return false;
+    auto prevHitLyr = m_comb_candidate->hot(popPrev).layer;
+    auto prevHitIdx = m_comb_candidate->hot(popPrev).index;
+    if (popHitLyr != prevHitLyr || prevHitIdx < 0)
+      return false;
+    lastHitIdx_ = popPrev;
+
+    --nFoundHits_;
+    chi2_ -= popChi2;
+    --nOverlapHits_;
+    return true;
+  }
   //==============================================================================
 
   class EventOfCombCandidates {

RecoTracker/MkFitCore/src/CandCloner.cc

@@ -24,11 +24,13 @@ namespace mkfit {
 
   void CandCloner::begin_eta_bin(EventOfCombCandidates *e_o_ccs,
                                  std::vector<UpdateIndices> *update_list,
+                                 std::vector<UpdateIndices> *overlap_list,
                                  std::vector<std::vector<TrackCand>> *extra_cands,
                                  int start_seed,
                                  int n_seeds) {
     mp_event_of_comb_candidates = e_o_ccs;
     mp_kalman_update_list = update_list;
+    mp_kalman_overlap_list = overlap_list;
     mp_extra_cands = extra_cands;
     m_start_seed = start_seed;
     m_n_seeds = n_seeds;
@@ -50,6 +52,7 @@ namespace mkfit {
     m_idx_max_prev = 0;
 
     mp_kalman_update_list->clear();
+    mp_kalman_overlap_list->clear();
 
 #ifdef CC_TIME_LAYER
     t_lay = dtime();
@@ -193,14 +196,21 @@ namespace mkfit {
             break;
 
           if (h2a.hitIdx >= 0) {
-            mp_kalman_update_list->emplace_back(UpdateIndices(m_start_seed + is, n_pushed, h2a.hitIdx));
-
-            // set the overlap if we have it and and pT > pTCutOverlap
+            // set the overlap if we have it and pT > pTCutOverlap
             HitMatch *hm;
             if (tc.pT() > mp_iteration_params->pTCutOverlap &&
                 (hm = ccand[h2a.trkIdx].findOverlap(h2a.hitIdx, h2a.module))) {
-              tc.addHitIdx(hm->m_hit_idx, m_layer, 0);
-              tc.incOverlapCount();
+              if (mp_iteration_params->recheckOverlap) {
+                // Special overlap_list if the overlap hit needs to be re-checked after primary update.
+                mp_kalman_overlap_list->emplace_back(
+                    UpdateIndices(m_start_seed + is, n_pushed, h2a.hitIdx, hm->m_hit_idx));
+              } else {
+                tc.addHitIdx(hm->m_hit_idx, m_layer, 0);
+                tc.incOverlapCount();
+                mp_kalman_update_list->emplace_back(UpdateIndices(m_start_seed + is, n_pushed, h2a.hitIdx, -1));
+              }
+            } else {
+              mp_kalman_update_list->emplace_back(UpdateIndices(m_start_seed + is, n_pushed, h2a.hitIdx, -1));
             }
           }
 

RecoTracker/MkFitCore/src/CandCloner.h

@@ -25,6 +25,7 @@ namespace mkfit {
 
     void begin_eta_bin(EventOfCombCandidates *e_o_ccs,
                        std::vector<UpdateIndices> *update_list,
+                       std::vector<UpdateIndices> *overlap_list,
                        std::vector<std::vector<TrackCand>> *extra_cands,
                        int start_seed,
                        int n_seeds);
@@ -56,7 +57,7 @@ namespace mkfit {
 
     const IterationParams *mp_iteration_params = nullptr;
     EventOfCombCandidates *mp_event_of_comb_candidates;
-    std::vector<UpdateIndices> *mp_kalman_update_list;
+    std::vector<UpdateIndices> *mp_kalman_update_list, *mp_kalman_overlap_list;
     std::vector<std::vector<TrackCand>> *mp_extra_cands;
 
 #if defined(CC_TIME_ETA) or defined(CC_TIME_LAYER)

RecoTracker/MkFitCore/src/FindingFoos.h

@@ -16,6 +16,10 @@ namespace mkfit {
   const MPlexLS &, const MPlexLV &, MPlexQI &, const MPlexHS &, const MPlexHV &, MPlexLS &, MPlexLV &, MPlexQI &, \
       const int, const PropagationFlags &, const bool
 
+#define COMPUTE_CHI2_AND_UPDATE_ARGS                                                                              \
+  const MPlexLS &, const MPlexLV &, MPlexQI &, const MPlexHS &, const MPlexHV &, MPlexQF &, MPlexLS &, MPlexLV &, \
+      MPlexQI &, const int, const PropagationFlags, const bool
+
   class FindingFoos {
   public:
     void (*m_compute_chi2_foo)(COMPUTE_CHI2_ARGS);

RecoTracker/MkFitCore/src/HitStructures.cc

@@ -91,6 +91,9 @@ namespace mkfit {
       m_hit_infos.reserve(m_n_hits);
     }
 
+    // Factor to get from hit sigma to half-length in q direction.
+    const float hl_fac = is_pixel() ? 3.0f : std::sqrt(3.0f);
+
     for (unsigned int i = 0; i < m_n_hits; ++i) {
       const Hit &h = hitv[i];
 
@@ -100,13 +103,12 @@ namespace mkfit {
       m_binnor.register_entry_safe(phi, q);
 
       if (Config::usePhiQArrays) {
-        const float sqrt3 = std::sqrt(3);
         float half_length, qbar;
         if (m_is_barrel) {
-          half_length = sqrt3 * std::sqrt(h.ezz());
+          half_length = hl_fac * std::sqrt(h.ezz());
           qbar = h.r();
         } else {
-          half_length = sqrt3 * std::sqrt(h.exx() + h.eyy());
+          half_length = hl_fac * std::sqrt(h.exx() + h.eyy());
           qbar = h.z();
         }
         hinfos.emplace_back(HitInfo({phi, q, half_length, qbar}));
@@ -168,13 +170,14 @@ namespace mkfit {
     m_binnor.register_entry_safe(phi, q);
 
     if (Config::usePhiQArrays) {
-      const float sqrt3 = std::sqrt(3);
+      // Factor to get from hit sigma to half-length in q direction.
+      const float hl_fac = is_pixel() ? 3.0f : std::sqrt(3.0f);
       float half_length, qbar;
       if (m_is_barrel) {
-        half_length = sqrt3 * std::sqrt(h.ezz());
+        half_length = hl_fac * std::sqrt(h.ezz());
         qbar = h.r();
       } else {
-        half_length = sqrt3 * std::sqrt(h.exx() + h.eyy());
+        half_length = hl_fac * std::sqrt(h.exx() + h.eyy());
         qbar = h.z();
       }
       m_hit_infos.emplace_back(HitInfo({phi, q, half_length, qbar}));

RecoTracker/MkFitCore/src/IterationConfig.cc

@@ -61,6 +61,8 @@ namespace mkfit {
                                    /* float */ chi2Cut_min,
                                    /* float */ chi2CutOverlap,
                                    /* float */ pTCutOverlap,
+                                   /* bool */ recheckOverlap,
+                                   /* bool */ useHitSelectionV2,
                                    /* int */ minHitsQF,
                                    /* float */ minPtCut,
                                    /* unsigned int */ maxClusterSize)

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/MkBuilder.cc

@@ -979,15 +979,16 @@ namespace mkfit {
     const int n_seeds = end_seed - start_seed;
 
     std::vector<std::pair<int, int>> seed_cand_idx;
-    std::vector<UpdateIndices> seed_cand_update_idx;
+    std::vector<UpdateIndices> seed_cand_update_idx, seed_cand_overlap_idx;
     seed_cand_idx.reserve(n_seeds * params.maxCandsPerSeed);
     seed_cand_update_idx.reserve(n_seeds * params.maxCandsPerSeed);
+    seed_cand_overlap_idx.reserve(n_seeds * params.maxCandsPerSeed);
 
     std::vector<std::vector<TrackCand>> extra_cands(n_seeds);
     for (int ii = 0; ii < n_seeds; ++ii)
       extra_cands[ii].reserve(params.maxCandsPerSeed);
 
-    cloner.begin_eta_bin(&eoccs, &seed_cand_update_idx, &extra_cands, start_seed, n_seeds);
+    cloner.begin_eta_bin(&eoccs, &seed_cand_update_idx, &seed_cand_overlap_idx, &extra_cands, start_seed, n_seeds);
 
     // Loop over layers, starting from after the seed.
 
@@ -1087,7 +1088,10 @@ namespace mkfit {
 
         dprint("now get hit range");
 
-        mkfndr->selectHitIndices(layer_of_hits, end - itrack);
+        if (iter_params.useHitSelectionV2)
+          mkfndr->selectHitIndicesV2(layer_of_hits, end - itrack);
+        else
+          mkfndr->selectHitIndices(layer_of_hits, end - itrack);
 
         find_tracks_handle_missed_layers(
             mkfndr, layer_info, extra_cands, seed_cand_idx, region, start_seed, itrack, end);
@@ -1115,6 +1119,9 @@ namespace mkfit {
 
       // Update loop of best candidates. CandCloner prepares the list of those
       // that need update (excluding all those with negative last hit index).
+      // This is split into two sections - candidates without overlaps and with overlaps.
+      // On CMS PU-50 the ratio of those is ~ 65 : 35 over all iterations.
+      // Note, overlap recheck is only enabled for some iterations, e.g. pixelLess.
 
       const int theEndUpdater = seed_cand_update_idx.size();
 
@@ -1129,6 +1136,43 @@ namespace mkfit {
         mkfndr->copyOutParErr(eoccs.refCandidates_nc(), end - itrack, false);
       }
 
+      const int theEndOverlapper = seed_cand_overlap_idx.size();
+
+      for (int itrack = 0; itrack < theEndOverlapper; itrack += NN) {
+        const int end = std::min(itrack + NN, theEndOverlapper);
+
+        mkfndr->inputTracksAndHits(eoccs.refCandidates(), layer_of_hits, seed_cand_overlap_idx, itrack, end, true);
+
+        mkfndr->updateWithLoadedHit(end - itrack, layer_of_hits, fnd_foos);
+
+        mkfndr->copyOutParErr(eoccs.refCandidates_nc(), end - itrack, false);
+
+        mkfndr->inputOverlapHits(layer_of_hits, seed_cand_overlap_idx, itrack, end);
+
+        // XXXX Could also be calcChi2AndUpdate(), then copy-out would have to be done
+        // below, choosing appropriate slot (with or without the overlap hit).
+        // Probably in a dedicated MkFinder copyOutXyzz function.
+        mkfndr->chi2OfLoadedHit(end - itrack, fnd_foos);
+
+        for (int ii = itrack; ii < end; ++ii) {
+          const int fi = ii - itrack;
+          TrackCand &tc = eoccs[seed_cand_overlap_idx[ii].seed_idx][seed_cand_overlap_idx[ii].cand_idx];
+
+          // XXXX For now we DO NOT use chi2 as this was how things were done before the post-update
+          // chi2 check. To use it we should retune scoring function (might be even simpler).
+          auto chi2Ovlp = mkfndr->m_Chi2[fi];
+          if (mkfndr->m_FailFlag[fi] == 0 && chi2Ovlp >= 0.0f && chi2Ovlp <= 60.0f) {
+            auto scoreCand =
+                getScoreCand(st_par.m_track_scorer, tc, true /*penalizeTailMissHits*/, true /*inFindCandidates*/);
+            tc.addHitIdx(seed_cand_overlap_idx[ii].ovlp_idx, curr_layer, chi2Ovlp);
+            tc.incOverlapCount();
+            auto scoreCandOvlp = getScoreCand(st_par.m_track_scorer, tc, true, true);
+            if (scoreCand > scoreCandOvlp)
+              tc.popOverlap();
+          }
+        }
+      }
+
       // Check if cands are sorted, as expected.
 #ifdef DEBUG
       for (int iseed = start_seed; iseed < end_seed; ++iseed) {