PrimaryVertexProducer/src/AdaptiveChisquarePrimaryVertexFitter.cc

0001 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0002 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0003 #include "RecoVertex/PrimaryVertexProducer/interface/AdaptiveChisquarePrimaryVertexFitter.h"
0004 #include "vdt/vdtMath.h"
0005
0006 //#define PVTX_DEBUG
0007
0008 AdaptiveChisquarePrimaryVertexFitter::AdaptiveChisquarePrimaryVertexFitter(const double chicutoff,
0009                                                                            const double zcutoff,
0010                                                                            const double mintrkweight,
0011                                                                            const bool multivertexfit)
0012     : chi_cutoff_(chicutoff), z_cutoff_(zcutoff), min_trackweight_(mintrkweight), multivertexfit_(multivertexfit) {
0013 #ifdef PVTX_DEBUG
0014   edm::LogInfo("AdaptiveChisquarePrimaryVertexFitter")
0015       << "instantiating AdaptiveChisquarePrimaryVertexFitter with "
0016       << " chi cutoff =" << chi_cutoff_ << " z cutoff =" << z_cutoff_ << " min_trackweight=" << min_trackweight_
0017       << " multivertex mode " << multivertexfit_ << std::endl;
0018 #endif
0019 }
0020
0021 double AdaptiveChisquarePrimaryVertexFitter::track_in_vertex_chsq(const TrackInfo &ti,
0022                                                                   const double xvtx,
0023                                                                   const double yvtx,
0024                                                                   const double zvtx) {
0025   double F1 = ti.b1 + xvtx * ti.H1[0] + yvtx * ti.H1[1];         // using H1[2]=0
0026   double F2 = ti.b2 + xvtx * ti.H2[0] + yvtx * ti.H2[1] - zvtx;  // using H2[2]=-1
0027   double chsq = F1 * F1 * ti.S11 + F2 * F2 * ti.S22 + 2. * F1 * F2 * ti.S12;
0028 #ifdef PVTX_DEBUG
0029   assert((chsq >= 0) && " negative chi**2");
0030 #endif
0031   return chsq;
0032 }
0033
0034 void AdaptiveChisquarePrimaryVertexFitter::fill_trackinfo(const std::vector<reco::TransientTrack> &tracks,
0035                                                           const reco::BeamSpot &beamSpot) {
0036   /* fill track information used during fits into arrays, parallell to the list of input tracks */
0037
0038   trackinfo_.clear();
0039   trackinfo_.reserve(tracks.size());
0040
0041   for (auto &trk : tracks) {
0042     TrackInfo ti;
0043     // F1,F2 are the perigee parameters (3,4)
0044     const auto tspca = trk.stateAtBeamLine().trackStateAtPCA();  // freeTrajectoryState
0045     const auto tspca_pe = PerigeeConversions::ftsToPerigeeError(tspca);
0046     const auto momentum = tspca.momentum();
0047     auto const cos_phi = momentum.x() / momentum.perp();
0048     auto const sin_phi = momentum.y() / momentum.perp();
0049     auto const tan_lambda = momentum.z() / momentum.perp();
0050
0051     // covariance matrix of (F1,F2)
0052     double cov11 = tspca_pe.covarianceMatrix()(3, 3);
0053     double cov22 = tspca_pe.covarianceMatrix()(4, 4);
0054     double cov12 = tspca_pe.covarianceMatrix()(3, 4);
0055
0056     // S = cov^{-1}
0057     double DetV = cov11 * cov22 - cov12 * cov12;
0058     if (fabs(DetV) < 1.e-16) {
0059       edm::LogWarning("AdaptiveChisquarePrimaryVertexFitter")
0060           << "Warning, det(V) almost vanishes : " << DetV << " !! This should not happen!" << std::endl;
0061       ti.S11 = 0;
0062       ti.S22 = 0;
0063       ti.S12 = 0;
0064     } else {
0065       ti.S11 = cov22 / DetV;
0066       ti.S22 = cov11 / DetV;
0067       ti.S12 = -cov12 / DetV;
0068     }
0069     ti.b1 = tspca.position().x() * sin_phi - tspca.position().y() * cos_phi;
0070     ti.H1[0] = -sin_phi;
0071     ti.H1[1] = cos_phi;
0072     ti.H1[2] = 0;
0073     ti.b2 = tspca.position().z() - (tspca.position().x() * cos_phi + tspca.position().y() * sin_phi) * tan_lambda;
0074     ti.H2[0] = cos_phi * tan_lambda;
0075     ti.H2[1] = sin_phi * tan_lambda;
0076     ti.H2[2] = -1.;
0077
0078     for (int k = 0; k < 3; k++) {
0079       double SH1k = (ti.S11 * ti.H1[k] + ti.S12 * ti.H2[k]);
0080       double SH2k = (ti.S12 * ti.H1[k] + ti.S22 * ti.H2[k]);
0081       ti.g[k] = ti.b1 * SH1k + ti.b2 * SH2k;
0082       for (int l = 0; l < 3; l++) {
0083         ti.C(l, k) = ti.H1[l] * SH1k + ti.H2[l] * SH2k;
0084       }
0085     }
0086
0087     ti.zpca = tspca.position().z();
0088     ti.dzError = trk.track().dzError();
0089     trackinfo_.push_back(ti);
0090   }
0091 }
0092
0093 void AdaptiveChisquarePrimaryVertexFitter::make_vtx_trk_map(double zrange_scale) {
0094   unsigned const int nv = xv_.size();
0095   unsigned const int nt = trackinfo_.size();
0096
0097 #ifdef PVTX_DEBUG
0098   if (nv < 1) {
0099     edm::LogWarning("AdaptiveChisquarePrimaryVertexFit") << " empty vertex list with " << nt << " tracks" << std::endl;
0100     return;
0101   }
0102 #endif
0103
0104   // parallel lists for track to vertex mapping, tracks are not sorted
0105   tkmap_.clear();     // index in trackinfo_
0106   tkweight_.clear();  // weight in vertex
0107   tkfirstv_.clear();  // each vertex k owns a section of those list : tkfirstv_[k] .. tkfirstv_[k+1]-1
0108
0109   if (nv == 1) {
0110     // always accept all tracks for a single vertex fit
0111     tkfirstv_.push_back(0);
0112     tkfirstv_.push_back(nt);
0113     tkweight_.assign(nt, 0.);
0114     tkmap_.reserve(nt);
0115     for (unsigned int i = 0; i < nt; i++) {
0116       tkmap_.emplace_back(i);
0117     }
0118     return;
0119   }
0120
0121   // n > 1
0122   tkmap_.reserve(nv * 100);
0123   tkweight_.reserve(nv * 100);
0124   for (unsigned int k = 0; k < nv; k++) {
0125     tkfirstv_.emplace_back(tkmap_.size());
0126     for (unsigned int i = 0; i < nt; i++) {
0127       auto &ti = trackinfo_[i];
0128       const double zrange = zrange_scale * ti.dzError;
0129       if (std::abs(zv_[k] - ti.zpca) < z_cutoff_) {
0130         const double dztrk = ti.b2 + xv_[k] * ti.H2[0] + yv_[k] * ti.H2[1] - zv_[k];
0131         if (std::abs(dztrk) < zrange) {
0132           tkmap_.emplace_back(i);
0133           tkweight_.emplace_back(0.);
0134         }
0135       }
0136     }
0137   }
0138   tkfirstv_.emplace_back(tkmap_.size());  // extra entry, simplifies loops, every vertex has a "successor" now
0139 }
0140
0141 void AdaptiveChisquarePrimaryVertexFitter::fill_weights(const reco::BeamSpot &beamspot, double beta) {
0142   // multi-vertex version
0143   unsigned const int nt = trackinfo_.size();
0144   unsigned const int nv = xv_.size();
0145   const double beta_over_2 = 0.5 * beta;
0146   const double argmax = beta_over_2 * chi_cutoff_ * chi_cutoff_ * 5;
0147   const double Z_cutoff = vdt::fast_exp(-beta_over_2 * chi_cutoff_ * chi_cutoff_);
0148
0149   std::vector<double> Z_track(nt, Z_cutoff);
0150
0151   // evaluate and cache track-vertex assignment chi**2 for all clusters and sum up Z
0152   for (unsigned int k = 0; k < nv; k++) {
0153     for (unsigned int j = tkfirstv_[k]; j < tkfirstv_[k + 1]; j++) {
0154       const unsigned int i = tkmap_[j];
0155       double arg = beta_over_2 * track_in_vertex_chsq(trackinfo_[i], xv_[k], yv_[k], zv_[k]);
0156       if (arg < argmax) {
0157         const double e = vdt::fast_exp(-arg);
0158         tkweight_[j] = e;  // must later be normalized by the proper Z_track[i]
0159         Z_track[i] += e;   // sum up exponentials to normalize
0160       } else {
0161         tkweight_[j] = 0.;
0162       }
0163     }
0164   }
0165
0166   // now we have the partition function, Z_i and can evaluate assignment probabilities (aka weights)
0167   for (unsigned int j = 0; j < tkmap_.size(); j++) {
0168     const unsigned int i = tkmap_[j];
0169 #ifdef PVT_DEBUG
0170     assert((i < nt) && "tkmap out of range");
0171     assert((tkmap_.size() == tkweight_.size()) && "map and list not aliged");
0172 #endif
0173     tkweight_[j] /= Z_track[i];
0174   }
0175 }
0176
0177 bool AdaptiveChisquarePrimaryVertexFitter::clean() {
0178   /* in multi-vertex fitting, nearby vertices can fall on top of each other,
0179      even when the initial seeds don't, some kind of duplicate removal is required
0180      the approach in this method is similar to the method applied in clustering:
0181      at least two tracks with a weight above a threshold (trkweight_threshold) are required.
0182      vertices that don't fulfill this are either insignficant or very close
0183      to another vertex
0184    */
0185   const double trkweight_threshold = 0.7;
0186   unsigned int nv = xv_.size();
0187   if (nv < 2)
0188     return false;
0189
0190   // sum of weights per vertex
0191   std::vector<double> wsumhi(nv, 0);
0192   for (unsigned int k = 0; k < nv; k++) {
0193     for (unsigned int j = tkfirstv_[k]; j < tkfirstv_[k + 1]; j++) {
0194       if (tkweight_[j] > trkweight_threshold)
0195         wsumhi[k] += tkweight_[j];
0196     }
0197   }
0198
0199   double dzmin = 0;
0200   unsigned int k_dzmin = 0;
0201   for (unsigned int k = 0; k < nv - 1; k++) {
0202     double dz = std::abs(zv_[k + 1] - zv_[k]);
0203     if ((k == 0) || (dz < dzmin)) {
0204       dzmin = dz;
0205       k_dzmin = k;
0206     }
0207   }
0208
0209   if ((std::abs(dzmin) < 0.0200) && (std::min(wsumhi[k_dzmin], wsumhi[k_dzmin + 1]) < 0.5)) {
0210     if (wsumhi[k_dzmin] < wsumhi[k_dzmin + 1]) {
0211       remove_vertex(k_dzmin);
0212     } else {
0213       remove_vertex(k_dzmin + 1);
0214     }
0215   }
0216
0217   return true;
0218 }
0219
0220 void AdaptiveChisquarePrimaryVertexFitter::remove_vertex(unsigned int k) {
0221   // remove a vertex or rather merge it with it's neighbour
0222   // used for multi-vertex fits only
0223   unsigned int nv = xv_.size();
0224   if (nv < 2)
0225     return;
0226
0227   // 1) remove the vertex from the vertex list
0228   xv_.erase(xv_.begin() + k);
0229   yv_.erase(yv_.begin() + k);
0230   zv_.erase(zv_.begin() + k);
0231   covv_.erase(covv_.begin() + k);
0232
0233   // 2) adjust the track-map map
0234   // 2a) remove the map entries that belong the deleted vertex
0235   const unsigned int num_erased_map_entries = tkfirstv_[k + 1] - tkfirstv_[k];
0236   tkmap_.erase(tkmap_.begin() + tkfirstv_[k], tkmap_.begin() + tkfirstv_[k + 1]);
0237   tkweight_.erase(tkweight_.begin() + tkfirstv_[k], tkweight_.begin() + tkfirstv_[k + 1]);
0238   // 2b) adjust pointers for the following vertices, including the dummy entry behind the last (now [nv-1])
0239   for (unsigned int k1 = k + 1; k1 < nv + 1; k1++) {
0240     tkfirstv_[k1] -= num_erased_map_entries;
0241   }
0242   // 2c) erase the pointer of the removed vertex
0243   tkfirstv_.erase(tkfirstv_.begin() + k);
0244 }
0245
0246 double AdaptiveChisquarePrimaryVertexFitter::update(const reco::BeamSpot &beamspot,
0247                                                     const float beam_weight,
0248                                                     const bool fill_covariances) {
0249   double rho_vtx = 0;
0250   double delta_z = 0;
0251   double delta_x = 0;
0252   double delta_y = 0;
0253   unsigned const int nt = trackinfo_.size();
0254   unsigned const int nv = xv_.size();
0255   if (fill_covariances) {
0256     covv_.clear();
0257   }
0258
0259   // initial value for S, 0 or inverse of the beamspot covariance matrix
0260   Error3 S0;
0261   double c_beam[3] = {0, 0, 0};
0262   if (beam_weight > 0) {
0263     S0 = get_inverse_beam_covariance(beamspot);
0264     for (unsigned int j = 0; j < 3; j++) {
0265       c_beam[j] = -(S0(j, 0) * beamspot.x0() + S0(j, 1) * beamspot.y0() + S0(j, 2) * beamspot.z0());
0266     }
0267   }
0268
0269   for (unsigned int k = 0; k < nv; k++) {
0270     rho_vtx = 0;
0271     Error3 S(S0);
0272     // sum track contributions
0273     double c[3] = {c_beam[0], c_beam[1], c_beam[2]};
0274     for (unsigned int j = tkfirstv_[k]; j < tkfirstv_[k + 1]; j++) {
0275       const unsigned int i = tkmap_[j];
0276       const auto w = tkweight_[j];
0277       rho_vtx += w;
0278       S += w * trackinfo_[i].C;
0279       for (unsigned int l = 0; l < 3; l++) {
0280         c[l] += w * trackinfo_[i].g[l];
0281       }
0282     }
0283
0284 #ifdef PVTX_DEBUG
0285     if ((fabs(S(1, 2) - S(2, 1)) > 1e-3) || (fabs(S(0, 2) - S(2, 0)) > 1e-3) || (fabs(S(0, 1) - S(1, 0)) > 1e-3) ||
0286         (S(0, 0) <= 0) || (S(0, 0) <= 0) || (S(0, 0) <= 0)) {
0287       edm::LogWarning("AdaptiveChisquarePrimaryVertexFitter") << "update()  bad S-matrix   S=" << std::endl
0288                                                               << S << std::endl;
0289       edm::LogWarning("AdaptiveChisquarePrimaryVertexFitter")
0290           << "n-vertex = " << nv << "  n-track = " << nt << std::endl;
0291     }
0292 #endif
0293
0294     const auto xold = xv_[k];
0295     const auto yold = yv_[k];
0296     const auto zold = zv_[k];
0297
0298     if (S.Invert()) {
0299       xv_[k] = -(S(0, 0) * c[0] + S(0, 1) * c[1] + S(0, 2) * c[2]);
0300       yv_[k] = -(S(1, 0) * c[0] + S(1, 1) * c[1] + S(1, 2) * c[2]);
0301       zv_[k] = -(S(2, 0) * c[0] + S(2, 1) * c[1] + S(2, 2) * c[2]);
0302       if (fill_covariances) {
0303         covv_.emplace_back(S);
0304       }
0305     } else {
0306       edm::LogWarning("AdaptiveChisquarePrimaryVertexFitter") << "update()   Matrix inversion failed" << S << std::endl;
0307       if (fill_covariances) {
0308         Error3 covv_dummy;
0309         covv_dummy(0, 0) = 100.;
0310         covv_dummy(1, 1) = 100.;
0311         covv_dummy(2, 2) = 100.;
0312         covv_.emplace_back(covv_dummy);
0313       }
0314     }
0315
0316     if ((nt > 1) && (rho_vtx > 1.0)) {
0317       delta_x = std::max(delta_x, std::abs(xv_[k] - xold));
0318       delta_y = std::max(delta_y, std::abs(yv_[k] - yold));
0319       delta_z = std::max(delta_z, std::abs(zv_[k] - zold));
0320     }
0321
0322   }  // vertex loop
0323
0324   return std::max(delta_z, std::max(delta_x, delta_y));
0325 }
0326
0327 AdaptiveChisquarePrimaryVertexFitter::Error3 AdaptiveChisquarePrimaryVertexFitter::get_inverse_beam_covariance(
0328     const reco::BeamSpot &beamspot) {
0329   auto SBeam = beamspot.rotatedCovariance3D();
0330   if (SBeam.Invert()) {
0331     return SBeam;
0332   } else {
0333     edm::LogWarning("AdaptiveChisquarePrimaryVertexFitter")
0334         << "Warning, beam-spot covariance matrix inversion failed " << std::endl;
0335     Error3 S0;
0336     S0(0, 0) = 1. / pow(beamspot.BeamWidthX(), 2);
0337     S0(1, 1) = 1. / pow(beamspot.BeamWidthY(), 2);
0338     S0(2, 2) = 1. / pow(beamspot.sigmaZ(), 2);
0339     return S0;
0340   }
0341 }
0342
0343 TransientVertex AdaptiveChisquarePrimaryVertexFitter::get_TransientVertex(
0344     const unsigned int k,
0345     const std::vector<std::pair<unsigned int, float>> &vertex_track_weights,
0346     const std::vector<reco::TransientTrack> &tracks,
0347     const float beam_weight,
0348     const reco::BeamSpot &beamspot) {
0349   const GlobalPoint pos(xv_[k], yv_[k], zv_[k]);
0350   const GlobalError posError(
0351       covv_[k](0, 0), covv_[k](1, 0), covv_[k](1, 1), covv_[k](2, 0), covv_[k](2, 1), covv_[k](2, 2));
0352   float chi2 = 0.;
0353   float vtx_ndof = -3.;
0354   if (beam_weight > 0) {
0355     // add beam-spot chi**2 and degrees of freedom
0356     vtx_ndof = 3 * beam_weight;
0357     const auto S = get_inverse_beam_covariance(beamspot);
0358     const double dx = xv_[k] - beamspot.x0();
0359     const double dy = yv_[k] - beamspot.y0();
0360     const double dz = zv_[k] - beamspot.z0();
0361     chi2 = beam_weight * (S(0, 0) * dx * dx + S(1, 1) * dy * dy + 2 * S(0, 1) * dx * dy + S(2, 2) * dz * dz +
0362                           2 * S(0, 2) * dx * dz + 2 * S(1, 2) * dy * dz);
0363   }
0364
0365   std::vector<reco::TransientTrack> vertex_tracks;
0366   TransientVertex::TransientTrackToFloatMap trkWeightMap;
0367   for (const auto &tk : vertex_track_weights) {
0368     const unsigned int i = tk.first;
0369     const float track_weight = tk.second;
0370     if (track_weight >= min_trackweight_) {
0371       vertex_tracks.emplace_back(tracks[i]);
0372       trkWeightMap[tracks[i]] = track_weight;
0373       vtx_ndof += 2 * track_weight;
0374       chi2 += track_weight * track_in_vertex_chsq(trackinfo_[i], xv_[k], yv_[k], zv_[k]);
0375     }
0376   }
0377
0378   auto vtx = TransientVertex(pos, posError, vertex_tracks, chi2, vtx_ndof);
0379   vtx.weightMap(trkWeightMap);
0380
0381   return vtx;
0382 }
0383
0384 std::vector<TransientVertex> AdaptiveChisquarePrimaryVertexFitter::vertices(
0385     const std::vector<reco::TransientTrack> &tracks,
0386     const std::vector<TransientVertex> &clusters,
0387     const reco::BeamSpot &beamspot,
0388     const bool useBeamConstraint) {
0389   // simultaneous fit of all vertices in the input list
0390
0391   const int max_iterations = 50;
0392
0393   // initialize the vertices
0394   const unsigned int nv = clusters.size();
0395   xv_.clear();
0396   xv_.reserve(nv);
0397   yv_.clear();
0398   yv_.reserve(nv);
0399   zv_.clear();
0400   zv_.reserve(nv);
0401   tkfirstv_.clear();
0402   tkfirstv_.reserve(nv + 1);
0403   covv_.clear();
0404   covv_.reserve(nv);
0405
0406   // seeds
0407   for (auto &clu : clusters) {
0408     const double zclu = clu.position().z();
0409     xv_.emplace_back(beamspot.x(zclu));
0410     yv_.emplace_back(beamspot.y(zclu));
0411     zv_.emplace_back(zclu);
0412   }
0413
0414   fill_trackinfo(tracks, beamspot);
0415
0416   make_vtx_trk_map(5.);  // use tracks within 5 sigma windows (if that is less than z_cutoff_)
0417
0418   float beam_weight = useBeamConstraint ? 1. : 0.;
0419
0420   double delta = 0;
0421   unsigned int nit = 0;
0422   while ((nit == 0) || ((delta > 0.0001) && (nit < max_iterations))) {
0423     fill_weights(beamspot);
0424     delta = update(beamspot, beam_weight, false);
0425     nit++;
0426   }
0427   if ((nit >= max_iterations) && (delta > 0.01)) {
0428     edm::LogWarning("AdaptiveChisquarePrimaryVertexFitter")
0429         << "iteration limit reached " << nit << "  last delta = " << delta << std::endl
0430         << " nv = " << nv << "    nt = " << tracks.size() << std::endl;
0431   }
0432
0433   // may need to remove collapsed vertices
0434   nit = 0;
0435   while ((xv_.size() > 1) && (nit < max_iterations) && (clean())) {
0436     fill_weights(beamspot);
0437     update(beamspot, beam_weight, false);
0438     nit++;
0439   }
0440
0441   // fill the covariance matrices
0442   update(beamspot, beam_weight, true);
0443
0444   // assign tracks to vertices
0445   std::vector<unsigned int> track_to_vertex(trackinfo_.size(), nv);
0446   // for each track identify the vertex that wants it most
0447   std::vector<double> maxweight(trackinfo_.size(), -1.);
0448   for (unsigned int k = 0; k < nv; k++) {
0449     for (unsigned int j = tkfirstv_[k]; j < tkfirstv_[k + 1]; j++) {
0450       const unsigned int i = tkmap_[j];
0451       if (tkweight_[j] > maxweight[i]) {
0452         maxweight[i] = tkweight_[j];
0453         track_to_vertex[i] = k;
0454       }
0455     }
0456   }
0457
0458   // fill the fit result into transient vertices
0459   std::vector<TransientVertex> pvs;
0460   for (unsigned int k = 0; k < xv_.size(); k++) {
0461     std::vector<std::pair<unsigned int, float>> vertex_tracks_weights;
0462     for (unsigned int j = tkfirstv_[k]; j < tkfirstv_[k + 1]; j++) {
0463       unsigned int i = tkmap_[j];
0464       if (track_to_vertex[i] == k) {
0465         vertex_tracks_weights.emplace_back(tkmap_[j], tkweight_[j]);
0466       }
0467     }
0468     pvs.emplace_back(get_TransientVertex(k, vertex_tracks_weights, tracks, beam_weight, beamspot));
0469   }
0470
0471   return pvs;
0472 }
0473
0474 TransientVertex AdaptiveChisquarePrimaryVertexFitter::refit(const TransientVertex &cluster,
0475                                                             const reco::BeamSpot &beamspot,
0476                                                             const bool useBeamConstraint) {
0477   // fit a single vertex using all tracks in the tracklist
0478   const unsigned int nt = cluster.originalTracks().size();
0479   const int max_iterations = 50;
0480
0481   // initialize, vectors with size=1 here to avoid code duplication from the multivertex case in update()
0482   const double zclu = cluster.position().z();
0483   xv_ = {beamspot.x(zclu)};
0484   yv_ = {beamspot.y(zclu)};
0485   zv_ = {zclu};
0486   tkfirstv_ = {0, nt};
0487   covv_.clear();
0488
0489   fill_trackinfo(cluster.originalTracks(), beamspot);
0490   tkweight_.assign(nt, 0.);
0491   tkmap_.clear();
0492   tkmap_.reserve(nt);
0493   for (unsigned int i = 0; i < nt; i++) {
0494     tkmap_.emplace_back(i);  // trivial map for single vertex fits
0495   }
0496
0497   float beam_weight = useBeamConstraint ? 1. : 0.;
0498
0499   double delta = 0;
0500   unsigned int nit = 0;
0501   while ((nit == 0) || ((delta > 0.0001) && (nit < max_iterations))) {
0502     fill_weights(beamspot);
0503     delta = update(beamspot, beam_weight, false);
0504     nit++;
0505   }
0506
0507   if ((nit >= max_iterations) && (delta > 0.1)) {
0508     edm::LogWarning("AdaptiveChisquarePrimaryVertexFitter")
0509         << "single vertex fit, iteration limit reached " << nit << "  last delta = " << delta << std::endl
0510         << "    nt = " << cluster.originalTracks().size() << std::endl;
0511   }
0512
0513   // fill the covariance matrices
0514   update(beamspot, beam_weight, true);
0515
0516   // put the result into a transient vertex
0517   std::vector<std::pair<unsigned int, float>> vertex_track_weights;
0518   for (unsigned int i = 0; i < nt; i++) {
0519     vertex_track_weights.emplace_back(i, tkweight_[i]);
0520   }
0521
0522   return get_TransientVertex(0, vertex_track_weights, cluster.originalTracks(), beam_weight, beamspot);
0523 }
0524
0525 //
0526 std::vector<TransientVertex> AdaptiveChisquarePrimaryVertexFitter::fit(const std::vector<reco::TransientTrack> &tracks,
0527                                                                        const std::vector<TransientVertex> &clusters,
0528                                                                        const reco::BeamSpot &beamspot,
0529                                                                        const bool useBeamConstraint) {
0530   if (multivertexfit_) {
0531     return vertices(tracks, clusters, beamspot, useBeamConstraint);
0532
0533   } else {
0534     // fit the clusters one-by-one using the tracklist of the clusters (ignores the "tracks" argument)
0535     std::vector<TransientVertex> pvs;
0536     pvs.reserve(clusters.size());
0537     for (auto &cluster : clusters) {
0538       if (cluster.originalTracks().size() > (useBeamConstraint ? 0 : 1)) {
0539         auto pv = refit(cluster, beamspot, useBeamConstraint);
0540         if (pv.isValid()) {
0541           pvs.emplace_back(pv);
0542         }
0543       }
0544     }
0545     return pvs;
0546   }
0547 }