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

Jump to bottom

modified merging scheme and tighter outlier suppression in primary vertex reconstruction #33341

Merged
merged 2 commits into from Apr 10, 2021
Merged

modified merging scheme and tighter outlier suppression in primary vertex reconstruction #33341

Changes from 1 commit
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
71b1f6f
remove Tc control for merging and make the outlier suppression term p…
werdmann Apr 6, 2021
1107cea
fix DEBUG mode code that involves Tc
werdmann Apr 9, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
56 changes: 25 additions & 31 deletions RecoVertex/PrimaryVertexProducer/src/DAClusterizerInZ_vect.cc
View file
Open in desktop
Expand Up @@ -483,31 +483,26 @@ bool DAClusterizerInZ_vect::merge(vertex_t& y, track_t& tks, double& beta) const
for (unsigned int ik = 0; ik < critical.size(); ik++) {
unsigned int k = critical[ik].second;
double rho = y.rho[k] + y.rho[k + 1];
double swE = y.swE[k] + y.swE[k + 1] - y.rho[k] * y.rho[k + 1] / rho * std::pow(y.zvtx[k + 1] - y.zvtx[k], 2);
double Tc = 2 * swE / (y.sw[k] + y.sw[k + 1]);
Copy link
Contributor

Choose a reason for hiding this comment

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

Since you are removing the local variable Tc, please also remove it from the cout at L490

Copy link
Contributor Author

Choose a reason for hiding this comment

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

thanks for catching. I'll also modify the "dump" routine, which is only used when DEBUG is enabled, to print meaningful Tc.


if (Tc * beta < 1) {
#ifdef DEBUG
assert((k + 1) < nv);
if (DEBUGLEVEL > 1) {
std::cout << "merging " << fixed << setprecision(4) << y.zvtx[k + 1] << " and " << y.zvtx[k] << " Tc = " << Tc
<< " sw = " << y.sw[k] + y.sw[k + 1] << std::endl;
}
assert((k + 1) < nv);
if (DEBUGLEVEL > 1) {
std::cout << "merging " << fixed << setprecision(4) << y.zvtx[k + 1] << " and " << y.zvtx[k] << " Tc = " << Tc
<< " sw = " << y.sw[k] + y.sw[k + 1] << std::endl;
}
#endif

if (rho > 0) {
y.zvtx[k] = (y.rho[k] * y.zvtx[k] + y.rho[k + 1] * y.zvtx[k + 1]) / rho;
} else {
y.zvtx[k] = 0.5 * (y.zvtx[k] + y.zvtx[k + 1]);
}
y.rho[k] = rho;
y.sw[k] += y.sw[k + 1];
y.swE[k] = swE;
y.removeItem(k + 1, tks);
set_vtx_range(beta, tks, y);
y.extractRaw();
return true;
if (rho > 0) {
y.zvtx[k] = (y.rho[k] * y.zvtx[k] + y.rho[k + 1] * y.zvtx[k + 1]) / rho;
} else {
y.zvtx[k] = 0.5 * (y.zvtx[k] + y.zvtx[k + 1]);
}
y.rho[k] = rho;
y.sw[k] += y.sw[k + 1];
y.removeItem(k + 1, tks);
set_vtx_range(beta, tks, y);
y.extractRaw();
return true;
}

return false;
Expand Down Expand Up @@ -645,9 +640,11 @@ double DAClusterizerInZ_vect::beta0(double betamax, track_t const& tks, vertex_t

bool DAClusterizerInZ_vect::split(const double beta, track_t& tks, vertex_t& y, double threshold) const {
// split only critical vertices (Tc >~ T=1/beta <==> beta*Tc>~1)
// an update must have been made just before doing this (same beta, no merging)
// returns true if at least one cluster was split

// update the sums needed for Tc, rho0 is never non-zero while splitting is still active
update(beta, tks, y, 0., true); // the "true" flag enables Tc evaluation

constexpr double epsilon = 1e-3; // minimum split size
unsigned int nv = y.getSize();

Expand Down Expand Up @@ -794,9 +791,8 @@ vector<TransientVertex> DAClusterizerInZ_vect::vertices(const vector<reco::Trans
double betafreeze = betamax_ * sqrt(coolingFactor_);

while (beta < betafreeze) {
update(beta, tks, y, rho0, true);
while (merge(y, tks, beta)) {
update(beta, tks, y, rho0, true);
update(beta, tks, y, rho0, false);
}
split(beta, tks, y);

Expand All @@ -814,20 +810,19 @@ vector<TransientVertex> DAClusterizerInZ_vect::vertices(const vector<reco::Trans
#endif

set_vtx_range(beta, tks, y);
update(beta, tks, y, rho0, true); // merge() and split() use Tc
update(beta, tks, y, rho0, false);

while (merge(y, tks, beta)) {
set_vtx_range(beta, tks, y);
update(beta, tks, y, rho0, true);
update(beta, tks, y, rho0, false);
}

unsigned int ntry = 0;
double threshold = 1.0;
while (split(beta, tks, y, threshold) && (ntry++ < 10)) {
thermalize(beta, tks, y, delta_highT_, rho0); // rho0 = 0. here
update(beta, tks, y, rho0, true);
while (merge(y, tks, beta)) {
update(beta, tks, y, rho0, true);
update(beta, tks, y, rho0, false);
}

// relax splitting a bit to reduce multiple split-merge cycles of the same cluster
Expand All @@ -842,16 +837,15 @@ vector<TransientVertex> DAClusterizerInZ_vect::vertices(const vector<reco::Trans
}
#endif

// switch on outlier rejection at T=Tmin, doesn't do much at high PU
// switch on outlier rejection at T=Tmin
if (dzCutOff_ > 0) {
rho0 = 1. / nt;
rho0 = y.getSize() > 1 ? 1. / y.getSize() : 1.;
for (unsigned int a = 0; a < 5; a++) {
update(beta, tks, y, a * rho0 / 5.); // adiabatic turn-on
}
}

thermalize(beta, tks, y, delta_lowT_, rho0);
update(beta, tks, y, rho0, true);

#ifdef DEBUG
verify(y, tks);
Expand All @@ -866,7 +860,7 @@ vector<TransientVertex> DAClusterizerInZ_vect::vertices(const vector<reco::Trans
// merge again (some cluster split by outliers collapse here)
while (merge(y, tks, beta)) {
set_vtx_range(beta, tks, y);
update(beta, tks, y, rho0, true);
update(beta, tks, y, rho0, false);
}

#ifdef DEBUG
Expand Down