Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-10-25 10:02:29

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 #include "oneapi/tbb/parallel_for.h"
0012 
0013 namespace mkfit {
0014 
0015   namespace StdSeq {
0016 
0017     //=========================================================================
0018     // Hit processing
0019     //=========================================================================
0020 
0021     void loadDeads(EventOfHits &eoh, const std::vector<DeadVec> &deadvectors) {
0022       for (size_t il = 0; il < deadvectors.size(); il++) {
0023         eoh.suckInDeads(int(il), deadvectors[il]);
0024       }
0025     }
0026 
0027     // Loading hits in CMSSW from two "large multi-layer vectors".
0028     // orig_hitvectors[0] - pixels,
0029     // orig_hitvectors[1] - strips.
0030 
0031     void cmssw_LoadHits_Begin(EventOfHits &eoh, const std::vector<const HitVec *> &orig_hitvectors) {
0032       eoh.reset();
0033       for (int i = 0; i < eoh.nLayers(); ++i) {
0034         auto &&l = eoh[i];
0035         l.beginRegistrationOfHits(*orig_hitvectors[l.is_pixel() ? 0 : 1]);
0036       }
0037     }
0038 
0039     // Loop with LayerOfHits::registerHit(int idx) - it takes Hit out of original HitVec to
0040     // extract phi, r/z, and calculate qphifines
0041     //
0042     // Something like what is done in MkFitInputConverter::convertHits
0043     //
0044     // Problem is I don't know layers for each large-vector;
0045     // Also, layer is calculated for each detset when looping over the HitCollection
0046 
0047     void cmssw_LoadHits_End(EventOfHits &eoh) {
0048       for (int i = 0; i < eoh.nLayers(); ++i) {
0049         auto &&l = eoh[i];
0050         l.endRegistrationOfHits(false);
0051       }
0052     }
0053 
0054     //=========================================================================
0055     // Hit-index mapping / remapping
0056     //=========================================================================
0057 
0058     void cmssw_Map_TrackHitIndices(const EventOfHits &eoh, TrackVec &seeds) {
0059       for (auto &&track : seeds) {
0060         for (int i = 0; i < track.nTotalHits(); ++i) {
0061           const int hitidx = track.getHitIdx(i);
0062           const int hitlyr = track.getHitLyr(i);
0063           if (hitidx >= 0) {
0064             const auto &loh = eoh[hitlyr];
0065             track.setHitIdx(i, loh.getHitIndexFromOriginal(hitidx));
0066           }
0067         }
0068       }
0069     }
0070 
0071     void cmssw_ReMap_TrackHitIndices(const EventOfHits &eoh, TrackVec &out_tracks) {
0072       for (auto &&track : out_tracks) {
0073         for (int i = 0; i < track.nTotalHits(); ++i) {
0074           const int hitidx = track.getHitIdx(i);
0075           const int hitlyr = track.getHitLyr(i);
0076           if (hitidx >= 0) {
0077             const auto &loh = eoh[hitlyr];
0078             track.setHitIdx(i, loh.getOriginalHitIndex(hitidx));
0079           }
0080         }
0081       }
0082     }
0083 
0084     //=========================================================================
0085     // Seed cleaning (multi-iter)
0086     //=========================================================================
0087     int clean_cms_seedtracks_iter(TrackVec &seeds, const IterationConfig &itrcfg, const BeamSpot &bspot) {
0088       using Algo = TrackBase::TrackAlgorithm;
0089 
0090       const float etamax_brl = Config::c_etamax_brl;
0091       const float dpt_common = Config::c_dpt_common;
0092 
0093       const float dzmax_bh = itrcfg.sc_dzmax_bh;
0094       const float drmax_bh = itrcfg.sc_drmax_bh;
0095       const float dzmax_eh = itrcfg.sc_dzmax_eh;
0096       const float drmax_eh = itrcfg.sc_drmax_eh;
0097       const float dzmax_bl = itrcfg.sc_dzmax_bl;
0098       const float drmax_bl = itrcfg.sc_drmax_bl;
0099       const float dzmax_el = itrcfg.sc_dzmax_el;
0100       const float drmax_el = itrcfg.sc_drmax_el;
0101 
0102       const float ptmin_hpt = itrcfg.sc_ptthr_hpt;
0103 
0104       const float dzmax2_inv_bh = 1.f / (dzmax_bh * dzmax_bh);
0105       const float drmax2_inv_bh = 1.f / (drmax_bh * drmax_bh);
0106       const float dzmax2_inv_eh = 1.f / (dzmax_eh * dzmax_eh);
0107       const float drmax2_inv_eh = 1.f / (drmax_eh * drmax_eh);
0108       const float dzmax2_inv_bl = 1.f / (dzmax_bl * dzmax_bl);
0109       const float drmax2_inv_bl = 1.f / (drmax_bl * drmax_bl);
0110       const float dzmax2_inv_el = 1.f / (dzmax_el * dzmax_el);
0111       const float drmax2_inv_el = 1.f / (drmax_el * drmax_el);
0112 
0113       // Merge hits from overlapping seeds?
0114       // For now always true, we require extra hits after seed,
0115       // except for lowPtQuadStep, where we only merge hits for seeds at low pT and large pseudo-rapidity
0116       const bool merge_hits = true;  // itrcfg.merge_seed_hits_during_cleaning();
0117       const float ptmax_merge_lowPtQuad = 0.2;
0118       const float etamin_merge_lowPtQuad = 1.5;
0119 
0120       if (seeds.empty())
0121         return 0;
0122 
0123       const int ns = seeds.size();
0124 #ifdef DEBUG
0125       std::cout << "before seed cleaning " << seeds.size() << std::endl;
0126 #endif
0127       TrackVec cleanSeedTracks;
0128       cleanSeedTracks.reserve(ns);
0129       std::vector<bool> writetrack(ns, true);
0130 
0131       const float invR1GeV = 1.f / Config::track1GeVradius;
0132 
0133       std::vector<int> nHits(ns);
0134       std::vector<int> charge(ns);
0135       std::vector<float> oldPhi(ns);
0136       std::vector<float> pos2(ns);
0137       std::vector<float> eta(ns);
0138       std::vector<float> ctheta(ns);
0139       std::vector<float> invptq(ns);
0140       std::vector<float> pt(ns);
0141       std::vector<float> x(ns);
0142       std::vector<float> y(ns);
0143       std::vector<float> z(ns);
0144       std::vector<float> d0(ns);
0145       int i1, i2;  //for the sorting
0146 
0147       axis_pow2_u1<float, unsigned short, 16, 8> ax_phi(-Const::PI, Const::PI);
0148       axis<float, unsigned short, 8, 8> ax_eta(-3.0, 3.0, 30u);
0149       binnor<unsigned int, decltype(ax_phi), decltype(ax_eta), 24, 8> phi_eta_binnor(ax_phi, ax_eta);
0150 
0151       phi_eta_binnor.begin_registration(ns);
0152 
0153       for (int ts = 0; ts < ns; ts++) {
0154         const Track &tk = seeds[ts];
0155         nHits[ts] = tk.nFoundHits();
0156         charge[ts] = tk.charge();
0157         oldPhi[ts] = tk.momPhi();
0158         pos2[ts] = std::pow(tk.x(), 2) + std::pow(tk.y(), 2);
0159         eta[ts] = tk.momEta();
0160         ctheta[ts] = 1.f / std::tan(tk.theta());
0161         invptq[ts] = tk.charge() * tk.invpT();
0162         pt[ts] = tk.pT();
0163         x[ts] = tk.x();
0164         y[ts] = tk.y();
0165         z[ts] = tk.z();
0166         d0[ts] = tk.d0BeamSpot(bspot.x, bspot.y);
0167 
0168         phi_eta_binnor.register_entry_safe(oldPhi[ts], eta[ts]);
0169         // If one is sure values are *within* axis ranges: b.register_entry(oldPhi[ts], eta[ts]);
0170       }
0171 
0172       phi_eta_binnor.finalize_registration();
0173 
0174       for (int sorted_ts = 0; sorted_ts < ns; sorted_ts++) {
0175         int ts = phi_eta_binnor.m_ranks[sorted_ts];
0176 
0177         if (not writetrack[ts])
0178           continue;  // Note: this speed up prevents transitive masking (possibly marginal gain).
0179 
0180         const float oldPhi1 = oldPhi[ts];
0181         const float pos2_first = pos2[ts];
0182         const float eta1 = eta[ts];
0183         const float pt1 = pt[ts];
0184         const float invptq_first = invptq[ts];
0185 
0186         // To study some more details -- need EventOfHits for this
0187         int n_ovlp_hits_added = 0;
0188 
0189         auto phi_rng = ax_phi.from_R_rdr_to_N_bins(oldPhi[ts], 0.08f);
0190         auto eta_rng = ax_eta.from_R_rdr_to_N_bins(eta[ts], .1f);
0191 
0192         for (auto i_phi = phi_rng.begin; i_phi != phi_rng.end; i_phi = ax_phi.next_N_bin(i_phi)) {
0193           for (auto i_eta = eta_rng.begin; i_eta != eta_rng.end; i_eta = ax_eta.next_N_bin(i_eta)) {
0194             const auto cbin = phi_eta_binnor.get_content(i_phi, i_eta);
0195             for (auto i = cbin.first; i < cbin.end(); ++i) {
0196               int tss = phi_eta_binnor.m_ranks[i];
0197 
0198               if (not writetrack[ts])
0199                 continue;
0200               if (not writetrack[tss])
0201                 continue;
0202               if (tss == ts)
0203                 continue;
0204 
0205               const float pt2 = pt[tss];
0206 
0207               // Always require charge consistency. If different charge is assigned, do not remove seed-track
0208               if (charge[tss] != charge[ts])
0209                 continue;
0210 
0211               const float thisDPt = std::abs(pt2 - pt1);
0212               // Require pT consistency between seeds. If dpT is large, do not remove seed-track.
0213               if (thisDPt > dpt_common * pt1)
0214                 continue;
0215 
0216               const float eta2 = eta[tss];
0217               const float deta2 = std::pow(eta1 - eta2, 2);
0218 
0219               const float oldPhi2 = oldPhi[tss];
0220 
0221               const float pos2_second = pos2[tss];
0222               const float thisDXYSign05 = pos2_second > pos2_first ? -0.5f : 0.5f;
0223 
0224               const float thisDXY = thisDXYSign05 * sqrt(std::pow(x[ts] - x[tss], 2) + std::pow(y[ts] - y[tss], 2));
0225 
0226               const float invptq_second = invptq[tss];
0227 
0228               const float newPhi1 = oldPhi1 - thisDXY * invR1GeV * invptq_first;
0229               const float newPhi2 = oldPhi2 + thisDXY * invR1GeV * invptq_second;
0230 
0231               const float dphi = cdist(std::abs(newPhi1 - newPhi2));
0232 
0233               const float dr2 = deta2 + dphi * dphi;
0234 
0235               const float thisDZ = z[ts] - z[tss] - thisDXY * (ctheta[ts] + ctheta[tss]);
0236               const float dz2 = thisDZ * thisDZ;
0237 
0238               // Reject tracks within dR-dz elliptical window.
0239               // Adaptive thresholds, based on observation that duplicates are more abundant at large pseudo-rapidity and low track pT
0240               bool overlapping = false;
0241               if (std::abs(eta1) < etamax_brl) {
0242                 if (pt1 > ptmin_hpt) {
0243                   if (dz2 * dzmax2_inv_bh + dr2 * drmax2_inv_bh < 1.0f)
0244                     overlapping = true;
0245                 } else {
0246                   if (dz2 * dzmax2_inv_bl + dr2 * drmax2_inv_bl < 1.0f)
0247                     overlapping = true;
0248                 }
0249               } else {
0250                 if (pt1 > ptmin_hpt) {
0251                   if (dz2 * dzmax2_inv_eh + dr2 * drmax2_inv_eh < 1.0f)
0252                     overlapping = true;
0253                 } else {
0254                   if (dz2 * dzmax2_inv_el + dr2 * drmax2_inv_el < 1.0f)
0255                     overlapping = true;
0256                 }
0257               }
0258 
0259               if (overlapping) {
0260                 //Mark tss as a duplicate
0261                 i1 = ts;
0262                 i2 = tss;
0263                 if (d0[tss] > d0[ts])
0264                   writetrack[tss] = false;
0265                 else {
0266                   writetrack[ts] = false;
0267                   i2 = ts;
0268                   i1 = tss;
0269                 }
0270                 // Add hits from tk2 to the seed we are keeping.
0271                 // NOTE: We have a limit in Track::Status for the number of seed hits.
0272                 //       There is a check at entry and after adding of a new hit.
0273                 Track &tk = seeds[i1];
0274                 if (merge_hits && tk.nTotalHits() < Track::Status::kMaxSeedHits &&
0275                     (Algo(itrcfg.m_track_algorithm) != Algo::lowPtQuadStep ||
0276                      (pt1 < ptmax_merge_lowPtQuad && std::abs(eta1) > etamin_merge_lowPtQuad))) {
0277                   const Track &tk2 = seeds[i2];
0278                   //We are not actually fitting to the extra hits; use chi2 of 0
0279                   float fakeChi2 = 0.0;
0280 
0281                   for (int j = 0; j < tk2.nTotalHits(); ++j) {
0282                     int hitidx = tk2.getHitIdx(j);
0283                     int hitlyr = tk2.getHitLyr(j);
0284                     if (hitidx >= 0) {
0285                       bool unique = true;
0286                       for (int i = 0; i < tk.nTotalHits(); ++i) {
0287                         if ((hitidx == tk.getHitIdx(i)) && (hitlyr == tk.getHitLyr(i))) {
0288                           unique = false;
0289                           break;
0290                         }
0291                       }
0292                       if (unique) {
0293                         tk.addHitIdx(tk2.getHitIdx(j), tk2.getHitLyr(j), fakeChi2);
0294                         ++n_ovlp_hits_added;
0295                         if (tk.nTotalHits() >= Track::Status::kMaxSeedHits)
0296                           break;
0297                       }
0298                     }
0299                   }
0300                 }
0301                 if (n_ovlp_hits_added > 0)
0302                   tk.sortHitsByLayer();
0303               }
0304             }  //end of inner loop over tss
0305           }    //eta bin
0306         }      //phi bin
0307 
0308         if (writetrack[ts]) {
0309           cleanSeedTracks.emplace_back(seeds[ts]);
0310         }
0311       }
0312 
0313       seeds.swap(cleanSeedTracks);
0314 
0315 #ifdef DEBUG
0316       {
0317         const int ns2 = seeds.size();
0318         printf("Number of CMS seeds before %d --> after %d cleaning\n", ns, ns2);
0319 
0320         for (int it = 0; it < ns2; it++) {
0321           const Track &ss = seeds[it];
0322           printf("  %3i q=%+i pT=%7.3f eta=% 7.3f nHits=%i label=% i\n",
0323                  it,
0324                  ss.charge(),
0325                  ss.pT(),
0326                  ss.momEta(),
0327                  ss.nFoundHits(),
0328                  ss.label());
0329         }
0330       }
0331 #endif
0332 
0333 #ifdef DEBUG
0334       std::cout << "AFTER seed cleaning " << seeds.size() << std::endl;
0335 #endif
0336 
0337       return seeds.size();
0338     }
0339 
0340     namespace {
0341       CMS_SA_ALLOW struct register_seed_cleaners {
0342         register_seed_cleaners() {
0343           IterationConfig::register_seed_cleaner("phase1:default", clean_cms_seedtracks_iter);
0344         }
0345       } rsc_instance;
0346     }  // namespace
0347 
0348     //=========================================================================
0349     // Duplicate cleaning
0350     //=========================================================================
0351 
0352     void remove_duplicates(TrackVec &tracks) {
0353       tracks.erase(std::remove_if(tracks.begin(), tracks.end(), [](auto track) { return track.getDuplicateValue(); }),
0354                    tracks.end());
0355     }
0356 
0357     void clean_duplicates(TrackVec &tracks, const IterationConfig &) {
0358       const auto ntracks = tracks.size();
0359       float eta1, phi1, pt1, deta, dphi, dr2;
0360 
0361       if (ntracks == 0) {
0362         return;
0363       }
0364       for (auto itrack = 0U; itrack < ntracks - 1; itrack++) {
0365         auto &track = tracks[itrack];
0366         eta1 = track.momEta();
0367         phi1 = track.momPhi();
0368         pt1 = track.pT();
0369         for (auto jtrack = itrack + 1; jtrack < ntracks; jtrack++) {
0370           auto &track2 = tracks[jtrack];
0371           if (track.label() == track2.label())
0372             continue;
0373           if (track.algoint() != track2.algoint())
0374             continue;
0375 
0376           deta = std::abs(track2.momEta() - eta1);
0377           if (deta > Config::maxdEta)
0378             continue;
0379 
0380           dphi = std::abs(squashPhiMinimal(phi1 - track2.momPhi()));
0381           if (dphi > Config::maxdPhi)
0382             continue;
0383 
0384           float maxdR = Config::maxdR;
0385           float maxdRSquared = maxdR * maxdR;
0386           if (std::abs(eta1) > 2.5f)
0387             maxdRSquared *= 16.0f;
0388           else if (std::abs(eta1) > 1.44f)
0389             maxdRSquared *= 9.0f;
0390           dr2 = dphi * dphi + deta * deta;
0391           if (dr2 < maxdRSquared) {
0392             //Keep track with best score
0393             if (track.score() > track2.score())
0394               track2.setDuplicateValue(true);
0395             else
0396               track.setDuplicateValue(true);
0397             continue;
0398           } else {
0399             if (pt1 == 0)
0400               continue;
0401             if (track2.pT() == 0)
0402               continue;
0403 
0404             if (std::abs((1 / track2.pT()) - (1 / pt1)) < Config::maxdPt) {
0405               if (Config::useHitsForDuplicates) {
0406                 float numHitsShared = 0;
0407                 for (int ihit2 = 0; ihit2 < track2.nTotalHits(); ihit2++) {
0408                   const int hitidx2 = track2.getHitIdx(ihit2);
0409                   const int hitlyr2 = track2.getHitLyr(ihit2);
0410                   if (hitidx2 >= 0) {
0411                     auto const it = std::find_if(track.beginHitsOnTrack(),
0412                                                  track.endHitsOnTrack(),
0413                                                  [&hitidx2, &hitlyr2](const HitOnTrack &element) {
0414                                                    return (element.index == hitidx2 && element.layer == hitlyr2);
0415                                                  });
0416                     if (it != track.endHitsOnTrack())
0417                       numHitsShared++;
0418                   }
0419                 }
0420 
0421                 float fracHitsShared = numHitsShared / std::min(track.nFoundHits(), track2.nFoundHits());
0422                 //Only remove one of the tracks if they share at least X% of the hits (denominator is the shorter track)
0423                 if (fracHitsShared < Config::minFracHitsShared)
0424                   continue;
0425               }
0426               //Keep track with best score
0427               if (track.score() > track2.score())
0428                 track2.setDuplicateValue(true);
0429               else
0430                 track.setDuplicateValue(true);
0431             }  //end of if dPt
0432           }    //end of else
0433         }      //end of loop over track2
0434       }        //end of loop over track1
0435 
0436       remove_duplicates(tracks);
0437     }
0438 
0439     //=========================================================================
0440     // SHARED HITS DUPLICATE CLEANING
0441     //=========================================================================
0442 
0443     void clean_duplicates_sharedhits(TrackVec &tracks, const IterationConfig &itconf) {
0444       const float fraction = itconf.dc_fracSharedHits;
0445       const auto ntracks = tracks.size();
0446 
0447       std::vector<float> ctheta(ntracks);
0448       std::multimap<int, int> hitMap;
0449 
0450       for (auto itrack = 0U; itrack < ntracks; itrack++) {
0451         auto &trk = tracks[itrack];
0452         ctheta[itrack] = 1.f / std::tan(trk.theta());
0453         for (int i = 0; i < trk.nTotalHits(); ++i) {
0454           if (trk.getHitIdx(i) < 0)
0455             continue;
0456           int a = trk.getHitLyr(i);
0457           int b = trk.getHitIdx(i);
0458           hitMap.insert(std::make_pair(b * 1000 + a, i > 0 ? itrack : -itrack));  //negative for first hit in trk
0459         }
0460       }
0461 
0462       for (auto itrack = 0U; itrack < ntracks; itrack++) {
0463         auto &trk = tracks[itrack];
0464         auto phi1 = trk.momPhi();
0465         auto ctheta1 = ctheta[itrack];
0466 
0467         std::map<int, int> sharingMap;
0468         for (int i = 0; i < trk.nTotalHits(); ++i) {
0469           if (trk.getHitIdx(i) < 0)
0470             continue;
0471           int a = trk.getHitLyr(i);
0472           int b = trk.getHitIdx(i);
0473           auto range = hitMap.equal_range(b * 1000 + a);
0474           for (auto it = range.first; it != range.second; ++it) {
0475             if (std::abs(it->second) >= (int)itrack)
0476               continue;  // don't check your own hits (==) nor sym. checks (>)
0477             if (i == 0 && it->second < 0)
0478               continue;  // shared first - first is not counted
0479             sharingMap[std::abs(it->second)]++;
0480           }
0481         }
0482 
0483         for (const auto &elem : sharingMap) {
0484           auto &track2 = tracks[elem.first];
0485 
0486           // broad dctheta-dphi compatibility checks; keep mostly to preserve consistency with old results
0487           auto dctheta = std::abs(ctheta[elem.first] - ctheta1);
0488           if (dctheta > 1.)
0489             continue;
0490 
0491           auto dphi = std::abs(squashPhiMinimal(phi1 - track2.momPhi()));
0492           if (dphi > 1.)
0493             continue;
0494 
0495           if (elem.second >= std::min(trk.nFoundHits(), track2.nFoundHits()) * fraction) {
0496             if (trk.score() > track2.score())
0497               track2.setDuplicateValue(true);
0498             else
0499               trk.setDuplicateValue(true);
0500           }
0501         }  // end sharing hits loop
0502       }    // end trk loop
0503 
0504       remove_duplicates(tracks);
0505     }
0506 
0507     void clean_duplicates_sharedhits_pixelseed(TrackVec &tracks, const IterationConfig &itconf) {
0508       const float fraction = itconf.dc_fracSharedHits;
0509       const float drth_central = itconf.dc_drth_central;
0510       const float drth_obarrel = itconf.dc_drth_obarrel;
0511       const float drth_forward = itconf.dc_drth_forward;
0512       const auto ntracks = tracks.size();
0513 
0514       std::vector<float> ctheta(ntracks);
0515       for (auto itrack = 0U; itrack < ntracks; itrack++) {
0516         auto &trk = tracks[itrack];
0517         ctheta[itrack] = 1.f / std::tan(trk.theta());
0518       }
0519 
0520       float phi1, invpt1, dctheta, ctheta1, dphi, dr2;
0521       for (auto itrack = 0U; itrack < ntracks; itrack++) {
0522         auto &trk = tracks[itrack];
0523         phi1 = trk.momPhi();
0524         invpt1 = trk.invpT();
0525         ctheta1 = ctheta[itrack];
0526         for (auto jtrack = itrack + 1; jtrack < ntracks; jtrack++) {
0527           auto &track2 = tracks[jtrack];
0528           if (trk.label() == track2.label())
0529             continue;
0530 
0531           dctheta = std::abs(ctheta[jtrack] - ctheta1);
0532 
0533           if (dctheta > Config::maxdcth)
0534             continue;
0535 
0536           dphi = std::abs(squashPhiMinimal(phi1 - track2.momPhi()));
0537 
0538           if (dphi > Config::maxdphi)
0539             continue;
0540 
0541           float maxdRSquared = drth_central * drth_central;
0542           if (std::abs(ctheta1) > Config::maxcth_fw)
0543             maxdRSquared = drth_forward * drth_forward;
0544           else if (std::abs(ctheta1) > Config::maxcth_ob)
0545             maxdRSquared = drth_obarrel * drth_obarrel;
0546           dr2 = dphi * dphi + dctheta * dctheta;
0547           if (dr2 < maxdRSquared) {
0548             //Keep track with best score
0549             if (trk.score() > track2.score())
0550               track2.setDuplicateValue(true);
0551             else
0552               trk.setDuplicateValue(true);
0553             continue;
0554           }
0555 
0556           if (std::abs(track2.invpT() - invpt1) > Config::maxd1pt)
0557             continue;
0558 
0559           auto sharedCount = 0;
0560           auto sharedFirst = 0;
0561           const auto minFoundHits = std::min(trk.nFoundHits(), track2.nFoundHits());
0562 
0563           for (int i = 0; i < trk.nTotalHits(); ++i) {
0564             if (trk.getHitIdx(i) < 0)
0565               continue;
0566             const int a = trk.getHitLyr(i);
0567             const int b = trk.getHitIdx(i);
0568             for (int j = 0; j < track2.nTotalHits(); ++j) {
0569               if (track2.getHitIdx(j) < 0)
0570                 continue;
0571               const int c = track2.getHitLyr(j);
0572               const int d = track2.getHitIdx(j);
0573 
0574               //this is to count once shared matched hits (may be done more properly...)
0575               if (a == c && b == d)
0576                 sharedCount += 1;
0577               if (j == 0 && i == 0 && a == c && b == d)
0578                 sharedFirst += 1;
0579 
0580               if ((sharedCount - sharedFirst) >= ((minFoundHits - sharedFirst) * fraction))
0581                 continue;
0582             }
0583             if ((sharedCount - sharedFirst) >= ((minFoundHits - sharedFirst) * fraction))
0584               continue;
0585           }
0586 
0587           //selection here - 11percent fraction of shared hits to label a duplicate
0588           if ((sharedCount - sharedFirst) >= ((minFoundHits - sharedFirst) * fraction)) {
0589             if (trk.score() > track2.score())
0590               track2.setDuplicateValue(true);
0591             else
0592               trk.setDuplicateValue(true);
0593           }
0594         }
0595       }  //end loop one over tracks
0596 
0597       remove_duplicates(tracks);
0598     }
0599 
0600     namespace {
0601       CMS_SA_ALLOW struct register_duplicate_cleaners {
0602         register_duplicate_cleaners() {
0603           IterationConfig::register_duplicate_cleaner("phase1:clean_duplicates", clean_duplicates);
0604           IterationConfig::register_duplicate_cleaner("phase1:clean_duplicates_sharedhits",
0605                                                       clean_duplicates_sharedhits);
0606           IterationConfig::register_duplicate_cleaner("phase1:clean_duplicates_sharedhits_pixelseed",
0607                                                       clean_duplicates_sharedhits_pixelseed);
0608         }
0609       } rdc_instance;
0610     }  // namespace
0611 
0612     //=========================================================================
0613     // Quality filters
0614     //=========================================================================
0615 
0616     // quality filter for n hits with seed hit "penalty" for strip-based seeds
0617     //   this implicitly separates triplets and doublet seeds with glued layers
0618     template <class TRACK>
0619     bool qfilter_n_hits(const TRACK &t, const MkJob &j) {
0620       int seedHits = t.getNSeedHits();
0621       int seedReduction = (seedHits <= 5) ? 2 : 3;
0622       return t.nFoundHits() - seedReduction >= j.params_cur().minHitsQF;
0623     }
0624 
0625     // simple hit-count quality filter; used with pixel-based seeds
0626     template <class TRACK>
0627     bool qfilter_n_hits_pixseed(const TRACK &t, const MkJob &j) {
0628       return t.nFoundHits() >= j.params_cur().minHitsQF;
0629     }
0630 
0631     // layer-dependent quality filter
0632     // includes ad hoc tuning for phase-1
0633     template <class TRACK>
0634     bool qfilter_n_layers(const TRACK &t, const MkJob &j) {
0635       const BeamSpot &bspot = j.m_beam_spot;
0636       const TrackerInfo &trk_inf = j.m_trk_info;
0637       int enhits = t.nHitsByTypeEncoded(trk_inf);
0638       int npixhits = t.nPixelDecoded(enhits);
0639       int enlyrs = t.nLayersByTypeEncoded(trk_inf);
0640       int npixlyrs = t.nPixelDecoded(enlyrs);
0641       int nmatlyrs = t.nTotMatchDecoded(enlyrs);
0642       int llyr = t.getLastFoundHitLyr();
0643       int lplyr = t.getLastFoundPixelHitLyr();
0644       float invpt = t.invpT();
0645 
0646       // based on fr and eff vs pt (convert to native invpt)
0647       float invptmin = 1.43;  // min 1/pT (=1/0.7) for full filter on (npixhits<=3 .or. npixlyrs<=3)
0648       float d0BS = t.d0BeamSpot(bspot.x, bspot.y);
0649       float d0_max = 0.1;  // 1 mm, max for somewhat prompt
0650 
0651       // next-to-outermost pixel layers (almost): BPIX3 or FPIX1
0652       bool endsInsidePix = (llyr == 2 || llyr == 18 || llyr == 45);
0653       // not last pixel layers: BPIX[123] or FPIX[12]
0654       bool lastInsidePix = ((0 <= lplyr && lplyr < 3) || (18 <= lplyr && lplyr < 20) || (45 <= lplyr && lplyr < 47));
0655       // reject short tracks missing last pixel layer except for prompt-looking
0656       return !(((npixhits <= 3 || npixlyrs <= 3) && endsInsidePix &&
0657                 (invpt < invptmin || (invpt >= invptmin && std::abs(d0BS) > d0_max))) ||
0658                ((npixlyrs <= 3 && nmatlyrs <= 6) && lastInsidePix && llyr != lplyr && std::abs(d0BS) > d0_max));
0659     }
0660 
0661     /// quality filter tuned for pixelLess iteration during forward search
0662     // includes ad hoc tuning for phase-1
0663     template <class TRACK>
0664     bool qfilter_pixelLessFwd(const TRACK &t, const MkJob &j) {
0665       const BeamSpot &bspot = j.m_beam_spot;
0666       const TrackerInfo &tk_info = j.m_trk_info;
0667       float d0BS = t.d0BeamSpot(bspot.x, bspot.y);
0668       float d0_max = 0.05;  // 0.5 mm, max for somewhat prompt
0669 
0670       int encoded;
0671       encoded = t.nLayersByTypeEncoded(tk_info);
0672       int nLyrs = t.nTotMatchDecoded(encoded);
0673       encoded = t.nHitsByTypeEncoded(tk_info);
0674       int nHits = t.nTotMatchDecoded(encoded);
0675 
0676       // to subtract stereo seed layers to count just r-phi seed layers (better pt err)
0677       int seedReduction = (t.getNSeedHits() <= 5) ? 2 : 3;
0678 
0679       // based on fr and eff vs pt and eta (convert to native invpt and theta)
0680       float invpt = t.invpT();
0681       float invptmin = 1.11;  // =1/0.9
0682 
0683       float thetasym = std::abs(t.theta() - Const::PIOver2);
0684       float thetasymmin = 1.11;  // -> |eta|=1.45
0685 
0686       // accept longer tracks, reject too short and displaced
0687       return (((t.nFoundHits() - seedReduction >= 4 && invpt < invptmin) ||
0688                (t.nFoundHits() - seedReduction >= 3 && invpt > invptmin && thetasym <= thetasymmin) ||
0689                (t.nFoundHits() - seedReduction >= 4 && invpt > invptmin && thetasym > thetasymmin)) &&
0690               !((nLyrs <= 4 || nHits <= 4) && std::abs(d0BS) > d0_max && invpt < invptmin));
0691     }
0692 
0693     /// quality filter tuned for pixelLess iteration during backward search
0694     // includes ad hoc tuning for phase-1
0695     template <class TRACK>
0696     bool qfilter_pixelLessBkwd(const TRACK &t, const MkJob &j) {
0697       const BeamSpot &bspot = j.m_beam_spot;
0698       const TrackerInfo &tk_info = j.m_trk_info;
0699       float d0BS = t.d0BeamSpot(bspot.x, bspot.y);
0700       float d0_max = 0.1;  // 1 mm
0701 
0702       int encoded;
0703       encoded = t.nLayersByTypeEncoded(tk_info);
0704       int nLyrs = t.nTotMatchDecoded(encoded);
0705       encoded = t.nHitsByTypeEncoded(tk_info);
0706       int nHits = t.nTotMatchDecoded(encoded);
0707 
0708       // based on fr and eff vs pt and eta (convert to native invpt and theta)
0709       float invpt = t.invpT();
0710       float invptmin = 1.11;  // =1/0.9
0711 
0712       float thetasym = std::abs(t.theta() - Const::PIOver2);
0713       float thetasymmin_l = 0.80;  // -> |eta|=0.9
0714       float thetasymmin_h = 1.11;  // -> |eta|=1.45
0715 
0716       // reject too short or too displaced tracks
0717       return !(
0718           ((nLyrs <= 3 || nHits <= 3)) ||
0719           ((nLyrs <= 4 || nHits <= 4) && (invpt < invptmin || (thetasym > thetasymmin_l && std::abs(d0BS) > d0_max))) ||
0720           ((nLyrs <= 5 || nHits <= 5) && (invpt > invptmin && thetasym > thetasymmin_h && std::abs(d0BS) > d0_max)));
0721     }
0722 
0723     namespace {
0724       CMS_SA_ALLOW struct register_quality_filters {
0725         register_quality_filters() {
0726           IterationConfig::register_candidate_filter("phase1:qfilter_n_hits", qfilter_n_hits<TrackCand>);
0727           IterationConfig::register_candidate_filter("phase1:qfilter_n_hits_pixseed",
0728                                                      qfilter_n_hits_pixseed<TrackCand>);
0729           IterationConfig::register_candidate_filter("phase1:qfilter_n_layers", qfilter_n_layers<TrackCand>);
0730           IterationConfig::register_candidate_filter("phase1:qfilter_pixelLessFwd", qfilter_pixelLessFwd<TrackCand>);
0731           IterationConfig::register_candidate_filter("phase1:qfilter_pixelLessBkwd", qfilter_pixelLessBkwd<TrackCand>);
0732         }
0733       } rqf_instance;
0734     }  // namespace
0735 
0736     //=========================================================================
0737     // Track scoring
0738     //=========================================================================
0739 
0740     float trackScoreDefault(const int nfoundhits,
0741                             const int ntailholes,
0742                             const int noverlaphits,
0743                             const int nmisshits,
0744                             const float chi2,
0745                             const float pt,
0746                             const bool inFindCandidates) {
0747       float maxBonus = 8.0;
0748       float bonus = Config::validHitSlope_ * nfoundhits + Config::validHitBonus_;
0749       float penalty = Config::missingHitPenalty_;
0750       float tailPenalty = Config::tailMissingHitPenalty_;
0751       float overlapBonus = Config::overlapHitBonus_;
0752       if (pt < 0.9) {
0753         penalty *= inFindCandidates ? 1.7f : 1.5f;
0754         bonus = std::min(bonus * (inFindCandidates ? 0.9f : 1.0f), maxBonus);
0755       }
0756       float score =
0757           bonus * nfoundhits + overlapBonus * noverlaphits - penalty * nmisshits - tailPenalty * ntailholes - chi2;
0758       return score;
0759     }
0760 
0761     namespace {
0762       CMS_SA_ALLOW struct register_track_scorers {
0763         register_track_scorers() {
0764           IterationConfig::register_track_scorer("default", trackScoreDefault);
0765           IterationConfig::register_track_scorer("phase1:default", trackScoreDefault);
0766         }
0767       } rts_instance;
0768     }  // namespace
0769 
0770   }  // namespace StdSeq
0771 }  // namespace mkfit