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