Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-09-07 04:37:09 

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() + 1);
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     static constexpr int kZ0Size = 12,  // Width of z-position (40cm / 0.1)
0765         kPtSize = 14,                   // Width of pt
0766         kPtMagSize = 9,                 // Width of pt magnitude (unsigned)
0767         kReducedPrecisionPt = 7         // Width of the reduced precision, integer only, pt
0768         ;
0769 
0770     static constexpr int
0771         kBinSize = 8,                     // Width of a single bin in z
0772         kBinFixedSize = 8,                // Width of a single z0 bin in fixed point representation
0773         kBinFixedMagSize = 5,             // Width (magnitude) of a single z0 bin in fixed point representation
0774         kSlidingSumSize = 11,             // Width of the sum of a window of bins
0775         kInverseSize = 14,                // Width of the inverse sum
0776         kInverseMagSize = 1,              // Width of the inverse sum magnitude (unsigned)
0777         kWeightedSlidingSumSize = 20,     // Width of the pT weighted sliding sum
0778         kWeightedSlidingSumMagSize = 10,  // Width of the pT weighted sliding sum magnitude (signed)
0779         kWindowSize = 3,                  // Number of bins in the window used to sum histogram bins
0780         kSumPtLinkSize = 9,  // Number of bits used to represent the sum of track pts in a single bin from a single link
0781 
0782         kSumPtWindowBits = BitsToRepresent(kWindowSize * (1 << kSumPtLinkSize)),
0783         // Number of bits to represent the untruncated sum of track pts in a single bin from a single link
0784         kSumPtUntruncatedLinkSize = kPtSize + 2, kSumPtUntruncatedLinkMagSize = kPtMagSize + 2;
0785 
0786     static constexpr unsigned int kTableSize = ((1 << kSumPtLinkSize) - 1) * kWindowSize;
0787 
0788     typedef ap_ufixed<kPtSize, kPtMagSize, AP_RND_CONV, AP_SAT> pt_t;
0789     // Same size as TTTrack_TrackWord::z0_t, but now taking into account the sign bit (i.e. 2's complement)
0790     typedef ap_int<kZ0Size> z0_t;
0791     // 7 bits chosen to represent values between [0,127]
0792     // This is the next highest power of 2 value to our chosen track pt saturation value (100)
0793     typedef ap_ufixed<kReducedPrecisionPt, kReducedPrecisionPt, AP_RND_INF, AP_SAT> track_pt_fixed_t;
0794     // Histogram bin index
0795     typedef ap_uint<kBinSize> histbin_t;
0796     // Histogram bin in fixed point representation, before truncation
0797     typedef ap_ufixed<kBinFixedSize, kBinFixedMagSize, AP_RND_INF, AP_SAT> histbin_fixed_t;
0798     // This type is slightly arbitrary, but 2 bits larger than untruncated track pt to store sums in histogram bins
0799     // with truncation just before vertex-finding
0800     typedef ap_ufixed<kSumPtUntruncatedLinkSize, kSumPtUntruncatedLinkMagSize, AP_RND_INF, AP_SAT>
0801         histbin_pt_sum_fixed_t;
0802     // This value is slightly arbitrary, but small enough that the windows sums aren't too big.
0803     typedef ap_ufixed<kSumPtLinkSize, kSumPtLinkSize, AP_RND_INF, AP_SAT> link_pt_sum_fixed_t;
0804     // Enough bits to store kWindowSize * (2**kSumPtLinkSize)
0805     typedef ap_ufixed<kSumPtWindowBits, kSumPtWindowBits, AP_RND_INF, AP_SAT> window_pt_sum_fixed_t;
0806     // pt weighted sum of bins in window
0807     typedef ap_fixed<kWeightedSlidingSumSize, kWeightedSlidingSumMagSize, AP_RND_INF, AP_SAT> zsliding_t;
0808     // Sum of histogram bins in window
0809     typedef ap_uint<kSlidingSumSize> slidingsum_t;
0810     // Inverse of sum of bins in a given window
0811     typedef ap_ufixed<kInverseSize, kInverseMagSize, AP_RND_INF, AP_SAT> inverse_t;
0812 
0813     auto track_quality_check = [&](const track_pt_fixed_t& pt) -> bool {
0814       // Track quality cuts
0815       if (pt.to_double() < settings_->vx_TrackMinPt())
0816         return false;
0817       return true;
0818     };
0819 
0820     auto fetch_bin = [&](const z0_t& z0, int nbins) -> std::pair<histbin_t, bool> {
0821       // Increase the the number of bits in the word to allow for additional dynamic range
0822       ap_int<kZ0Size + 1> z0_13 = z0;
0823       // Add a number equal to half of the range in z0, meaning that the range is now [0, 2*z0_max]
0824       ap_int<kZ0Size + 1> absz0_13 = z0_13 + (1 << (kZ0Size - 1));
0825       // Shift the bits down to truncate the dynamic range to the most significant kBinFixedSize bits
0826       ap_int<kZ0Size + 1> absz0_13_reduced = absz0_13 >> (kZ0Size - kBinFixedSize);
0827       // Put the relevant bits into the histbin_t container
0828       histbin_t bin = absz0_13_reduced.range(kBinFixedSize - 1, 0);
0829 
0830       if (settings_->debug() > 2) {
0831         edm::LogInfo("VertexProducer")
0832             << "fastHistoEmulation::fetchBin() Checking the mapping from z0 to bin index ... \n"
0833             << "histbin_fixed_t(1.0 / settings_->vx_histogram_binwidth()) = "
0834             << histbin_fixed_t(1.0 / settings_->vx_histogram_binwidth()) << "\n"
0835             << "histbin_t(std::floor(nbins / 2) = " << histbin_t(std::floor(nbins / 2.)) << "\n"
0836             << "z0 = " << z0 << "\n"
0837             << "bin = " << bin;
0838       }
0839       bool valid = true;
0840       if (bin < 0) {
0841         return std::make_pair(0, false);
0842       } else if (bin > (nbins - 1)) {
0843         return std::make_pair(0, false);
0844       }
0845       return std::make_pair(bin, valid);
0846     };
0847 
0848     // Replace with https://stackoverflow.com/questions/13313980/populate-an-array-using-constexpr-at-compile-time ?
0849     auto init_inversion_table = [&]() -> std::vector<inverse_t> {
0850       std::vector<inverse_t> table_out(kTableSize, 0.);
0851       for (unsigned int ii = 0; ii < kTableSize; ii++) {
0852         // Compute lookup table function. This matches the format of the GTT HLS code.
0853         // Biased generation f(x) = 1 / (x + 1) is inverted by g(y) = inversion(x - 1) = 1 / (x - 1 + 1) = 1 / y
0854         table_out.at(ii) = (1.0 / (ii + 1));
0855       }
0856       return table_out;
0857     };
0858 
0859     auto inversion = [&](slidingsum_t& data_den) -> inverse_t {
0860       std::vector<inverse_t> inversion_table = init_inversion_table();
0861 
0862       // Index into the lookup table based on data
0863       int index;
0864       if (data_den < 0)
0865         data_den = 0;
0866       if (data_den > (kTableSize - 1))
0867         data_den = kTableSize - 1;
0868       index = data_den;
0869       return inversion_table.at(index);
0870     };
0871 
0872     auto bin_center = [&](zsliding_t iz, int nbins) -> l1t::VertexWord::vtxz0_t {
0873       zsliding_t z = iz - histbin_t(std::floor(nbins / 2.));
0874       std::unique_ptr<edm::LogInfo> log;
0875       if (settings_->debug() >= 1) {
0876         log = std::make_unique<edm::LogInfo>("VertexProducer");
0877         *log << "bin_center information ...\n"
0878              << "iz = " << iz << "\n"
0879              << "histbin_t(std::floor(nbins / 2.)) = " << histbin_t(std::floor(nbins / 2.)) << "\n"
0880              << "binwidth = " << zsliding_t(settings_->vx_histogram_binwidth()) << "\n"
0881              << "z = " << z << "\n"
0882              << "zsliding_t(z * zsliding_t(binwidth)) = " << std::setprecision(7)
0883              << l1t::VertexWord::vtxz0_t(z * zsliding_t(settings_->vx_histogram_binwidth()));
0884       }
0885       return l1t::VertexWord::vtxz0_t(z * zsliding_t(settings_->vx_histogram_binwidth()));
0886     };
0887 
0888     auto weighted_position = [&](histbin_t b_max,
0889                                  const std::vector<link_pt_sum_fixed_t>& binpt,
0890                                  slidingsum_t maximums,
0891                                  int nbins) -> zsliding_t {
0892       zsliding_t zvtx_sliding = 0;
0893       slidingsum_t zvtx_sliding_sum = 0;
0894       inverse_t inv = 0;
0895 
0896       std::unique_ptr<edm::LogInfo> log;
0897       if (settings_->debug() >= 1) {
0898         log = std::make_unique<edm::LogInfo>("VertexProducer");
0899         *log << "Progression of weighted_position() ...\n"
0900              << "zvtx_sliding_sum = ";
0901       }
0902 
0903       // Find the weighted position within the window in index space (width = 1)
0904       for (ap_uint<BitsToRepresent(kWindowSize)> w = 0; w < kWindowSize; ++w) {
0905         zvtx_sliding_sum += (binpt.at(w) * w);
0906         if (settings_->debug() >= 1) {
0907           *log << "(" << w << " * " << binpt.at(w) << ")";
0908           if (w < kWindowSize - 1) {
0909             *log << " + ";
0910           }
0911         }
0912       }
0913 
0914       if (settings_->debug() >= 1) {
0915         *log << " = " << zvtx_sliding_sum << "\n";
0916       }
0917 
0918       if (maximums != 0) {
0919         //match F/W inversion_lut offset (inversion[x] = 1 / (x + 1); inversion[x - 1] = 1 / x;), for consistency
0920         slidingsum_t offsetmaximums = maximums - 1;
0921         inv = inversion(offsetmaximums);
0922         zvtx_sliding = zvtx_sliding_sum * inv;
0923       } else {
0924         zvtx_sliding = (settings_->vx_windowSize() / 2.0) + (((int(settings_->vx_windowSize()) % 2) != 0) ? 0.5 : 0.0);
0925       }
0926       if (settings_->debug() >= 1) {
0927         *log << "inversion(" << maximums << ") = " << inv << "\nzvtx_sliding = " << zvtx_sliding << "\n";
0928       }
0929 
0930       // 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
0931       zvtx_sliding += b_max;
0932       zvtx_sliding += ap_ufixed<1, 0>(0.5);
0933       if (settings_->debug() >= 1) {
0934         *log << "b_max = " << b_max << "\n";
0935         *log << "zvtx_sliding + b_max + 0.5 = " << zvtx_sliding << "\n";
0936       }
0937 
0938       // Shift the z position from index space into z [cm] space
0939       zvtx_sliding = bin_center(zvtx_sliding, nbins);
0940       if (settings_->debug() >= 1) {
0941         *log << "bin_center(zvtx_sliding + b_max + 0.5, nbins) = " << std::setprecision(7) << zvtx_sliding;
0942         log.reset();
0943       }
0944       return zvtx_sliding;
0945     };
0946 
0947     // Create the histogram
0948     unsigned int nbins = std::round((settings_->vx_histogram_max() - settings_->vx_histogram_min()) /
0949                                     settings_->vx_histogram_binwidth());
0950     unsigned int nsums = nbins - settings_->vx_windowSize() + 1;
0951     std::vector<link_pt_sum_fixed_t> hist(nbins, 0);
0952     std::vector<histbin_pt_sum_fixed_t> hist_untruncated(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 saturation taken care of by the data type specification, now delayed to end of histogramming
0962       pt_t tkpt = 0;
0963       tkpt.V = track.getTTTrackPtr()->getTrackWord()(TTTrack_TrackWord::TrackBitLocations::kRinvMSB - 1,
0964                                                      TTTrack_TrackWord::TrackBitLocations::kRinvLSB);
0965       z0_t tkZ0 = track.getTTTrackPtr()->getZ0Word();
0966 
0967       if ((settings_->vx_DoQualityCuts() && track_quality_check(tkpt)) || (!settings_->vx_DoQualityCuts())) {
0968         //
0969         // Check bin validity of bin found for the current track
0970         //
0971         std::pair<histbin_t, bool> bin = fetch_bin(tkZ0, nbins);
0972         assert(bin.first >= 0 && bin.first < nbins);
0973 
0974         //
0975         // If the bin is valid then sum the tracks
0976         //
0977         if (settings_->debug() > 2) {
0978           edm::LogInfo("VertexProducer") << "fastHistoEmulation::Checking the track word ... \n"
0979                                          << "track word = " << track.getTTTrackPtr()->getTrackWord().to_string(2)
0980                                          << "\n"
0981                                          << "tkZ0 = " << tkZ0.to_double() << "(" << tkZ0.to_string(2)
0982                                          << ")\ttkpt = " << tkpt.to_double() << "(" << tkpt.to_string(2)
0983                                          << ")\tbin = " << bin.first.to_int() << "\n"
0984                                          << "pt sum in bin " << bin.first.to_int()
0985                                          << " BEFORE adding track = " << hist_untruncated.at(bin.first).to_double();
0986         }
0987         if (bin.second) {
0988           hist_untruncated.at(bin.first) = hist_untruncated.at(bin.first) + tkpt;
0989         }
0990         if (settings_->debug() > 2) {
0991           edm::LogInfo("VertexProducer") << "fastHistoEmulation::\npt sum in bin " << bin.first.to_int()
0992                                          << " AFTER adding track = " << hist_untruncated.at(bin.first).to_double();
0993         }
0994       } else {
0995         if (settings_->debug() > 2) {
0996           edm::LogInfo("VertexProducer") << "fastHistoEmulation::Did not add the following track ... \n"
0997                                          << "track word = " << track.getTTTrackPtr()->getTrackWord().to_string(2)
0998                                          << "\n"
0999                                          << "tkZ0 = " << tkZ0.to_double() << "(" << tkZ0.to_string(2)
1000                                          << ")\ttkpt = " << tkpt.to_double() << "(" << tkpt.to_string(2) << ")";
1001         }
1002       }
1003     }  // end loop over tracks
1004 
1005     // HLS histogramming used to truncate track pt before adding, using
1006     // track_pt_fixed_t pt_tmp = tkpt;
1007     // Now, truncation should happen after histograms are filled but prior to the vertex-finding part of the algo
1008     for (unsigned int hb = 0; hb < hist.size(); ++hb) {
1009       link_pt_sum_fixed_t bin_trunc = hist_untruncated.at(hb).range(
1010           kSumPtUntruncatedLinkSize - 1, kSumPtUntruncatedLinkSize - kSumPtUntruncatedLinkMagSize);
1011       hist.at(hb) = bin_trunc;
1012       if (settings_->debug() > 2) {
1013         edm::LogInfo("VertexProducer") << "fastHistoEmulation::truncating histogram bin pt once filling is complete \n"
1014                                        << "hist_untruncated.at(" << hb << ") = " << hist_untruncated.at(hb).to_double()
1015                                        << "(" << hist_untruncated.at(hb).to_string(2)
1016                                        << ")\tbin_trunc = " << bin_trunc.to_double() << "(" << bin_trunc.to_string(2)
1017                                        << ")\n\thist.at(" << hb << ") = " << hist.at(hb).to_double() << "("
1018                                        << hist.at(hb).to_string(2) << ")";
1019       }
1020     }
1021 
1022     // Loop through all bins, taking into account the fact that the last bin is nbins-window_width+1,
1023     // and compute the sums using sliding windows ... sum_i_i+(w-1) where i in (0,nbins-w) and w is the window size
1024     std::vector<window_pt_sum_fixed_t> hist_window_sums(nsums, 0);
1025     for (unsigned int b = 0; b < nsums; ++b) {
1026       for (unsigned int w = 0; w < kWindowSize; ++w) {
1027         unsigned int index = b + w;
1028         hist_window_sums.at(b) += hist.at(index);
1029       }
1030     }
1031 
1032     // Find the top N vertices
1033     std::vector<int> found;
1034     found.reserve(settings_->vx_nvtx());
1035     for (unsigned int ivtx = 0; ivtx < settings_->vx_nvtx(); ivtx++) {
1036       histbin_t b_max = 0;
1037       window_pt_sum_fixed_t max_pt = 0;
1038       zsliding_t zvtx_sliding = -999;
1039       std::vector<link_pt_sum_fixed_t> binpt_max(kWindowSize, 0);
1040 
1041       // Find the maxima of the sums
1042       for (unsigned int i = 0; i < hist_window_sums.size(); i++) {
1043         // Skip this window if it will already be returned
1044         if (find(found.begin(), found.end(), i) != found.end())
1045           continue;
1046         if (hist_window_sums.at(i) > max_pt) {
1047           b_max = i;
1048           max_pt = hist_window_sums.at(b_max);
1049           std::copy(std::begin(hist) + b_max, std::begin(hist) + b_max + kWindowSize, std::begin(binpt_max));
1050 
1051           // Find the weighted position only for the highest sum pt window
1052           zvtx_sliding = weighted_position(b_max, binpt_max, max_pt, nbins);
1053         }
1054       }
1055       if (settings_->debug() >= 1) {
1056         edm::LogInfo log("VertexProducer");
1057         log << "fastHistoEmulation::Checking the output parameters ... \n";
1058         if (found.empty()) {
1059           printHistogram<link_pt_sum_fixed_t, edm::LogInfo>(log, hist, 80, 0, -1, "fastHistoEmulation::hist", "\e[92m");
1060           printHistogram<window_pt_sum_fixed_t, edm::LogInfo>(
1061               log, hist_window_sums, 80, 0, -1, "fastHistoEmulation::hist_window_sums", "\e[92m");
1062         }
1063         printHistogram<link_pt_sum_fixed_t, edm::LogInfo>(
1064             log, binpt_max, 80, 0, -1, "fastHistoEmulation::binpt_max", "\e[92m");
1065         log << "bin index (not a VertexWord parameter) = " << b_max << "\n"
1066             << "sumPt = " << max_pt.to_double() << "\n"
1067             << "z0 = " << zvtx_sliding.to_double();
1068       }
1069       found.push_back(b_max);
1070       verticesEmulation_.emplace_back(l1t::VertexWord::vtxvalid_t(1),
1071                                       l1t::VertexWord::vtxz0_t(zvtx_sliding),
1072                                       l1t::VertexWord::vtxmultiplicity_t(0),
1073                                       l1t::VertexWord::vtxsumpt_t(max_pt),
1074                                       l1t::VertexWord::vtxquality_t(0),
1075                                       l1t::VertexWord::vtxinversemult_t(0),
1076                                       l1t::VertexWord::vtxunassigned_t(0));
1077     }
1078     pv_index_ = 0;
1079   }  // end of fastHistoEmulation
1080 
1081   void VertexFinder::NNVtxEmulation(tensorflow::Session* TrackWeightSesh,
1082                                     tensorflow::Session* PatternRecSesh,
1083                                     tensorflow::Session* AssociationSesh) {
1084     // #### Weight Tracks: ####
1085     // Loop over tracks -> weight the network -> set track weights
1086     tensorflow::Tensor inputTrkWeight(tensorflow::DT_FLOAT, {1, 3});  //Single batch of 3 values
1087     uint counter = 0;
1088 
1089     for (auto& track : fitTracks_) {
1090       // Quantised Network: Use values from L1GTTInputProducer pT, MVA1, eta
1091       auto& gttTrack = fitTracks_.at(counter);
1092 
1093       TTTrack_TrackWord::tanl_t etaEmulationBits = gttTrack.getTTTrackPtr()->getTanlWord();
1094       ap_fixed<16, 3> etaEmulation;
1095       etaEmulation.V = (etaEmulationBits.range());
1096 
1097       ap_uint<14> ptEmulationBits = gttTrack.getTTTrackPtr()->getTrackWord()(
1098           TTTrack_TrackWord::TrackBitLocations::kRinvMSB - 1, TTTrack_TrackWord::TrackBitLocations::kRinvLSB);
1099       ap_ufixed<14, 9> ptEmulation;
1100       ptEmulation.V = (ptEmulationBits.range());
1101 
1102       ap_ufixed<22, 9> ptEmulation_rescale;
1103       ptEmulation_rescale = ptEmulation.to_double();
1104 
1105       ap_ufixed<22, 9> etaEmulation_rescale;
1106       etaEmulation_rescale = abs(etaEmulation.to_double());
1107 
1108       ap_ufixed<22, 9> MVAEmulation_rescale;
1109       MVAEmulation_rescale = gttTrack.getTTTrackPtr()->getMVAQualityBits();
1110 
1111       inputTrkWeight.tensor<float, 2>()(0, 0) = ptEmulation_rescale.to_double();
1112       inputTrkWeight.tensor<float, 2>()(0, 1) = MVAEmulation_rescale.to_double();
1113       inputTrkWeight.tensor<float, 2>()(0, 2) = etaEmulation_rescale.to_double();
1114 
1115       // CNN output: track weight
1116       std::vector<tensorflow::Tensor> outputTrkWeight;
1117       tensorflow::run(TrackWeightSesh, {{"weight:0", inputTrkWeight}}, {"Identity:0"}, &outputTrkWeight);
1118       // Set track weight pack into tracks:
1119 
1120       ap_ufixed<16, 5> NNOutput;
1121       NNOutput = (double)outputTrkWeight[0].tensor<float, 2>()(0, 0);
1122 
1123       //std::cout<<"NNOutput_weight_network = "<< NNOutput <<std::endl;
1124 
1125       track.setWeight(NNOutput.to_double());
1126 
1127       ++counter;
1128     }
1129     // #### Find Vertices: ####
1130     tensorflow::Tensor inputPV(tensorflow::DT_FLOAT,
1131                                {1, settings_->vx_histogram_numbins(), 1});  //Single batch with 256 bins and 1 weight
1132     std::vector<tensorflow::Tensor> outputPV;
1133     RecoVertexCollection vertices(settings_->vx_histogram_numbins());
1134     std::map<float, int> vertexMap;
1135     std::map<int, float> histogram;
1136     std::map<int, float> nnOutput;
1137 
1138     float binWidth = settings_->vx_histogram_binwidth();
1139 
1140     int track_z = 0;
1141 
1142     for (int z = 0; z < settings_->vx_histogram_numbins(); z += 1) {
1143       counter = 0;
1144       double vxWeight = 0;
1145 
1146       for (const L1Track& track : fitTracks_) {
1147         auto& gttTrack = fitTracks_.at(counter);
1148         double temp_z0 = gttTrack.getTTTrackPtr()->z0();
1149 
1150         track_z = std::floor((temp_z0 + settings_->vx_histogram_max()) / binWidth);
1151 
1152         if (track_z >= z && track_z < (z + 1)) {
1153           vertices.at(z).insert(&track);
1154           vxWeight += track.weight();
1155         }
1156         ++counter;
1157       }
1158       // Get centre of bin before setting z0
1159       vertices.at(z).setZ0(((z + 0.5) * binWidth) - settings_->vx_histogram_max());
1160 
1161       vertexMap[vxWeight] = z;
1162       inputPV.tensor<float, 3>()(0, z, 0) = vxWeight;
1163       //Fill histogram for 3 bin sliding window:
1164       histogram[z] = vxWeight;
1165     }
1166 
1167     // Run PV Network:
1168     tensorflow::run(PatternRecSesh, {{"hist:0", inputPV}}, {"Identity:0"}, &outputPV);
1169     // Threshold needed due to rounding differences in internal CNN layer emulation versus firmware
1170     const float histogrammingThreshold_ = 0.0;
1171     for (int i(0); i < settings_->vx_histogram_numbins(); ++i) {
1172       if (outputPV[0].tensor<float, 3>()(0, i, 0) >= histogrammingThreshold_) {
1173         nnOutput[i] = outputPV[0].tensor<float, 3>()(0, i, 0);
1174       }
1175     }
1176 
1177     //Find max then find all occurances of it in histogram and average their position -> python argmax layer
1178     //Performance is not optimised for multiple peaks in histogram or spread peaks these are edge cases, need to revisit
1179     int max_index = 0;
1180     int num_maxes = 0;
1181     float max_element = 0.0;
1182     for (int i(0); i < settings_->vx_histogram_numbins(); ++i) {
1183       if (nnOutput[i] > max_element) {
1184         max_element = nnOutput[i];
1185       }
1186     }
1187 
1188     for (int i(0); i < settings_->vx_histogram_numbins(); ++i) {
1189       if (nnOutput[i] == max_element) {
1190         num_maxes++;
1191         max_index += i;
1192       }
1193     }
1194     int PV_index = ceil((float)max_index / (float)num_maxes);
1195 
1196     edm::LogVerbatim("VertexFinder") << " NNVtxEmulation Chosen PV: prob: " << nnOutput[PV_index]
1197                                      << " bin = " << PV_index << " z0 = " << vertices.at(PV_index).z0() << '\n';
1198 
1199     verticesEmulation_.emplace_back(1, vertices.at(PV_index).z0(), 0, vertices.at(PV_index).pt(), 0, 0, 0);
1200 
1201   }  // end of NNVtx Algorithm
1202 
1203 }  // namespace l1tVertexFinder
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!