Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-03-17 11:23:30 

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