From b87477e6a0bcfb5a923b340e4ae7f973a05f4f62 Mon Sep 17 00:00:00 2001
From: leonardo <leonardo.giannini@cern.ch>
Date: Wed, 12 Jan 2022 15:16:27 -0800
Subject: [PATCH 1/2] seed cleaning in eta-phi bins

---
 mkFit/MkStdSeqs.cc | 35 +++++++++++++++++++++++++++++++----
 1 file changed, 31 insertions(+), 4 deletions(-)

diff --git a/mkFit/MkStdSeqs.cc b/mkFit/MkStdSeqs.cc
index 14064b0e..d6e32144 100644
--- a/mkFit/MkStdSeqs.cc
+++ b/mkFit/MkStdSeqs.cc
@@ -168,6 +168,13 @@ int clean_cms_seedtracks_iter(TrackVec *seed_ptr, const IterationConfig& itrcfg,
   std::vector<float>  d0(ns);
   int i1,i2; //for the sorting
 
+  axis_pow2_u1<float, unsigned short, 16, 8> ax_phi(-Config::PI, Config::PI);
+  axis<float, unsigned short, 8, 8> ax_eta(-3.0, 3.0, 30u);
+  binnor<unsigned int, decltype(ax_phi), decltype(ax_eta), 24, 8> b(ax_phi, ax_eta);
+
+  b.begin_registration(ns);
+
+
   for(int ts=0; ts<ns; ts++){
     const Track & tk = seeds[ts];
     nHits[ts] = tk.nFoundHits();
@@ -182,9 +189,15 @@ int clean_cms_seedtracks_iter(TrackVec *seed_ptr, const IterationConfig& itrcfg,
     y[ts] = tk.y();
     z[ts] = tk.z();
     d0[ts] = tk.d0BeamSpot(bspot.x,bspot.y);
+
+    b.register_entry_safe(oldPhi[ts], eta[ts]);     // If one is sure values are *within* axis ranges: b.register_entry(oldPhi[ts], eta[ts]);
+
   }
 
-  for(int ts=0; ts<ns; ts++){
+  b.finalize_registration();
+
+  for(int sorted_ts=0; sorted_ts<ns; sorted_ts++){
+    int ts = b.m_ranks[sorted_ts];
 
     if (not writetrack[ts]) continue;//FIXME: this speed up prevents transitive masking; check build cost!
 
@@ -200,9 +213,21 @@ int clean_cms_seedtracks_iter(TrackVec *seed_ptr, const IterationConfig& itrcfg,
     // int  n_ovlp_hits_shared = 0;
     // int  n_ovlp_tracks = 0;
 
-    //#pragma simd /* Vectorization via simd had issues with icc */
-    for (int tss= ts+1; tss<ns; tss++)
-    {
+    //these two loops may be put inside the class
+    auto phi_rng = ax_phi.Rrdr_to_N_bins(oldPhi[ts], 0.08);
+    auto eta_rng = ax_eta.Rrdr_to_N_bins(eta[ts], .1);
+
+    for (auto i_phi = phi_rng.begin; i_phi != phi_rng.end; i_phi = ax_phi.next_N_bin(i_phi)){
+    for (auto i_eta = eta_rng.begin; i_eta != eta_rng.end; i_eta = ax_eta.next_N_bin(i_eta)){
+
+    const auto cbin = b.get_content(i_phi, i_eta);
+    for (auto i = cbin.first; i < cbin.end(); ++i){
+      int tss = b.m_ranks[i];
+
+      if (not writetrack[ts]) continue;
+      if (not writetrack[tss]) continue;
+      if (tss==ts) continue;
+
       const float Pt2 = pt[tss];
 
       ////// Always require charge consistency. If different charge is assigned, do not remove seed-track
@@ -296,6 +321,8 @@ int clean_cms_seedtracks_iter(TrackVec *seed_ptr, const IterationConfig& itrcfg,
            tk.sortHitsByLayer();
       }
     } //end of inner loop over tss
+    } //eta bin
+    } //phi bin
 
     if (writetrack[ts])
     {

From 3357b74a62dae6c3f016764d2bc67a844cc675bc Mon Sep 17 00:00:00 2001
From: leonardo <leonardo.giannini@cern.ch>
Date: Mon, 24 Jan 2022 11:39:51 -0800
Subject: [PATCH 2/2] temporary

---
 mkFit/MkStdSeqs.cc      | 44 +++++++++++++++++++++++++++++++++++++----
 mkFit/buildtestMPlex.cc | 13 ++++++++++--
 2 files changed, 51 insertions(+), 6 deletions(-)

diff --git a/mkFit/MkStdSeqs.cc b/mkFit/MkStdSeqs.cc
index d6e32144..ca4ce4a3 100644
--- a/mkFit/MkStdSeqs.cc
+++ b/mkFit/MkStdSeqs.cc
@@ -482,32 +482,64 @@ void handle_duplicates(Event *m_event)
 void find_duplicates_sharedhits(TrackVec &tracks, const float fraction)
 {
   const auto ntracks = tracks.size();
+  int ncfr_h=0;
+  int ncfr_k=0;
+
+  axis_pow2_u1<float, unsigned short, 16, 8> ax_phi(-Config::PI, Config::PI);
+  axis<float, unsigned short, 8, 8> ax_ctheta(-10.0, 10.0, 30u);
+  binnor<unsigned int, decltype(ax_phi), decltype(ax_ctheta), 24, 8> b(ax_phi, ax_ctheta);
+
+  b.begin_registration(ntracks);
 
   std::vector<float>  ctheta(ntracks);  
   for (auto itrack = 0U; itrack < ntracks; itrack++)
   {
     auto &trk = tracks[itrack];
     ctheta[itrack] = 1.f/std::tan(trk.theta());
+    b.register_entry_safe(trk.momPhi(), ctheta[itrack]);     // If one is sure values are *within* axis ranges: b.register_entry(oldPhi[ts], eta[ts]);
+
   }
 
-  for (auto itrack = 0U; itrack < ntracks; itrack++)
+  b.finalize_registration(); 
+
+  //for (auto itrack = 0U; itrack < ntracks; itrack++)
+  for (auto sorted_itrack = 0U; sorted_itrack < ntracks; sorted_itrack++)
   {
+   
+    int itrack = b.m_ranks[sorted_itrack];
     auto &trk = tracks[itrack];
     auto phi1 = trk.momPhi();
     auto ctheta1 = ctheta[itrack];
+    
+    auto phi_rng = ax_phi.Rrdr_to_N_bins(phi1, 1.0);
+    auto ctheta_rng = ax_ctheta.Rrdr_to_N_bins(ctheta1, 0.24);
+    
+    for (auto i_phi = phi_rng.begin; i_phi != phi_rng.end; i_phi = ax_phi.next_N_bin(i_phi)){
+    for (auto i_ctheta = ctheta_rng.begin; i_ctheta != ctheta_rng.end; i_ctheta = ax_ctheta.next_N_bin(i_ctheta)){
 
-    for (auto jtrack = itrack + 1; jtrack < ntracks; jtrack++)
+    const auto cbin = b.get_content(i_phi, i_ctheta);
+    for (auto i = cbin.first; i < cbin.end(); ++i)
     {
+      int jtrack = b.m_ranks[i];
+    //for (auto jtrack = itrack + 1; jtrack < ntracks; jtrack++)
+    //{
+      ncfr_k++;
+      if(itrack<=jtrack) continue;
+
       auto &track2 = tracks[jtrack];
       auto sharedCount=0;
       auto sharedFirst=0;
-
+      
+      //ncfr_k++;
+      
       auto dctheta = std::abs(ctheta[jtrack] - ctheta1);
       if (dctheta > 1.) continue;
 
       auto dphi = std::abs(squashPhiMinimal(phi1 - track2.momPhi()));
       if (dphi > 1.) continue;
-
+      
+      ncfr_h++;
+      
       for (int i = 0; i < trk.nTotalHits(); ++i)
       {
 	if (trk.getHitIdx(i)<0) continue;
@@ -524,13 +556,17 @@ void find_duplicates_sharedhits(TrackVec &tracks, const float fraction)
       }
       if  ((sharedCount - sharedFirst) >= ((std::min(trk.nFoundHits(), track2.nFoundHits()) - sharedFirst) * (fraction)) )
       {
+        std::cout << " "<< (squashPhiMinimal(phi1 - track2.momPhi())) << ", ";
 	if (trk.score() > track2.score())
 	  track2.setDuplicateValue(true);
 	else
 	  trk.setDuplicateValue(true);
       }
     }
+    }
+    }
   }
+  std::cout << "compare kin "<< ncfr_k << " compare hits "<< ncfr_h << " ---- NNN t ---- " <<ntracks<< std::endl;
   tracks.erase(std::remove_if(tracks.begin(),tracks.end(),[](auto track){return track.getDuplicateValue();}),tracks.end());
 }
 
diff --git a/mkFit/buildtestMPlex.cc b/mkFit/buildtestMPlex.cc
index 4f49d182..5f61aebe 100644
--- a/mkFit/buildtestMPlex.cc
+++ b/mkFit/buildtestMPlex.cc
@@ -584,8 +584,17 @@ std::vector<double> runBtpCe_MultiIter(Event& ev, const EventOfHits &eoh, MkBuil
        { return StdSeq::qfilter_nan_n_silly(t); });
 
       builder.select_best_comb_cands(true); // true -> clear m_tracks as they were already filled once above
-
-      StdSeq::find_and_remove_duplicates(builder.ref_tracks_nc(), itconf);
+      double t0 = dtime();
+      int isize = builder.ref_tracks_nc().size();
+      TrackVec xxx; 
+      for (int i=0;i<1;++i) {
+      xxx = builder.ref_tracks_nc();
+      StdSeq::find_and_remove_duplicates(xxx, itconf);
+      }
+      builder.ref_tracks_nc() = xxx;
+      //StdSeq::find_and_remove_duplicates(builder.ref_tracks_nc(), itconf);
+      printf("Dup clean of %d tracks, out_tracks %d, 1000 times, N_comparisons=, took %.5fs, iter. %d\n",
+             isize, builder.ref_tracks_nc().size(), dtime() - t0, itconf.m_track_algorithm);
       builder.export_tracks(ev.fitTracks_);
     }