Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-08-17 23:57:47

0001 #include "RecoTracker/MkFitCMS/interface/MkStdSeqs.h"
0002 
0003 #include "RecoTracker/MkFitCore/interface/HitStructures.h"
0004 #include "RecoTracker/MkFitCore/interface/IterationConfig.h"
0005 
0006 #include "RecoTracker/MkFitCore/interface/binnor.h"
0007 
0008 #include "oneapi/tbb/parallel_for.h"
0009 
0010 namespace mkfit {
0011 
0012   namespace StdSeq {
0013 
0014     //=========================================================================
0015     // Hit processing
0016     //=========================================================================
0017 
0018     void loadDeads(EventOfHits &eoh, const std::vector<DeadVec> &deadvectors) {
0019       for (size_t il = 0; il < deadvectors.size(); il++) {
0020         eoh.suckInDeads(int(il), deadvectors[il]);
0021       }
0022     }
0023 
0024     // Loading hits in CMSSW from two "large multi-layer vectors".
0025     // orig_hitvectors[0] - pixels,
0026     // orig_hitvectors[1] - strips.
0027 
0028     void cmssw_LoadHits_Begin(EventOfHits &eoh, const std::vector<const HitVec *> &orig_hitvectors) {
0029       eoh.reset();
0030       for (int i = 0; i < eoh.nLayers(); ++i) {
0031         auto &&l = eoh[i];
0032         l.beginRegistrationOfHits(*orig_hitvectors[l.is_pixel() ? 0 : 1]);
0033       }
0034     }
0035 
0036     // Loop with LayerOfHits::registerHit(int idx) - it takes Hit out of original HitVec to
0037     // extract phi, r/z, and calculate qphifines
0038     //
0039     // Something like what is done in MkFitInputConverter::convertHits
0040     //
0041     // Problem is I don't know layers for each large-vector;
0042     // Also, layer is calculated for each detset when looping over the HitCollection
0043 
0044     void cmssw_LoadHits_End(EventOfHits &eoh) {
0045       for (int i = 0; i < eoh.nLayers(); ++i) {
0046         auto &&l = eoh[i];
0047         l.endRegistrationOfHits(false);
0048       }
0049     }
0050 
0051     //=========================================================================
0052     // Hit-index mapping / remapping
0053     //=========================================================================
0054 
0055     void cmssw_Map_TrackHitIndices(const EventOfHits &eoh, TrackVec &seeds) {
0056       for (auto &&track : seeds) {
0057         for (int i = 0; i < track.nTotalHits(); ++i) {
0058           const int hitidx = track.getHitIdx(i);
0059           const int hitlyr = track.getHitLyr(i);
0060           if (hitidx >= 0) {
0061             const auto &loh = eoh[hitlyr];
0062             track.setHitIdx(i, loh.getHitIndexFromOriginal(hitidx));
0063           }
0064         }
0065       }
0066     }
0067 
0068     void cmssw_ReMap_TrackHitIndices(const EventOfHits &eoh, TrackVec &out_tracks) {
0069       for (auto &&track : out_tracks) {
0070         for (int i = 0; i < track.nTotalHits(); ++i) {
0071           const int hitidx = track.getHitIdx(i);
0072           const int hitlyr = track.getHitLyr(i);
0073           if (hitidx >= 0) {
0074             const auto &loh = eoh[hitlyr];
0075             track.setHitIdx(i, loh.getOriginalHitIndex(hitidx));
0076           }
0077         }
0078       }
0079     }
0080 
0081     //=========================================================================
0082     // Seed cleaning (multi-iter)
0083     //=========================================================================
0084     int clean_cms_seedtracks_iter(TrackVec *seed_ptr, const IterationConfig &itrcfg, const BeamSpot &bspot) {
0085       using Algo = TrackBase::TrackAlgorithm;
0086 
0087       const float etamax_brl = Config::c_etamax_brl;
0088       const float dpt_common = Config::c_dpt_common;
0089 
0090       const float dzmax_bh = itrcfg.m_params.c_dzmax_bh;
0091       const float drmax_bh = itrcfg.m_params.c_drmax_bh;
0092       const float dzmax_eh = itrcfg.m_params.c_dzmax_eh;
0093       const float drmax_eh = itrcfg.m_params.c_drmax_eh;
0094       const float dzmax_bl = itrcfg.m_params.c_dzmax_bl;
0095       const float drmax_bl = itrcfg.m_params.c_drmax_bl;
0096       const float dzmax_el = itrcfg.m_params.c_dzmax_el;
0097       const float drmax_el = itrcfg.m_params.c_drmax_el;
0098 
0099       const float ptmin_hpt = itrcfg.m_params.c_ptthr_hpt;
0100 
0101       const float dzmax2_inv_bh = 1.f / (dzmax_bh * dzmax_bh);
0102       const float drmax2_inv_bh = 1.f / (drmax_bh * drmax_bh);
0103       const float dzmax2_inv_eh = 1.f / (dzmax_eh * dzmax_eh);
0104       const float drmax2_inv_eh = 1.f / (drmax_eh * drmax_eh);
0105       const float dzmax2_inv_bl = 1.f / (dzmax_bl * dzmax_bl);
0106       const float drmax2_inv_bl = 1.f / (drmax_bl * drmax_bl);
0107       const float dzmax2_inv_el = 1.f / (dzmax_el * dzmax_el);
0108       const float drmax2_inv_el = 1.f / (drmax_el * drmax_el);
0109 
0110       // Merge hits from overlapping seeds?
0111       // For now always true, we require extra hits after seed,
0112       // except for lowPtQuadStep, where we only merge hits for seeds at low pT and large pseudo-rapidity
0113       const bool merge_hits = true;  // itrcfg.merge_seed_hits_during_cleaning();
0114       const float ptmax_merge_lowPtQuad = 0.2;
0115       const float etamin_merge_lowPtQuad = 1.5;
0116 
0117       if (seed_ptr == nullptr)
0118         return 0;
0119       TrackVec &seeds = *seed_ptr;
0120 
0121       const int ns = seeds.size();
0122 #ifdef DEBUG
0123       std::cout << "before seed cleaning " << seeds.size() << std::endl;
0124 #endif
0125       TrackVec cleanSeedTracks;
0126       cleanSeedTracks.reserve(ns);
0127       std::vector<bool> writetrack(ns, true);
0128 
0129       const float invR1GeV = 1.f / Config::track1GeVradius;
0130 
0131       std::vector<int> nHits(ns);
0132       std::vector<int> charge(ns);
0133       std::vector<float> oldPhi(ns);
0134       std::vector<float> pos2(ns);
0135       std::vector<float> eta(ns);
0136       std::vector<float> ctheta(ns);
0137       std::vector<float> invptq(ns);
0138       std::vector<float> pt(ns);
0139       std::vector<float> x(ns);
0140       std::vector<float> y(ns);
0141       std::vector<float> z(ns);
0142       std::vector<float> d0(ns);
0143       int i1, i2;  //for the sorting
0144 
0145       axis_pow2_u1<float, unsigned short, 16, 8> ax_phi(-Const::PI, Const::PI);
0146       axis<float, unsigned short, 8, 8> ax_eta(-3.0, 3.0, 30u);
0147       binnor<unsigned int, decltype(ax_phi), decltype(ax_eta), 24, 8> phi_eta_binnor(ax_phi, ax_eta);
0148 
0149       phi_eta_binnor.begin_registration(ns);
0150 
0151       for (int ts = 0; ts < ns; ts++) {
0152         const Track &tk = seeds[ts];
0153         nHits[ts] = tk.nFoundHits();
0154         charge[ts] = tk.charge();
0155         oldPhi[ts] = tk.momPhi();
0156         pos2[ts] = std::pow(tk.x(), 2) + std::pow(tk.y(), 2);
0157         eta[ts] = tk.momEta();
0158         ctheta[ts] = 1.f / std::tan(tk.theta());
0159         invptq[ts] = tk.charge() * tk.invpT();
0160         pt[ts] = tk.pT();
0161         x[ts] = tk.x();
0162         y[ts] = tk.y();
0163         z[ts] = tk.z();
0164         d0[ts] = tk.d0BeamSpot(bspot.x, bspot.y);
0165 
0166         phi_eta_binnor.register_entry_safe(oldPhi[ts], eta[ts]);
0167         // If one is sure values are *within* axis ranges: b.register_entry(oldPhi[ts], eta[ts]);
0168       }
0169 
0170       phi_eta_binnor.finalize_registration();
0171 
0172       for (int sorted_ts = 0; sorted_ts < ns; sorted_ts++) {
0173         int ts = phi_eta_binnor.m_ranks[sorted_ts];
0174 
0175         if (not writetrack[ts])
0176           continue;  // Note: this speed up prevents transitive masking (possibly marginal gain).
0177 
0178         const float oldPhi1 = oldPhi[ts];
0179         const float pos2_first = pos2[ts];
0180         const float eta1 = eta[ts];
0181         const float pt1 = pt[ts];
0182         const float invptq_first = invptq[ts];
0183 
0184         // To study some more details -- need EventOfHits for this
0185         int n_ovlp_hits_added = 0;
0186 
0187         auto phi_rng = ax_phi.from_R_rdr_to_N_bins(oldPhi[ts], 0.08f);
0188         auto eta_rng = ax_eta.from_R_rdr_to_N_bins(eta[ts], .1f);
0189 
0190         for (auto i_phi = phi_rng.begin; i_phi != phi_rng.end; i_phi = ax_phi.next_N_bin(i_phi)) {
0191           for (auto i_eta = eta_rng.begin; i_eta != eta_rng.end; i_eta = ax_eta.next_N_bin(i_eta)) {
0192             const auto cbin = phi_eta_binnor.get_content(i_phi, i_eta);
0193             for (auto i = cbin.first; i < cbin.end(); ++i) {
0194               int tss = phi_eta_binnor.m_ranks[i];
0195 
0196               if (not writetrack[ts])
0197                 continue;
0198               if (not writetrack[tss])
0199                 continue;
0200               if (tss == ts)
0201                 continue;
0202 
0203               const float pt2 = pt[tss];
0204 
0205               // Always require charge consistency. If different charge is assigned, do not remove seed-track
0206               if (charge[tss] != charge[ts])
0207                 continue;
0208 
0209               const float thisDPt = std::abs(pt2 - pt1);
0210               // Require pT consistency between seeds. If dpT is large, do not remove seed-track.
0211               if (thisDPt > dpt_common * pt1)
0212                 continue;
0213 
0214               const float eta2 = eta[tss];
0215               const float deta2 = std::pow(eta1 - eta2, 2);
0216 
0217               const float oldPhi2 = oldPhi[tss];
0218 
0219               const float pos2_second = pos2[tss];
0220               const float thisDXYSign05 = pos2_second > pos2_first ? -0.5f : 0.5f;
0221 
0222               const float thisDXY = thisDXYSign05 * sqrt(std::pow(x[ts] - x[tss], 2) + std::pow(y[ts] - y[tss], 2));
0223 
0224               const float invptq_second = invptq[tss];
0225 
0226               const float newPhi1 = oldPhi1 - thisDXY * invR1GeV * invptq_first;
0227               const float newPhi2 = oldPhi2 + thisDXY * invR1GeV * invptq_second;
0228 
0229               const float dphi = cdist(std::abs(newPhi1 - newPhi2));
0230 
0231               const float dr2 = deta2 + dphi * dphi;
0232 
0233               const float thisDZ = z[ts] - z[tss] - thisDXY * (ctheta[ts] + ctheta[tss]);
0234               const float dz2 = thisDZ * thisDZ;
0235 
0236               // Reject tracks within dR-dz elliptical window.
0237               // Adaptive thresholds, based on observation that duplicates are more abundant at large pseudo-rapidity and low track pT
0238               bool overlapping = false;
0239               if (std::abs(eta1) < etamax_brl) {
0240                 if (pt1 > ptmin_hpt) {
0241                   if (dz2 * dzmax2_inv_bh + dr2 * drmax2_inv_bh < 1.0f)
0242                     overlapping = true;
0243                 } else {
0244                   if (dz2 * dzmax2_inv_bl + dr2 * drmax2_inv_bl < 1.0f)
0245                     overlapping = true;
0246                 }
0247               } else {
0248                 if (pt1 > ptmin_hpt) {
0249                   if (dz2 * dzmax2_inv_eh + dr2 * drmax2_inv_eh < 1.0f)
0250                     overlapping = true;
0251                 } else {
0252                   if (dz2 * dzmax2_inv_el + dr2 * drmax2_inv_el < 1.0f)
0253                     overlapping = true;
0254                 }
0255               }
0256 
0257               if (overlapping) {
0258                 //Mark tss as a duplicate
0259                 i1 = ts;
0260                 i2 = tss;
0261                 if (d0[tss] > d0[ts])
0262                   writetrack[tss] = false;
0263                 else {
0264                   writetrack[ts] = false;
0265                   i2 = ts;
0266                   i1 = tss;
0267                 }
0268                 // Add hits from tk2 to the seed we are keeping.
0269                 // NOTE: We have a limit in Track::Status for the number of seed hits.
0270                 //       There is a check at entry and after adding of a new hit.
0271                 Track &tk = seeds[i1];
0272                 if (merge_hits && tk.nTotalHits() < Track::Status::kMaxSeedHits &&
0273                     (Algo(itrcfg.m_track_algorithm) != Algo::lowPtQuadStep ||
0274                      (pt1 < ptmax_merge_lowPtQuad && std::abs(eta1) > etamin_merge_lowPtQuad))) {
0275                   const Track &tk2 = seeds[i2];
0276                   //We are not actually fitting to the extra hits; use chi2 of 0
0277                   float fakeChi2 = 0.0;
0278 
0279                   for (int j = 0; j < tk2.nTotalHits(); ++j) {
0280                     int hitidx = tk2.getHitIdx(j);
0281                     int hitlyr = tk2.getHitLyr(j);
0282                     if (hitidx >= 0) {
0283                       bool unique = true;
0284                       for (int i = 0; i < tk.nTotalHits(); ++i) {
0285                         if ((hitidx == tk.getHitIdx(i)) && (hitlyr == tk.getHitLyr(i))) {
0286                           unique = false;
0287                           break;
0288                         }
0289                       }
0290                       if (unique) {
0291                         tk.addHitIdx(tk2.getHitIdx(j), tk2.getHitLyr(j), fakeChi2);
0292                         ++n_ovlp_hits_added;
0293                         if (tk.nTotalHits() >= Track::Status::kMaxSeedHits)
0294                           break;
0295                       }
0296                     }
0297                   }
0298                 }
0299                 if (n_ovlp_hits_added > 0)
0300                   tk.sortHitsByLayer();
0301               }
0302             }  //end of inner loop over tss
0303           }    //eta bin
0304         }      //phi bin
0305 
0306         if (writetrack[ts]) {
0307           cleanSeedTracks.emplace_back(seeds[ts]);
0308         }
0309       }
0310 
0311       seeds.swap(cleanSeedTracks);
0312 
0313 #ifdef DEBUG
0314       {
0315         const int ns2 = seeds.size();
0316         printf("Number of CMS seeds before %d --> after %d cleaning\n", ns, ns2);
0317 
0318         for (int it = 0; it < ns2; it++) {
0319           const Track &ss = seeds[it];
0320           printf("  %3i q=%+i pT=%7.3f eta=% 7.3f nHits=%i label=% i\n",
0321                  it,
0322                  ss.charge(),
0323                  ss.pT(),
0324                  ss.momEta(),
0325                  ss.nFoundHits(),
0326                  ss.label());
0327         }
0328       }
0329 #endif
0330 
0331 #ifdef DEBUG
0332       std::cout << "AFTER seed cleaning " << seeds.size() << std::endl;
0333 #endif
0334 
0335       return seeds.size();
0336     }
0337 
0338     //=========================================================================
0339     // Duplicate cleaning
0340     //=========================================================================
0341 
0342     void find_duplicates(TrackVec &tracks) {
0343       const auto ntracks = tracks.size();
0344       float eta1, phi1, pt1, deta, dphi, dr2;
0345 
0346       if (ntracks == 0) {
0347         return;
0348       }
0349       for (auto itrack = 0U; itrack < ntracks - 1; itrack++) {
0350         auto &track = tracks[itrack];
0351         using Algo = TrackBase::TrackAlgorithm;
0352         auto const algo = track.algorithm();
0353         if (algo == Algo::pixelLessStep || algo == Algo::tobTecStep)
0354           continue;
0355         eta1 = track.momEta();
0356         phi1 = track.momPhi();
0357         pt1 = track.pT();
0358         for (auto jtrack = itrack + 1; jtrack < ntracks; jtrack++) {
0359           auto &track2 = tracks[jtrack];
0360           if (track.label() == track2.label())
0361             continue;
0362           if (track.algoint() != track2.algoint())
0363             continue;
0364 
0365           deta = std::abs(track2.momEta() - eta1);
0366           if (deta > Config::maxdEta)
0367             continue;
0368 
0369           dphi = std::abs(squashPhiMinimal(phi1 - track2.momPhi()));
0370           if (dphi > Config::maxdPhi)
0371             continue;
0372 
0373           float maxdR = Config::maxdR;
0374           float maxdRSquared = maxdR * maxdR;
0375           if (std::abs(eta1) > 2.5f)
0376             maxdRSquared *= 16.0f;
0377           else if (std::abs(eta1) > 1.44f)
0378             maxdRSquared *= 9.0f;
0379           dr2 = dphi * dphi + deta * deta;
0380           if (dr2 < maxdRSquared) {
0381             //Keep track with best score
0382             if (track.score() > track2.score())
0383               track2.setDuplicateValue(true);
0384             else
0385               track.setDuplicateValue(true);
0386             continue;
0387           } else {
0388             if (pt1 == 0)
0389               continue;
0390             if (track2.pT() == 0)
0391               continue;
0392 
0393             if (std::abs((1 / track2.pT()) - (1 / pt1)) < Config::maxdPt) {
0394               if (Config::useHitsForDuplicates) {
0395                 float numHitsShared = 0;
0396                 for (int ihit2 = 0; ihit2 < track2.nTotalHits(); ihit2++) {
0397                   const int hitidx2 = track2.getHitIdx(ihit2);
0398                   const int hitlyr2 = track2.getHitLyr(ihit2);
0399                   if (hitidx2 >= 0) {
0400                     auto const it = std::find_if(track.beginHitsOnTrack(),
0401                                                  track.endHitsOnTrack(),
0402                                                  [&hitidx2, &hitlyr2](const HitOnTrack &element) {
0403                                                    return (element.index == hitidx2 && element.layer == hitlyr2);
0404                                                  });
0405                     if (it != track.endHitsOnTrack())
0406                       numHitsShared++;
0407                   }
0408                 }
0409 
0410                 float fracHitsShared = numHitsShared / std::min(track.nFoundHits(), track2.nFoundHits());
0411                 //Only remove one of the tracks if they share at least X% of the hits (denominator is the shorter track)
0412                 if (fracHitsShared < Config::minFracHitsShared)
0413                   continue;
0414               }
0415               //Keep track with best score
0416               if (track.score() > track2.score())
0417                 track2.setDuplicateValue(true);
0418               else
0419                 track.setDuplicateValue(true);
0420             }  //end of if dPt
0421           }    //end of else
0422         }      //end of loop over track2
0423       }        //end of loop over track1
0424     }
0425 
0426     void remove_duplicates(TrackVec &tracks) {
0427       tracks.erase(std::remove_if(tracks.begin(), tracks.end(), [](auto track) { return track.getDuplicateValue(); }),
0428                    tracks.end());
0429     }
0430 
0431     //=========================================================================
0432     // SHARED HITS DUPLICATE CLEANING
0433     //=========================================================================
0434 
0435     void find_duplicates_sharedhits(TrackVec &tracks, const float fraction) {
0436       const auto ntracks = tracks.size();
0437 
0438       std::vector<float> ctheta(ntracks);
0439       std::multimap<int, int> hitMap;
0440 
0441       for (auto itrack = 0U; itrack < ntracks; itrack++) {
0442         auto &trk = tracks[itrack];
0443         ctheta[itrack] = 1.f / std::tan(trk.theta());
0444         for (int i = 0; i < trk.nTotalHits(); ++i) {
0445           if (trk.getHitIdx(i) < 0)
0446             continue;
0447           int a = trk.getHitLyr(i);
0448           int b = trk.getHitIdx(i);
0449           hitMap.insert(std::make_pair(b * 1000 + a, i > 0 ? itrack : -itrack));  //negative for first hit in trk
0450         }
0451       }
0452 
0453       for (auto itrack = 0U; itrack < ntracks; itrack++) {
0454         auto &trk = tracks[itrack];
0455         auto phi1 = trk.momPhi();
0456         auto ctheta1 = ctheta[itrack];
0457 
0458         std::map<int, int> sharingMap;
0459         for (int i = 0; i < trk.nTotalHits(); ++i) {
0460           if (trk.getHitIdx(i) < 0)
0461             continue;
0462           int a = trk.getHitLyr(i);
0463           int b = trk.getHitIdx(i);
0464           auto range = hitMap.equal_range(b * 1000 + a);
0465           for (auto it = range.first; it != range.second; ++it) {
0466             if (std::abs(it->second) >= (int)itrack)
0467               continue;  // don't check your own hits (==) nor sym. checks (>)
0468             if (i == 0 && it->second < 0)
0469               continue;  // shared first - first is not counted
0470             sharingMap[std::abs(it->second)]++;
0471           }
0472         }
0473 
0474         for (const auto &elem : sharingMap) {
0475           auto &track2 = tracks[elem.first];
0476 
0477           // broad dctheta-dphi compatibility checks; keep mostly to preserve consistency with old results
0478           auto dctheta = std::abs(ctheta[elem.first] - ctheta1);
0479           if (dctheta > 1.)
0480             continue;
0481 
0482           auto dphi = std::abs(squashPhiMinimal(phi1 - track2.momPhi()));
0483           if (dphi > 1.)
0484             continue;
0485 
0486           if (elem.second >= std::min(trk.nFoundHits(), track2.nFoundHits()) * fraction) {
0487             if (trk.score() > track2.score())
0488               track2.setDuplicateValue(true);
0489             else
0490               trk.setDuplicateValue(true);
0491           }
0492         }  // end sharing hits loop
0493       }    // end trk loop
0494 
0495       tracks.erase(std::remove_if(tracks.begin(), tracks.end(), [](auto track) { return track.getDuplicateValue(); }),
0496                    tracks.end());
0497     }
0498 
0499     void find_duplicates_sharedhits_pixelseed(TrackVec &tracks,
0500                                               const float fraction,
0501                                               const float drth_central,
0502                                               const float drth_obarrel,
0503                                               const float drth_forward) {
0504       const auto ntracks = tracks.size();
0505 
0506       std::vector<float> ctheta(ntracks);
0507       for (auto itrack = 0U; itrack < ntracks; itrack++) {
0508         auto &trk = tracks[itrack];
0509         ctheta[itrack] = 1.f / std::tan(trk.theta());
0510       }
0511 
0512       float phi1, invpt1, dctheta, ctheta1, dphi, dr2;
0513       for (auto itrack = 0U; itrack < ntracks; itrack++) {
0514         auto &trk = tracks[itrack];
0515         phi1 = trk.momPhi();
0516         invpt1 = trk.invpT();
0517         ctheta1 = ctheta[itrack];
0518         for (auto jtrack = itrack + 1; jtrack < ntracks; jtrack++) {
0519           auto &track2 = tracks[jtrack];
0520           if (trk.label() == track2.label())
0521             continue;
0522 
0523           dctheta = std::abs(ctheta[jtrack] - ctheta1);
0524 
0525           if (dctheta > Config::maxdcth)
0526             continue;
0527 
0528           dphi = std::abs(squashPhiMinimal(phi1 - track2.momPhi()));
0529 
0530           if (dphi > Config::maxdphi)
0531             continue;
0532 
0533           float maxdRSquared = drth_central * drth_central;
0534           if (std::abs(ctheta1) > Config::maxcth_fw)
0535             maxdRSquared = drth_forward * drth_forward;
0536           else if (std::abs(ctheta1) > Config::maxcth_ob)
0537             maxdRSquared = drth_obarrel * drth_obarrel;
0538           dr2 = dphi * dphi + dctheta * dctheta;
0539           if (dr2 < maxdRSquared) {
0540             //Keep track with best score
0541             if (trk.score() > track2.score())
0542               track2.setDuplicateValue(true);
0543             else
0544               trk.setDuplicateValue(true);
0545             continue;
0546           }
0547 
0548           if (std::abs(track2.invpT() - invpt1) > Config::maxd1pt)
0549             continue;
0550 
0551           auto sharedCount = 0;
0552           auto sharedFirst = 0;
0553           const auto minFoundHits = std::min(trk.nFoundHits(), track2.nFoundHits());
0554 
0555           for (int i = 0; i < trk.nTotalHits(); ++i) {
0556             if (trk.getHitIdx(i) < 0)
0557               continue;
0558             const int a = trk.getHitLyr(i);
0559             const int b = trk.getHitIdx(i);
0560             for (int j = 0; j < track2.nTotalHits(); ++j) {
0561               if (track2.getHitIdx(j) < 0)
0562                 continue;
0563               const int c = track2.getHitLyr(j);
0564               const int d = track2.getHitIdx(j);
0565 
0566               //this is to count once shared matched hits (may be done more properly...)
0567               if (a == c && b == d)
0568                 sharedCount += 1;
0569               if (j == 0 && i == 0 && a == c && b == d)
0570                 sharedFirst += 1;
0571 
0572               if ((sharedCount - sharedFirst) >= ((minFoundHits - sharedFirst) * fraction))
0573                 continue;
0574             }
0575             if ((sharedCount - sharedFirst) >= ((minFoundHits - sharedFirst) * fraction))
0576               continue;
0577           }
0578 
0579           //selection here - 11percent fraction of shared hits to label a duplicate
0580           if ((sharedCount - sharedFirst) >= ((minFoundHits - sharedFirst) * fraction)) {
0581             if (trk.score() > track2.score())
0582               track2.setDuplicateValue(true);
0583             else
0584               trk.setDuplicateValue(true);
0585           }
0586         }
0587       }  //end loop one over tracks
0588 
0589       //removal here
0590       tracks.erase(std::remove_if(tracks.begin(), tracks.end(), [](auto track) { return track.getDuplicateValue(); }),
0591                    tracks.end());
0592     }
0593 
0594     //=========================================================================
0595     //
0596     //=========================================================================
0597 
0598     void find_and_remove_duplicates(TrackVec &tracks, const IterationConfig &itconf) {
0599 #ifdef DEBUG
0600       std::cout << " find_and_remove_duplicates: input track size " << tracks.size() << std::endl;
0601 #endif
0602       if (itconf.m_requires_quality_filter && !(itconf.m_requires_dupclean_tight)) {
0603         find_duplicates_sharedhits(tracks, itconf.m_params.fracSharedHits);
0604       } else if (itconf.m_requires_dupclean_tight) {
0605         find_duplicates_sharedhits_pixelseed(tracks,
0606                                              itconf.m_params.fracSharedHits,
0607                                              itconf.m_params.drth_central,
0608                                              itconf.m_params.drth_obarrel,
0609                                              itconf.m_params.drth_forward);
0610       } else {
0611         find_duplicates(tracks);
0612         remove_duplicates(tracks);
0613       }
0614 
0615 #ifdef DEBUG
0616       std::cout << " find_and_remove_duplicates: output track size " << tracks.size() << std::endl;
0617       for (auto const &tk : tracks) {
0618         std::cout << tk.parameters() << std::endl;
0619       }
0620 #endif
0621     }
0622 
0623   }  // namespace StdSeq
0624 }  // namespace mkfit