diff --git a/RecoTracker/MkFitCMS/src/MkStdSeqs.cc b/RecoTracker/MkFitCMS/src/MkStdSeqs.cc index 3a6ad20c5dfdb..25181dd42e216 100644 --- a/RecoTracker/MkFitCMS/src/MkStdSeqs.cc +++ b/RecoTracker/MkFitCMS/src/MkStdSeqs.cc @@ -137,6 +137,12 @@ namespace mkfit { std::vector d0(ns); int i1, i2; //for the sorting + axis_pow2_u1 ax_phi(-Const::PI, Const::PI); + axis ax_eta(-3.0, 3.0, 30u); + binnor phi_eta_binnor(ax_phi, ax_eta); + + phi_eta_binnor.begin_registration(ns); + for (int ts = 0; ts < ns; ts++) { const Track &tk = seeds[ts]; nHits[ts] = tk.nFoundHits(); @@ -151,9 +157,16 @@ namespace mkfit { y[ts] = tk.y(); z[ts] = tk.z(); d0[ts] = tk.d0BeamSpot(bspot.x, bspot.y); + + phi_eta_binnor.register_entry_safe(oldPhi[ts], eta[ts]); + // If one is sure values are *within* axis ranges: b.register_entry(oldPhi[ts], eta[ts]); } - for (int ts = 0; ts < ns; ts++) { + phi_eta_binnor.finalize_registration(); + + for (int sorted_ts = 0; sorted_ts < ns; sorted_ts++) { + int ts = phi_eta_binnor.m_ranks[sorted_ts]; + if (not writetrack[ts]) continue; // Note: this speed up prevents transitive masking (possibly marginal gain). @@ -166,105 +179,122 @@ namespace mkfit { // To study some more details -- need EventOfHits for this int n_ovlp_hits_added = 0; - for (int tss = ts + 1; tss < ns; tss++) { - const float pt2 = pt[tss]; + auto phi_rng = ax_phi.from_R_rdr_to_N_bins(oldPhi[ts], 0.08f); + auto eta_rng = ax_eta.from_R_rdr_to_N_bins(eta[ts], .1f); - ////// Always require charge consistency. If different charge is assigned, do not remove seed-track - if (charge[tss] != charge[ts]) - continue; + for (auto i_phi = phi_rng.begin; i_phi != phi_rng.end; i_phi = ax_phi.next_N_bin(i_phi)) { + for (auto i_eta = eta_rng.begin; i_eta != eta_rng.end; i_eta = ax_eta.next_N_bin(i_eta)) { + const auto cbin = phi_eta_binnor.get_content(i_phi, i_eta); + for (auto i = cbin.first; i < cbin.end(); ++i) { + int tss = phi_eta_binnor.m_ranks[i]; - const float thisDPt = std::abs(pt2 - pt1); - ////// Require pT consistency between seeds. If dpT is large, do not remove seed-track. - if (thisDPt > dpt_common * (pt1)) - continue; + if (not writetrack[ts]) + continue; + if (not writetrack[tss]) + continue; + if (tss == ts) + continue; - const float eta2 = eta[tss]; - const float deta2 = std::pow(eta1 - eta2, 2); + const float pt2 = pt[tss]; - const float oldPhi2 = oldPhi[tss]; + // Always require charge consistency. If different charge is assigned, do not remove seed-track + if (charge[tss] != charge[ts]) + continue; - const float pos2_second = pos2[tss]; - const float thisDXYSign05 = pos2_second > pos2_first ? -0.5f : 0.5f; + const float thisDPt = std::abs(pt2 - pt1); + // Require pT consistency between seeds. If dpT is large, do not remove seed-track. + if (thisDPt > dpt_common * pt1) + continue; - const float thisDXY = thisDXYSign05 * sqrt(std::pow(x[ts] - x[tss], 2) + std::pow(y[ts] - y[tss], 2)); + const float eta2 = eta[tss]; + const float deta2 = std::pow(eta1 - eta2, 2); - const float invptq_second = invptq[tss]; + const float oldPhi2 = oldPhi[tss]; - const float newPhi1 = oldPhi1 - thisDXY * invR1GeV * invptq_first; - const float newPhi2 = oldPhi2 + thisDXY * invR1GeV * invptq_second; + const float pos2_second = pos2[tss]; + const float thisDXYSign05 = pos2_second > pos2_first ? -0.5f : 0.5f; - const float dphi = cdist(std::abs(newPhi1 - newPhi2)); + const float thisDXY = thisDXYSign05 * sqrt(std::pow(x[ts] - x[tss], 2) + std::pow(y[ts] - y[tss], 2)); - const float dr2 = deta2 + dphi * dphi; + const float invptq_second = invptq[tss]; - const float thisDZ = z[ts] - z[tss] - thisDXY * (ctheta[ts] + ctheta[tss]); - const float dz2 = thisDZ * thisDZ; + const float newPhi1 = oldPhi1 - thisDXY * invR1GeV * invptq_first; + const float newPhi2 = oldPhi2 + thisDXY * invR1GeV * invptq_second; - ////// Reject tracks within dR-dz elliptical window. - ////// Adaptive thresholds, based on observation that duplicates are more abundant at large pseudo-rapidity and low track pT - bool overlapping = false; - if (std::abs(eta1) < etamax_brl) { - if (pt1 > ptmin_hpt) { - if (dz2 * dzmax2_inv_bh + dr2 * drmax2_inv_bh < 1.0f) - overlapping = true; - } else { - if (dz2 * dzmax2_inv_bl + dr2 * drmax2_inv_bl < 1.0f) - overlapping = true; - } - } else { - if (pt1 > ptmin_hpt) { - if (dz2 * dzmax2_inv_eh + dr2 * drmax2_inv_eh < 1.0f) - overlapping = true; - } else { - if (dz2 * dzmax2_inv_el + dr2 * drmax2_inv_el < 1.0f) - overlapping = true; - } - } + const float dphi = cdist(std::abs(newPhi1 - newPhi2)); - if (overlapping) { - //Mark tss as a duplicate - i1 = ts; - i2 = tss; - if (d0[tss] > d0[ts]) - writetrack[tss] = false; - else { - writetrack[ts] = false; - i2 = ts; - i1 = tss; - } - // Add hits from tk2 to the seed we are keeping. - // NOTE: We have a limit in Track::Status for the number of seed hits. - // There is a check at entry and after adding of a new hit. - Track &tk = seeds[i1]; - if (merge_hits && tk.nTotalHits() < Track::Status::kMaxSeedHits) { - const Track &tk2 = seeds[i2]; - //We are not actually fitting to the extra hits; use chi2 of 0 - float fakeChi2 = 0.0; - - for (int j = 0; j < tk2.nTotalHits(); ++j) { - int hitidx = tk2.getHitIdx(j); - int hitlyr = tk2.getHitLyr(j); - if (hitidx >= 0) { - bool unique = true; - for (int i = 0; i < tk.nTotalHits(); ++i) { - if ((hitidx == tk.getHitIdx(i)) && (hitlyr == tk.getHitLyr(i))) { - unique = false; - break; + const float dr2 = deta2 + dphi * dphi; + + const float thisDZ = z[ts] - z[tss] - thisDXY * (ctheta[ts] + ctheta[tss]); + const float dz2 = thisDZ * thisDZ; + + // Reject tracks within dR-dz elliptical window. + // Adaptive thresholds, based on observation that duplicates are more abundant at large pseudo-rapidity and low track pT + bool overlapping = false; + if (std::abs(eta1) < etamax_brl) { + if (pt1 > ptmin_hpt) { + if (dz2 * dzmax2_inv_bh + dr2 * drmax2_inv_bh < 1.0f) + overlapping = true; + } else { + if (dz2 * dzmax2_inv_bl + dr2 * drmax2_inv_bl < 1.0f) + overlapping = true; + } + } else { + if (pt1 > ptmin_hpt) { + if (dz2 * dzmax2_inv_eh + dr2 * drmax2_inv_eh < 1.0f) + overlapping = true; + } else { + if (dz2 * dzmax2_inv_el + dr2 * drmax2_inv_el < 1.0f) + overlapping = true; + } + } + + if (overlapping) { + //Mark tss as a duplicate + i1 = ts; + i2 = tss; + if (d0[tss] > d0[ts]) + writetrack[tss] = false; + else { + writetrack[ts] = false; + i2 = ts; + i1 = tss; + } + // Add hits from tk2 to the seed we are keeping. + // NOTE: We have a limit in Track::Status for the number of seed hits. + // There is a check at entry and after adding of a new hit. + Track &tk = seeds[i1]; + if (merge_hits && tk.nTotalHits() < Track::Status::kMaxSeedHits) { + const Track &tk2 = seeds[i2]; + //We are not actually fitting to the extra hits; use chi2 of 0 + float fakeChi2 = 0.0; + + for (int j = 0; j < tk2.nTotalHits(); ++j) { + int hitidx = tk2.getHitIdx(j); + int hitlyr = tk2.getHitLyr(j); + if (hitidx >= 0) { + bool unique = true; + for (int i = 0; i < tk.nTotalHits(); ++i) { + if ((hitidx == tk.getHitIdx(i)) && (hitlyr == tk.getHitLyr(i))) { + unique = false; + break; + } + } + if (unique) { + tk.addHitIdx(tk2.getHitIdx(j), tk2.getHitLyr(j), fakeChi2); + ++n_ovlp_hits_added; + if (tk.nTotalHits() >= Track::Status::kMaxSeedHits) + break; + } } } - if (unique) { - tk.addHitIdx(tk2.getHitIdx(j), tk2.getHitLyr(j), fakeChi2); - ++n_ovlp_hits_added; - if (tk.nTotalHits() >= Track::Status::kMaxSeedHits) - break; - } } + if (n_ovlp_hits_added > 0) + tk.sortHitsByLayer(); } - } - if (n_ovlp_hits_added > 0) - tk.sortHitsByLayer(); - } - } //end of inner loop over tss + } //end of inner loop over tss + } //eta bin + } //phi bin if (writetrack[ts]) { cleanSeedTracks.emplace_back(seeds[ts]); @@ -399,9 +429,18 @@ namespace mkfit { const auto ntracks = tracks.size(); std::vector ctheta(ntracks); + std::multimap hitMap; + for (auto itrack = 0U; itrack < ntracks; itrack++) { auto &trk = tracks[itrack]; ctheta[itrack] = 1.f / std::tan(trk.theta()); + for (int i = 0; i < trk.nTotalHits(); ++i) { + if (trk.getHitIdx(i) < 0) + continue; + int a = trk.getHitLyr(i); + int b = trk.getHitIdx(i); + hitMap.insert(std::make_pair(b * 1000 + a, i > 0 ? itrack : -itrack)); //negative for first hit in trk + } } for (auto itrack = 0U; itrack < ntracks; itrack++) { @@ -409,12 +448,27 @@ namespace mkfit { auto phi1 = trk.momPhi(); auto ctheta1 = ctheta[itrack]; - for (auto jtrack = itrack + 1; jtrack < ntracks; jtrack++) { - auto &track2 = tracks[jtrack]; - auto sharedCount = 0; - auto sharedFirst = 0; + std::map sharingMap; + for (int i = 0; i < trk.nTotalHits(); ++i) { + if (trk.getHitIdx(i) < 0) + continue; + int a = trk.getHitLyr(i); + int b = trk.getHitIdx(i); + auto range = hitMap.equal_range(b * 1000 + a); + for (auto it = range.first; it != range.second; ++it) { + if (std::abs(it->second) >= (int)itrack) + continue; // don't check your own hits (==) nor sym. checks (>) + if (i == 0 && it->second < 0) + continue; // shared first - first is not counted + sharingMap[std::abs(it->second)]++; + } + } - auto dctheta = std::abs(ctheta[jtrack] - ctheta1); + for (const auto &elem : sharingMap) { + auto &track2 = tracks[elem.first]; + + // broad dctheta-dphi compatibility checks; keep mostly to preserve consistency with old results + auto dctheta = std::abs(ctheta[elem.first] - ctheta1); if (dctheta > 1.) continue; @@ -422,31 +476,15 @@ namespace mkfit { if (dphi > 1.) continue; - for (int i = 0; i < trk.nTotalHits(); ++i) { - if (trk.getHitIdx(i) < 0) - continue; - const int a = trk.getHitLyr(i); - const int b = trk.getHitIdx(i); - for (int j = 0; j < track2.nTotalHits(); ++j) { - if (track2.getHitIdx(j) < 0) - continue; - const int c = track2.getHitLyr(j); - const int d = track2.getHitIdx(j); - if (a == c && b == d) - sharedCount += 1; - if (a == c && b == d && j == 0 && i == 0) - sharedFirst += 1; - } - } - if ((sharedCount - sharedFirst) >= - ((std::min(trk.nFoundHits(), track2.nFoundHits()) - sharedFirst) * (fraction))) { + if (elem.second >= std::min(trk.nFoundHits(), track2.nFoundHits()) * fraction) { if (trk.score() > track2.score()) track2.setDuplicateValue(true); else trk.setDuplicateValue(true); } - } - } + } // end sharing hits loop + } // end trk loop + tracks.erase(std::remove_if(tracks.begin(), tracks.end(), [](auto track) { return track.getDuplicateValue(); }), tracks.end()); } diff --git a/RecoTracker/MkFitCore/interface/binnor.h b/RecoTracker/MkFitCore/interface/binnor.h index cac5aea4ddd65..aa4cf28b12598 100644 --- a/RecoTracker/MkFitCore/interface/binnor.h +++ b/RecoTracker/MkFitCore/interface/binnor.h @@ -5,7 +5,7 @@ #include #include #include - +#include #include namespace mkfit { @@ -47,26 +47,30 @@ namespace mkfit { m_M_fac(M_size / (max - min)), m_N_fac(N_size / (max - min)), m_last_M_bin(M_size - 1), - m_last_N_bin(N_size - 1) {} + m_last_N_bin(N_size - 1) { + // Requested number of bins must fit within the intended bit-field (declared by binnor, later). + assert(N_size <= (1 << N)); + } - I R_to_M_bin(R r) const { return (r - m_R_min) * m_M_fac; } - I R_to_N_bin(R r) const { return (r - m_R_min) * m_N_fac; } + I from_R_to_M_bin(R r) const { return (r - m_R_min) * m_M_fac; } + I from_R_to_N_bin(R r) const { return (r - m_R_min) * m_N_fac; } - I R_to_M_bin_safe(R r) const { return r <= m_R_min ? 0 : (r >= m_R_max ? m_last_M_bin : R_to_M_bin(r)); } - I R_to_N_bin_safe(R r) const { return r <= m_R_min ? 0 : (r >= m_R_max ? m_last_N_bin : R_to_N_bin(r)); } + I from_R_to_M_bin_safe(R r) const { return r <= m_R_min ? 0 : (r >= m_R_max ? m_last_M_bin : from_R_to_M_bin(r)); } + I from_R_to_N_bin_safe(R r) const { return r <= m_R_min ? 0 : (r >= m_R_max ? m_last_N_bin : from_R_to_N_bin(r)); } - I M_bin_to_N_bin(I m) const { return m >> c_M2N_shift; } + I from_M_bin_to_N_bin(I m) const { return m >> c_M2N_shift; } - I_pair Rminmax_to_N_bins(R rmin, R rmax) const { - return I_pair(R_to_N_bin_safe(rmin), R_to_N_bin_safe(rmax) + I{1}); + I_pair from_R_minmax_to_N_bins(R rmin, R rmax) const { + return I_pair(from_R_to_N_bin_safe(rmin), from_R_to_N_bin_safe(rmax) + I{1}); } - I_pair Rrdr_to_N_bins(R r, R dr) const { return Rminmax_to_N_bins(r - dr, r + dr); } + I_pair from_R_rdr_to_N_bins(R r, R dr) const { return from_R_minmax_to_N_bins(r - dr, r + dr); } I next_N_bin(I bin) const { return bin + 1; } }; // axis_pow2_base //--------------- + // Assumes the numbers of fine/normal bins are powers of 2 that are inferred directly from the number of bits. template struct axis_pow2_base : public axis_base { static constexpr unsigned c_M_end = 1 << M; @@ -80,6 +84,8 @@ namespace mkfit { // axis_pow2_u1 //------------- + // Specialization of axis_pow2 for the "U(1)" case where the coordinate is periodic with period (Rmax - Rmin). + // In the "safe" methods below, bit masking serves as the modulo operator for out-of-range bin numbers. template struct axis_pow2_u1 : public axis_pow2_base { static constexpr I c_M_mask = (1 << M) - 1; @@ -87,14 +93,17 @@ namespace mkfit { axis_pow2_u1(R min, R max) : axis_pow2_base(min, max) {} - I R_to_M_bin_safe(R r) const { return this->R_to_M_bin(r) & c_M_mask; } - I R_to_N_bin_safe(R r) const { return this->R_to_N_bin(r) & c_N_mask; } + I from_R_to_M_bin_safe(R r) const { return this->from_R_to_M_bin(r) & c_M_mask; } + I from_R_to_N_bin_safe(R r) const { return this->from_R_to_N_bin(r) & c_N_mask; } - typename axis_base::I_pair Rminmax_to_N_bins(R rmin, R rmax) const { - return typename axis_base::I_pair(R_to_N_bin_safe(rmin), (this->R_to_N_bin(rmax) + I{1}) & c_N_mask); + typename axis_base::I_pair from_R_minmax_to_N_bins(R rmin, R rmax) const { + return typename axis_base::I_pair(from_R_to_N_bin_safe(rmin), + (this->from_R_to_N_bin(rmax) + I{1}) & c_N_mask); } - typename axis_base::I_pair Rrdr_to_N_bins(R r, R dr) const { return Rminmax_to_N_bins(r - dr, r + dr); } + typename axis_base::I_pair from_R_rdr_to_N_bins(R r, R dr) const { + return from_R_minmax_to_N_bins(r - dr, r + dr); + } I next_N_bin(I bin) const { return (bin + 1) & c_N_mask; } }; @@ -130,7 +139,18 @@ namespace mkfit { // binnor //--------------- - // C - bin content type + // To build and populate bins, do the following: + // 1. Construct two axis objects, giving numbers of bits and bins, and extents. + // 2. Construct a binnor from the axis objects, and begin_registration on it. + // 3. Loop register_entry (optional: _safe) over pairs of coordinates, to fill + // m_cons with the corresponding pairs of bin indices (B_pairs). + // 4. Call finalize_registration, which sorts m_ranks based on m_cons, making + // m_ranks into an in-order map into m_cons (as well as the inputs that were + // used to fill it). Final counts for all the bins, as well as starting + // indices for the bins (within m_ranks), are computed and stored in packed + // form (i.e., bit-fields) in m_bins. + // + // C - bin content type, to hold "bin population coordinates" in packed form (bit-fields) // A1, A2 - axis types // NB_first, NB_count - number of bits for storage of { first, count } pairs @@ -139,12 +159,12 @@ namespace mkfit { static_assert(std::is_same()); static_assert(A1::c_M + A2::c_M <= 32); - static constexpr unsigned c_A1_mask = (1 << A1::c_M) - 1; - static constexpr unsigned c_A2_Mout_mask = ~(((1 << A2::c_M2N_shift) - 1) << A1::c_M); + static constexpr unsigned int c_A1_mask = (1 << A1::c_M) - 1; + static constexpr unsigned int c_A2_Mout_mask = ~(((1 << A2::c_M2N_shift) - 1) << A1::c_M); // Pair of axis bin indices packed into unsigned. struct B_pair { - unsigned int packed_value; // bin1 in A1::c_M lower bits, bin2 above + unsigned int packed_value; // bin1 in A1::c_M lower bits, bin2 above B_pair() : packed_value(0) {} B_pair(typename A1::index_t i1, typename A2::index_t i2) : packed_value(i2 << A1::c_M | i1) {} @@ -155,7 +175,7 @@ namespace mkfit { unsigned int mask_A2_M_bins() const { return packed_value & c_A2_Mout_mask; } }; - // Bin content pair. + // Bin content pair (bit-fields). struct C_pair { C first : NB_first; C count : NB_count; @@ -176,12 +196,14 @@ namespace mkfit { // Access - B_pair m_bin_to_n_bin(B_pair m_bin) { return {m_a1.M_bin_to_N_bin(m_bin.bin1()), m_a2.M_bin_to_N_bin(m_bin.bin2())}; } + B_pair m_bin_to_n_bin(B_pair m_bin) { + return {m_a1.from_M_bin_to_N_bin(m_bin.bin1()), m_a2.from_M_bin_to_N_bin(m_bin.bin2())}; + } B_pair get_n_bin(typename A1::index_t n1, typename A2::index_t n2) const { return {n1, n2}; } B_pair get_n_bin(typename A1::real_t r1, typename A2::real_t r2) const { - return {m_a1.R_to_N_bin(r1), m_a2.R_to_N_bin(r2)}; + return {m_a1.from_R_to_N_bin(r1), m_a2.from_R_to_N_bin(r2)}; } C_pair &ref_content(B_pair n_bin) { return m_bins[n_bin.bin2() * m_a1.size_of_N() + n_bin.bin1()]; } @@ -193,7 +215,7 @@ namespace mkfit { } C_pair get_content(typename A1::real_t r1, typename A2::real_t r2) const { - return get_content(m_a1.R_to_N_bin(r1), m_a2.R_to_N_bin(r2)); + return get_content(m_a1.from_R_to_N_bin(r1), m_a2.from_R_to_N_bin(r2)); } // Filling @@ -207,11 +229,11 @@ namespace mkfit { void begin_registration(C n_items) { m_cons.reserve(n_items); } void register_entry(typename A1::real_t r1, typename A2::real_t r2) { - m_cons.push_back({m_a1.R_to_M_bin(r1), m_a2.R_to_M_bin(r2)}); + m_cons.push_back({m_a1.from_R_to_M_bin(r1), m_a2.from_R_to_M_bin(r2)}); } void register_entry_safe(typename A1::real_t r1, typename A2::real_t r2) { - m_cons.push_back({m_a1.R_to_M_bin_safe(r1), m_a2.R_to_M_bin_safe(r2)}); + m_cons.push_back({m_a1.from_R_to_M_bin_safe(r1), m_a2.from_R_to_M_bin_safe(r2)}); } // Do M-binning outside, potentially using R_to_M_bin_safe(). diff --git a/RecoTracker/MkFitCore/standalone/test/binnor_demo.cxx b/RecoTracker/MkFitCore/standalone/test/binnor_demo.cxx index 1919c22b6055f..c7a165674d2af 100644 --- a/RecoTracker/MkFitCore/standalone/test/binnor_demo.cxx +++ b/RecoTracker/MkFitCore/standalone/test/binnor_demo.cxx @@ -1,5 +1,6 @@ #include "../../interface/binnor.h" #include +#include // build as: // c++ -o binnor_demo -std=c++17 binnor_demo.cxx @@ -21,9 +22,9 @@ int main() /* for (float p = -TwoPI; p < TwoPI; p += TwoPI / 15.4f) { printf(" phi=%-9f m=%5d n=%3d m2n=%3d n_safe=%3d\n", p, - phi.R_to_M_bin(p), phi.R_to_N_bin(p), - phi.M_bin_to_N_bin( phi.R_to_M_bin(p) ), - phi.R_to_N_bin_safe(p) ); + phi.from_R_to_M_bin(p), phi.from_R_to_N_bin(p), + phi.from_M_bin_to_N_bin( phi.from_R_to_M_bin(p) ), + phi.from_R_to_N_bin_safe(p) ); } */ @@ -48,6 +49,8 @@ int main() std::vector tracks; tracks.reserve(NN); + auto start = std::chrono::high_resolution_clock::now(); + b.begin_registration(NN); // optional, reserves construction vector for (int i = 0; i < NN; ++i) @@ -59,6 +62,9 @@ int main() b.finalize_registration(); + auto stop = std::chrono::high_resolution_clock::now(); + auto duration = std::chrono::duration_cast(stop - start); + // for (int i = 0; i < NN; ++i) // { // const track &t = tracks[ b.m_ranks[i] ]; @@ -75,8 +81,8 @@ int main() } printf("\n\n--- Range access for phi=[(-PI+0.02 +- 0.1], eta=[1.3 +- .2]:\n\n"); - auto phi_rng = phi.Rrdr_to_N_bins(-PI+0.02, 0.1); - auto eta_rng = eta.Rrdr_to_N_bins(1.3, .2); + auto phi_rng = phi.from_R_rdr_to_N_bins(-PI+0.02, 0.1); + auto eta_rng = eta.from_R_rdr_to_N_bins(1.3, .2); printf("phi bin range: %u, %u; eta %u, %u\n", phi_rng.begin, phi_rng.end, eta_rng.begin, eta_rng.end); for (auto i_phi = phi_rng.begin; i_phi != phi_rng.end; i_phi = phi.next_N_bin(i_phi)) { @@ -91,291 +97,9 @@ int main() } } + printf("\nBinning time for %d points: %f sec\n", NN, 1.e-6*duration.count()); b.reset_contents(); return 0; } - - - -// buildtestMPlex.cc::runBtpCe_MultiIter(), loop over seed cleaning multiple times to measure time: -/* - if ( itconf.m_requires_dupclean_tight ) { - double t0 = dtime(); - TrackVec xxx; int n_comp; - for (int i=0;i<1000;++i) { - xxx = seeds; - n_comp = StdSeq::clean_cms_seedtracks_iter(&xxx, itconf, eoh.m_beam_spot); - } - printf("Seedacleena of %d seeds, out_seeds %d, 1000 times, N_comparisons=%d, took %.5fs\n", - (int)seeds.size(), (int)xxx.size(), n_comp, dtime() - t0); - seeds = xxx; - } -*/ - -// Example clean seeds using binnor. -// Further enhancements possible by moving the iteration into binnor class (+iterator), -// doing pre-selection on m_cons fine m-bins. -// Perf notes: https://gist.github.com/osschar/2dcd2b01e7c15cc25aa6489f3b242ccb -/* -//========================================================================= -// Seed cleaning (multi-iter) -//========================================================================= -int clean_cms_seedtracks_iter(TrackVec *seed_ptr, const IterationConfig& itrcfg, const BeamSpot &bspot) -{ - const float etamax_brl = Config::c_etamax_brl; - const float dpt_common = Config::c_dpt_common; - - const float dzmax_bh = itrcfg.m_params.c_dzmax_bh; - const float drmax_bh = itrcfg.m_params.c_drmax_bh; - const float dzmax_eh = itrcfg.m_params.c_dzmax_eh; - const float drmax_eh = itrcfg.m_params.c_drmax_eh; - const float dzmax_bl = itrcfg.m_params.c_dzmax_bl; - const float drmax_bl = itrcfg.m_params.c_drmax_bl; - const float dzmax_el = itrcfg.m_params.c_dzmax_el; - const float drmax_el = itrcfg.m_params.c_drmax_el; - - const float ptmin_hpt = itrcfg.m_params.c_ptthr_hpt; - - const float dzmax2_inv_bh = 1.f/(dzmax_bh*dzmax_bh); - const float drmax2_inv_bh = 1.f/(drmax_bh*drmax_bh); - const float dzmax2_inv_eh = 1.f/(dzmax_eh*dzmax_eh); - const float drmax2_inv_eh = 1.f/(drmax_eh*drmax_eh); - const float dzmax2_inv_bl = 1.f/(dzmax_bl*dzmax_bl); - const float drmax2_inv_bl = 1.f/(drmax_bl*drmax_bl); - const float dzmax2_inv_el = 1.f/(dzmax_el*dzmax_el); - const float drmax2_inv_el = 1.f/(drmax_el*drmax_el); - - // Merge hits from overlapping seeds? - // For now always true, we require extra hits after seed. - const bool merge_hits = true; // itrcfg.merge_seed_hits_during_cleaning(); - - if (seed_ptr == nullptr) return 0; - TrackVec &seeds = *seed_ptr; - - const int ns = seeds.size(); - #ifdef DEBUG - std::cout << "before seed cleaning "<< seeds.size()< writetrack(ns, true); - - const float invR1GeV = 1.f/Config::track1GeVradius; - - std::vector nHits(ns); - std::vector charge(ns); - std::vector oldPhi(ns); - std::vector pos2(ns); - std::vector eta(ns); - std::vector ctheta(ns); - std::vector invptq(ns); - std::vector pt(ns); - std::vector x(ns); - std::vector y(ns); - std::vector z(ns); - std::vector d0(ns); - int i1,i2; //for the sorting - - axis_pow2_u1 ax_phi(-Config::PI, Config::PI); - axis ax_eta(-2.6, 2.6, 30u); - - binnor b(ax_phi, ax_eta); - b.begin_registration(ns); - - for(int ts=0; ts dpt_common*(Pt1) ) - // continue; - break; // following seeds will only be farther away in pT - - ++n_comparisons; - - const float Eta2 = eta[tss]; - const float deta2 = std::pow(Eta1-Eta2, 2); - - const float oldPhi2 = oldPhi[tss]; - - const float pos2_second = pos2[tss]; - const float thisDXYSign05 = pos2_second > pos2_first ? -0.5f : 0.5f; - - const float thisDXY = thisDXYSign05*sqrt( std::pow(x[ts]-x[tss], 2) + std::pow(y[ts]-y[tss], 2) ); - - const float invptq_second = invptq[tss]; - - const float newPhi1 = oldPhi1-thisDXY*invR1GeV*invptq_first; - const float newPhi2 = oldPhi2+thisDXY*invR1GeV*invptq_second; - - const float dphi = cdist(std::abs(newPhi1-newPhi2)); - - const float dr2 = deta2+dphi*dphi; - - const float thisDZ = z[ts]-z[tss]-thisDXY*(ctheta[ts]+ctheta[tss]); - const float dz2 = thisDZ*thisDZ; - - ////// Reject tracks within dR-dz elliptical window. - ////// Adaptive thresholds, based on observation that duplicates are more abundant at large pseudo-rapidity and low track pT - bool overlapping = false; - if(std::abs(Eta1)ptmin_hpt){if(dz2*dzmax2_inv_bh+dr2*drmax2_inv_bh<1.0f) overlapping=true; } - else{if(dz2*dzmax2_inv_bl+dr2*drmax2_inv_bl<1.0f) overlapping=true; } - } - else { - if(Pt1>ptmin_hpt){if(dz2*dzmax2_inv_eh+dr2*drmax2_inv_eh<1.0f) overlapping=true; } - else{if(dz2*dzmax2_inv_el+dr2*drmax2_inv_el<1.0f) overlapping=true; } - } - - if(overlapping){ - //Mark tss as a duplicate - i1=ts; - i2=tss; - if (d0[tss]>d0[ts]) - writetrack[tss] = false; - else { - writetrack[ts] = false; - i2 = ts; - i1 = tss; - } - // Add hits from tk2 to the seed we are keeping. - // NOTE: We only have 3 bits in Track::Status for number of seed hits. - // There is a check at entry and after adding of a new hit. - Track &tk = seeds[i1]; - if (merge_hits && tk.nTotalHits() < 15) - { - const Track &tk2 = seeds[i2]; - //We are not actually fitting to the extra hits; use chi2 of 0 - float fakeChi2 = 0.0; - - for (int j = 0; j < tk2.nTotalHits(); ++j) - { - int hitidx = tk2.getHitIdx(j); - int hitlyr = tk2.getHitLyr(j); - if (hitidx >= 0) - { - bool unique = true; - for (int i = 0; i < tk.nTotalHits(); ++i) - { - if ((hitidx == tk.getHitIdx(i)) && (hitlyr == tk.getHitLyr(i))) { - unique = false; - break; - } - } - if (unique) { - tk.addHitIdx(tk2.getHitIdx(j), tk2.getHitLyr(j), fakeChi2); - ++n_ovlp_hits_added; - if (tk.nTotalHits() >= 15) - break; - } - } - } - } - if (n_ovlp_hits_added > 0) { - tk.sortHitsByLayer(); - n_ovlp_hits_added = 0; - } - - if ( ! writetrack[ts]) goto end_ts_loop; - } - } //end of inner loop over tss - } - } - - if (writetrack[ts]) - { - cleanSeedTracks.emplace_back(seeds[ts]); - } -end_ts_loop: ; - } - - seeds.swap(cleanSeedTracks); - -#ifdef DEBUG - { - const int ns2 = seeds.size(); - printf("Number of CMS seeds before %d --> after %d cleaning\n", ns, ns2); - - for (int it = 0; it < ns2; it++) - { - const Track& ss = seeds[it]; - printf(" %3i q=%+i pT=%7.3f eta=% 7.3f nHits=%i label=% i\n", - it,ss.charge(),ss.pT(),ss.momEta(),ss.nFoundHits(),ss.label()); - } - } -#endif - -#ifdef DEBUG - std::cout << "AFTER seed cleaning "<< seeds.size()<