File indexing completed on 2023-03-17 11:23:30
0001 #include "RecoVertex/PrimaryVertexProducer/interface/DAClusterizerInZ_vect.h"
0002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0003 #include "DataFormats/GeometryCommonDetAlgo/interface/Measurement1D.h"
0004 #include "RecoVertex/VertexPrimitives/interface/VertexException.h"
0005
0006 #include <cmath>
0007 #include <cassert>
0008 #include <limits>
0009 #include <iomanip>
0010 #include "FWCore/Utilities/interface/isFinite.h"
0011 #include "vdt/vdtMath.h"
0012
0013 using namespace std;
0014
0015
0016 #ifdef DEBUG
0017 #define DEBUGLEVEL 0
0018 #endif
0019
0020 DAClusterizerInZ_vect::DAClusterizerInZ_vect(const edm::ParameterSet& conf) {
0021
0022 maxIterations_ = 1000;
0023 mintrkweight_ = 0.5;
0024
0025
0026 #ifdef DEBUG
0027 zdumpcenter_ = conf.getUntrackedParameter<double>("zdumpcenter", 0.);
0028 zdumpwidth_ = conf.getUntrackedParameter<double>("zdumpwidth", 20.);
0029 #endif
0030
0031
0032 double Tmin = conf.getParameter<double>("Tmin");
0033 double Tpurge = conf.getParameter<double>("Tpurge");
0034 double Tstop = conf.getParameter<double>("Tstop");
0035 vertexSize_ = conf.getParameter<double>("vertexSize");
0036 coolingFactor_ = conf.getParameter<double>("coolingFactor");
0037 d0CutOff_ = conf.getParameter<double>("d0CutOff");
0038 dzCutOff_ = conf.getParameter<double>("dzCutOff");
0039 uniquetrkweight_ = conf.getParameter<double>("uniquetrkweight");
0040 uniquetrkminp_ = conf.getParameter<double>("uniquetrkminp");
0041 zmerge_ = conf.getParameter<double>("zmerge");
0042 sel_zrange_ = conf.getParameter<double>("zrange");
0043 convergence_mode_ = conf.getParameter<int>("convergence_mode");
0044 delta_lowT_ = conf.getParameter<double>("delta_lowT");
0045 delta_highT_ = conf.getParameter<double>("delta_highT");
0046 runInBlocks_ = conf.getParameter<bool>("runInBlocks");
0047 block_size_ = conf.getParameter<unsigned int>("block_size");
0048 overlap_frac_ = conf.getParameter<double>("overlap_frac");
0049
0050 #ifdef DEBUG
0051 std::cout << "DAClusterizerinZ_vect: mintrkweight = " << mintrkweight << std::endl;
0052 std::cout << "DAClusterizerinZ_vect: uniquetrkweight = " << uniquetrkweight_ << std::endl;
0053 std::cout << "DAClusterizerInZ_vect: uniquetrkminp = " << uniquetrkminp_ << std::endl;
0054 std::cout << "DAClusterizerinZ_vect: zmerge = " << zmerge_ << std::endl;
0055 std::cout << "DAClusterizerinZ_vect: Tmin = " << Tmin << std::endl;
0056 std::cout << "DAClusterizerinZ_vect: Tpurge = " << Tpurge << std::endl;
0057 std::cout << "DAClusterizerinZ_vect: Tstop = " << Tstop << std::endl;
0058 std::cout << "DAClusterizerinZ_vect: vertexSize = " << vertexSize_ << std::endl;
0059 std::cout << "DAClusterizerinZ_vect: coolingFactor = " << coolingFactor_ << std::endl;
0060 std::cout << "DAClusterizerinZ_vect: d0CutOff = " << d0CutOff_ << std::endl;
0061 std::cout << "DAClusterizerinZ_vect: dzCutOff = " << dzCutOff_ << std::endl;
0062 std::cout << "DAClusterizerInZ_vect: zrange = " << sel_zrange_ << std::endl;
0063 std::cout << "DAClusterizerinZ_vect: convergence mode = " << convergence_mode_ << std::endl;
0064 std::cout << "DAClusterizerinZ_vect: delta_highT = " << delta_highT_ << std::endl;
0065 std::cout << "DAClusterizerinZ_vect: delta_lowT = " << delta_lowT_ << std::endl;
0066 std::cout << "DAClusterizerinZ_vect: DEBUGLEVEL " << DEBUGLEVEL << std::endl;
0067 #endif
0068
0069 if (convergence_mode_ > 1) {
0070 edm::LogWarning("DAClusterizerinZ_vect")
0071 << "DAClusterizerInZ_vect: invalid convergence_mode " << convergence_mode_ << " reset to default " << 0;
0072 convergence_mode_ = 0;
0073 }
0074
0075 if (Tmin == 0) {
0076 betamax_ = 1.0;
0077 edm::LogWarning("DAClusterizerinZ_vect")
0078 << "DAClusterizerInZ_vect: invalid Tmin " << Tmin << " reset to default " << 1. / betamax_;
0079 } else {
0080 betamax_ = 1. / Tmin;
0081 }
0082
0083 if ((Tpurge > Tmin) || (Tpurge == 0)) {
0084 edm::LogWarning("DAClusterizerinZ_vect")
0085 << "DAClusterizerInZ_vect: invalid Tpurge " << Tpurge << " set to " << Tmin;
0086 Tpurge = Tmin;
0087 }
0088 betapurge_ = 1. / Tpurge;
0089
0090 if ((Tstop > Tpurge) || (Tstop == 0)) {
0091 edm::LogWarning("DAClusterizerinZ_vect")
0092 << "DAClusterizerInZ_vect: invalid Tstop " << Tstop << " set to " << max(1., Tpurge);
0093 Tstop = max(1., Tpurge);
0094 }
0095 betastop_ = 1. / Tstop;
0096 }
0097
0098 namespace {
0099 inline double local_exp(double const& inp) { return vdt::fast_exp(inp); }
0100
0101 inline void local_exp_list_range(double const* __restrict__ arg_inp,
0102 double* __restrict__ arg_out,
0103 const int kmin,
0104 const int kmax) {
0105 for (auto i = kmin; i != kmax; ++i)
0106 arg_out[i] = vdt::fast_exp(arg_inp[i]);
0107 }
0108
0109 }
0110
0111 void DAClusterizerInZ_vect::verify(const vertex_t& v, const track_t& tks, unsigned int nv, unsigned int nt) const {
0112 if (!(nv == 999999)) {
0113 assert(nv == v.getSize());
0114 } else {
0115 nv = v.getSize();
0116 }
0117
0118 if (!(nt == 999999)) {
0119 assert(nt == tks.getSize());
0120 } else {
0121 nt = tks.getSize();
0122 }
0123
0124 assert(v.zvtx_vec.size() == nv);
0125 assert(v.rho_vec.size() == nv);
0126 assert(v.swz_vec.size() == nv);
0127 assert(v.exp_arg_vec.size() == nv);
0128 assert(v.exp_vec.size() == nv);
0129 assert(v.se_vec.size() == nv);
0130 assert(v.swz_vec.size() == nv);
0131 assert(v.swE_vec.size() == nv);
0132
0133 assert(v.zvtx == &v.zvtx_vec.front());
0134 assert(v.rho == &v.rho_vec.front());
0135 assert(v.exp_arg == &v.exp_arg_vec.front());
0136 assert(v.sw == &v.sw_vec.front());
0137 assert(v.swz == &v.swz_vec.front());
0138 assert(v.se == &v.se_vec.front());
0139 assert(v.swE == &v.swE_vec.front());
0140
0141 for (unsigned int k = 0; k < nv - 1; k++) {
0142 if (v.zvtx_vec[k] <= v.zvtx_vec[k + 1])
0143 continue;
0144 cout << " Z, cluster z-ordering assertion failure z[" << k << "] =" << v.zvtx_vec[k] << " z[" << k + 1
0145 << "] =" << v.zvtx_vec[k + 1] << endl;
0146 }
0147
0148 assert(nt == tks.zpca_vec.size());
0149 assert(nt == tks.dz2_vec.size());
0150 assert(nt == tks.tt.size());
0151 assert(nt == tks.tkwt_vec.size());
0152 assert(nt == tks.sum_Z_vec.size());
0153 assert(nt == tks.kmin.size());
0154 assert(nt == tks.kmax.size());
0155
0156 assert(tks.zpca == &tks.zpca_vec.front());
0157 assert(tks.dz2 == &tks.dz2_vec.front());
0158 assert(tks.tkwt == &tks.tkwt_vec.front());
0159 assert(tks.sum_Z == &tks.sum_Z_vec.front());
0160
0161 for (unsigned int i = 0; i < nt; i++) {
0162 if ((tks.kmin[i] < tks.kmax[i]) && (tks.kmax[i] <= nv))
0163 continue;
0164 cout << "track vertex range assertion failure" << i << "/" << nt << " kmin,kmax=" << tks.kmin[i] << ", "
0165 << tks.kmax[i] << " nv=" << nv << endl;
0166 }
0167
0168 for (unsigned int i = 0; i < nt; i++) {
0169 assert((tks.kmin[i] < tks.kmax[i]) && (tks.kmax[i] <= nv));
0170 }
0171 }
0172
0173 DAClusterizerInZ_vect::track_t DAClusterizerInZ_vect::fill(const vector<reco::TransientTrack>& tracks) const {
0174
0175 track_t tks;
0176 double sumtkwt = 0.;
0177 for (auto it = tracks.begin(); it != tracks.end(); it++) {
0178 if (!(*it).isValid())
0179 continue;
0180 double t_tkwt = 1.;
0181 double t_z = ((*it).stateAtBeamLine().trackStateAtPCA()).position().z();
0182 if (std::fabs(t_z) > 1000.)
0183 continue;
0184 auto const& t_mom = (*it).stateAtBeamLine().trackStateAtPCA().momentum();
0185
0186 reco::BeamSpot beamspot = (it->stateAtBeamLine()).beamSpot();
0187 double t_dz2 = std::pow((*it).track().dzError(), 2)
0188 + (std::pow(beamspot.BeamWidthX() * t_mom.x(), 2) + std::pow(beamspot.BeamWidthY() * t_mom.y(), 2)) *
0189 std::pow(t_mom.z(), 2) / std::pow(t_mom.perp2(), 2)
0190 + std::pow(vertexSize_, 2);
0191 t_dz2 = 1. / t_dz2;
0192 if (edm::isNotFinite(t_dz2) || t_dz2 < std::numeric_limits<double>::min())
0193 continue;
0194 if (d0CutOff_ > 0) {
0195 Measurement1D atIP = (*it).stateAtBeamLine().transverseImpactParameter();
0196 t_tkwt = 1. / (1. + local_exp(std::pow(atIP.value() / atIP.error(), 2) -
0197 std::pow(d0CutOff_, 2)));
0198 if (edm::isNotFinite(t_tkwt) || t_tkwt < std::numeric_limits<double>::epsilon()) {
0199 edm::LogWarning("DAClusterizerinZ_vect") << "rejected track t_tkwt " << t_tkwt;
0200 continue;
0201 }
0202 }
0203 tks.addItemSorted(t_z, t_dz2, &(*it), t_tkwt);
0204 sumtkwt += t_tkwt;
0205 }
0206
0207 tks.extractRaw();
0208 tks.osumtkwt = sumtkwt > 0 ? 1. / sumtkwt : 0.;
0209
0210 #ifdef DEBUG
0211 if (DEBUGLEVEL > 0) {
0212 std::cout << "Track count (Z) " << tks.getSize() << std::endl;
0213 }
0214 #endif
0215
0216 return tks;
0217 }
0218
0219 namespace {
0220 inline double Eik(double t_z, double k_z, double t_dz2) { return std::pow(t_z - k_z, 2) * t_dz2; }
0221 }
0222
0223 void DAClusterizerInZ_vect::set_vtx_range(double beta, track_t& gtracks, vertex_t& gvertices) const {
0224 const unsigned int nv = gvertices.getSize();
0225 const unsigned int nt = gtracks.getSize();
0226
0227 if (nv == 0) {
0228 edm::LogWarning("DAClusterizerinZ_vect") << "empty cluster list in set_vtx_range";
0229 return;
0230 }
0231
0232 for (auto itrack = 0U; itrack < nt; ++itrack) {
0233 double zrange = max(sel_zrange_ / sqrt(beta * gtracks.dz2[itrack]), zrange_min_);
0234
0235 double zmin = gtracks.zpca[itrack] - zrange;
0236 unsigned int kmin = min(nv - 1, gtracks.kmin[itrack]);
0237
0238 if (gvertices.zvtx[kmin] > zmin) {
0239 while ((kmin > 0) && (gvertices.zvtx[kmin - 1] > zmin)) {
0240 kmin--;
0241 }
0242 } else {
0243 while ((kmin < (nv - 1)) && (gvertices.zvtx[kmin] < zmin)) {
0244 kmin++;
0245 }
0246 }
0247
0248 double zmax = gtracks.zpca[itrack] + zrange;
0249 unsigned int kmax = min(nv - 1, gtracks.kmax[itrack] - 1);
0250
0251
0252 if (gvertices.zvtx[kmax] < zmax) {
0253 while ((kmax < (nv - 1)) && (gvertices.zvtx[kmax + 1] < zmax)) {
0254 kmax++;
0255 }
0256 } else {
0257 while ((kmax > 0) && (gvertices.zvtx[kmax] > zmax)) {
0258 kmax--;
0259 }
0260 }
0261
0262 if (kmin <= kmax) {
0263 gtracks.kmin[itrack] = kmin;
0264 gtracks.kmax[itrack] = kmax + 1;
0265 } else {
0266 gtracks.kmin[itrack] = max(0U, min(kmin, kmax));
0267 gtracks.kmax[itrack] = min(nv, max(kmin, kmax) + 1);
0268 }
0269 }
0270 }
0271
0272 void DAClusterizerInZ_vect::clear_vtx_range(track_t& gtracks, vertex_t& gvertices) const {
0273 const unsigned int nt = gtracks.getSize();
0274 const unsigned int nv = gvertices.getSize();
0275 for (auto itrack = 0U; itrack < nt; ++itrack) {
0276 gtracks.kmin[itrack] = 0;
0277 gtracks.kmax[itrack] = nv;
0278 }
0279 }
0280
0281 double DAClusterizerInZ_vect::update(
0282 double beta, track_t& gtracks, vertex_t& gvertices, const double rho0, const bool updateTc) const {
0283
0284
0285
0286
0287 const unsigned int nt = gtracks.getSize();
0288 const unsigned int nv = gvertices.getSize();
0289 auto osumtkwt = gtracks.osumtkwt;
0290
0291 double Z_init = 0;
0292 if (rho0 > 0) {
0293 Z_init = rho0 * local_exp(-beta * dzCutOff_ * dzCutOff_);
0294 }
0295
0296
0297 auto kernel_calc_exp_arg_range = [beta](const unsigned int itrack,
0298 track_t const& tracks,
0299 vertex_t const& vertices,
0300 const unsigned int kmin,
0301 const unsigned int kmax) {
0302 const double track_z = tracks.zpca[itrack];
0303 const double botrack_dz2 = -beta * tracks.dz2[itrack];
0304
0305
0306 for (unsigned int ivertex = kmin; ivertex < kmax; ++ivertex) {
0307 auto mult_res = track_z - vertices.zvtx[ivertex];
0308 vertices.exp_arg[ivertex] = botrack_dz2 * (mult_res * mult_res);
0309 }
0310 };
0311
0312 auto kernel_add_Z_range = [Z_init](
0313 vertex_t const& vertices, const unsigned int kmin, const unsigned int kmax) -> double {
0314 double ZTemp = Z_init;
0315 for (unsigned int ivertex = kmin; ivertex < kmax; ++ivertex) {
0316 ZTemp += vertices.rho[ivertex] * vertices.exp[ivertex];
0317 }
0318 return ZTemp;
0319 };
0320
0321 auto kernel_calc_normalization_range = [updateTc](const unsigned int track_num,
0322 track_t& tracks,
0323 vertex_t& vertices,
0324 const unsigned int kmin,
0325 const unsigned int kmax) {
0326 auto o_trk_sum_Z = tracks.tkwt[track_num] / tracks.sum_Z[track_num];
0327 auto o_trk_dz2 = tracks.dz2[track_num];
0328 auto tmp_trk_z = tracks.zpca[track_num];
0329
0330
0331 if (updateTc) {
0332 #pragma GCC ivdep
0333 for (unsigned int k = kmin; k < kmax; ++k) {
0334 vertices.se[k] += vertices.exp[k] * o_trk_sum_Z;
0335 auto w = vertices.rho[k] * vertices.exp[k] * (o_trk_sum_Z * o_trk_dz2);
0336 vertices.sw[k] += w;
0337 vertices.swz[k] += w * tmp_trk_z;
0338 vertices.swE[k] += w * vertices.exp_arg[k];
0339 }
0340 } else {
0341
0342 #pragma GCC ivdep
0343 for (unsigned int k = kmin; k < kmax; ++k) {
0344 vertices.se[k] += vertices.exp[k] * o_trk_sum_Z;
0345 auto w = vertices.rho[k] * vertices.exp[k] * (o_trk_sum_Z * o_trk_dz2);
0346 vertices.sw[k] += w;
0347 vertices.swz[k] += w * tmp_trk_z;
0348 }
0349 }
0350 };
0351
0352 if (updateTc) {
0353 for (auto ivertex = 0U; ivertex < nv; ++ivertex) {
0354 gvertices.se[ivertex] = 0.0;
0355 gvertices.sw[ivertex] = 0.0;
0356 gvertices.swz[ivertex] = 0.0;
0357 gvertices.swE[ivertex] = 0.0;
0358 }
0359 } else {
0360 for (auto ivertex = 0U; ivertex < nv; ++ivertex) {
0361 gvertices.se[ivertex] = 0.0;
0362 gvertices.sw[ivertex] = 0.0;
0363 gvertices.swz[ivertex] = 0.0;
0364 }
0365 }
0366
0367
0368 for (auto itrack = 0U; itrack < nt; ++itrack) {
0369 const unsigned int kmin = gtracks.kmin[itrack];
0370 const unsigned int kmax = gtracks.kmax[itrack];
0371
0372 kernel_calc_exp_arg_range(itrack, gtracks, gvertices, kmin, kmax);
0373 local_exp_list_range(gvertices.exp_arg, gvertices.exp, kmin, kmax);
0374 gtracks.sum_Z[itrack] = kernel_add_Z_range(gvertices, kmin, kmax);
0375
0376 if (edm::isNotFinite(gtracks.sum_Z[itrack]))
0377 gtracks.sum_Z[itrack] = 0.0;
0378
0379 if (gtracks.sum_Z[itrack] > 1.e-100) {
0380 kernel_calc_normalization_range(itrack, gtracks, gvertices, kmin, kmax);
0381 }
0382 }
0383
0384
0385 if (updateTc) {
0386 auto obeta = -1. / beta;
0387 for (auto ivertex = 0U; ivertex < nv; ++ivertex) {
0388 gvertices.swE[ivertex] *= obeta;
0389 }
0390 }
0391
0392
0393 auto kernel_calc_z = [osumtkwt, nv](vertex_t& vertices) -> double {
0394 double delta = 0;
0395
0396 for (unsigned int ivertex = 0; ivertex < nv; ++ivertex) {
0397 if (vertices.sw[ivertex] > 0) {
0398 auto znew = vertices.swz[ivertex] / vertices.sw[ivertex];
0399
0400 delta = max(std::abs(vertices.zvtx[ivertex] - znew), delta);
0401 vertices.zvtx[ivertex] = znew;
0402 }
0403 }
0404
0405 for (unsigned int ivertex = 0; ivertex < nv; ++ivertex)
0406 vertices.rho[ivertex] = vertices.rho[ivertex] * vertices.se[ivertex] * osumtkwt;
0407
0408 return delta;
0409 };
0410
0411 double delta = kernel_calc_z(gvertices);
0412
0413
0414 return delta;
0415 }
0416
0417 unsigned int DAClusterizerInZ_vect::thermalize(
0418 double beta, track_t& tks, vertex_t& v, const double delta_max0, const double rho0) const {
0419 unsigned int niter = 0;
0420 double delta = 0;
0421 double delta_max = delta_lowT_;
0422
0423 if (convergence_mode_ == 0) {
0424 delta_max = delta_max0;
0425 } else if (convergence_mode_ == 1) {
0426 delta_max = delta_lowT_ / sqrt(std::max(beta, 1.0));
0427 }
0428
0429 set_vtx_range(beta, tks, v);
0430 double delta_sum_range = 0;
0431 std::vector<double> z0 = v.zvtx_vec;
0432
0433 while (niter++ < maxIterations_) {
0434 delta = update(beta, tks, v, rho0);
0435 delta_sum_range += delta;
0436
0437 if (delta_sum_range > zrange_min_) {
0438 for (unsigned int k = 0; k < v.getSize(); k++) {
0439 if (std::abs(v.zvtx_vec[k] - z0[k]) > zrange_min_) {
0440 set_vtx_range(beta, tks, v);
0441 delta_sum_range = 0;
0442 z0 = v.zvtx_vec;
0443 break;
0444 }
0445 }
0446 }
0447
0448 if (delta < delta_max) {
0449 break;
0450 }
0451 }
0452
0453 #ifdef DEBUG
0454 if (DEBUGLEVEL > 0) {
0455 std::cout << "DAClusterizerInZ_vect.thermalize niter = " << niter << " at T = " << 1 / beta
0456 << " nv = " << v.getSize() << std::endl;
0457 if (DEBUGLEVEL > 2)
0458 dump(beta, v, tks, 0, rho0);
0459 }
0460 #endif
0461
0462 return niter;
0463 }
0464
0465 bool DAClusterizerInZ_vect::merge(vertex_t& y, track_t& tks, double& beta) const {
0466
0467
0468
0469 const unsigned int nv = y.getSize();
0470
0471 if (nv < 2)
0472 return false;
0473
0474
0475 std::vector<std::pair<double, unsigned int>> critical;
0476 for (unsigned int k = 0; (k + 1) < nv; k++) {
0477 if (std::fabs(y.zvtx[k + 1] - y.zvtx[k]) < zmerge_) {
0478 critical.push_back(make_pair(std::fabs(y.zvtx[k + 1] - y.zvtx[k]), k));
0479 }
0480 }
0481 if (critical.empty())
0482 return false;
0483
0484 std::stable_sort(critical.begin(), critical.end(), std::less<std::pair<double, unsigned int>>());
0485
0486 for (unsigned int ik = 0; ik < critical.size(); ik++) {
0487 unsigned int k = critical[ik].second;
0488 double rho = y.rho[k] + y.rho[k + 1];
0489
0490 #ifdef DEBUG
0491 assert((k + 1) < nv);
0492 if (DEBUGLEVEL > 1) {
0493 std::cout << "merging " << fixed << setprecision(4) << y.zvtx[k + 1] << " and " << y.zvtx[k]
0494 << " sw = " << y.sw[k] + y.sw[k + 1] << std::endl;
0495 }
0496 #endif
0497
0498 if (rho > 0) {
0499 y.zvtx[k] = (y.rho[k] * y.zvtx[k] + y.rho[k + 1] * y.zvtx[k + 1]) / rho;
0500 } else {
0501 y.zvtx[k] = 0.5 * (y.zvtx[k] + y.zvtx[k + 1]);
0502 }
0503 y.rho[k] = rho;
0504 y.sw[k] += y.sw[k + 1];
0505 y.removeItem(k + 1, tks);
0506 set_vtx_range(beta, tks, y);
0507 y.extractRaw();
0508 return true;
0509 }
0510
0511 return false;
0512 }
0513
0514 bool DAClusterizerInZ_vect::purge(vertex_t& y, track_t& tks, double& rho0, const double beta) const {
0515 constexpr double eps = 1.e-100;
0516
0517 const unsigned int nv = y.getSize();
0518 const unsigned int nt = tks.getSize();
0519
0520 if (nv < 2)
0521 return false;
0522
0523 std::vector<double> sump_v(nv), arg_cache_v(nv), exp_cache_v(nv), pcut_cache_v(nv);
0524 std::vector<int> nUnique_v(nv);
0525 double* __restrict__ parg_cache;
0526 double* __restrict__ pexp_cache;
0527 double* __restrict__ ppcut_cache;
0528 double* __restrict__ psump;
0529 int* __restrict__ pnUnique;
0530 int constexpr nunique_min_ = 2;
0531
0532 set_vtx_range(beta, tks, y);
0533
0534 parg_cache = arg_cache_v.data();
0535 pexp_cache = exp_cache_v.data();
0536 ppcut_cache = pcut_cache_v.data();
0537 psump = sump_v.data();
0538 pnUnique = nUnique_v.data();
0539
0540 const auto rhoconst = rho0 * local_exp(-beta * dzCutOff_ * dzCutOff_);
0541 for (unsigned int k = 0; k < nv; k++) {
0542 const double pmax = y.rho[k] / (y.rho[k] + rhoconst);
0543 ppcut_cache[k] = uniquetrkweight_ * pmax;
0544 }
0545
0546 for (unsigned int i = 0; i < nt; i++) {
0547 const auto invZ = ((tks.sum_Z[i] > eps) && (tks.tkwt[i] > uniquetrkminp_)) ? 1. / tks.sum_Z[i] : 0.;
0548 const auto track_z = tks.zpca[i];
0549 const auto botrack_dz2 = -beta * tks.dz2[i];
0550 const auto kmin = tks.kmin[i];
0551 const auto kmax = tks.kmax[i];
0552
0553 for (unsigned int k = kmin; k < kmax; k++) {
0554 const auto mult_resz = track_z - y.zvtx[k];
0555 parg_cache[k] = botrack_dz2 * (mult_resz * mult_resz);
0556 }
0557
0558 local_exp_list_range(parg_cache, pexp_cache, kmin, kmax);
0559
0560 for (unsigned int k = kmin; k < kmax; k++) {
0561 const double p = y.rho[k] * pexp_cache[k] * invZ;
0562 psump[k] += p;
0563 pnUnique[k] += (p > ppcut_cache[k]) ? 1 : 0;
0564 }
0565 }
0566
0567 double sumpmin = nt;
0568 unsigned int k0 = nv;
0569 for (unsigned k = 0; k < nv; k++) {
0570 if ((pnUnique[k] < nunique_min_) && (psump[k] < sumpmin)) {
0571 sumpmin = psump[k];
0572 k0 = k;
0573 }
0574 }
0575
0576 if (k0 != nv) {
0577 #ifdef DEBUG
0578 assert(k0 < y.getSize());
0579 if (DEBUGLEVEL > 1) {
0580 std::cout << "eliminating prototype at " << std::setw(10) << std::setprecision(4) << y.zvtx[k0]
0581 << " with sump=" << sumpmin << " rho*nt =" << y.rho[k0] * nt << endl;
0582 }
0583 #endif
0584
0585 y.removeItem(k0, tks);
0586 set_vtx_range(beta, tks, y);
0587 return true;
0588 } else {
0589 return false;
0590 }
0591 }
0592
0593 double DAClusterizerInZ_vect::beta0(double betamax, track_t const& tks, vertex_t const& y) const {
0594 double T0 = 0;
0595
0596 const unsigned int nt = tks.getSize();
0597 const unsigned int nv = y.getSize();
0598
0599 for (unsigned int k = 0; k < nv; k++) {
0600
0601 double sumwz = 0;
0602 double sumw = 0;
0603 for (unsigned int i = 0; i < nt; i++) {
0604 double w = tks.tkwt[i] * tks.dz2[i];
0605 sumwz += w * tks.zpca[i];
0606 sumw += w;
0607 }
0608
0609 y.zvtx[k] = sumwz / sumw;
0610
0611
0612 double a = 0, b = 0;
0613 for (unsigned int i = 0; i < nt; i++) {
0614 double dx = tks.zpca[i] - y.zvtx[k];
0615 double w = tks.tkwt[i] * tks.dz2[i];
0616 a += w * std::pow(dx, 2) * tks.dz2[i];
0617 b += w;
0618 }
0619 double Tc = 2. * a / b;
0620
0621 if (Tc > T0)
0622 T0 = Tc;
0623
0624 }
0625
0626 #ifdef DEBUG
0627 if (DEBUGLEVEL > 0) {
0628 std::cout << "DAClusterizerInZ_vect.beta0: Tc = " << T0 << std::endl;
0629 int coolingsteps = 1 - int(std::log(T0 * betamax) / std::log(coolingFactor_));
0630 std::cout << "DAClusterizerInZ_vect.beta0: nstep = " << coolingsteps << std::endl;
0631 }
0632 #endif
0633
0634 if (T0 > 1. / betamax) {
0635 int coolingsteps = 1 - int(std::log(T0 * betamax) / std::log(coolingFactor_));
0636
0637 return betamax * std::pow(coolingFactor_, coolingsteps);
0638 } else {
0639
0640 return betamax * coolingFactor_;
0641 }
0642 }
0643
0644 bool DAClusterizerInZ_vect::split(const double beta, track_t& tks, vertex_t& y, double threshold) const {
0645
0646
0647
0648
0649 update(beta, tks, y, 0., true);
0650
0651 constexpr double epsilon = 1e-3;
0652 unsigned int nv = y.getSize();
0653
0654
0655
0656 std::vector<std::pair<double, unsigned int>> critical;
0657 for (unsigned int k = 0; k < nv; k++) {
0658 double Tc = 2 * y.swE[k] / y.sw[k];
0659 if (beta * Tc > threshold) {
0660 critical.push_back(make_pair(Tc, k));
0661 }
0662 }
0663 if (critical.empty())
0664 return false;
0665
0666 std::stable_sort(critical.begin(), critical.end(), std::greater<std::pair<double, unsigned int>>());
0667
0668 bool split = false;
0669 const unsigned int nt = tks.getSize();
0670
0671 for (unsigned int ic = 0; ic < critical.size(); ic++) {
0672 unsigned int k = critical[ic].second;
0673
0674
0675 double p1 = 0, z1 = 0, w1 = 0;
0676 double p2 = 0, z2 = 0, w2 = 0;
0677 for (unsigned int i = 0; i < nt; i++) {
0678 if (tks.sum_Z[i] > 1.e-100) {
0679
0680 double tl = tks.zpca[i] < y.zvtx[k] ? 1. : 0.;
0681 double tr = 1. - tl;
0682
0683
0684 double arg = (tks.zpca[i] - y.zvtx[k]) * sqrt(beta * tks.dz2[i]);
0685 if (std::fabs(arg) < 20) {
0686 double t = local_exp(-arg);
0687 tl = t / (t + 1.);
0688 tr = 1 / (t + 1.);
0689 }
0690
0691 double p = y.rho[k] * tks.tkwt[i] * local_exp(-beta * Eik(tks.zpca[i], y.zvtx[k], tks.dz2[i])) / tks.sum_Z[i];
0692 double w = p * tks.dz2[i];
0693 p1 += p * tl;
0694 z1 += w * tl * tks.zpca[i];
0695 w1 += w * tl;
0696 p2 += p * tr;
0697 z2 += w * tr * tks.zpca[i];
0698 w2 += w * tr;
0699 }
0700 }
0701
0702 if (w1 > 0) {
0703 z1 = z1 / w1;
0704 } else {
0705 z1 = y.zvtx[k] - epsilon;
0706 }
0707 if (w2 > 0) {
0708 z2 = z2 / w2;
0709 } else {
0710 z2 = y.zvtx[k] + epsilon;
0711 }
0712
0713
0714 if ((k > 0) && (z1 < (0.6 * y.zvtx[k] + 0.4 * y.zvtx[k - 1]))) {
0715 z1 = 0.6 * y.zvtx[k] + 0.4 * y.zvtx[k - 1];
0716 }
0717 if ((k + 1 < nv) && (z2 > (0.6 * y.zvtx[k] + 0.4 * y.zvtx[k + 1]))) {
0718 z2 = 0.6 * y.zvtx[k] + 0.4 * y.zvtx[k + 1];
0719 }
0720
0721 #ifdef DEBUG
0722 assert(k < nv);
0723 if (DEBUGLEVEL > 1) {
0724 if (std::fabs(y.zvtx[k] - zdumpcenter_) < zdumpwidth_) {
0725 std::cout << " T= " << std::setw(8) << 1. / beta << " Tc= " << critical[ic].first << " splitting "
0726 << std::fixed << std::setprecision(4) << y.zvtx[k] << " --> " << z1 << "," << z2 << " [" << p1
0727 << "," << p2 << "]";
0728 if (std::fabs(z2 - z1) > epsilon) {
0729 std::cout << std::endl;
0730 } else {
0731 std::cout << " rejected " << std::endl;
0732 }
0733 }
0734 }
0735 #endif
0736
0737
0738 if ((z2 - z1) > epsilon) {
0739 split = true;
0740 double pk1 = p1 * y.rho[k] / (p1 + p2);
0741 double pk2 = p2 * y.rho[k] / (p1 + p2);
0742 y.zvtx[k] = z2;
0743 y.rho[k] = pk2;
0744 y.insertItem(k, z1, pk1, tks);
0745 if (k == 0)
0746 y.extractRaw();
0747
0748 nv++;
0749
0750
0751 for (unsigned int jc = ic; jc < critical.size(); jc++) {
0752 if (critical[jc].second >= k) {
0753 critical[jc].second++;
0754 }
0755 }
0756 } else {
0757 #ifdef DEBUG
0758 std::cout << "warning ! split rejected, too small." << endl;
0759 #endif
0760 }
0761 }
0762
0763 return split;
0764 }
0765
0766 vector<TransientVertex> DAClusterizerInZ_vect::vertices(const vector<reco::TransientTrack>& tracks) const {
0767 track_t&& tks = fill(tracks);
0768 tks.extractRaw();
0769
0770 unsigned int nt = tks.getSize();
0771 double rho0 = 0.0;
0772
0773 vector<TransientVertex> clusters;
0774 if (tks.getSize() == 0)
0775 return clusters;
0776
0777 vertex_t y;
0778
0779
0780 y.addItem(0, 1.0);
0781 clear_vtx_range(tks, y);
0782
0783
0784 double beta = beta0(betamax_, tks, y);
0785 #ifdef DEBUG
0786 if (DEBUGLEVEL > 0)
0787 std::cout << "Beta0 is " << beta << std::endl;
0788 #endif
0789
0790 thermalize(beta, tks, y, delta_highT_);
0791
0792
0793
0794 double betafreeze = betamax_ * sqrt(coolingFactor_);
0795
0796 while (beta < betafreeze) {
0797 while (merge(y, tks, beta)) {
0798 update(beta, tks, y, rho0, false);
0799 }
0800 split(beta, tks, y);
0801
0802 beta = beta / coolingFactor_;
0803 thermalize(beta, tks, y, delta_highT_);
0804 }
0805
0806 #ifdef DEBUG
0807 verify(y, tks);
0808
0809 if (DEBUGLEVEL > 0) {
0810 std::cout << "DAClusterizerInZ_vect::vertices :"
0811 << "last round of splitting" << std::endl;
0812 }
0813 #endif
0814
0815 set_vtx_range(beta, tks, y);
0816 update(beta, tks, y, rho0, false);
0817
0818 while (merge(y, tks, beta)) {
0819 set_vtx_range(beta, tks, y);
0820 update(beta, tks, y, rho0, false);
0821 }
0822
0823 unsigned int ntry = 0;
0824 double threshold = 1.0;
0825 while (split(beta, tks, y, threshold) && (ntry++ < 10)) {
0826 thermalize(beta, tks, y, delta_highT_, rho0);
0827 while (merge(y, tks, beta)) {
0828 update(beta, tks, y, rho0, false);
0829 }
0830
0831
0832 threshold *= 1.1;
0833 }
0834
0835 #ifdef DEBUG
0836 verify(y, tks);
0837 if (DEBUGLEVEL > 0) {
0838 std::cout << "DAClusterizerInZ_vect::vertices :"
0839 << "turning on outlier rejection at T=" << 1 / beta << std::endl;
0840 }
0841 #endif
0842
0843
0844 if (dzCutOff_ > 0) {
0845 rho0 = y.getSize() > 1 ? 1. / y.getSize() : 1.;
0846 for (unsigned int a = 0; a < 5; a++) {
0847 update(beta, tks, y, a * rho0 / 5.);
0848 }
0849 }
0850
0851 thermalize(beta, tks, y, delta_lowT_, rho0);
0852
0853 #ifdef DEBUG
0854 verify(y, tks);
0855 if (DEBUGLEVEL > 0) {
0856 std::cout << "DAClusterizerInZ_vect::vertices :"
0857 << "merging with outlier rejection at T=" << 1 / beta << std::endl;
0858 }
0859 if (DEBUGLEVEL > 2)
0860 dump(beta, y, tks, 2, rho0);
0861 #endif
0862
0863
0864 while (merge(y, tks, beta)) {
0865 set_vtx_range(beta, tks, y);
0866 update(beta, tks, y, rho0, false);
0867 }
0868
0869 #ifdef DEBUG
0870 verify(y, tks);
0871 if (DEBUGLEVEL > 0) {
0872 std::cout << "DAClusterizerInZ_vect::vertices :"
0873 << "after merging with outlier rejection at T=" << 1 / beta << std::endl;
0874 }
0875 if (DEBUGLEVEL > 2)
0876 dump(beta, y, tks, 2, rho0);
0877 #endif
0878
0879
0880 while (beta < betapurge_) {
0881 beta = min(beta / coolingFactor_, betapurge_);
0882 thermalize(beta, tks, y, delta_lowT_, rho0);
0883 }
0884
0885 #ifdef DEBUG
0886 verify(y, tks);
0887 if (DEBUGLEVEL > 0) {
0888 std::cout << "DAClusterizerInZ_vect::vertices :"
0889 << "purging at T=" << 1 / beta << std::endl;
0890 }
0891 #endif
0892
0893
0894 while (purge(y, tks, rho0, beta)) {
0895 thermalize(beta, tks, y, delta_lowT_, rho0);
0896 }
0897
0898 #ifdef DEBUG
0899 verify(y, tks);
0900 if (DEBUGLEVEL > 0) {
0901 std::cout << "DAClusterizerInZ_vect::vertices :"
0902 << "last cooling T=" << 1 / beta << std::endl;
0903 }
0904 #endif
0905
0906
0907 while (beta < betastop_) {
0908 beta = min(beta / coolingFactor_, betastop_);
0909 thermalize(beta, tks, y, delta_lowT_, rho0);
0910 }
0911
0912 #ifdef DEBUG
0913 verify(y, tks);
0914 if (DEBUGLEVEL > 0) {
0915 std::cout << "DAClusterizerInZ_vect::vertices :"
0916 << "stop cooling at T=" << 1 / beta << std::endl;
0917 }
0918 if (DEBUGLEVEL > 2)
0919 dump(beta, y, tks, 2, rho0);
0920 #endif
0921
0922
0923
0924 set_vtx_range(beta, tks, y);
0925 const unsigned int nv = y.getSize();
0926 for (unsigned int k = 0; k < nv; k++) {
0927 if (edm::isNotFinite(y.rho[k]) || edm::isNotFinite(y.zvtx[k])) {
0928 y.rho[k] = 0;
0929 y.zvtx[k] = 0;
0930 }
0931 }
0932
0933 const auto z_sum_init = rho0 * local_exp(-beta * dzCutOff_ * dzCutOff_);
0934 std::vector<std::vector<unsigned int>> vtx_track_indices(nv);
0935 for (unsigned int i = 0; i < nt; i++) {
0936 const auto kmin = tks.kmin[i];
0937 const auto kmax = tks.kmax[i];
0938 for (auto k = kmin; k < kmax; k++) {
0939 y.exp_arg[k] = -beta * Eik(tks.zpca[i], y.zvtx[k], tks.dz2[i]);
0940 }
0941
0942 local_exp_list_range(y.exp_arg, y.exp, kmin, kmax);
0943
0944 tks.sum_Z[i] = z_sum_init;
0945 for (auto k = kmin; k < kmax; k++) {
0946 tks.sum_Z[i] += y.rho[k] * y.exp[k];
0947 }
0948 const double invZ = tks.sum_Z[i] > 1e-100 ? 1. / tks.sum_Z[i] : 0.0;
0949
0950 for (auto k = kmin; k < kmax; k++) {
0951 double p = y.rho[k] * y.exp[k] * invZ;
0952 if (p > mintrkweight_) {
0953
0954 vtx_track_indices[k].push_back(i);
0955 break;
0956 }
0957 }
0958
0959 }
0960
0961 GlobalError dummyError(0.01, 0, 0.01, 0., 0., 0.01);
0962 for (unsigned int k = 0; k < nv; k++) {
0963 if (!vtx_track_indices[k].empty()) {
0964 GlobalPoint pos(0, 0, y.zvtx[k]);
0965 vector<reco::TransientTrack> vertexTracks;
0966 for (auto i : vtx_track_indices[k]) {
0967 vertexTracks.push_back(*(tks.tt[i]));
0968 }
0969 TransientVertex v(pos, dummyError, vertexTracks, 0);
0970 clusters.push_back(v);
0971 }
0972 }
0973
0974 return clusters;
0975 }
0976
0977 vector<TransientVertex> DAClusterizerInZ_vect::vertices_in_blocks(const vector<reco::TransientTrack>& tracks) const {
0978 vector<reco::TransientTrack> sorted_tracks;
0979 vector<pair<float, float>> vertices_tot;
0980 for (unsigned int i = 0; i < tracks.size(); i++) {
0981 sorted_tracks.push_back(tracks[i]);
0982 }
0983 double rho0, beta;
0984 std::sort(sorted_tracks.begin(),
0985 sorted_tracks.end(),
0986 [](const reco::TransientTrack& a, const reco::TransientTrack& b) -> bool {
0987 return (a.stateAtBeamLine().trackStateAtPCA()).position().z() <
0988 (b.stateAtBeamLine().trackStateAtPCA()).position().z();
0989 });
0990
0991 unsigned int nBlocks = (unsigned int)std::floor(sorted_tracks.size() / (block_size_ * (1 - overlap_frac_)));
0992 if (nBlocks < 1) {
0993 nBlocks = 1;
0994 edm::LogWarning("DAClusterizerinZ_vect")
0995 << "Warning nBlocks was 0 with ntracks = " << sorted_tracks.size() << " block_size = " << block_size_
0996 << " and overlap fraction = " << overlap_frac_ << ". Setting nBlocks = 1";
0997 }
0998 for (unsigned int block = 0; block < nBlocks; block++) {
0999 vector<reco::TransientTrack> block_tracks;
1000 unsigned int begin = (unsigned int)(block * block_size_ * (1 - overlap_frac_));
1001 unsigned int end = (unsigned int)std::min(begin + block_size_, (unsigned int)sorted_tracks.size());
1002 for (unsigned int i = begin; i < end; i++) {
1003 block_tracks.push_back(sorted_tracks[i]);
1004 }
1005 if (block_tracks.empty()) {
1006 continue;
1007 }
1008
1009 #ifdef DEBUG
1010 std::cout << "Running vertices_in_blocks on" << std::endl;
1011 std::cout << "- block no." << block << " on " << nBlocks << " blocks " << std::endl;
1012 std::cout << "- block track size: " << sorted_tracks.size() << " - block size: " << block_size_ << std::endl;
1013 #endif
1014 track_t&& tks = fill(block_tracks);
1015 tks.extractRaw();
1016
1017 rho0 = 0.0;
1018
1019 vertex_t y;
1020
1021
1022 y.addItem(0, 1.0);
1023 clear_vtx_range(tks, y);
1024
1025
1026 beta = beta0(betamax_, tks, y);
1027 #ifdef DEBUG
1028 if (DEBUGLEVEL > 0)
1029 std::cout << "Beta0 is " << beta << std::endl;
1030 #endif
1031
1032 thermalize(beta, tks, y, delta_highT_);
1033
1034
1035
1036 double betafreeze = betamax_ * sqrt(coolingFactor_);
1037 while (beta < betafreeze) {
1038 while (merge(y, tks, beta)) {
1039 update(beta, tks, y, rho0, false);
1040 }
1041 split(beta, tks, y);
1042
1043 beta = beta / coolingFactor_;
1044 thermalize(beta, tks, y, delta_highT_);
1045 }
1046
1047 #ifdef DEBUG
1048 verify(y, tks);
1049
1050 if (DEBUGLEVEL > 0) {
1051 std::cout << "DAClusterizerInZSubCluster_vect::vertices :"
1052 << "last round of splitting" << std::endl;
1053 }
1054 #endif
1055
1056 set_vtx_range(beta, tks, y);
1057 update(beta, tks, y, rho0, false);
1058
1059 while (merge(y, tks, beta)) {
1060 set_vtx_range(beta, tks, y);
1061 update(beta, tks, y, rho0, false);
1062 }
1063
1064 unsigned int ntry = 0;
1065 double threshold = 1.0;
1066 while (split(beta, tks, y, threshold) && (ntry++ < 10)) {
1067 thermalize(beta, tks, y, delta_highT_, rho0);
1068 while (merge(y, tks, beta)) {
1069 update(beta, tks, y, rho0, false);
1070 }
1071
1072
1073 threshold *= 1.1;
1074 }
1075
1076 #ifdef DEBUG
1077 verify(y, tks);
1078 if (DEBUGLEVEL > 0) {
1079 std::cout << "DAClusterizerInZSubCluster_vect::vertices :"
1080 << "turning on outlier rejection at T=" << 1 / beta << std::endl;
1081 }
1082 #endif
1083
1084
1085 if (dzCutOff_ > 0) {
1086 rho0 = y.getSize() > 1 ? 1. / y.getSize() : 1.;
1087 for (unsigned int a = 0; a < 5; a++) {
1088 update(beta, tks, y, a * rho0 / 5.);
1089 }
1090 }
1091
1092 thermalize(beta, tks, y, delta_lowT_, rho0);
1093
1094 #ifdef DEBUG
1095 verify(y, tks);
1096 if (DEBUGLEVEL > 0) {
1097 std::cout << "DAClusterizerInZSubCluster_vect::vertices :"
1098 << "merging with outlier rejection at T=" << 1 / beta << std::endl;
1099 }
1100 if (DEBUGLEVEL > 2)
1101 dump(beta, y, tks, 2, rho0);
1102 #endif
1103
1104
1105 while (merge(y, tks, beta)) {
1106 set_vtx_range(beta, tks, y);
1107 update(beta, tks, y, rho0, false);
1108 }
1109
1110 #ifdef DEBUG
1111 verify(y, tks);
1112 if (DEBUGLEVEL > 0) {
1113 std::cout << "DAClusterizerInZSubCluster_vect::vertices :"
1114 << "after merging with outlier rejection at T=" << 1 / beta << std::endl;
1115 }
1116 if (DEBUGLEVEL > 2)
1117 dump(beta, y, tks, 2, rho0);
1118 #endif
1119
1120
1121 while (beta < betapurge_) {
1122 beta = min(beta / coolingFactor_, betapurge_);
1123 thermalize(beta, tks, y, delta_lowT_, rho0);
1124 }
1125
1126 #ifdef DEBUG
1127 verify(y, tks);
1128 if (DEBUGLEVEL > 0) {
1129 std::cout << "DAClusterizerInZSubCluster_vect::vertices :"
1130 << "purging at T=" << 1 / beta << std::endl;
1131 }
1132 #endif
1133
1134
1135 while (purge(y, tks, rho0, beta)) {
1136 thermalize(beta, tks, y, delta_lowT_, rho0);
1137 }
1138
1139 #ifdef DEBUG
1140 verify(y, tks);
1141 if (DEBUGLEVEL > 0) {
1142 std::cout << "DAClusterizerInZSubCluster_vect::vertices :"
1143 << "last cooling T=" << 1 / beta << std::endl;
1144 }
1145 #endif
1146
1147
1148 while (beta < betastop_) {
1149 beta = min(beta / coolingFactor_, betastop_);
1150 thermalize(beta, tks, y, delta_lowT_, rho0);
1151 }
1152
1153 #ifdef DEBUG
1154 verify(y, tks);
1155 if (DEBUGLEVEL > 0) {
1156 std::cout << "DAClusterizerInZSubCluster_vect::vertices :"
1157 << "stop cooling at T=" << 1 / beta << std::endl;
1158 }
1159 if (DEBUGLEVEL > 2)
1160 dump(beta, y, tks, 2, rho0);
1161 #endif
1162
1163 for (unsigned int ivertex = 0; ivertex < y.getSize(); ivertex++) {
1164 if (y.zvtx_vec[ivertex] != 0 && y.rho_vec[ivertex] != 0) {
1165 vertices_tot.push_back(pair(y.zvtx_vec[ivertex], y.rho_vec[ivertex]));
1166 #ifdef DEBUG
1167 std::cout << "Found new vertex " << y.zvtx_vec[ivertex] << " , " << y.rho_vec[ivertex] << std::endl;
1168 #endif
1169 }
1170 }
1171 }
1172
1173 std::sort(vertices_tot.begin(),
1174 vertices_tot.end(),
1175 [](const pair<float, float>& a, const pair<float, float>& b) -> bool { return a.first < b.first; });
1176
1177
1178 track_t&& tracks_tot = fill(tracks);
1179 const unsigned int nv = vertices_tot.size();
1180 const unsigned int nt = tracks_tot.getSize();
1181
1182 for (auto itrack = 0U; itrack < nt; ++itrack) {
1183 double zrange = max(sel_zrange_ / sqrt(beta * tracks_tot.dz2[itrack]), zrange_min_);
1184
1185 double zmin = tracks_tot.zpca[itrack] - zrange;
1186 unsigned int kmin = min(nv - 1, tracks_tot.kmin[itrack]);
1187
1188 if (vertices_tot[kmin].first > zmin) {
1189 while ((kmin > 0) && (vertices_tot[kmin - 1].first > zmin)) {
1190 kmin--;
1191 }
1192 } else {
1193 while ((kmin < (nv - 1)) && (vertices_tot[kmin].first < zmin)) {
1194 kmin++;
1195 }
1196 }
1197
1198 double zmax = tracks_tot.zpca[itrack] + zrange;
1199 unsigned int kmax = min(nv - 1, tracks_tot.kmax[itrack] - 1);
1200
1201
1202 if (vertices_tot[kmax].first < zmax) {
1203 while ((kmax < (nv - 1)) && (vertices_tot[kmax + 1].first < zmax)) {
1204 kmax++;
1205 }
1206 } else {
1207 while ((kmax > 0) && (vertices_tot[kmax].first > zmax)) {
1208 kmax--;
1209 }
1210 }
1211
1212 if (kmin <= kmax) {
1213 tracks_tot.kmin[itrack] = kmin;
1214 tracks_tot.kmax[itrack] = kmax + 1;
1215 } else {
1216 tracks_tot.kmin[itrack] = max(0U, min(kmin, kmax));
1217 tracks_tot.kmax[itrack] = min(nv, max(kmin, kmax) + 1);
1218 }
1219 }
1220
1221 rho0 = nv > 1 ? 1. / nv : 1.;
1222 const auto z_sum_init = rho0 * local_exp(-beta * dzCutOff_ * dzCutOff_);
1223
1224 std::vector<std::vector<unsigned int>> vtx_track_indices(nv);
1225 for (unsigned int i = 0; i < nt; i++) {
1226 const auto kmin = tracks_tot.kmin[i];
1227 const auto kmax = tracks_tot.kmax[i];
1228 double p_max = -1;
1229 unsigned int iMax = 10000;
1230 float sum_Z = z_sum_init;
1231 for (auto k = kmin; k < kmax; k++) {
1232 float v_exp = local_exp(-beta * Eik(tracks_tot.zpca[i], vertices_tot[k].first, tracks_tot.dz2[i]));
1233 sum_Z += vertices_tot[k].second * v_exp;
1234 }
1235 double invZ = sum_Z > 1e-100 ? 1. / sum_Z : 0.0;
1236 for (auto k = kmin; k < kmax && invZ != 0.0; k++) {
1237 float v_exp = local_exp(-beta * Eik(tracks_tot.zpca[i], vertices_tot[k].first, tracks_tot.dz2[i]));
1238 double p = vertices_tot[k].second * v_exp * invZ;
1239 if (p > p_max && p > mintrkweight_) {
1240 p_max = p;
1241 iMax = k;
1242 }
1243 }
1244 if (iMax < vtx_track_indices.size())
1245 vtx_track_indices[iMax].push_back(i);
1246 }
1247 #ifdef DEBUG
1248 for (auto itrack = 0U; itrack < nt; ++itrack) {
1249 std::cout << "itrack " << itrack << " , " << tracks_tot.kmin[itrack] << " , " << tracks_tot.kmax[itrack]
1250 << std::endl;
1251 }
1252 #endif
1253
1254 vector<TransientVertex> clusters;
1255
1256 GlobalError dummyError(0.01, 0, 0.01, 0., 0., 0.01);
1257 for (unsigned int k = 0; k < nv; k++) {
1258 if (!vtx_track_indices[k].empty()) {
1259 GlobalPoint pos(0, 0, vertices_tot[k].first);
1260 vector<reco::TransientTrack> vertexTracks;
1261 for (auto i : vtx_track_indices[k]) {
1262 vertexTracks.push_back(*(tracks_tot.tt[i]));
1263 #ifdef DEBUG
1264 std::cout << y.zvtx[k] << "," << (*tks.tt[i]).stateAtBeamLine().trackStateAtPCA().position().z() << std::endl;
1265 #endif
1266 }
1267 TransientVertex v(pos, dummyError, vertexTracks, 0);
1268 clusters.push_back(v);
1269 }
1270 }
1271
1272 return clusters;
1273 }
1274
1275 vector<vector<reco::TransientTrack>> DAClusterizerInZ_vect::clusterize(
1276 const vector<reco::TransientTrack>& tracks) const {
1277 vector<vector<reco::TransientTrack>> clusters;
1278
1279 vector<TransientVertex> pv;
1280 if (runInBlocks_ and (block_size_ < tracks.size()))
1281 pv = vertices_in_blocks(tracks);
1282 else
1283 pv = vertices(tracks);
1284
1285 #ifdef DEBUG
1286 if (DEBUGLEVEL > 0) {
1287 std::cout << "###################################################" << endl;
1288 std::cout << "# vectorized DAClusterizerInZ_vect::clusterize nt=" << tracks.size() << endl;
1289 std::cout << "# DAClusterizerInZ_vect::clusterize pv.size=" << pv.size() << endl;
1290 std::cout << "###################################################" << endl;
1291 }
1292 #endif
1293
1294 if (pv.empty()) {
1295 return clusters;
1296 }
1297
1298
1299 vector<reco::TransientTrack> aCluster = pv.begin()->originalTracks();
1300
1301 for (auto k = pv.begin() + 1; k != pv.end(); k++) {
1302 if (std::abs(k->position().z() - (k - 1)->position().z()) > (2 * vertexSize_)) {
1303
1304 if (aCluster.size() > 1) {
1305 clusters.push_back(aCluster);
1306 }
1307 #ifdef DEBUG
1308 else {
1309 std::cout << " one track cluster at " << k->position().z() << " suppressed" << std::endl;
1310 }
1311 #endif
1312 aCluster.clear();
1313 }
1314 for (unsigned int i = 0; i < k->originalTracks().size(); i++) {
1315 aCluster.push_back(k->originalTracks()[i]);
1316 }
1317 }
1318 clusters.emplace_back(std::move(aCluster));
1319
1320 return clusters;
1321 }
1322
1323 void DAClusterizerInZ_vect::dump(
1324 const double beta, const vertex_t& y, const track_t& tks, const int verbosity, const double rho0) const {
1325 #ifdef DEBUG
1326 const unsigned int nv = y.getSize();
1327 const unsigned int nt = tks.getSize();
1328
1329 vertex_t y_local = y;
1330 track_t tks_local = tks;
1331 update(beta, tks_local, y_local, rho0, true);
1332
1333 std::vector<unsigned int> iz;
1334 for (unsigned int j = 0; j < nt; j++) {
1335 iz.push_back(j);
1336 }
1337 std::sort(iz.begin(), iz.end(), [tks](unsigned int a, unsigned int b) { return tks.zpca[a] < tks.zpca[b]; });
1338 std::cout << std::endl;
1339 std::cout << "-----DAClusterizerInZ::dump ----" << nv << " clusters " << std::endl;
1340 std::cout << " ";
1341 for (unsigned int ivertex = 0; ivertex < nv; ++ivertex) {
1342 if (std::fabs(y.zvtx[ivertex] - zdumpcenter_) < zdumpwidth_) {
1343 std::cout << " " << setw(3) << ivertex << " ";
1344 }
1345 }
1346 std::cout << endl;
1347 std::cout << " z= ";
1348 std::cout << setprecision(4);
1349 for (unsigned int ivertex = 0; ivertex < nv; ++ivertex) {
1350 if (std::fabs(y.zvtx[ivertex] - zdumpcenter_) < zdumpwidth_) {
1351 std::cout << setw(8) << fixed << y.zvtx[ivertex];
1352 }
1353 }
1354 std::cout << endl
1355 << "T=" << setw(15) << 1. / beta << " Tmin =" << setw(10) << 1. / betamax_
1356 << " Tc= ";
1357 for (unsigned int ivertex = 0; ivertex < nv; ++ivertex) {
1358 if (std::fabs(y.zvtx[ivertex] - zdumpcenter_) < zdumpwidth_) {
1359 double Tc = 2 * y_local.swE[ivertex] / y_local.sw[ivertex];
1360 std::cout << setw(8) << fixed << setprecision(1) << Tc;
1361 }
1362 }
1363 std::cout << endl;
1364
1365 std::cout << " pk= ";
1366 double sumpk = 0;
1367 for (unsigned int ivertex = 0; ivertex < nv; ++ivertex) {
1368 sumpk += y.rho[ivertex];
1369 if (std::fabs(y.zvtx[ivertex] - zdumpcenter_) > zdumpwidth_)
1370 continue;
1371 std::cout << setw(8) << setprecision(4) << fixed << y.rho[ivertex];
1372 }
1373 std::cout << endl;
1374
1375 std::cout << " nt= ";
1376 for (unsigned int ivertex = 0; ivertex < nv; ++ivertex) {
1377 if (std::fabs(y.zvtx[ivertex] - zdumpcenter_) > zdumpwidth_)
1378 continue;
1379 std::cout << setw(8) << setprecision(1) << fixed << y.rho[ivertex] * nt;
1380 }
1381 std::cout << endl;
1382
1383 if (verbosity > 0) {
1384 double E = 0, F = 0;
1385 std::cout << endl;
1386 std::cout << "---- z +/- dz ip +/-dip pt phi eta weights ----" << endl;
1387 std::cout << setprecision(4);
1388 for (unsigned int i0 = 0; i0 < nt; i0++) {
1389 unsigned int i = iz[i0];
1390 if (tks.sum_Z[i] > 0) {
1391 F -= std::log(tks.sum_Z[i]) / beta;
1392 }
1393 double tz = tks.zpca[i];
1394
1395 if (std::fabs(tz - zdumpcenter_) > zdumpwidth_)
1396 continue;
1397 std::cout << setw(4) << i << ")" << setw(8) << fixed << setprecision(4) << tz << " +/-" << setw(6)
1398 << sqrt(1. / tks.dz2[i]);
1399 if ((tks.tt[i] == nullptr)) {
1400 std::cout << " effective track ";
1401 } else {
1402 if (tks.tt[i]->track().quality(reco::TrackBase::highPurity)) {
1403 std::cout << " *";
1404 } else {
1405 std::cout << " ";
1406 }
1407 if (tks.tt[i]->track().hitPattern().hasValidHitInPixelLayer(PixelSubdetector::SubDetector::PixelBarrel, 1)) {
1408 std::cout << "+";
1409 } else {
1410 std::cout << "-";
1411 }
1412 std::cout << setw(1)
1413 << tks.tt[i]
1414 ->track()
1415 .hitPattern()
1416 .pixelBarrelLayersWithMeasurement();
1417 std::cout << setw(1) << tks.tt[i]->track().hitPattern().pixelEndcapLayersWithMeasurement();
1418 std::cout << setw(1) << hex
1419 << tks.tt[i]->track().hitPattern().trackerLayersWithMeasurement() -
1420 tks.tt[i]->track().hitPattern().pixelLayersWithMeasurement()
1421 << dec;
1422 std::cout << "=" << setw(1) << hex
1423 << tks.tt[i]->track().hitPattern().numberOfLostHits(reco::HitPattern::MISSING_OUTER_HITS) << dec;
1424
1425 Measurement1D IP = tks.tt[i]->stateAtBeamLine().transverseImpactParameter();
1426 std::cout << setw(8) << IP.value() << "+/-" << setw(6) << IP.error();
1427 std::cout << " " << setw(6) << setprecision(2) << tks.tt[i]->track().pt() * tks.tt[i]->track().charge();
1428 std::cout << " " << setw(5) << setprecision(2) << tks.tt[i]->track().phi() << " " << setw(5) << setprecision(2)
1429 << tks.tt[i]->track().eta();
1430 }
1431
1432 double sump = 0.;
1433 for (unsigned int ivertex = 0; ivertex < nv; ++ivertex) {
1434 if (std::fabs(y.zvtx[ivertex] - zdumpcenter_) > zdumpwidth_)
1435 continue;
1436
1437 if ((tks.tkwt[i] > 0) && (tks.sum_Z[i] > 0)) {
1438
1439 double p = y.rho[ivertex] * local_exp(-beta * Eik(tks.zpca[i], y.zvtx[ivertex], tks.dz2[i])) / tks.sum_Z[i];
1440 if (p > 0.0001) {
1441 std::cout << setw(8) << setprecision(3) << p;
1442 } else {
1443 std::cout << " . ";
1444 }
1445 E += p * Eik(tks.zpca[i], y.zvtx[ivertex], tks.dz2[i]);
1446 sump += p;
1447 } else {
1448 std::cout << " ";
1449 }
1450 }
1451 std::cout << " ( " << std::setw(3) << tks.kmin[i] << "," << std::setw(3) << tks.kmax[i] - 1 << " ) ";
1452 std::cout << endl;
1453 }
1454 std::cout << " ";
1455 for (unsigned int ivertex = 0; ivertex < nv; ++ivertex) {
1456 if (std::fabs(y.zvtx[ivertex] - zdumpcenter_) < zdumpwidth_) {
1457 std::cout << " " << setw(3) << ivertex << " ";
1458 }
1459 }
1460 std::cout << endl;
1461 std::cout << " z= ";
1462 std::cout << setprecision(4);
1463 for (unsigned int ivertex = 0; ivertex < nv; ++ivertex) {
1464 if (std::fabs(y.zvtx[ivertex] - zdumpcenter_) < zdumpwidth_) {
1465 std::cout << setw(8) << fixed << y.zvtx[ivertex];
1466 }
1467 }
1468 std::cout << endl;
1469 std::cout << endl
1470 << "T=" << 1 / beta << " E=" << E << " n=" << y.getSize() << " F= " << F << endl
1471 << "----------" << endl;
1472 }
1473 #endif
1474 }
1475
1476 void DAClusterizerInZ_vect::fillPSetDescription(edm::ParameterSetDescription& desc) {
1477 desc.addUntracked<double>("zdumpcenter", 0.);
1478 desc.addUntracked<double>("zdumpwidth", 20.);
1479 desc.add<double>("d0CutOff", 3.0);
1480 desc.add<double>("Tmin", 2.0);
1481 desc.add<double>("delta_lowT", 0.001);
1482 desc.add<double>("zmerge", 0.01);
1483 desc.add<double>("dzCutOff", 3.0);
1484 desc.add<double>("Tpurge", 2.0);
1485 desc.add<int>("convergence_mode", 0);
1486 desc.add<double>("delta_highT", 0.01);
1487 desc.add<double>("Tstop", 0.5);
1488 desc.add<double>("coolingFactor", 0.6);
1489 desc.add<double>("vertexSize", 0.006);
1490 desc.add<double>("uniquetrkweight", 0.8);
1491 desc.add<double>("uniquetrkminp", 0.0);
1492 desc.add<double>("zrange", 4.0);
1493 desc.add<bool>("runInBlocks", false);
1494 desc.add<unsigned int>("block_size", 512);
1495 desc.add<double>("overlap_frac", 0.5);
1496 }