Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-04-01 23:54:11

0001 #include "L1Trigger/VertexFinder/interface/VertexFinder.h"
0002 
0003 using namespace std;
0004 
0005 namespace l1tVertexFinder {
0006 
0007   void VertexFinder::computeAndSetVertexParameters(RecoVertex<>& vertex,
0008                                                    const std::vector<float>& bin_centers,
0009                                                    const std::vector<unsigned int>& counts) {
0010     double pt = 0.;
0011     double z0 = -999.;
0012     double z0width = 0.;
0013     bool highPt = false;
0014     double highestPt = 0.;
0015     unsigned int numHighPtTracks = 0;
0016 
0017     float SumZ = 0.;
0018     float z0square = 0.;
0019     float trackPt = 0.;
0020 
0021     std::vector<double> bin_pt(bin_centers.size(), 0.0);
0022     unsigned int ibin = 0;
0023     unsigned int itrack = 0;
0024 
0025     for (const L1Track* track : vertex.tracks()) {
0026       itrack++;
0027       trackPt = track->pt();
0028 
0029       // Skip the bins with no tracks
0030       while (ibin < counts.size() && counts[ibin] == 0)
0031         ibin++;
0032 
0033       if (trackPt > settings_->vx_TrackMaxPt()) {
0034         highPt = true;
0035         numHighPtTracks++;
0036         highestPt = (trackPt > highestPt) ? trackPt : highestPt;
0037         if (settings_->vx_TrackMaxPtBehavior() == 0)
0038           continue;  // ignore this track
0039         else if (settings_->vx_TrackMaxPtBehavior() == 1)
0040           trackPt = settings_->vx_TrackMaxPt();  // saturate
0041       }
0042 
0043       pt += std::pow(trackPt, settings_->vx_weightedmean());
0044       if (bin_centers.empty() && counts.empty()) {
0045         SumZ += track->z0() * std::pow(trackPt, settings_->vx_weightedmean());
0046         z0square += track->z0() * track->z0();
0047       } else {
0048         bin_pt[ibin] += std::pow(trackPt, settings_->vx_weightedmean());
0049         if (itrack == counts[ibin]) {
0050           SumZ += bin_centers[ibin] * bin_pt[ibin];
0051           z0square += bin_centers[ibin] * bin_centers[ibin];
0052           itrack = 0;
0053           ibin++;
0054         }
0055       }
0056     }
0057 
0058     z0 = SumZ / ((settings_->vx_weightedmean() > 0) ? pt : vertex.numTracks());
0059     z0square /= vertex.numTracks();
0060     z0width = sqrt(std::abs(z0 * z0 - z0square));
0061 
0062     vertex.setParameters(pt, z0, z0width, highPt, numHighPtTracks, highestPt);
0063   }
0064 
0065   void VertexFinder::GapClustering() {
0066     sort(fitTracks_.begin(), fitTracks_.end(), SortTracksByZ0());
0067     iterations_ = 0;
0068     RecoVertex Vertex;
0069     for (unsigned int i = 0; i < fitTracks_.size(); ++i) {
0070       Vertex.insert(&fitTracks_[i]);
0071       iterations_++;
0072       if ((i + 1 < fitTracks_.size() and fitTracks_[i + 1].z0() - fitTracks_[i].z0() > settings_->vx_distance()) or
0073           i == fitTracks_.size() - 1) {
0074         if (Vertex.numTracks() >= settings_->vx_minTracks()) {
0075           computeAndSetVertexParameters(Vertex, {}, {});
0076           vertices_.push_back(Vertex);
0077         }
0078         Vertex.clear();
0079       }
0080     }
0081   }
0082 
0083   float VertexFinder::maxDistance(RecoVertex<> cluster0, RecoVertex<> cluster1) {
0084     float distance = 0;
0085     for (const L1Track* track0 : cluster0.tracks()) {
0086       for (const L1Track* track1 : cluster1.tracks()) {
0087         if (std::abs(track0->z0() - track1->z0()) > distance) {
0088           distance = std::abs(track0->z0() - track1->z0());
0089         }
0090       }
0091     }
0092 
0093     return distance;
0094   }
0095 
0096   float VertexFinder::minDistance(RecoVertex<> cluster0, RecoVertex<> cluster1) {
0097     float distance = 9999;
0098     for (const L1Track* track0 : cluster0.tracks()) {
0099       for (const L1Track* track1 : cluster1.tracks()) {
0100         if (std::abs(track0->z0() - track1->z0()) < distance) {
0101           distance = std::abs(track0->z0() - track1->z0());
0102         }
0103       }
0104     }
0105 
0106     return distance;
0107   }
0108 
0109   float VertexFinder::meanDistance(RecoVertex<> cluster0, RecoVertex<> cluster1) {
0110     float distanceSum = 0;
0111 
0112     for (const L1Track* track0 : cluster0.tracks()) {
0113       for (const L1Track* track1 : cluster1.tracks()) {
0114         distanceSum += std::abs(track0->z0() - track1->z0());
0115       }
0116     }
0117 
0118     float distance = distanceSum / (cluster0.numTracks() * cluster1.numTracks());
0119     return distance;
0120   }
0121 
0122   float VertexFinder::centralDistance(RecoVertex<> cluster0, RecoVertex<> cluster1) {
0123     computeAndSetVertexParameters(cluster0, {}, {});
0124     computeAndSetVertexParameters(cluster1, {}, {});
0125     float distance = std::abs(cluster0.z0() - cluster1.z0());
0126     return distance;
0127   }
0128 
0129   void VertexFinder::agglomerativeHierarchicalClustering() {
0130     iterations_ = 0;
0131 
0132     sort(fitTracks_.begin(), fitTracks_.end(), SortTracksByZ0());
0133 
0134     RecoVertexCollection vClusters;
0135     vClusters.resize(fitTracks_.size());
0136 
0137     for (unsigned int i = 0; i < fitTracks_.size(); ++i) {
0138       vClusters[i].insert(&fitTracks_[i]);
0139       // iterations_++;
0140     }
0141 
0142     while (true) {
0143       float MinimumScore = 9999;
0144 
0145       unsigned int clusterId0 = 0;
0146       unsigned int clusterId1 = 0;
0147       for (unsigned int iClust = 0; iClust < vClusters.size() - 1; iClust++) {
0148         iterations_++;
0149 
0150         float M = 0;
0151         if (settings_->vx_distanceType() == 0)
0152           M = maxDistance(vClusters[iClust], vClusters[iClust + 1]);
0153         else if (settings_->vx_distanceType() == 1)
0154           M = minDistance(vClusters[iClust], vClusters[iClust + 1]);
0155         else if (settings_->vx_distanceType() == 2)
0156           M = meanDistance(vClusters[iClust], vClusters[iClust + 1]);
0157         else
0158           M = centralDistance(vClusters[iClust], vClusters[iClust + 1]);
0159 
0160         if (M < MinimumScore) {
0161           MinimumScore = M;
0162           clusterId0 = iClust;
0163           clusterId1 = iClust + 1;
0164         }
0165       }
0166       if (MinimumScore > settings_->vx_distance() or vClusters[clusterId1].tracks().empty())
0167         break;
0168       for (const L1Track* track : vClusters[clusterId0].tracks()) {
0169         vClusters[clusterId1].insert(track);
0170       }
0171       vClusters.erase(vClusters.begin() + clusterId0);
0172     }
0173 
0174     for (RecoVertex clust : vClusters) {
0175       if (clust.numTracks() >= settings_->vx_minTracks()) {
0176         computeAndSetVertexParameters(clust, {}, {});
0177         vertices_.push_back(clust);
0178       }
0179     }
0180   }
0181 
0182   void VertexFinder::DBSCAN() {
0183     // std::vector<RecoVertex> vClusters;
0184     std::vector<unsigned int> visited;
0185     std::vector<unsigned int> saved;
0186 
0187     sort(fitTracks_.begin(), fitTracks_.end(), SortTracksByPt());
0188     iterations_ = 0;
0189 
0190     for (unsigned int i = 0; i < fitTracks_.size(); ++i) {
0191       if (find(visited.begin(), visited.end(), i) != visited.end())
0192         continue;
0193 
0194       // if(fitTracks_[i]->pt()>10.){
0195       visited.push_back(i);
0196       std::set<unsigned int> neighbourTrackIds;
0197       unsigned int numDensityTracks = 0;
0198       if (fitTracks_[i].pt() > settings_->vx_dbscan_pt())
0199         numDensityTracks++;
0200       else
0201         continue;
0202       for (unsigned int k = 0; k < fitTracks_.size(); ++k) {
0203         iterations_++;
0204         if (k != i and std::abs(fitTracks_[k].z0() - fitTracks_[i].z0()) < settings_->vx_distance()) {
0205           neighbourTrackIds.insert(k);
0206           if (fitTracks_[k].pt() > settings_->vx_dbscan_pt()) {
0207             numDensityTracks++;
0208           }
0209         }
0210       }
0211 
0212       if (numDensityTracks < settings_->vx_dbscan_mintracks()) {
0213         // mark track as noise
0214       } else {
0215         RecoVertex vertex;
0216         vertex.insert(&fitTracks_[i]);
0217         saved.push_back(i);
0218         for (unsigned int id : neighbourTrackIds) {
0219           if (find(visited.begin(), visited.end(), id) == visited.end()) {
0220             visited.push_back(id);
0221             std::vector<unsigned int> neighbourTrackIds2;
0222             for (unsigned int k = 0; k < fitTracks_.size(); ++k) {
0223               iterations_++;
0224               if (std::abs(fitTracks_[k].z0() - fitTracks_[id].z0()) < settings_->vx_distance()) {
0225                 neighbourTrackIds2.push_back(k);
0226               }
0227             }
0228 
0229             // if (neighbourTrackIds2.size() >= settings_->vx_minTracks()) {
0230             for (unsigned int id2 : neighbourTrackIds2) {
0231               neighbourTrackIds.insert(id2);
0232             }
0233             // }
0234           }
0235           if (find(saved.begin(), saved.end(), id) == saved.end())
0236             vertex.insert(&fitTracks_[id]);
0237         }
0238         computeAndSetVertexParameters(vertex, {}, {});
0239         if (vertex.numTracks() >= settings_->vx_minTracks())
0240           vertices_.push_back(vertex);
0241       }
0242       // }
0243     }
0244   }
0245 
0246   void VertexFinder::PVR() {
0247     bool start = true;
0248     FitTrackCollection discardedTracks, acceptedTracks;
0249     iterations_ = 0;
0250     for (const L1Track& track : fitTracks_) {
0251       acceptedTracks.push_back(track);
0252     }
0253 
0254     while (discardedTracks.size() >= settings_->vx_minTracks() or start == true) {
0255       start = false;
0256       bool removing = true;
0257       discardedTracks.clear();
0258       while (removing) {
0259         float oldDistance = 0.;
0260 
0261         if (settings_->debug() > 2)
0262           edm::LogInfo("VertexFinder") << "PVR::AcceptedTracks " << acceptedTracks.size();
0263 
0264         float z0start = 0;
0265         for (const L1Track& track : acceptedTracks) {
0266           z0start += track.z0();
0267           iterations_++;
0268         }
0269 
0270         z0start /= acceptedTracks.size();
0271         if (settings_->debug() > 2)
0272           edm::LogInfo("VertexFinder") << "PVR::z0 vertex " << z0start;
0273         FitTrackCollection::iterator badTrackIt = acceptedTracks.end();
0274         removing = false;
0275 
0276         for (FitTrackCollection::iterator it = acceptedTracks.begin(); it < acceptedTracks.end(); ++it) {
0277           const L1Track* track = &*it;
0278           iterations_++;
0279           if (std::abs(track->z0() - z0start) > settings_->vx_distance() and
0280               std::abs(track->z0() - z0start) > oldDistance) {
0281             badTrackIt = it;
0282             oldDistance = std::abs(track->z0() - z0start);
0283             removing = true;
0284           }
0285         }
0286 
0287         if (removing) {
0288           const L1Track badTrack = *badTrackIt;
0289           if (settings_->debug() > 2)
0290             edm::LogInfo("VertexFinder") << "PVR::Removing track " << badTrack.z0() << " at distance " << oldDistance;
0291           discardedTracks.push_back(badTrack);
0292           acceptedTracks.erase(badTrackIt);
0293         }
0294       }
0295 
0296       if (acceptedTracks.size() >= settings_->vx_minTracks()) {
0297         RecoVertex vertex;
0298         for (const L1Track& track : acceptedTracks) {
0299           vertex.insert(&track);
0300         }
0301         computeAndSetVertexParameters(vertex, {}, {});
0302         vertices_.push_back(vertex);
0303       }
0304       if (settings_->debug() > 2)
0305         edm::LogInfo("VertexFinder") << "PVR::DiscardedTracks size " << discardedTracks.size();
0306       acceptedTracks.clear();
0307       acceptedTracks = discardedTracks;
0308     }
0309   }
0310 
0311   void VertexFinder::adaptiveVertexReconstruction() {
0312     bool start = true;
0313     iterations_ = 0;
0314     FitTrackCollection discardedTracks, acceptedTracks, discardedTracks2;
0315 
0316     for (const L1Track& track : fitTracks_) {
0317       discardedTracks.push_back(track);
0318     }
0319 
0320     while (discardedTracks.size() >= settings_->vx_minTracks() or start == true) {
0321       start = false;
0322       discardedTracks2.clear();
0323       FitTrackCollection::iterator it = discardedTracks.begin();
0324       const L1Track track = *it;
0325       acceptedTracks.push_back(track);
0326       float z0sum = track.z0();
0327 
0328       for (FitTrackCollection::iterator it2 = discardedTracks.begin(); it2 < discardedTracks.end(); ++it2) {
0329         if (it2 != it) {
0330           const L1Track secondTrack = *it2;
0331           // Calculate new vertex z0 adding this track
0332           z0sum += secondTrack.z0();
0333           float z0vertex = z0sum / (acceptedTracks.size() + 1);
0334           // Calculate chi2 of new vertex
0335           float chi2 = 0.;
0336           float dof = 0.;
0337           for (const L1Track& accTrack : acceptedTracks) {
0338             iterations_++;
0339             float Residual = accTrack.z0() - z0vertex;
0340             if (std::abs(accTrack.eta()) < 1.2)
0341               Residual /= 0.1812;  // Assumed z0 resolution
0342             else if (std::abs(accTrack.eta()) >= 1.2 && std::abs(accTrack.eta()) < 1.6)
0343               Residual /= 0.2912;
0344             else if (std::abs(accTrack.eta()) >= 1.6 && std::abs(accTrack.eta()) < 2.)
0345               Residual /= 0.4628;
0346             else
0347               Residual /= 0.65;
0348 
0349             chi2 += Residual * Residual;
0350             dof = (acceptedTracks.size() + 1) * 2 - 1;
0351           }
0352           if (chi2 / dof < settings_->vx_chi2cut()) {
0353             acceptedTracks.push_back(secondTrack);
0354           } else {
0355             discardedTracks2.push_back(secondTrack);
0356             z0sum -= secondTrack.z0();
0357           }
0358         }
0359       }
0360 
0361       if (acceptedTracks.size() >= settings_->vx_minTracks()) {
0362         RecoVertex vertex;
0363         for (const L1Track& track : acceptedTracks) {
0364           vertex.insert(&track);
0365         }
0366         computeAndSetVertexParameters(vertex, {}, {});
0367         vertices_.push_back(vertex);
0368       }
0369 
0370       acceptedTracks.clear();
0371       discardedTracks.clear();
0372       discardedTracks = discardedTracks2;
0373     }
0374   }
0375 
0376   void VertexFinder::HPV() {
0377     iterations_ = 0;
0378     sort(fitTracks_.begin(), fitTracks_.end(), SortTracksByPt());
0379 
0380     RecoVertex vertex;
0381     bool first = true;
0382     float z = 99.;
0383     for (const L1Track& track : fitTracks_) {
0384       if (track.pt() < 50.) {
0385         if (first) {
0386           first = false;
0387           z = track.z0();
0388           vertex.insert(&track);
0389         } else {
0390           if (std::abs(track.z0() - z) < settings_->vx_distance())
0391             vertex.insert(&track);
0392         }
0393       }
0394     }
0395 
0396     computeAndSetVertexParameters(vertex, {}, {});
0397     vertex.setZ0(z);
0398     vertices_.push_back(vertex);
0399   }
0400 
0401   void VertexFinder::Kmeans() {
0402     unsigned int NumberOfClusters = settings_->vx_kmeans_nclusters();
0403 
0404     vertices_.resize(NumberOfClusters);
0405     float ClusterSeparation = 30. / NumberOfClusters;
0406 
0407     for (unsigned int i = 0; i < NumberOfClusters; ++i) {
0408       float ClusterCentre = -15. + ClusterSeparation * (i + 0.5);
0409       vertices_[i].setZ0(ClusterCentre);
0410     }
0411     unsigned int iterations = 0;
0412     // Initialise Clusters
0413     while (iterations < settings_->vx_kmeans_iterations()) {
0414       for (unsigned int i = 0; i < NumberOfClusters; ++i) {
0415         vertices_[i].clear();
0416       }
0417 
0418       for (const L1Track& track : fitTracks_) {
0419         float distance = 9999;
0420         if (iterations == settings_->vx_kmeans_iterations() - 3)
0421           distance = settings_->vx_distance() * 2;
0422         if (iterations > settings_->vx_kmeans_iterations() - 3)
0423           distance = settings_->vx_distance();
0424         unsigned int ClusterId;
0425         bool NA = true;
0426         for (unsigned int id = 0; id < NumberOfClusters; ++id) {
0427           if (std::abs(track.z0() - vertices_[id].z0()) < distance) {
0428             distance = std::abs(track.z0() - vertices_[id].z0());
0429             ClusterId = id;
0430             NA = false;
0431           }
0432         }
0433         if (!NA) {
0434           vertices_[ClusterId].insert(&track);
0435         }
0436       }
0437       for (unsigned int i = 0; i < NumberOfClusters; ++i) {
0438         if (vertices_[i].numTracks() >= settings_->vx_minTracks())
0439           computeAndSetVertexParameters(vertices_[i], {}, {});
0440       }
0441       iterations++;
0442     }
0443   }
0444 
0445   void VertexFinder::findPrimaryVertex() {
0446     if (settings_->vx_precision() == Precision::Emulation) {
0447       pv_index_ = std::distance(verticesEmulation_.begin(),
0448                                 std::max_element(verticesEmulation_.begin(),
0449                                                  verticesEmulation_.end(),
0450                                                  [](const l1t::VertexWord& vertex0, const l1t::VertexWord& vertex1) {
0451                                                    return (vertex0.pt() < vertex1.pt());
0452                                                  }));
0453     } else {
0454       pv_index_ = std::distance(
0455           vertices_.begin(),
0456           std::max_element(
0457               vertices_.begin(), vertices_.end(), [](const RecoVertex<>& vertex0, const RecoVertex<>& vertex1) {
0458                 return (vertex0.pt() < vertex1.pt());
0459               }));
0460     }
0461   }
0462 
0463   // Possible Formatting Codes: https://misc.flogisoft.com/bash/tip_colors_and_formatting
0464   template <class data_type, typename stream_type>
0465   void VertexFinder::printHistogram(stream_type& stream,
0466                                     std::vector<data_type> data,
0467                                     int width,
0468                                     int minimum,
0469                                     int maximum,
0470                                     std::string title,
0471                                     std::string color) {
0472     int tableSize = data.size();
0473 
0474     if (maximum == -1) {
0475       maximum = float(*std::max_element(std::begin(data), std::end(data))) * 1.05;
0476     } else if (maximum <= minimum) {
0477       maximum = float(*std::max_element(std::begin(data), std::end(data))) * 1.05;
0478       minimum = float(*std::min_element(std::begin(data), std::end(data)));
0479     }
0480 
0481     if (minimum < 0) {
0482       minimum *= 1.05;
0483     } else {
0484       minimum = 0;
0485     }
0486 
0487     std::vector<std::string> intervals(tableSize, "");
0488     std::vector<std::string> values(tableSize, "");
0489     char buffer[128];
0490     int intervalswidth = 0, valueswidth = 0, tmpwidth = 0;
0491     for (int i = 0; i < tableSize; i++) {
0492       //Format the bin labels
0493       tmpwidth = sprintf(buffer, "[%-.5g, %-.5g)", float(i), float(i + 1));
0494       intervals[i] = buffer;
0495       if (i == (tableSize - 1)) {
0496         intervals[i][intervals[i].size() - 1] = ']';
0497       }
0498       if (tmpwidth > intervalswidth)
0499         intervalswidth = tmpwidth;
0500 
0501       //Format the values
0502       tmpwidth = sprintf(buffer, "%-.5g", float(data[i]));
0503       values[i] = buffer;
0504       if (tmpwidth > valueswidth)
0505         valueswidth = tmpwidth;
0506     }
0507 
0508     sprintf(buffer, "%-.5g", float(minimum));
0509     std::string minimumtext = buffer;
0510     sprintf(buffer, "%-.5g", float(maximum));
0511     std::string maximumtext = buffer;
0512 
0513     int plotwidth =
0514         std::max(int(minimumtext.size() + maximumtext.size()), width - (intervalswidth + 1 + valueswidth + 1 + 2));
0515     std::string scale =
0516         minimumtext + std::string(plotwidth + 2 - minimumtext.size() - maximumtext.size(), ' ') + maximumtext;
0517 
0518     float norm = float(plotwidth) / float(maximum - minimum);
0519     int zero = std::round((0.0 - minimum) * norm);
0520     std::vector<char> line(plotwidth, '-');
0521 
0522     if ((minimum != 0) && (0 <= zero) && (zero < plotwidth)) {
0523       line[zero] = '+';
0524     }
0525     std::string capstone =
0526         std::string(intervalswidth + 1 + valueswidth + 1, ' ') + "+" + std::string(line.begin(), line.end()) + "+";
0527 
0528     std::vector<std::string> out;
0529     if (!title.empty()) {
0530       out.push_back(title);
0531       out.push_back(std::string(title.size(), '='));
0532     }
0533     out.push_back(std::string(intervalswidth + valueswidth + 2, ' ') + scale);
0534     out.push_back(capstone);
0535     for (int i = 0; i < tableSize; i++) {
0536       std::string interval = intervals[i];
0537       std::string value = values[i];
0538       data_type x = data[i];
0539       std::fill_n(line.begin(), plotwidth, ' ');
0540 
0541       int pos = std::round((float(x) - minimum) * norm);
0542       if (x < 0) {
0543         std::fill_n(line.begin() + pos, zero - pos, '*');
0544       } else {
0545         std::fill_n(line.begin() + zero, pos - zero, '*');
0546       }
0547 
0548       if ((minimum != 0) && (0 <= zero) && (zero < plotwidth)) {
0549         line[zero] = '|';
0550       }
0551 
0552       sprintf(buffer,
0553               "%-*s %-*s |%s|",
0554               intervalswidth,
0555               interval.c_str(),
0556               valueswidth,
0557               value.c_str(),
0558               std::string(line.begin(), line.end()).c_str());
0559       out.push_back(buffer);
0560     }
0561     out.push_back(capstone);
0562     if (!color.empty())
0563       stream << color;
0564     for (const auto& o : out) {
0565       stream << o << "\n";
0566     }
0567     if (!color.empty())
0568       stream << "\e[0m";
0569     stream << "\n";
0570   }
0571 
0572   void VertexFinder::sortVerticesInPt() {
0573     if (settings_->vx_precision() == Precision::Emulation) {
0574       std::sort(
0575           verticesEmulation_.begin(),
0576           verticesEmulation_.end(),
0577           [](const l1t::VertexWord& vertex0, const l1t::VertexWord& vertex1) { return (vertex0.pt() > vertex1.pt()); });
0578     } else {
0579       std::sort(vertices_.begin(), vertices_.end(), [](const RecoVertex<>& vertex0, const RecoVertex<>& vertex1) {
0580         return (vertex0.pt() > vertex1.pt());
0581       });
0582     }
0583   }
0584 
0585   void VertexFinder::sortVerticesInZ0() {
0586     if (settings_->vx_precision() == Precision::Emulation) {
0587       std::sort(
0588           verticesEmulation_.begin(),
0589           verticesEmulation_.end(),
0590           [](const l1t::VertexWord& vertex0, const l1t::VertexWord& vertex1) { return (vertex0.z0() < vertex1.z0()); });
0591     } else {
0592       std::sort(vertices_.begin(), vertices_.end(), [](const RecoVertex<>& vertex0, const RecoVertex<>& vertex1) {
0593         return (vertex0.z0() < vertex1.z0());
0594       });
0595     }
0596   }
0597 
0598   void VertexFinder::associatePrimaryVertex(double trueZ0) {
0599     double distance = 999.;
0600     for (unsigned int id = 0; id < vertices_.size(); ++id) {
0601       if (std::abs(trueZ0 - vertices_[id].z0()) < distance) {
0602         distance = std::abs(trueZ0 - vertices_[id].z0());
0603         pv_index_ = id;
0604       }
0605     }
0606   }
0607 
0608   void VertexFinder::fastHistoLooseAssociation() {
0609     float vxPt = 0.;
0610     RecoVertex leading_vertex;
0611 
0612     for (float z = settings_->vx_histogram_min(); z < settings_->vx_histogram_max();
0613          z += settings_->vx_histogram_binwidth()) {
0614       RecoVertex vertex;
0615       for (const L1Track& track : fitTracks_) {
0616         if (std::abs(z - track.z0()) < settings_->vx_width()) {
0617           vertex.insert(&track);
0618         }
0619       }
0620       computeAndSetVertexParameters(vertex, {}, {});
0621       vertex.setZ0(z);
0622       if (vertex.pt() > vxPt) {
0623         leading_vertex = vertex;
0624         vxPt = vertex.pt();
0625       }
0626     }
0627 
0628     vertices_.emplace_back(leading_vertex);
0629     pv_index_ = 0;  // by default fastHistoLooseAssociation algorithm finds only hard PV
0630   }                 // end of fastHistoLooseAssociation
0631 
0632   void VertexFinder::fastHisto(const TrackerTopology* tTopo) {
0633     // Create the histogram
0634     int nbins =
0635         std::ceil((settings_->vx_histogram_max() - settings_->vx_histogram_min()) / settings_->vx_histogram_binwidth());
0636     std::vector<RecoVertex<>> hist(nbins);
0637     std::vector<RecoVertex<>> sums(nbins - settings_->vx_windowSize());
0638     std::vector<float> bounds(nbins + 1);
0639     strided_iota(std::begin(bounds),
0640                  std::next(std::begin(bounds), nbins + 1),
0641                  settings_->vx_histogram_min(),
0642                  settings_->vx_histogram_binwidth());
0643 
0644     // Loop over the tracks and fill the histogram
0645     for (const L1Track& track : fitTracks_) {
0646       if ((track.z0() < settings_->vx_histogram_min()) || (track.z0() > settings_->vx_histogram_max()))
0647         continue;
0648       if (track.getTTTrackPtr()->chi2() > settings_->vx_TrackMaxChi2())
0649         continue;
0650       if (track.pt() < settings_->vx_TrackMinPt())
0651         continue;
0652 
0653       // Get the number of stubs and the number of stubs in PS layers
0654       float nPS = 0., nstubs = 0;
0655 
0656       // Get pointers to stubs associated to the L1 track
0657       const auto& theStubs = track.getTTTrackPtr()->getStubRefs();
0658       if (theStubs.empty()) {
0659         edm::LogWarning("VertexFinder") << "fastHisto::Could not retrieve the vector of stubs.";
0660         continue;
0661       }
0662 
0663       // Loop over the stubs
0664       for (const auto& stub : theStubs) {
0665         nstubs++;
0666         bool isPS = false;
0667         DetId detId(stub->getDetId());
0668         if (detId.det() == DetId::Detector::Tracker) {
0669           if (detId.subdetId() == StripSubdetector::TOB && tTopo->tobLayer(detId) <= 3)
0670             isPS = true;
0671           else if (detId.subdetId() == StripSubdetector::TID && tTopo->tidRing(detId) <= 9)
0672             isPS = true;
0673         }
0674         if (isPS)
0675           nPS++;
0676       }  // End loop over stubs
0677       if (nstubs < settings_->vx_NStubMin())
0678         continue;
0679       if (nPS < settings_->vx_NStubPSMin())
0680         continue;
0681 
0682       // Quality cuts, may need to be re-optimized
0683       int trk_nstub = (int)track.getTTTrackPtr()->getStubRefs().size();
0684       float chi2dof = track.getTTTrackPtr()->chi2() / (2 * trk_nstub - 4);
0685 
0686       if (settings_->vx_DoPtComp()) {
0687         float trk_consistency = track.getTTTrackPtr()->stubPtConsistency();
0688         if (trk_nstub == 4) {
0689           if (std::abs(track.eta()) < 2.2 && trk_consistency > 10)
0690             continue;
0691           else if (std::abs(track.eta()) > 2.2 && chi2dof > 5.0)
0692             continue;
0693         }
0694       }
0695       if (settings_->vx_DoTightChi2()) {
0696         if (track.pt() > 10.0 && chi2dof > 5.0)
0697           continue;
0698       }
0699 
0700       // Assign the track to the correct vertex
0701       // The values are ordered with bounds [lower, upper)
0702       // Values below bounds.begin() return 0 as the index (underflow)
0703       // Values above bounds.end() will return the index of the last bin (overflow)
0704       auto upper_bound = std::upper_bound(bounds.begin(), bounds.end(), track.z0());
0705       int index = std::distance(bounds.begin(), upper_bound) - 1;
0706       if (index == -1)
0707         index = 0;
0708       hist.at(index).insert(&track);
0709     }  // end loop over tracks
0710 
0711     // Compute the sums
0712     // sliding windows ... sum_i_i+(w-1) where i in (0,nbins-w) and w is the window size
0713     std::vector<float> bin_centers(settings_->vx_windowSize(), 0.0);
0714     std::vector<unsigned int> counts(settings_->vx_windowSize(), 0);
0715     for (unsigned int i = 0; i < sums.size(); i++) {
0716       for (unsigned int j = 0; j < settings_->vx_windowSize(); j++) {
0717         bin_centers[j] = settings_->vx_histogram_min() + ((i + j) * settings_->vx_histogram_binwidth()) +
0718                          (0.5 * settings_->vx_histogram_binwidth());
0719         counts[j] = hist.at(i + j).numTracks();
0720         sums.at(i) += hist.at(i + j);
0721       }
0722       computeAndSetVertexParameters(sums.at(i), bin_centers, counts);
0723     }
0724 
0725     // Find the maxima of the sums
0726     float sigma_max = -999;
0727     int imax = -999;
0728     std::vector<int> found;
0729     found.reserve(settings_->vx_nvtx());
0730     for (unsigned int ivtx = 0; ivtx < settings_->vx_nvtx(); ivtx++) {
0731       sigma_max = -999;
0732       imax = -999;
0733       for (unsigned int i = 0; i < sums.size(); i++) {
0734         // Skip this window if it will already be returned
0735         if (find(found.begin(), found.end(), i) != found.end())
0736           continue;
0737         if (sums.at(i).pt() > sigma_max) {
0738           sigma_max = sums.at(i).pt();
0739           imax = i;
0740         }
0741       }
0742       found.push_back(imax);
0743       vertices_.emplace_back(sums.at(imax));
0744     }
0745     pv_index_ = 0;
0746 
0747     if (settings_->debug() >= 1) {
0748       edm::LogInfo log("VertexProducer");
0749       log << "fastHisto::Checking the output parameters ... \n";
0750       std::vector<double> tmp;
0751       std::transform(std::begin(sums), std::end(sums), std::back_inserter(tmp), [](const RecoVertex<>& v) -> double {
0752         return v.pt();
0753       });
0754       printHistogram<double, edm::LogInfo>(log, tmp, 80, 0, -1, "fastHisto::sums", "\e[92m");
0755       for (unsigned int i = 0; i < found.size(); i++) {
0756         log << "RecoVertex " << i << ": bin index = " << found[i] << "\tsumPt = " << sums.at(imax).pt()
0757             << "\tz0 = " << sums.at(imax).z0();
0758       }
0759     }
0760   }  // end of fastHisto
0761 
0762   void VertexFinder::fastHistoEmulation() {
0763     // Relevant constants for the track word
0764     enum TrackBitWidths {
0765       kZ0Size = 12,             // Width of z-position (40cm / 0.1)
0766       kZ0MagSize = 5,           // Width of z-position magnitude (signed)
0767       kPtSize = 14,             // Width of pt
0768       kPtMagSize = 9,           // Width of pt magnitude (unsigned)
0769       kReducedPrecisionPt = 7,  // Width of the reduced precision, integer only, pt
0770     };
0771 
0772     enum HistogramBitWidths {
0773       kBinSize = 10,                    // Width of a single bin in z
0774       kBinFixedSize = 7,                // Width of a single z0 bin in fixed point representation
0775       kBinFixedMagSize = 4,             // Width (magnitude) of a single z0 bin in fixed point representation
0776       kSlidingSumSize = 11,             // Width of the sum of a window of bins
0777       kInverseSize = 14,                // Width of the inverse sum
0778       kInverseMagSize = 1,              // Width of the inverse sum magnitude (unsigned)
0779       kWeightedSlidingSumSize = 20,     // Width of the pT weighted sliding sum
0780       kWeightedSlidingSumMagSize = 10,  // Width of the pT weighted sliding sum magnitude (signed)
0781       kWindowSize = 3,                  // Number of bins in the window used to sum histogram bins
0782       kSumPtLinkSize = 9,  // Number of bits used to represent the sum of track pts in a single bin from a single link
0783       kSumPtWindowBits = BitsToRepresent(HistogramBitWidths::kWindowSize * (1 << HistogramBitWidths::kSumPtLinkSize)),
0784     };
0785 
0786     static constexpr unsigned int kTableSize =
0787         ((1 << HistogramBitWidths::kSumPtLinkSize) - 1) * HistogramBitWidths::kWindowSize;
0788     static constexpr double kZ0Scale =
0789         (TTTrack_TrackWord::stepZ0 *
0790          (1 << (TrackBitWidths::kZ0Size - TrackBitWidths::kZ0MagSize)));  // scale = 1.27932032
0791 
0792     typedef ap_ufixed<TrackBitWidths::kPtSize, TrackBitWidths::kPtMagSize, AP_RND_CONV, AP_SAT> pt_t;
0793     typedef ap_fixed<TrackBitWidths::kZ0Size, TrackBitWidths::kZ0MagSize, AP_RND_CONV, AP_SAT> z0_t;
0794     // 7 bits chosen to represent values between [0,127]
0795     // This is the next highest power of 2 value to our chosen track pt saturation value (100)
0796     typedef ap_ufixed<TrackBitWidths::kReducedPrecisionPt, TrackBitWidths::kReducedPrecisionPt, AP_RND_INF, AP_SAT>
0797         track_pt_fixed_t;
0798     // Histogram bin index
0799     typedef ap_uint<HistogramBitWidths::kBinSize> histbin_t;
0800     // Histogram bin in fixed point representation, before truncation
0801     typedef ap_ufixed<HistogramBitWidths::kBinFixedSize, HistogramBitWidths::kBinFixedMagSize, AP_RND_INF, AP_SAT>
0802         histbin_fixed_t;
0803     // This value is slightly arbitrary, but small enough that the windows sums aren't too big.
0804     typedef ap_ufixed<HistogramBitWidths::kSumPtLinkSize, HistogramBitWidths::kSumPtLinkSize, AP_RND_INF, AP_SAT>
0805         link_pt_sum_fixed_t;
0806     // Enough bits to store HistogramBitWidths::kWindowSize * (2**HistogramBitWidths::kSumPtLinkSize)
0807     typedef ap_ufixed<HistogramBitWidths::kSumPtWindowBits, HistogramBitWidths::kSumPtWindowBits, AP_RND_INF, AP_SAT>
0808         window_pt_sum_fixed_t;
0809     // pt weighted sum of bins in window
0810     typedef ap_fixed<HistogramBitWidths::kWeightedSlidingSumSize,
0811                      HistogramBitWidths::kWeightedSlidingSumMagSize,
0812                      AP_RND_INF,
0813                      AP_SAT>
0814         zsliding_t;
0815     // Sum of histogram bins in window
0816     typedef ap_uint<HistogramBitWidths::kSlidingSumSize> slidingsum_t;
0817     // Inverse of sum of bins in a given window
0818     typedef ap_ufixed<HistogramBitWidths::kInverseSize, HistogramBitWidths::kInverseMagSize, AP_RND_INF, AP_SAT>
0819         inverse_t;
0820 
0821     auto track_quality_check = [&](const track_pt_fixed_t& pt) -> bool {
0822       // track quality cuts
0823       if (pt.to_double() < settings_->vx_TrackMinPt())
0824         return false;
0825       return true;
0826     };
0827 
0828     auto fetch_bin = [&](const z0_t& z0, int nbins) -> std::pair<histbin_t, bool> {
0829       histbin_t bin = (z0 * histbin_fixed_t(1.0 / settings_->vx_histogram_binwidth())) +
0830                       histbin_t(std::floor(
0831                           nbins / 2.));  // Rounding down (std::floor) taken care of by implicitly casting to ap_uint
0832       if (settings_->debug() > 2) {
0833         edm::LogInfo("VertexProducer")
0834             << "fastHistoEmulation::fetchBin() Checking the mapping from z0 to bin index ... \n"
0835             << "histbin_fixed_t(1.0 / settings_->vx_histogram_binwidth()) = "
0836             << histbin_fixed_t(1.0 / settings_->vx_histogram_binwidth()) << "\n"
0837             << "histbin_t(std::floor(nbins / 2) = " << histbin_t(std::floor(nbins / 2.)) << "\n"
0838             << "z0 = " << z0 << "\n"
0839             << "bin = " << bin;
0840       }
0841       bool valid = true;
0842       if (bin < 0) {
0843         return std::make_pair(0, false);
0844       } else if (bin > (nbins - 1)) {
0845         return std::make_pair(0, false);
0846       }
0847       return std::make_pair(bin, valid);
0848     };
0849 
0850     // Replace with https://stackoverflow.com/questions/13313980/populate-an-array-using-constexpr-at-compile-time ?
0851     auto init_inversion_table = [&]() -> std::vector<inverse_t> {
0852       std::vector<inverse_t> table_out(kTableSize, 0.);
0853       for (unsigned int ii = 0; ii < kTableSize; ii++) {
0854         // First, convert from table index to X-value (unsigned 8-bit, range 0 to +1533)
0855         float in_val = 1533.0 * (ii / float(kTableSize));
0856         // Next, compute lookup table function
0857         table_out.at(ii) = (in_val > 0) ? (1.0 / in_val) : 0.0;
0858       }
0859       return table_out;
0860     };
0861 
0862     auto inversion = [&](slidingsum_t& data_den) -> inverse_t {
0863       std::vector<inverse_t> inversion_table = init_inversion_table();
0864 
0865       // Index into the lookup table based on data
0866       int index;
0867       if (data_den < 0)
0868         data_den = 0;
0869       if (data_den > (kTableSize - 1))
0870         data_den = kTableSize - 1;
0871       index = data_den;
0872       return inversion_table.at(index);
0873     };
0874 
0875     auto bin_center = [&](zsliding_t iz, int nbins) -> z0_t {
0876       zsliding_t z = iz - histbin_t(std::floor(nbins / 2.));
0877       std::unique_ptr<edm::LogInfo> log;
0878       if (settings_->debug() >= 1) {
0879         log = std::make_unique<edm::LogInfo>("VertexProducer");
0880         *log << "bin_center information ...\n"
0881              << "iz = " << iz << "\n"
0882              << "histbin_t(std::floor(nbins / 2.)) = " << histbin_t(std::floor(nbins / 2.)) << "\n"
0883              << "binwidth = " << zsliding_t(settings_->vx_histogram_binwidth()) << "\n"
0884              << "z = " << z << "\n"
0885              << "zsliding_t(z * zsliding_t(binwidth)) = " << std::setprecision(7)
0886              << z0_t(z * zsliding_t(settings_->vx_histogram_binwidth()));
0887       }
0888       return z0_t(z * zsliding_t(settings_->vx_histogram_binwidth()));
0889     };
0890 
0891     auto weighted_position = [&](histbin_t b_max,
0892                                  const std::vector<link_pt_sum_fixed_t>& binpt,
0893                                  slidingsum_t maximums,
0894                                  int nbins) -> zsliding_t {
0895       zsliding_t zvtx_sliding = 0;
0896       slidingsum_t zvtx_sliding_sum = 0;
0897       inverse_t inv = 0;
0898 
0899       std::unique_ptr<edm::LogInfo> log;
0900       if (settings_->debug() >= 1) {
0901         log = std::make_unique<edm::LogInfo>("VertexProducer");
0902         *log << "Progression of weighted_position() ...\n"
0903              << "zvtx_sliding_sum = ";
0904       }
0905 
0906       // Find the weighted position within the window in index space (width = 1)
0907       for (ap_uint<BitsToRepresent(HistogramBitWidths::kWindowSize)> w = 0; w < HistogramBitWidths::kWindowSize; ++w) {
0908         zvtx_sliding_sum += (binpt.at(w) * w);
0909         if (settings_->debug() >= 1) {
0910           *log << "(" << w << " * " << binpt.at(w) << ")";
0911           if (w < HistogramBitWidths::kWindowSize - 1) {
0912             *log << " + ";
0913           }
0914         }
0915       }
0916 
0917       if (settings_->debug() >= 1) {
0918         *log << " = " << zvtx_sliding_sum << "\n";
0919       }
0920 
0921       if (maximums != 0) {
0922         inv = inversion(maximums);
0923         zvtx_sliding = zvtx_sliding_sum * inv;
0924       } else {
0925         zvtx_sliding = (settings_->vx_windowSize() / 2.0) + (((int(settings_->vx_windowSize()) % 2) != 0) ? 0.5 : 0.0);
0926       }
0927       if (settings_->debug() >= 1) {
0928         *log << "inversion(" << maximums << ") = " << inv << "\nzvtx_sliding = " << zvtx_sliding << "\n";
0929       }
0930 
0931       // Add the starting index plus half an index to shift the z position to its weighted position (still in inxex space) within all of the bins
0932       zvtx_sliding += b_max;
0933       zvtx_sliding += ap_ufixed<1, 0>(0.5);
0934       if (settings_->debug() >= 1) {
0935         *log << "b_max = " << b_max << "\n";
0936         *log << "zvtx_sliding + b_max + 0.5 = " << zvtx_sliding << "\n";
0937       }
0938 
0939       // Shift the z position from index space into z [cm] space
0940       zvtx_sliding = bin_center(zvtx_sliding, nbins);
0941       if (settings_->debug() >= 1) {
0942         *log << "bin_center(zvtx_sliding + b_max + 0.5, nbins) = " << std::setprecision(7) << zvtx_sliding;
0943         log.reset();
0944       }
0945       return zvtx_sliding;
0946     };
0947 
0948     // Create the histogram
0949     unsigned int nbins =
0950         std::ceil((settings_->vx_histogram_max() - settings_->vx_histogram_min()) / settings_->vx_histogram_binwidth());
0951     unsigned int nsums = nbins - settings_->vx_windowSize();
0952     std::vector<link_pt_sum_fixed_t> hist(nbins, 0);
0953 
0954     // Loop over the tracks and fill the histogram
0955     if (settings_->debug() > 2) {
0956       edm::LogInfo("VertexProducer") << "fastHistoEmulation::Processing " << fitTracks_.size() << " tracks";
0957     }
0958     for (const L1Track& track : fitTracks_) {
0959       // Get the track pt and z0
0960       // Convert them to an appropriate data format
0961       // Truncation and saturdation taken care of by the data type specification
0962       pt_t tkpt = 0;
0963       tkpt.V = track.getTTTrackPtr()->getTrackWord()(TTTrack_TrackWord::TrackBitLocations::kRinvMSB - 1,
0964                                                      TTTrack_TrackWord::TrackBitLocations::kRinvLSB);
0965       track_pt_fixed_t pt_tmp = tkpt;
0966       //z0_t tkZ0 = track.getTTTrackPtr()->getZ0();
0967       z0_t tkZ0 = 0;
0968       tkZ0.V = track.getTTTrackPtr()->getTrackWord()(TTTrack_TrackWord::TrackBitLocations::kZ0MSB,
0969                                                      TTTrack_TrackWord::TrackBitLocations::kZ0LSB);
0970       ap_ufixed<32, 1> kZ0Scale_fixed = kZ0Scale;
0971       tkZ0 *= kZ0Scale_fixed;
0972 
0973       if ((settings_->vx_DoQualityCuts() && track_quality_check(tkpt)) || (!settings_->vx_DoQualityCuts())) {
0974         //
0975         // Check bin validity of bin found for the current track
0976         //
0977         std::pair<histbin_t, bool> bin = fetch_bin(tkZ0, nbins);
0978         assert(bin.first >= 0 && bin.first < nbins);
0979 
0980         //
0981         // If the bin is valid then sum the tracks
0982         //
0983         if (settings_->debug() > 2) {
0984           edm::LogInfo("VertexProducer") << "fastHistoEmulation::Checking the track word ... \n"
0985                                          << "track word = " << track.getTTTrackPtr()->getTrackWord().to_string(2)
0986                                          << "\n"
0987                                          << "tkZ0 = " << tkZ0.to_double() << "(" << tkZ0.to_string(2)
0988                                          << ")\ttkpt = " << tkpt.to_double() << "(" << tkpt.to_string(2)
0989                                          << ")\tpt_tmp = " << pt_tmp << "\tbin = " << bin.first.to_int() << "\n"
0990                                          << "pt sum in bin " << bin.first.to_int()
0991                                          << " BEFORE adding track = " << hist.at(bin.first).to_double();
0992         }
0993         if (bin.second) {
0994           hist.at(bin.first) = hist.at(bin.first) + pt_tmp;
0995         }
0996         if (settings_->debug() > 2) {
0997           edm::LogInfo("VertexProducer") << "fastHistoEmulation::\npt sum in bin " << bin.first.to_int()
0998                                          << " AFTER adding track = " << hist.at(bin.first).to_double();
0999         }
1000       } else {
1001         if (settings_->debug() > 2) {
1002           edm::LogInfo("VertexProducer") << "fastHistoEmulation::Did not add the following track ... \n"
1003                                          << "track word = " << track.getTTTrackPtr()->getTrackWord().to_string(2)
1004                                          << "\n"
1005                                          << "tkZ0 = " << tkZ0.to_double() << "(" << tkZ0.to_string(2)
1006                                          << ")\ttkpt = " << tkpt.to_double() << "(" << tkpt.to_string(2)
1007                                          << ")\tpt_tmp = " << pt_tmp;
1008         }
1009       }
1010     }  // end loop over tracks
1011 
1012     // Loop through all bins, taking into account the fact that the last bin is nbins-window_width+1,
1013     // and compute the sums using sliding windows ... sum_i_i+(w-1) where i in (0,nbins-w) and w is the window size
1014     std::vector<window_pt_sum_fixed_t> hist_window_sums(nsums, 0);
1015     for (unsigned int b = 0; b < nsums; ++b) {
1016       for (unsigned int w = 0; w < HistogramBitWidths::kWindowSize; ++w) {
1017         unsigned int index = b + w;
1018         hist_window_sums.at(b) += hist.at(index);
1019       }
1020     }
1021 
1022     // Find the top N vertices
1023     std::vector<int> found;
1024     found.reserve(settings_->vx_nvtx());
1025     for (unsigned int ivtx = 0; ivtx < settings_->vx_nvtx(); ivtx++) {
1026       histbin_t b_max = 0;
1027       window_pt_sum_fixed_t max_pt = 0;
1028       zsliding_t zvtx_sliding = -999;
1029       std::vector<link_pt_sum_fixed_t> binpt_max(HistogramBitWidths::kWindowSize, 0);
1030 
1031       // Find the maxima of the sums
1032       for (unsigned int i = 0; i < hist_window_sums.size(); i++) {
1033         // Skip this window if it will already be returned
1034         if (find(found.begin(), found.end(), i) != found.end())
1035           continue;
1036         if (hist_window_sums.at(i) > max_pt) {
1037           b_max = i;
1038           max_pt = hist_window_sums.at(b_max);
1039           std::copy(std::begin(hist) + b_max,
1040                     std::begin(hist) + b_max + HistogramBitWidths::kWindowSize,
1041                     std::begin(binpt_max));
1042 
1043           // Find the weighted position only for the highest sum pt window
1044           zvtx_sliding = weighted_position(b_max, binpt_max, max_pt, nbins);
1045         }
1046       }
1047       if (settings_->debug() >= 1) {
1048         edm::LogInfo log("VertexProducer");
1049         log << "fastHistoEmulation::Checking the output parameters ... \n";
1050         if (found.empty()) {
1051           printHistogram<link_pt_sum_fixed_t, edm::LogInfo>(log, hist, 80, 0, -1, "fastHistoEmulation::hist", "\e[92m");
1052           printHistogram<window_pt_sum_fixed_t, edm::LogInfo>(
1053               log, hist_window_sums, 80, 0, -1, "fastHistoEmulation::hist_window_sums", "\e[92m");
1054         }
1055         printHistogram<link_pt_sum_fixed_t, edm::LogInfo>(
1056             log, binpt_max, 80, 0, -1, "fastHistoEmulation::binpt_max", "\e[92m");
1057         log << "bin index (not a VertexWord parameter) = " << b_max << "\n"
1058             << "sumPt = " << max_pt.to_double() << "\n"
1059             << "z0 = " << zvtx_sliding.to_double();
1060       }
1061       found.push_back(b_max);
1062       verticesEmulation_.emplace_back(l1t::VertexWord::vtxvalid_t(1),
1063                                       l1t::VertexWord::vtxz0_t(zvtx_sliding),
1064                                       l1t::VertexWord::vtxmultiplicity_t(0),
1065                                       l1t::VertexWord::vtxsumpt_t(max_pt),
1066                                       l1t::VertexWord::vtxquality_t(0),
1067                                       l1t::VertexWord::vtxinversemult_t(0),
1068                                       l1t::VertexWord::vtxunassigned_t(0));
1069     }
1070     pv_index_ = 0;
1071   }  // end of fastHistoEmulation
1072 
1073 }  // namespace l1tVertexFinder