Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-04-11 22:44:13

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