Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

speed up duplicate and seed cleaning #82

Closed
wants to merge 5 commits into from
Closed

speed up duplicate and seed cleaning #82

Show file tree
Hide file tree
Changes from 2 commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
e25ef16
speed up duplicate and seed cleaning
leonardogiannini Feb 27, 2022
527b471
improved code style and sync with cmssw master
leonardogiannini Feb 28, 2022
414e807
improved code (except for binnor)
leonardogiannini Mar 1, 2022
f67a0ac
Rename value/pair conversions functions starting with capital letter …
osschar Mar 1, 2022
c4d7142
Incorporate comments from @srlantz.
osschar Mar 1, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
257 changes: 149 additions & 108 deletions RecoTracker/MkFitCMS/src/MkStdSeqs.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -137,6 +137,12 @@ namespace mkfit {
std::vector<float> d0(ns);
int i1, i2; //for the sorting

axis_pow2_u1<float, unsigned short, 16, 8> ax_phi(-Const::PI, Const::PI);
axis<float, unsigned short, 8, 8> ax_eta(-3.0, 3.0, 30u);
binnor<unsigned int, decltype(ax_phi), decltype(ax_eta), 24, 8> 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();
Expand All @@ -151,9 +157,17 @@ 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).

Expand All @@ -166,105 +180,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.Rrdr_to_N_bins(oldPhi[ts], 0.08f);
auto eta_rng = ax_eta.Rrdr_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]);
Expand Down Expand Up @@ -399,54 +430,64 @@ namespace mkfit {
const auto ntracks = tracks.size();

std::vector<float> ctheta(ntracks);
std::multimap<int, int> theHitMap;

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);
theHitMap.insert(std::make_pair(b * 1000 + a, i > 0 ? itrack : -itrack)); //negative for first hit in trk
}
}

for (auto itrack = 0U; itrack < ntracks; itrack++) {
auto &trk = tracks[itrack];
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<int, int> theSharingMap;
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 = theHitMap.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
theSharingMap[std::abs(it->second)]++;
}
}

auto dctheta = std::abs(ctheta[jtrack] - ctheta1);
for (const auto &elem : theSharingMap) {
auto &track2 = tracks[elem.first];

//these dctheta dphi are wasting time here -> need to remove these checks and revisit the fraction WP
auto dctheta = std::abs(ctheta[elem.first] - ctheta1);
if (dctheta > 1.)
continue;
if (dctheta > 1.)
continue;

auto dphi = std::abs(squashPhiMinimal(phi1 - track2.momPhi()));
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());
}
Expand Down
6 changes: 4 additions & 2 deletions RecoTracker/MkFitCore/interface/binnor.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -144,7 +144,7 @@ namespace mkfit {

// 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) {}
Expand Down Expand Up @@ -176,7 +176,9 @@ 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.M_bin_to_N_bin(m_bin.bin1()), m_a2.M_bin_to_N_bin(m_bin.bin2())};
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

another place to lower the case

Suggested change
return {m_a1.M_bin_to_N_bin(m_bin.bin1()), m_a2.M_bin_to_N_bin(m_bin.bin2())};
return {m_a1.m_bin_to_N_bin(m_bin.bin1()), m_a2.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}; }

Expand Down