Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-10-17 22:59:00

0001 #include "RecoTracker/MkFitCore/interface/cms_common_macros.h"
0002 #include "RecoTracker/MkFitCMS/interface/MkStdSeqs.h"
0003 
0004 #include "RecoTracker/MkFitCore/interface/HitStructures.h"
0005 #include "RecoTracker/MkFitCore/interface/IterationConfig.h"
0006 #include "RecoTracker/MkFitCore/interface/MkJob.h"
0007 #include "RecoTracker/MkFitCore/interface/TrackStructures.h"
0008 
0009 #include "RecoTracker/MkFitCore/interface/binnor.h"
0010 
0011 namespace mkfit {
0012 
0013   namespace StdSeq {
0014 
0015     //=========================================================================
0016     // Hit processing
0017     //=========================================================================
0018 
0019     void loadDeads(EventOfHits &eoh, const std::vector<DeadVec> &deadvectors) {
0020       for (size_t il = 0; il < deadvectors.size(); il++) {
0021         eoh.suckInDeads(int(il), deadvectors[il]);
0022       }
0023     }
0024 
0025     // Loading hits in CMSSW from two "large multi-layer vectors".
0026     // orig_hitvectors[0] - pixels,
0027     // orig_hitvectors[1] - strips.
0028 
0029     void cmssw_LoadHits_Begin(EventOfHits &eoh, const std::vector<const HitVec *> &orig_hitvectors) {
0030       eoh.reset();
0031       for (int i = 0; i < eoh.nLayers(); ++i) {
0032         auto &&l = eoh[i];
0033         l.beginRegistrationOfHits(*orig_hitvectors[l.is_pixel() ? 0 : 1]);
0034       }
0035     }
0036 
0037     // Loop with LayerOfHits::registerHit(int idx) - it takes Hit out of original HitVec to
0038     // extract phi, r/z, and calculate qphifines
0039     //
0040     // Something like what is done in MkFitInputConverter::convertHits
0041     //
0042     // Problem is I don't know layers for each large-vector;
0043     // Also, layer is calculated for each detset when looping over the HitCollection
0044 
0045     void cmssw_LoadHits_End(EventOfHits &eoh) {
0046       for (int i = 0; i < eoh.nLayers(); ++i) {
0047         auto &&l = eoh[i];
0048         l.endRegistrationOfHits(false);
0049       }
0050     }
0051 
0052     //=========================================================================
0053     // Hit-index mapping / remapping
0054     //=========================================================================
0055 
0056     void cmssw_Map_TrackHitIndices(const EventOfHits &eoh, TrackVec &seeds) {
0057       for (auto &&track : seeds) {
0058         for (int i = 0; i < track.nTotalHits(); ++i) {
0059           const int hitidx = track.getHitIdx(i);
0060           const int hitlyr = track.getHitLyr(i);
0061           if (hitidx >= 0) {
0062             const auto &loh = eoh[hitlyr];
0063             track.setHitIdx(i, loh.getHitIndexFromOriginal(hitidx));
0064           }
0065         }
0066       }
0067     }
0068 
0069     void cmssw_ReMap_TrackHitIndices(const EventOfHits &eoh, TrackVec &out_tracks) {
0070       for (auto &&track : out_tracks) {
0071         for (int i = 0; i < track.nTotalHits(); ++i) {
0072           const int hitidx = track.getHitIdx(i);
0073           const int hitlyr = track.getHitLyr(i);
0074           if (hitidx >= 0) {
0075             const auto &loh = eoh[hitlyr];
0076             track.setHitIdx(i, loh.getOriginalHitIndex(hitidx));
0077           }
0078         }
0079       }
0080     }
0081 
0082     //=========================================================================
0083     // Seed cleaning (multi-iter)
0084     //=========================================================================
0085     int clean_cms_seedtracks_iter(TrackVec &seeds, const IterationConfig &itrcfg, const BeamSpot &bspot) {
0086       using Algo = TrackBase::TrackAlgorithm;
0087 
0088       const float etamax_brl = Config::c_etamax_brl;
0089       const float dpt_common = Config::c_dpt_common;
0090 
0091       const float dzmax_bh = itrcfg.sc_dzmax_bh;
0092       const float drmax_bh = itrcfg.sc_drmax_bh;
0093       const float dzmax_eh = itrcfg.sc_dzmax_eh;
0094       const float drmax_eh = itrcfg.sc_drmax_eh;
0095       const float dzmax_bl = itrcfg.sc_dzmax_bl;
0096       const float drmax_bl = itrcfg.sc_drmax_bl;
0097       const float dzmax_el = itrcfg.sc_dzmax_el;
0098       const float drmax_el = itrcfg.sc_drmax_el;
0099 
0100       const float ptmin_hpt = itrcfg.sc_ptthr_hpt;
0101 
0102       const float dzmax2_inv_bh = 1.f / (dzmax_bh * dzmax_bh);
0103       const float drmax2_inv_bh = 1.f / (drmax_bh * drmax_bh);
0104       const float dzmax2_inv_eh = 1.f / (dzmax_eh * dzmax_eh);
0105       const float drmax2_inv_eh = 1.f / (drmax_eh * drmax_eh);
0106       const float dzmax2_inv_bl = 1.f / (dzmax_bl * dzmax_bl);
0107       const float drmax2_inv_bl = 1.f / (drmax_bl * drmax_bl);
0108       const float dzmax2_inv_el = 1.f / (dzmax_el * dzmax_el);
0109       const float drmax2_inv_el = 1.f / (drmax_el * drmax_el);
0110 
0111       // Merge hits from overlapping seeds?
0112       // For now always true, we require extra hits after seed,
0113       // except for lowPtQuadStep, where we only merge hits for seeds at low pT and large pseudo-rapidity
0114       const bool merge_hits = true;  // itrcfg.merge_seed_hits_during_cleaning();
0115       const float ptmax_merge_lowPtQuad = 0.2;
0116       const float etamin_merge_lowPtQuad = 1.5;
0117 
0118       if (seeds.empty())
0119         return 0;
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     namespace {
0339       CMS_SA_ALLOW struct register_seed_cleaners {
0340         register_seed_cleaners() {
0341           IterationConfig::register_seed_cleaner("phase1:default", clean_cms_seedtracks_iter);
0342         }
0343       } rsc_instance;
0344     }  // namespace
0345 
0346     //=========================================================================
0347     // Duplicate cleaning
0348     //=========================================================================
0349 
0350     void remove_duplicates(TrackVec &tracks) {
0351       tracks.erase(std::remove_if(tracks.begin(), tracks.end(), [](auto track) { return track.getDuplicateValue(); }),
0352                    tracks.end());
0353     }
0354 
0355     void clean_duplicates(TrackVec &tracks, const IterationConfig &) {
0356       const auto ntracks = tracks.size();
0357       float eta1, phi1, pt1, deta, dphi, dr2;
0358 
0359       if (ntracks == 0) {
0360         return;
0361       }
0362       for (auto itrack = 0U; itrack < ntracks - 1; itrack++) {
0363         auto &track = tracks[itrack];
0364         eta1 = track.momEta();
0365         phi1 = track.momPhi();
0366         pt1 = track.pT();
0367         for (auto jtrack = itrack + 1; jtrack < ntracks; jtrack++) {
0368           auto &track2 = tracks[jtrack];
0369           if (track.label() == track2.label())
0370             continue;
0371           if (track.algoint() != track2.algoint())
0372             continue;
0373 
0374           deta = std::abs(track2.momEta() - eta1);
0375           if (deta > Config::maxdEta)
0376             continue;
0377 
0378           dphi = std::abs(squashPhiMinimal(phi1 - track2.momPhi()));
0379           if (dphi > Config::maxdPhi)
0380             continue;
0381 
0382           float maxdR = Config::maxdR;
0383           float maxdRSquared = maxdR * maxdR;
0384           if (std::abs(eta1) > 2.5f)
0385             maxdRSquared *= 16.0f;
0386           else if (std::abs(eta1) > 1.44f)
0387             maxdRSquared *= 9.0f;
0388           dr2 = dphi * dphi + deta * deta;
0389           if (dr2 < maxdRSquared) {
0390             //Keep track with best score
0391             if (track.score() > track2.score())
0392               track2.setDuplicateValue(true);
0393             else
0394               track.setDuplicateValue(true);
0395             continue;
0396           } else {
0397             if (pt1 == 0)
0398               continue;
0399             if (track2.pT() == 0)
0400               continue;
0401 
0402             if (std::abs((1 / track2.pT()) - (1 / pt1)) < Config::maxdPt) {
0403               if (Config::useHitsForDuplicates) {
0404                 float numHitsShared = 0;
0405                 for (int ihit2 = 0; ihit2 < track2.nTotalHits(); ihit2++) {
0406                   const int hitidx2 = track2.getHitIdx(ihit2);
0407                   const int hitlyr2 = track2.getHitLyr(ihit2);
0408                   if (hitidx2 >= 0) {
0409                     auto const it = std::find_if(track.beginHitsOnTrack(),
0410                                                  track.endHitsOnTrack(),
0411                                                  [&hitidx2, &hitlyr2](const HitOnTrack &element) {
0412                                                    return (element.index == hitidx2 && element.layer == hitlyr2);
0413                                                  });
0414                     if (it != track.endHitsOnTrack())
0415                       numHitsShared++;
0416                   }
0417                 }
0418 
0419                 float fracHitsShared = numHitsShared / std::min(track.nFoundHits(), track2.nFoundHits());
0420                 //Only remove one of the tracks if they share at least X% of the hits (denominator is the shorter track)
0421                 if (fracHitsShared < Config::minFracHitsShared)
0422                   continue;
0423               }
0424               //Keep track with best score
0425               if (track.score() > track2.score())
0426                 track2.setDuplicateValue(true);
0427               else
0428                 track.setDuplicateValue(true);
0429             }  //end of if dPt
0430           }  //end of else
0431         }  //end of loop over track2
0432       }  //end of loop over track1
0433 
0434       remove_duplicates(tracks);
0435     }
0436 
0437     //=========================================================================
0438     // SHARED HITS DUPLICATE CLEANING
0439     //=========================================================================
0440 
0441     void clean_duplicates_sharedhits(TrackVec &tracks, const IterationConfig &itconf) {
0442       const float fraction = itconf.dc_fracSharedHits;
0443       const auto ntracks = tracks.size();
0444 
0445       std::vector<float> ctheta(ntracks);
0446       std::multimap<int, int> hitMap;
0447 
0448       for (auto itrack = 0U; itrack < ntracks; itrack++) {
0449         auto &trk = tracks[itrack];
0450         ctheta[itrack] = 1.f / std::tan(trk.theta());
0451         for (int i = 0; i < trk.nTotalHits(); ++i) {
0452           if (trk.getHitIdx(i) < 0)
0453             continue;
0454           int a = trk.getHitLyr(i);
0455           int b = trk.getHitIdx(i);
0456           hitMap.insert(std::make_pair(b * 1000 + a, i > 0 ? itrack : -itrack));  //negative for first hit in trk
0457         }
0458       }
0459 
0460       for (auto itrack = 0U; itrack < ntracks; itrack++) {
0461         auto &trk = tracks[itrack];
0462         auto phi1 = trk.momPhi();
0463         auto ctheta1 = ctheta[itrack];
0464 
0465         std::map<int, int> sharingMap;
0466         for (int i = 0; i < trk.nTotalHits(); ++i) {
0467           if (trk.getHitIdx(i) < 0)
0468             continue;
0469           int a = trk.getHitLyr(i);
0470           int b = trk.getHitIdx(i);
0471           auto range = hitMap.equal_range(b * 1000 + a);
0472           for (auto it = range.first; it != range.second; ++it) {
0473             if (std::abs(it->second) >= (int)itrack)
0474               continue;  // don't check your own hits (==) nor sym. checks (>)
0475             if (i == 0 && it->second < 0)
0476               continue;  // shared first - first is not counted
0477             sharingMap[std::abs(it->second)]++;
0478           }
0479         }
0480 
0481         for (const auto &elem : sharingMap) {
0482           auto &track2 = tracks[elem.first];
0483 
0484           // broad dctheta-dphi compatibility checks; keep mostly to preserve consistency with old results
0485           auto dctheta = std::abs(ctheta[elem.first] - ctheta1);
0486           if (dctheta > 1.)
0487             continue;
0488 
0489           auto dphi = std::abs(squashPhiMinimal(phi1 - track2.momPhi()));
0490           if (dphi > 1.)
0491             continue;
0492 
0493           if (elem.second >= std::min(trk.nFoundHits(), track2.nFoundHits()) * fraction) {
0494             if (trk.score() > track2.score())
0495               track2.setDuplicateValue(true);
0496             else
0497               trk.setDuplicateValue(true);
0498           }
0499         }  // end sharing hits loop
0500       }  // end trk loop
0501 
0502       remove_duplicates(tracks);
0503     }
0504 
0505     void clean_duplicates_sharedhits_pixelseed(TrackVec &tracks, const IterationConfig &itconf) {
0506       const float fraction = itconf.dc_fracSharedHits;
0507       const float drth_central = itconf.dc_drth_central;
0508       const float drth_obarrel = itconf.dc_drth_obarrel;
0509       const float drth_forward = itconf.dc_drth_forward;
0510       const auto ntracks = tracks.size();
0511 
0512       std::vector<float> ctheta(ntracks);
0513       for (auto itrack = 0U; itrack < ntracks; itrack++) {
0514         auto &trk = tracks[itrack];
0515         ctheta[itrack] = 1.f / std::tan(trk.theta());
0516       }
0517 
0518       float phi1, invpt1, dctheta, ctheta1, dphi, dr2;
0519       for (auto itrack = 0U; itrack < ntracks; itrack++) {
0520         auto &trk = tracks[itrack];
0521         phi1 = trk.momPhi();
0522         invpt1 = trk.invpT();
0523         ctheta1 = ctheta[itrack];
0524         for (auto jtrack = itrack + 1; jtrack < ntracks; jtrack++) {
0525           auto &track2 = tracks[jtrack];
0526           if (trk.label() == track2.label())
0527             continue;
0528 
0529           dctheta = std::abs(ctheta[jtrack] - ctheta1);
0530 
0531           if (dctheta > Config::maxdcth)
0532             continue;
0533 
0534           dphi = std::abs(squashPhiMinimal(phi1 - track2.momPhi()));
0535 
0536           if (dphi > Config::maxdphi)
0537             continue;
0538 
0539           float maxdRSquared = drth_central * drth_central;
0540           if (std::abs(ctheta1) > Config::maxcth_fw)
0541             maxdRSquared = drth_forward * drth_forward;
0542           else if (std::abs(ctheta1) > Config::maxcth_ob)
0543             maxdRSquared = drth_obarrel * drth_obarrel;
0544           dr2 = dphi * dphi + dctheta * dctheta;
0545           if (dr2 < maxdRSquared) {
0546             //Keep track with best score
0547             if (trk.score() > track2.score())
0548               track2.setDuplicateValue(true);
0549             else
0550               trk.setDuplicateValue(true);
0551             continue;
0552           }
0553 
0554           if (std::abs(track2.invpT() - invpt1) > Config::maxd1pt)
0555             continue;
0556 
0557           auto sharedCount = 0;
0558           auto sharedFirst = 0;
0559           const auto minFoundHits = std::min(trk.nFoundHits(), track2.nFoundHits());
0560 
0561           for (int i = 0; i < trk.nTotalHits(); ++i) {
0562             if (trk.getHitIdx(i) < 0)
0563               continue;
0564             const int a = trk.getHitLyr(i);
0565             const int b = trk.getHitIdx(i);
0566             for (int j = 0; j < track2.nTotalHits(); ++j) {
0567               if (track2.getHitIdx(j) < 0)
0568                 continue;
0569               const int c = track2.getHitLyr(j);
0570               const int d = track2.getHitIdx(j);
0571 
0572               //this is to count once shared matched hits (may be done more properly...)
0573               if (a == c && b == d)
0574                 sharedCount += 1;
0575               if (j == 0 && i == 0 && a == c && b == d)
0576                 sharedFirst += 1;
0577 
0578               if ((sharedCount - sharedFirst) >= ((minFoundHits - sharedFirst) * fraction))
0579                 continue;
0580             }
0581             if ((sharedCount - sharedFirst) >= ((minFoundHits - sharedFirst) * fraction))
0582               continue;
0583           }
0584 
0585           //selection here - 11percent fraction of shared hits to label a duplicate
0586           if ((sharedCount - sharedFirst) >= ((minFoundHits - sharedFirst) * fraction)) {
0587             if (trk.score() > track2.score())
0588               track2.setDuplicateValue(true);
0589             else
0590               trk.setDuplicateValue(true);
0591           }
0592         }
0593       }  //end loop one over tracks
0594 
0595       remove_duplicates(tracks);
0596     }
0597 
0598     namespace {
0599       CMS_SA_ALLOW struct register_duplicate_cleaners {
0600         register_duplicate_cleaners() {
0601           IterationConfig::register_duplicate_cleaner("phase1:clean_duplicates", clean_duplicates);
0602           IterationConfig::register_duplicate_cleaner("phase1:clean_duplicates_sharedhits",
0603                                                       clean_duplicates_sharedhits);
0604           IterationConfig::register_duplicate_cleaner("phase1:clean_duplicates_sharedhits_pixelseed",
0605                                                       clean_duplicates_sharedhits_pixelseed);
0606         }
0607       } rdc_instance;
0608     }  // namespace
0609 
0610     //=========================================================================
0611     // Quality filters
0612     //=========================================================================
0613 
0614     // quality filter for n hits with seed hit "penalty" for strip-based seeds
0615     //   this implicitly separates triplets and doublet seeds with glued layers
0616     template <class TRACK>
0617     bool qfilter_n_hits(const TRACK &t, const MkJob &j) {
0618       int seedHits = t.getNSeedHits();
0619       int seedReduction = (seedHits <= 5) ? 2 : 3;
0620       return t.nFoundHits() - seedReduction >= j.params_cur().minHitsQF;
0621     }
0622 
0623     // simple hit-count quality filter; used with pixel-based seeds
0624     template <class TRACK>
0625     bool qfilter_n_hits_pixseed(const TRACK &t, const MkJob &j) {
0626       return t.nFoundHits() >= j.params_cur().minHitsQF;
0627     }
0628 
0629     // layer-dependent quality filter
0630     // includes ad hoc tuning for phase-1
0631     template <class TRACK>
0632     bool qfilter_n_layers(const TRACK &t, const MkJob &j) {
0633       const BeamSpot &bspot = j.m_beam_spot;
0634       const TrackerInfo &trk_inf = j.m_trk_info;
0635       int enhits = t.nHitsByTypeEncoded(trk_inf);
0636       int npixhits = t.nPixelDecoded(enhits);
0637       int enlyrs = t.nLayersByTypeEncoded(trk_inf);
0638       int npixlyrs = t.nPixelDecoded(enlyrs);
0639       int nmatlyrs = t.nTotMatchDecoded(enlyrs);
0640       int llyr = t.getLastFoundHitLyr();
0641       int lplyr = t.getLastFoundPixelHitLyr();
0642       float invpt = t.invpT();
0643 
0644       // based on fr and eff vs pt (convert to native invpt)
0645       float invptmin = 1.43;  // min 1/pT (=1/0.7) for full filter on (npixhits<=3 .or. npixlyrs<=3)
0646       float d0BS = t.d0BeamSpot(bspot.x, bspot.y);
0647       float d0_max = 0.1;  // 1 mm, max for somewhat prompt
0648 
0649       // next-to-outermost pixel layers (almost): BPIX3 or FPIX1
0650       bool endsInsidePix = (llyr == 2 || llyr == 18 || llyr == 45);
0651       // not last pixel layers: BPIX[123] or FPIX[12]
0652       bool lastInsidePix = ((0 <= lplyr && lplyr < 3) || (18 <= lplyr && lplyr < 20) || (45 <= lplyr && lplyr < 47));
0653       // reject short tracks missing last pixel layer except for prompt-looking
0654       return !(((npixhits <= 3 || npixlyrs <= 3) && endsInsidePix &&
0655                 (invpt < invptmin || (invpt >= invptmin && std::abs(d0BS) > d0_max))) ||
0656                ((npixlyrs <= 3 && nmatlyrs <= 6) && lastInsidePix && llyr != lplyr && std::abs(d0BS) > d0_max));
0657     }
0658 
0659     /// quality filter tuned for pixelLess iteration during forward search
0660     // includes ad hoc tuning for phase-1
0661     template <class TRACK>
0662     bool qfilter_pixelLessFwd(const TRACK &t, const MkJob &j) {
0663       const BeamSpot &bspot = j.m_beam_spot;
0664       const TrackerInfo &tk_info = j.m_trk_info;
0665       float d0BS = t.d0BeamSpot(bspot.x, bspot.y);
0666       float d0_max = 0.05;  // 0.5 mm, max for somewhat prompt
0667 
0668       int encoded;
0669       encoded = t.nLayersByTypeEncoded(tk_info);
0670       int nLyrs = t.nTotMatchDecoded(encoded);
0671       encoded = t.nHitsByTypeEncoded(tk_info);
0672       int nHits = t.nTotMatchDecoded(encoded);
0673 
0674       // to subtract stereo seed layers to count just r-phi seed layers (better pt err)
0675       int seedReduction = (t.getNSeedHits() <= 5) ? 2 : 3;
0676 
0677       // based on fr and eff vs pt and eta (convert to native invpt and theta)
0678       float invpt = t.invpT();
0679       float invptmin = 1.11;  // =1/0.9
0680 
0681       float thetasym = std::abs(t.theta() - Const::PIOver2);
0682       float thetasymmin = 1.11;  // -> |eta|=1.45
0683 
0684       // accept longer tracks, reject too short and displaced
0685       return (((t.nFoundHits() - seedReduction >= 4 && invpt < invptmin) ||
0686                (t.nFoundHits() - seedReduction >= 3 && invpt > invptmin && thetasym <= thetasymmin) ||
0687                (t.nFoundHits() - seedReduction >= 4 && invpt > invptmin && thetasym > thetasymmin)) &&
0688               !((nLyrs <= 4 || nHits <= 4) && std::abs(d0BS) > d0_max && invpt < invptmin));
0689     }
0690 
0691     /// quality filter tuned for pixelLess iteration during backward search
0692     // includes ad hoc tuning for phase-1
0693     template <class TRACK>
0694     bool qfilter_pixelLessBkwd(const TRACK &t, const MkJob &j) {
0695       const BeamSpot &bspot = j.m_beam_spot;
0696       const TrackerInfo &tk_info = j.m_trk_info;
0697       float d0BS = t.d0BeamSpot(bspot.x, bspot.y);
0698       float d0_max = 0.1;  // 1 mm
0699 
0700       int encoded;
0701       encoded = t.nLayersByTypeEncoded(tk_info);
0702       int nLyrs = t.nTotMatchDecoded(encoded);
0703       encoded = t.nHitsByTypeEncoded(tk_info);
0704       int nHits = t.nTotMatchDecoded(encoded);
0705 
0706       // based on fr and eff vs pt and eta (convert to native invpt and theta)
0707       float invpt = t.invpT();
0708       float invptmin = 1.11;  // =1/0.9
0709 
0710       float thetasym = std::abs(t.theta() - Const::PIOver2);
0711       float thetasymmin_l = 0.80;  // -> |eta|=0.9
0712       float thetasymmin_h = 1.11;  // -> |eta|=1.45
0713 
0714       // reject too short or too displaced tracks
0715       return !(
0716           ((nLyrs <= 3 || nHits <= 3)) ||
0717           ((nLyrs <= 4 || nHits <= 4) && (invpt < invptmin || (thetasym > thetasymmin_l && std::abs(d0BS) > d0_max))) ||
0718           ((nLyrs <= 5 || nHits <= 5) && (invpt > invptmin && thetasym > thetasymmin_h && std::abs(d0BS) > d0_max)));
0719     }
0720 
0721     namespace {
0722       CMS_SA_ALLOW struct register_quality_filters {
0723         register_quality_filters() {
0724           IterationConfig::register_candidate_filter("phase1:qfilter_n_hits", qfilter_n_hits<TrackCand>);
0725           IterationConfig::register_candidate_filter("phase1:qfilter_n_hits_pixseed",
0726                                                      qfilter_n_hits_pixseed<TrackCand>);
0727           IterationConfig::register_candidate_filter("phase1:qfilter_n_layers", qfilter_n_layers<TrackCand>);
0728           IterationConfig::register_candidate_filter("phase1:qfilter_pixelLessFwd", qfilter_pixelLessFwd<TrackCand>);
0729           IterationConfig::register_candidate_filter("phase1:qfilter_pixelLessBkwd", qfilter_pixelLessBkwd<TrackCand>);
0730         }
0731       } rqf_instance;
0732     }  // namespace
0733 
0734     //=========================================================================
0735     // Track scoring
0736     //=========================================================================
0737 
0738     float trackScoreDefault(const int nfoundhits,
0739                             const int ntailholes,
0740                             const int noverlaphits,
0741                             const int nmisshits,
0742                             const float chi2,
0743                             const float pt,
0744                             const bool inFindCandidates) {
0745       float maxBonus = 8.0;
0746       float bonus = Config::validHitSlope_ * nfoundhits + Config::validHitBonus_;
0747       float penalty = Config::missingHitPenalty_;
0748       float tailPenalty = Config::tailMissingHitPenalty_;
0749       float overlapBonus = Config::overlapHitBonus_;
0750       if (pt < 0.9) {
0751         penalty *= inFindCandidates ? 1.7f : 1.5f;
0752         bonus = std::min(bonus * (inFindCandidates ? 0.9f : 1.0f), maxBonus);
0753       }
0754       float score =
0755           bonus * nfoundhits + overlapBonus * noverlaphits - penalty * nmisshits - tailPenalty * ntailholes - chi2;
0756       return score;
0757     }
0758 
0759     namespace {
0760       CMS_SA_ALLOW struct register_track_scorers {
0761         register_track_scorers() {
0762           IterationConfig::register_track_scorer("default", trackScoreDefault);
0763           IterationConfig::register_track_scorer("phase1:default", trackScoreDefault);
0764         }
0765       } rts_instance;
0766     }  // namespace
0767 
0768   }  // namespace StdSeq
0769 }  // namespace mkfit