Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-08-13 22:17:07 

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

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!