l1tVertexFinder::VertexFinder Class Reference

#include <VertexFinder.h>

Classes
struct	SortTracksByPt
struct	SortTracksByZ0
	Helper structs/classes. More...

Public Member Functions
void	adaptiveVertexReconstruction ()
	Adaptive Vertex Reconstruction algorithm. More...
void	agglomerativeHierarchicalClustering ()
	Simple Merge Algorithm. More...
void	associatePrimaryVertex (double trueZ0)
	Associate the primary vertex with the real one. More...
float	centralDistance (RecoVertex<> cluster0, RecoVertex<> cluster1)
	Find distance between centres of two clusters. More...
void	computeAndSetVertexParameters (RecoVertex<> &vertex, const std::vector< float > &bin_centers, const std::vector< unsigned int > &counts)
	Vertexing algorithms. More...
void	DBSCAN ()
	DBSCAN algorithm. More...
void	fastHisto (const TrackerTopology *tTopo)
	Histogramming algorithm. More...
void	fastHistoEmulation ()
	Histogramming algorithm (emulation) More...
void	fastHistoLooseAssociation ()
	High pT Vertex Algorithm. More...
void	findPrimaryVertex ()
	Find the primary vertex. More...
const FitTrackCollection &	fitTracks () const
	Accessors. More...
void	GapClustering ()
	Gap Clustering Algorithm. More...
void	HPV ()
	High pT Vertex Algorithm. More...
unsigned int	iterationsPerTrack () const
	Number of iterations. More...
void	Kmeans ()
	Kmeans Algorithm. More...
float	maxDistance (RecoVertex<> cluster0, RecoVertex<> cluster1)
	Find maximum distance in two clusters of tracks. More...
float	meanDistance (RecoVertex<> cluster0, RecoVertex<> cluster1)
	Find average distance in two clusters of tracks. More...
float	minDistance (RecoVertex<> cluster0, RecoVertex<> cluster1)
	Find minimum distance in two clusters of tracks. More...
void	NNVtxEmulation (tensorflow::Session *TrackWeightSesh=nullptr, tensorflow::Session *PatternRecSesh=nullptr, tensorflow::Session *AssociationSesh=nullptr)
	NNVtx algorithm. More...
unsigned int	numInputTracks () const
	Storage for tracks out of the L1 Track finder. More...
unsigned int	numIterations () const
	Number of iterations. More...
unsigned int	numVertices () const
	Number of reconstructed vertices. More...
unsigned int	numVerticesEmulation () const
	Number of emulation vertices. More...
template<typename T >
T	PrimaryVertex () const
	Reconstructed primary vertex. More...
unsigned int	primaryVertexId () const
	Reconstructed Primary Vertex Id. More...
template<class data_type , typename stream_type = std::ostream>
void	printHistogram (stream_type &stream, std::vector< data_type > data, int width=80, int minimum=0, int maximum=-1, std::string title="", std::string color="")
	Print an ASCII histogram. More...
void	PVR ()
	Find maximum distance in two clusters of tracks. More...
void	sortVerticesInPt ()
	Sort vertices in pT. More...
void	sortVerticesInZ0 ()
	Sort vertices in z. More...
template<typename ForwardIterator , typename T >
void	strided_iota (ForwardIterator first, ForwardIterator last, T value, T stride)
	VertexFinder (FitTrackCollection &fitTracks, const AlgoSettings &settings)
	Constructor and destructor. More...
const std::vector< RecoVertex<> > &	vertices () const
	Returns the z positions of the reconstructed primary vertices. More...
const l1t::VertexWordCollection &	verticesEmulation () const
	Returns the emulation primary vertices. More...
	~VertexFinder ()

Private Attributes
FitTrackCollection	fitTracks_
unsigned int	iterations_
unsigned int	numMatchedVertices_
unsigned int	pv_index_
const AlgoSettings *	settings_
RecoVertexCollection	vertices_
l1t::VertexWordCollection	verticesEmulation_

Detailed Description

Definition at line 27 of file VertexFinder.h.

Constructor & Destructor Documentation

VertexFinder()

l1tVertexFinder::VertexFinder::VertexFinder ( FitTrackCollection &  fitTracks, const AlgoSettings &  settings  )

inline

Constructor and destructor.

Definition at line 30 of file VertexFinder.h.

References fitTracks(), fitTracks_, and settings_.

                                                                              {
      fitTracks_ = fitTracks;
      settings_ = &settings;
    }

~VertexFinder()

l1tVertexFinder::VertexFinder::~VertexFinder ( )

inline

Definition at line 34 of file VertexFinder.h.

{}

Member Function Documentation

adaptiveVertexReconstruction()

void l1tVertexFinder::VertexFinder::adaptiveVertexReconstruction

(

)

Adaptive Vertex Reconstruction algorithm.

Definition at line 311 of file VertexFinder.cc.

References funct::abs(), nano_mu_local_reco_cff::chi2, ALPAKA_ACCELERATOR_NAMESPACE::vertexFinder::it, or, HLT_2024v13_cff::track, bphysicsOniaDQM_cfi::vertex, and l1tVertexFinder::L1Track::z0().

Referenced by VertexProducer::produce().

                                                  {
    bool start = true;
    iterations_ = 0;
    FitTrackCollection discardedTracks, acceptedTracks, discardedTracks2;

    for (const L1Track& track : fitTracks_) {
      discardedTracks.push_back(track);
    }

    while (discardedTracks.size() >= settings_->vx_minTracks() or start == true) {
      start = false;
      discardedTracks2.clear();
      FitTrackCollection::iterator it = discardedTracks.begin();
      const L1Track track = *it;
      acceptedTracks.push_back(track);
      float z0sum = track.z0();

      for (FitTrackCollection::iterator it2 = discardedTracks.begin(); it2 < discardedTracks.end(); ++it2) {
        if (it2 != it) {
          const L1Track secondTrack = *it2;
          // Calculate new vertex z0 adding this track
          z0sum += secondTrack.z0();
          float z0vertex = z0sum / (acceptedTracks.size() + 1);
          // Calculate chi2 of new vertex
          float chi2 = 0.;
          float dof = 0.;
          for (const L1Track& accTrack : acceptedTracks) {
            iterations_++;
            float Residual = accTrack.z0() - z0vertex;
            if (std::abs(accTrack.eta()) < 1.2)
              Residual /= 0.1812;  // Assumed z0 resolution
            else if (std::abs(accTrack.eta()) >= 1.2 && std::abs(accTrack.eta()) < 1.6)
              Residual /= 0.2912;
            else if (std::abs(accTrack.eta()) >= 1.6 && std::abs(accTrack.eta()) < 2.)
              Residual /= 0.4628;
            else
              Residual /= 0.65;

            chi2 += Residual * Residual;
            dof = (acceptedTracks.size() + 1) * 2 - 1;
          }
          if (chi2 / dof < settings_->vx_chi2cut()) {
            acceptedTracks.push_back(secondTrack);
          } else {
            discardedTracks2.push_back(secondTrack);
            z0sum -= secondTrack.z0();
          }
        }
      }

      if (acceptedTracks.size() >= settings_->vx_minTracks()) {
        RecoVertex vertex;
        for (const L1Track& track : acceptedTracks) {
          vertex.insert(&track);
        }
        computeAndSetVertexParameters(vertex, {}, {});
        vertices_.push_back(vertex);
      }

      acceptedTracks.clear();
      discardedTracks.clear();
      discardedTracks = discardedTracks2;
    }
  }

agglomerativeHierarchicalClustering()

void l1tVertexFinder::VertexFinder::agglomerativeHierarchicalClustering

(

)

Simple Merge Algorithm.

Definition at line 129 of file VertexFinder.cc.

References mps_fire::i, particleFlowClusterHGC_cfi::maxDistance, or, jetUpdater_cfi::sort, HLT_2024v13_cff::track, and DiMuonV_cfg::tracks.

Referenced by VertexProducer::produce().

                                                         {
    iterations_ = 0;

    sort(fitTracks_.begin(), fitTracks_.end(), SortTracksByZ0());

    RecoVertexCollection vClusters;
    vClusters.resize(fitTracks_.size());

    for (unsigned int i = 0; i < fitTracks_.size(); ++i) {
      vClusters[i].insert(&fitTracks_[i]);
      // iterations_++;
    }

    while (true) {
      float MinimumScore = 9999;

      unsigned int clusterId0 = 0;
      unsigned int clusterId1 = 0;
      for (unsigned int iClust = 0; iClust < vClusters.size() - 1; iClust++) {
        iterations_++;

        float M = 0;
        if (settings_->vx_distanceType() == 0)
          M = maxDistance(vClusters[iClust], vClusters[iClust + 1]);
        else if (settings_->vx_distanceType() == 1)
          M = minDistance(vClusters[iClust], vClusters[iClust + 1]);
        else if (settings_->vx_distanceType() == 2)
          M = meanDistance(vClusters[iClust], vClusters[iClust + 1]);
        else
          M = centralDistance(vClusters[iClust], vClusters[iClust + 1]);

        if (M < MinimumScore) {
          MinimumScore = M;
          clusterId0 = iClust;
          clusterId1 = iClust + 1;
        }
      }
      if (MinimumScore > settings_->vx_distance() or vClusters[clusterId1].tracks().empty())
        break;
      for (const L1Track* track : vClusters[clusterId0].tracks()) {
        vClusters[clusterId1].insert(track);
      }
      vClusters.erase(vClusters.begin() + clusterId0);
    }

    for (RecoVertex clust : vClusters) {
      if (clust.numTracks() >= settings_->vx_minTracks()) {
        computeAndSetVertexParameters(clust, {}, {});
        vertices_.push_back(clust);
      }
    }
  }

associatePrimaryVertex()

void l1tVertexFinder::VertexFinder::associatePrimaryVertex

(

double

trueZ0

)

Associate the primary vertex with the real one.

Definition at line 598 of file VertexFinder.cc.

References funct::abs(), HLT_2024v13_cff::distance, and l1ctLayer2EG_cff::id.

                                                         {
    double distance = 999.;
    for (unsigned int id = 0; id < vertices_.size(); ++id) {
      if (std::abs(trueZ0 - vertices_[id].z0()) < distance) {
        distance = std::abs(trueZ0 - vertices_[id].z0());
        pv_index_ = id;
      }
    }
  }

centralDistance()

float l1tVertexFinder::VertexFinder::centralDistance	(	RecoVertex<>	cluster0,
		RecoVertex<>	cluster1
	)

Find distance between centres of two clusters.

Definition at line 122 of file VertexFinder.cc.

References funct::abs(), HLT_2024v13_cff::distance, and l1tVertexFinder::RecoVertex< T >::z0().

                                                                                  {
    computeAndSetVertexParameters(cluster0, {}, {});
    computeAndSetVertexParameters(cluster1, {}, {});
    float distance = std::abs(cluster0.z0() - cluster1.z0());
    return distance;
  }

computeAndSetVertexParameters()

void l1tVertexFinder::VertexFinder::computeAndSetVertexParameters	(	RecoVertex<> &	vertex,
		const std::vector< float > &	bin_centers,
		const std::vector< unsigned int > &	counts
	)

Vertexing algorithms.

Compute the vertex parameters

Definition at line 7 of file VertexFinder.cc.

References funct::abs(), dqmiodumpmetadata::counts, funct::pow(), DiDispStaMuonMonitor_cfi::pt, mathSSE::sqrt(), HLT_2024v13_cff::track, listHistos::trackPt, and bphysicsOniaDQM_cfi::vertex.

                                                                                          {
    double pt = 0.;
    double z0 = -999.;
    double z0width = 0.;
    bool highPt = false;
    double highestPt = 0.;
    unsigned int numHighPtTracks = 0;

    float SumZ = 0.;
    float z0square = 0.;
    float trackPt = 0.;

    std::vector<double> bin_pt(bin_centers.size(), 0.0);
    unsigned int ibin = 0;
    unsigned int itrack = 0;

    for (const L1Track* track : vertex.tracks()) {
      itrack++;
      trackPt = track->pt();

      // Skip the bins with no tracks
      while (ibin < counts.size() && counts[ibin] == 0)
        ibin++;

      if (trackPt > settings_->vx_TrackMaxPt()) {
        highPt = true;
        numHighPtTracks++;
        highestPt = (trackPt > highestPt) ? trackPt : highestPt;
        if (settings_->vx_TrackMaxPtBehavior() == 0)
          continue;  // ignore this track
        else if (settings_->vx_TrackMaxPtBehavior() == 1)
          trackPt = settings_->vx_TrackMaxPt();  // saturate
      }

      pt += std::pow(trackPt, settings_->vx_weightedmean());
      if (bin_centers.empty() && counts.empty()) {
        SumZ += track->z0() * std::pow(trackPt, settings_->vx_weightedmean());
        z0square += track->z0() * track->z0();
      } else {
        bin_pt[ibin] += std::pow(trackPt, settings_->vx_weightedmean());
        if (itrack == counts[ibin]) {
          SumZ += bin_centers[ibin] * bin_pt[ibin];
          z0square += bin_centers[ibin] * bin_centers[ibin];
          itrack = 0;
          ibin++;
        }
      }
    }

    z0 = SumZ / ((settings_->vx_weightedmean() > 0) ? pt : vertex.numTracks());
    z0square /= vertex.numTracks();
    z0width = sqrt(std::abs(z0 * z0 - z0square));

    vertex.setParameters(pt, z0, z0width, highPt, numHighPtTracks, highestPt);
  }

DBSCAN()

void l1tVertexFinder::VertexFinder::DBSCAN

(

)

DBSCAN algorithm.

Definition at line 182 of file VertexFinder.cc.

References funct::abs(), spr::find(), mps_fire::i, l1ctLayer2EG_cff::id, globals_cff::id2, dqmdumpme::k, DiDispStaMuonMonitor_cfi::pt, jetUpdater_cfi::sort, bphysicsOniaDQM_cfi::vertex, and class-composition::visited.

Referenced by VertexProducer::produce().

                            {
    // std::vector<RecoVertex> vClusters;
    std::vector<unsigned int> visited;
    std::vector<unsigned int> saved;

    sort(fitTracks_.begin(), fitTracks_.end(), SortTracksByPt());
    iterations_ = 0;

    for (unsigned int i = 0; i < fitTracks_.size(); ++i) {
      if (find(visited.begin(), visited.end(), i) != visited.end())
        continue;

      // if(fitTracks_[i]->pt()>10.){
      visited.push_back(i);
      std::set<unsigned int> neighbourTrackIds;
      unsigned int numDensityTracks = 0;
      if (fitTracks_[i].pt() > settings_->vx_dbscan_pt())
        numDensityTracks++;
      else
        continue;
      for (unsigned int k = 0; k < fitTracks_.size(); ++k) {
        iterations_++;
        if (k != i and std::abs(fitTracks_[k].z0() - fitTracks_[i].z0()) < settings_->vx_distance()) {
          neighbourTrackIds.insert(k);
          if (fitTracks_[k].pt() > settings_->vx_dbscan_pt()) {
            numDensityTracks++;
          }
        }
      }

      if (numDensityTracks < settings_->vx_dbscan_mintracks()) {
        // mark track as noise
      } else {
        RecoVertex vertex;
        vertex.insert(&fitTracks_[i]);
        saved.push_back(i);
        for (unsigned int id : neighbourTrackIds) {
          if (find(visited.begin(), visited.end(), id) == visited.end()) {
            visited.push_back(id);
            std::vector<unsigned int> neighbourTrackIds2;
            for (unsigned int k = 0; k < fitTracks_.size(); ++k) {
              iterations_++;
              if (std::abs(fitTracks_[k].z0() - fitTracks_[id].z0()) < settings_->vx_distance()) {
                neighbourTrackIds2.push_back(k);
              }
            }

            // if (neighbourTrackIds2.size() >= settings_->vx_minTracks()) {
            for (unsigned int id2 : neighbourTrackIds2) {
              neighbourTrackIds.insert(id2);
            }
            // }
          }
          if (find(saved.begin(), saved.end(), id) == saved.end())
            vertex.insert(&fitTracks_[id]);
        }
        computeAndSetVertexParameters(vertex, {}, {});
        if (vertex.numTracks() >= settings_->vx_minTracks())
          vertices_.push_back(vertex);
      }
      // }
    }
  }

fastHisto()

void l1tVertexFinder::VertexFinder::fastHisto

(

const TrackerTopology *

tTopo

)

Histogramming algorithm.

Definition at line 632 of file VertexFinder.cc.

References funct::abs(), reco::ceil(), dqmiodumpmetadata::counts, hcalRecHitTable_cff::detId, HLT_2024v13_cff::distance, spr::find(), newFWLiteAna::found, compareTotals::hist, mps_fire::i, createfilelist::int, dqmiolumiharvest::j, dqm-mbProfile::log, LaserClient_cfi::nbins, GetRecoTauVFromDQM_MC_cff::next, StripSubdetector::TID, TrackerTopology::tidRing(), createJobs::tmp, StripSubdetector::TOB, TrackerTopology::tobLayer(), HLT_2024v13_cff::track, align::Tracker, HcalDetIdTransform::transform(), pfDeepBoostedJetPreprocessParams_cfi::upper_bound, and findQualityFiles::v.

Referenced by VertexProducer::produce().

                                                           {
    // Create the histogram
    int nbins =
        std::ceil((settings_->vx_histogram_max() - settings_->vx_histogram_min()) / settings_->vx_histogram_binwidth());
    std::vector<RecoVertex<>> hist(nbins);
    std::vector<RecoVertex<>> sums(nbins - settings_->vx_windowSize() + 1);
    std::vector<float> bounds(nbins + 1);
    strided_iota(std::begin(bounds),
                 std::next(std::begin(bounds), nbins + 1),
                 settings_->vx_histogram_min(),
                 settings_->vx_histogram_binwidth());

    // Loop over the tracks and fill the histogram
    for (const L1Track& track : fitTracks_) {
      if ((track.z0() < settings_->vx_histogram_min()) || (track.z0() > settings_->vx_histogram_max()))
        continue;
      if (track.getTTTrackPtr()->chi2() > settings_->vx_TrackMaxChi2())
        continue;
      if (track.pt() < settings_->vx_TrackMinPt())
        continue;

      // Get the number of stubs and the number of stubs in PS layers
      float nPS = 0., nstubs = 0;

      // Get pointers to stubs associated to the L1 track
      const auto& theStubs = track.getTTTrackPtr()->getStubRefs();
      if (theStubs.empty()) {
        edm::LogWarning("VertexFinder") << "fastHisto::Could not retrieve the vector of stubs.";
        continue;
      }

      // Loop over the stubs
      for (const auto& stub : theStubs) {
        nstubs++;
        bool isPS = false;
        DetId detId(stub->getDetId());
        if (detId.det() == DetId::Detector::Tracker) {
          if (detId.subdetId() == StripSubdetector::TOB && tTopo->tobLayer(detId) <= 3)
            isPS = true;
          else if (detId.subdetId() == StripSubdetector::TID && tTopo->tidRing(detId) <= 9)
            isPS = true;
        }
        if (isPS)
          nPS++;
      }  // End loop over stubs
      if (nstubs < settings_->vx_NStubMin())
        continue;
      if (nPS < settings_->vx_NStubPSMin())
        continue;

      // Quality cuts, may need to be re-optimized
      int trk_nstub = (int)track.getTTTrackPtr()->getStubRefs().size();
      float chi2dof = track.getTTTrackPtr()->chi2() / (2 * trk_nstub - 4);

      if (settings_->vx_DoPtComp()) {
        float trk_consistency = track.getTTTrackPtr()->stubPtConsistency();
        if (trk_nstub == 4) {
          if (std::abs(track.eta()) < 2.2 && trk_consistency > 10)
            continue;
          else if (std::abs(track.eta()) > 2.2 && chi2dof > 5.0)
            continue;
        }
      }
      if (settings_->vx_DoTightChi2()) {
        if (track.pt() > 10.0 && chi2dof > 5.0)
          continue;
      }

      // Assign the track to the correct vertex
      // The values are ordered with bounds [lower, upper)
      // Values below bounds.begin() return 0 as the index (underflow)
      // Values above bounds.end() will return the index of the last bin (overflow)
      auto upper_bound = std::upper_bound(bounds.begin(), bounds.end(), track.z0());
      int index = std::distance(bounds.begin(), upper_bound) - 1;
      if (index == -1)
        index = 0;
      hist.at(index).insert(&track);
    }  // end loop over tracks

    // Compute the sums
    // sliding windows ... sum_i_i+(w-1) where i in (0,nbins-w) and w is the window size
    std::vector<float> bin_centers(settings_->vx_windowSize(), 0.0);
    std::vector<unsigned int> counts(settings_->vx_windowSize(), 0);
    for (unsigned int i = 0; i < sums.size(); i++) {
      for (unsigned int j = 0; j < settings_->vx_windowSize(); j++) {
        bin_centers[j] = settings_->vx_histogram_min() + ((i + j) * settings_->vx_histogram_binwidth()) +
                         (0.5 * settings_->vx_histogram_binwidth());
        counts[j] = hist.at(i + j).numTracks();
        sums.at(i) += hist.at(i + j);
      }
      computeAndSetVertexParameters(sums.at(i), bin_centers, counts);
    }

    // Find the maxima of the sums
    float sigma_max = -999;
    int imax = -999;
    std::vector<int> found;
    found.reserve(settings_->vx_nvtx());
    for (unsigned int ivtx = 0; ivtx < settings_->vx_nvtx(); ivtx++) {
      sigma_max = -999;
      imax = -999;
      for (unsigned int i = 0; i < sums.size(); i++) {
        // Skip this window if it will already be returned
        if (find(found.begin(), found.end(), i) != found.end())
          continue;
        if (sums.at(i).pt() > sigma_max) {
          sigma_max = sums.at(i).pt();
          imax = i;
        }
      }
      found.push_back(imax);
      vertices_.emplace_back(sums.at(imax));
    }
    pv_index_ = 0;

    if (settings_->debug() >= 1) {
      edm::LogInfo log("VertexProducer");
      log << "fastHisto::Checking the output parameters ... \n";
      std::vector<double> tmp;
      std::transform(std::begin(sums), std::end(sums), std::back_inserter(tmp), [](const RecoVertex<>& v) -> double {
        return v.pt();
      });
      printHistogram<double, edm::LogInfo>(log, tmp, 80, 0, -1, "fastHisto::sums", "\e[92m");
      for (unsigned int i = 0; i < found.size(); i++) {
        log << "RecoVertex " << i << ": bin index = " << found[i] << "\tsumPt = " << sums.at(imax).pt()
            << "\tz0 = " << sums.at(imax).z0();
      }
    }
  }  // end of fastHisto

fastHistoEmulation()

void l1tVertexFinder::VertexFinder::fastHistoEmulation

(

)

Histogramming algorithm (emulation)

Definition at line 762 of file VertexFinder.cc.

References cms::cuda::assert(), b, newFWLiteAna::bin, l1tVertexFinder::BitsToRepresent(), ALPAKA_ACCELERATOR_NAMESPACE::brokenline::constexpr(), filterCSVwithJSON::copy, spr::find(), newFWLiteAna::found, hcalSimParameters_cfi::hb, compareTotals::hist, mps_fire::i, cuy::ii, createfilelist::int, l1tmetemu::kPtMagSize, dqm-mbProfile::log, LaserClient_cfi::nbins, DiDispStaMuonMonitor_cfi::pt, HLT_2024v13_cff::track, validateGeometry_cfg::valid, and w().

Referenced by VertexProducer::produce().

                                        {
    // Relevant constants for the track word
    static constexpr int kZ0Size = 12,  // Width of z-position (40cm / 0.1)
        kPtSize = 14,                   // Width of pt
        kPtMagSize = 9,                 // Width of pt magnitude (unsigned)
        kReducedPrecisionPt = 7         // Width of the reduced precision, integer only, pt
        ;

    static constexpr int
        kBinSize = 8,                     // Width of a single bin in z
        kBinFixedSize = 8,                // Width of a single z0 bin in fixed point representation
        kBinFixedMagSize = 5,             // Width (magnitude) of a single z0 bin in fixed point representation
        kSlidingSumSize = 11,             // Width of the sum of a window of bins
        kInverseSize = 14,                // Width of the inverse sum
        kInverseMagSize = 1,              // Width of the inverse sum magnitude (unsigned)
        kWeightedSlidingSumSize = 20,     // Width of the pT weighted sliding sum
        kWeightedSlidingSumMagSize = 10,  // Width of the pT weighted sliding sum magnitude (signed)
        kWindowSize = 3,                  // Number of bins in the window used to sum histogram bins
        kSumPtLinkSize = 9,  // Number of bits used to represent the sum of track pts in a single bin from a single link

        kSumPtWindowBits = BitsToRepresent(kWindowSize * (1 << kSumPtLinkSize)),
        // Number of bits to represent the untruncated sum of track pts in a single bin from a single link
        kSumPtUntruncatedLinkSize = kPtSize + 2, kSumPtUntruncatedLinkMagSize = kPtMagSize + 2;

    static constexpr unsigned int kTableSize = ((1 << kSumPtLinkSize) - 1) * kWindowSize;

    typedef ap_ufixed<kPtSize, kPtMagSize, AP_RND_CONV, AP_SAT> pt_t;
    // Same size as TTTrack_TrackWord::z0_t, but now taking into account the sign bit (i.e. 2's complement)
    typedef ap_int<kZ0Size> z0_t;
    // 7 bits chosen to represent values between [0,127]
    // This is the next highest power of 2 value to our chosen track pt saturation value (100)
    typedef ap_ufixed<kReducedPrecisionPt, kReducedPrecisionPt, AP_RND_INF, AP_SAT> track_pt_fixed_t;
    // Histogram bin index
    typedef ap_uint<kBinSize> histbin_t;
    // Histogram bin in fixed point representation, before truncation
    typedef ap_ufixed<kBinFixedSize, kBinFixedMagSize, AP_RND_INF, AP_SAT> histbin_fixed_t;
    // This type is slightly arbitrary, but 2 bits larger than untruncated track pt to store sums in histogram bins
    // with truncation just before vertex-finding
    typedef ap_ufixed<kSumPtUntruncatedLinkSize, kSumPtUntruncatedLinkMagSize, AP_RND_INF, AP_SAT>
        histbin_pt_sum_fixed_t;
    // This value is slightly arbitrary, but small enough that the windows sums aren't too big.
    typedef ap_ufixed<kSumPtLinkSize, kSumPtLinkSize, AP_RND_INF, AP_SAT> link_pt_sum_fixed_t;
    // Enough bits to store kWindowSize * (2**kSumPtLinkSize)
    typedef ap_ufixed<kSumPtWindowBits, kSumPtWindowBits, AP_RND_INF, AP_SAT> window_pt_sum_fixed_t;
    // pt weighted sum of bins in window
    typedef ap_fixed<kWeightedSlidingSumSize, kWeightedSlidingSumMagSize, AP_RND_INF, AP_SAT> zsliding_t;
    // Sum of histogram bins in window
    typedef ap_uint<kSlidingSumSize> slidingsum_t;
    // Inverse of sum of bins in a given window
    typedef ap_ufixed<kInverseSize, kInverseMagSize, AP_RND_INF, AP_SAT> inverse_t;

    auto track_quality_check = [&](const track_pt_fixed_t& pt) -> bool {
      // Track quality cuts
      if (pt.to_double() < settings_->vx_TrackMinPt())
        return false;
      return true;
    };

    auto fetch_bin = [&](const z0_t& z0, int nbins) -> std::pair<histbin_t, bool> {
      // Increase the the number of bits in the word to allow for additional dynamic range
      ap_int<kZ0Size + 1> z0_13 = z0;
      // Add a number equal to half of the range in z0, meaning that the range is now [0, 2*z0_max]
      ap_int<kZ0Size + 1> absz0_13 = z0_13 + (1 << (kZ0Size - 1));
      // Shift the bits down to truncate the dynamic range to the most significant kBinFixedSize bits
      ap_int<kZ0Size + 1> absz0_13_reduced = absz0_13 >> (kZ0Size - kBinFixedSize);
      // Put the relevant bits into the histbin_t container
      histbin_t bin = absz0_13_reduced.range(kBinFixedSize - 1, 0);

      if (settings_->debug() > 2) {
        edm::LogInfo("VertexProducer")
            << "fastHistoEmulation::fetchBin() Checking the mapping from z0 to bin index ... \n"
            << "histbin_fixed_t(1.0 / settings_->vx_histogram_binwidth()) = "
            << histbin_fixed_t(1.0 / settings_->vx_histogram_binwidth()) << "\n"
            << "histbin_t(std::floor(nbins / 2) = " << histbin_t(std::floor(nbins / 2.)) << "\n"
            << "z0 = " << z0 << "\n"
            << "bin = " << bin;
      }
      bool valid = true;
      if (bin < 0) {
        return std::make_pair(0, false);
      } else if (bin > (nbins - 1)) {
        return std::make_pair(0, false);
      }
      return std::make_pair(bin, valid);
    };

    // Replace with https://stackoverflow.com/questions/13313980/populate-an-array-using-constexpr-at-compile-time ?
    auto init_inversion_table = [&]() -> std::vector<inverse_t> {
      std::vector<inverse_t> table_out(kTableSize, 0.);
      for (unsigned int ii = 0; ii < kTableSize; ii++) {
        // Compute lookup table function. This matches the format of the GTT HLS code.
        // Biased generation f(x) = 1 / (x + 1) is inverted by g(y) = inversion(x - 1) = 1 / (x - 1 + 1) = 1 / y
        table_out.at(ii) = (1.0 / (ii + 1));
      }
      return table_out;
    };

    auto inversion = [&](slidingsum_t& data_den) -> inverse_t {
      std::vector<inverse_t> inversion_table = init_inversion_table();

      // Index into the lookup table based on data
      int index;
      if (data_den < 0)
        data_den = 0;
      if (data_den > (kTableSize - 1))
        data_den = kTableSize - 1;
      index = data_den;
      return inversion_table.at(index);
    };

    auto bin_center = [&](zsliding_t iz, int nbins) -> l1t::VertexWord::vtxz0_t {
      zsliding_t z = iz - histbin_t(std::floor(nbins / 2.));
      std::unique_ptr<edm::LogInfo> log;
      if (settings_->debug() >= 1) {
        log = std::make_unique<edm::LogInfo>("VertexProducer");
        *log << "bin_center information ...\n"
             << "iz = " << iz << "\n"
             << "histbin_t(std::floor(nbins / 2.)) = " << histbin_t(std::floor(nbins / 2.)) << "\n"
             << "binwidth = " << zsliding_t(settings_->vx_histogram_binwidth()) << "\n"
             << "z = " << z << "\n"
             << "zsliding_t(z * zsliding_t(binwidth)) = " << std::setprecision(7)
             << l1t::VertexWord::vtxz0_t(z * zsliding_t(settings_->vx_histogram_binwidth()));
      }
      return l1t::VertexWord::vtxz0_t(z * zsliding_t(settings_->vx_histogram_binwidth()));
    };

    auto weighted_position = [&](histbin_t b_max,
                                 const std::vector<link_pt_sum_fixed_t>& binpt,
                                 slidingsum_t maximums,
                                 int nbins) -> zsliding_t {
      zsliding_t zvtx_sliding = 0;
      slidingsum_t zvtx_sliding_sum = 0;
      inverse_t inv = 0;

      std::unique_ptr<edm::LogInfo> log;
      if (settings_->debug() >= 1) {
        log = std::make_unique<edm::LogInfo>("VertexProducer");
        *log << "Progression of weighted_position() ...\n"
             << "zvtx_sliding_sum = ";
      }

      // Find the weighted position within the window in index space (width = 1)
      for (ap_uint<BitsToRepresent(kWindowSize)> w = 0; w < kWindowSize; ++w) {
        zvtx_sliding_sum += (binpt.at(w) * w);
        if (settings_->debug() >= 1) {
          *log << "(" << w << " * " << binpt.at(w) << ")";
          if (w < kWindowSize - 1) {
            *log << " + ";
          }
        }
      }

      if (settings_->debug() >= 1) {
        *log << " = " << zvtx_sliding_sum << "\n";
      }

      if (maximums != 0) {
        //match F/W inversion_lut offset (inversion[x] = 1 / (x + 1); inversion[x - 1] = 1 / x;), for consistency
        slidingsum_t offsetmaximums = maximums - 1;
        inv = inversion(offsetmaximums);
        zvtx_sliding = zvtx_sliding_sum * inv;
      } else {
        zvtx_sliding = (settings_->vx_windowSize() / 2.0) + (((int(settings_->vx_windowSize()) % 2) != 0) ? 0.5 : 0.0);
      }
      if (settings_->debug() >= 1) {
        *log << "inversion(" << maximums << ") = " << inv << "\nzvtx_sliding = " << zvtx_sliding << "\n";
      }

      // 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
      zvtx_sliding += b_max;
      zvtx_sliding += ap_ufixed<1, 0>(0.5);
      if (settings_->debug() >= 1) {
        *log << "b_max = " << b_max << "\n";
        *log << "zvtx_sliding + b_max + 0.5 = " << zvtx_sliding << "\n";
      }

      // Shift the z position from index space into z [cm] space
      zvtx_sliding = bin_center(zvtx_sliding, nbins);
      if (settings_->debug() >= 1) {
        *log << "bin_center(zvtx_sliding + b_max + 0.5, nbins) = " << std::setprecision(7) << zvtx_sliding;
        log.reset();
      }
      return zvtx_sliding;
    };

    // Create the histogram
    unsigned int nbins = std::round((settings_->vx_histogram_max() - settings_->vx_histogram_min()) /
                                    settings_->vx_histogram_binwidth());
    unsigned int nsums = nbins - settings_->vx_windowSize() + 1;
    std::vector<link_pt_sum_fixed_t> hist(nbins, 0);
    std::vector<histbin_pt_sum_fixed_t> hist_untruncated(nbins, 0);

    // Loop over the tracks and fill the histogram
    if (settings_->debug() > 2) {
      edm::LogInfo("VertexProducer") << "fastHistoEmulation::Processing " << fitTracks_.size() << " tracks";
    }
    for (const L1Track& track : fitTracks_) {
      // Get the track pt and z0
      // Convert them to an appropriate data format
      // Truncation and saturation taken care of by the data type specification, now delayed to end of histogramming
      pt_t tkpt = 0;
      tkpt.V = track.getTTTrackPtr()->getTrackWord()(TTTrack_TrackWord::TrackBitLocations::kRinvMSB - 1,
                                                     TTTrack_TrackWord::TrackBitLocations::kRinvLSB);
      z0_t tkZ0 = track.getTTTrackPtr()->getZ0Word();

      if ((settings_->vx_DoQualityCuts() && track_quality_check(tkpt)) || (!settings_->vx_DoQualityCuts())) {
        //
        // Check bin validity of bin found for the current track
        //
        std::pair<histbin_t, bool> bin = fetch_bin(tkZ0, nbins);
        assert(bin.first >= 0 && bin.first < nbins);

        //
        // If the bin is valid then sum the tracks
        //
        if (settings_->debug() > 2) {
          edm::LogInfo("VertexProducer") << "fastHistoEmulation::Checking the track word ... \n"
                                         << "track word = " << track.getTTTrackPtr()->getTrackWord().to_string(2)
                                         << "\n"
                                         << "tkZ0 = " << tkZ0.to_double() << "(" << tkZ0.to_string(2)
                                         << ")\ttkpt = " << tkpt.to_double() << "(" << tkpt.to_string(2)
                                         << ")\tbin = " << bin.first.to_int() << "\n"
                                         << "pt sum in bin " << bin.first.to_int()
                                         << " BEFORE adding track = " << hist_untruncated.at(bin.first).to_double();
        }
        if (bin.second) {
          hist_untruncated.at(bin.first) = hist_untruncated.at(bin.first) + tkpt;
        }
        if (settings_->debug() > 2) {
          edm::LogInfo("VertexProducer") << "fastHistoEmulation::\npt sum in bin " << bin.first.to_int()
                                         << " AFTER adding track = " << hist_untruncated.at(bin.first).to_double();
        }
      } else {
        if (settings_->debug() > 2) {
          edm::LogInfo("VertexProducer") << "fastHistoEmulation::Did not add the following track ... \n"
                                         << "track word = " << track.getTTTrackPtr()->getTrackWord().to_string(2)
                                         << "\n"
                                         << "tkZ0 = " << tkZ0.to_double() << "(" << tkZ0.to_string(2)
                                         << ")\ttkpt = " << tkpt.to_double() << "(" << tkpt.to_string(2) << ")";
        }
      }
    }  // end loop over tracks

    // HLS histogramming used to truncate track pt before adding, using
    // track_pt_fixed_t pt_tmp = tkpt;
    // Now, truncation should happen after histograms are filled but prior to the vertex-finding part of the algo
    for (unsigned int hb = 0; hb < hist.size(); ++hb) {
      link_pt_sum_fixed_t bin_trunc = hist_untruncated.at(hb).range(
          kSumPtUntruncatedLinkSize - 1, kSumPtUntruncatedLinkSize - kSumPtUntruncatedLinkMagSize);
      hist.at(hb) = bin_trunc;
      if (settings_->debug() > 2) {
        edm::LogInfo("VertexProducer") << "fastHistoEmulation::truncating histogram bin pt once filling is complete \n"
                                       << "hist_untruncated.at(" << hb << ") = " << hist_untruncated.at(hb).to_double()
                                       << "(" << hist_untruncated.at(hb).to_string(2)
                                       << ")\tbin_trunc = " << bin_trunc.to_double() << "(" << bin_trunc.to_string(2)
                                       << ")\n\thist.at(" << hb << ") = " << hist.at(hb).to_double() << "("
                                       << hist.at(hb).to_string(2) << ")";
      }
    }

    // Loop through all bins, taking into account the fact that the last bin is nbins-window_width+1,
    // and compute the sums using sliding windows ... sum_i_i+(w-1) where i in (0,nbins-w) and w is the window size
    std::vector<window_pt_sum_fixed_t> hist_window_sums(nsums, 0);
    for (unsigned int b = 0; b < nsums; ++b) {
      for (unsigned int w = 0; w < kWindowSize; ++w) {
        unsigned int index = b + w;
        hist_window_sums.at(b) += hist.at(index);
      }
    }

    // Find the top N vertices
    std::vector<int> found;
    found.reserve(settings_->vx_nvtx());
    for (unsigned int ivtx = 0; ivtx < settings_->vx_nvtx(); ivtx++) {
      histbin_t b_max = 0;
      window_pt_sum_fixed_t max_pt = 0;
      zsliding_t zvtx_sliding = -999;
      std::vector<link_pt_sum_fixed_t> binpt_max(kWindowSize, 0);

      // Find the maxima of the sums
      for (unsigned int i = 0; i < hist_window_sums.size(); i++) {
        // Skip this window if it will already be returned
        if (find(found.begin(), found.end(), i) != found.end())
          continue;
        if (hist_window_sums.at(i) > max_pt) {
          b_max = i;
          max_pt = hist_window_sums.at(b_max);
          std::copy(std::begin(hist) + b_max, std::begin(hist) + b_max + kWindowSize, std::begin(binpt_max));

          // Find the weighted position only for the highest sum pt window
          zvtx_sliding = weighted_position(b_max, binpt_max, max_pt, nbins);
        }
      }
      if (settings_->debug() >= 1) {
        edm::LogInfo log("VertexProducer");
        log << "fastHistoEmulation::Checking the output parameters ... \n";
        if (found.empty()) {
          printHistogram<link_pt_sum_fixed_t, edm::LogInfo>(log, hist, 80, 0, -1, "fastHistoEmulation::hist", "\e[92m");
          printHistogram<window_pt_sum_fixed_t, edm::LogInfo>(
              log, hist_window_sums, 80, 0, -1, "fastHistoEmulation::hist_window_sums", "\e[92m");
        }
        printHistogram<link_pt_sum_fixed_t, edm::LogInfo>(
            log, binpt_max, 80, 0, -1, "fastHistoEmulation::binpt_max", "\e[92m");
        log << "bin index (not a VertexWord parameter) = " << b_max << "\n"
            << "sumPt = " << max_pt.to_double() << "\n"
            << "z0 = " << zvtx_sliding.to_double();
      }
      found.push_back(b_max);
      verticesEmulation_.emplace_back(l1t::VertexWord::vtxvalid_t(1),
                                      l1t::VertexWord::vtxz0_t(zvtx_sliding),
                                      l1t::VertexWord::vtxmultiplicity_t(0),
                                      l1t::VertexWord::vtxsumpt_t(max_pt),
                                      l1t::VertexWord::vtxquality_t(0),
                                      l1t::VertexWord::vtxinversemult_t(0),
                                      l1t::VertexWord::vtxunassigned_t(0));
    }
    pv_index_ = 0;
  }  // end of fastHistoEmulation

fastHistoLooseAssociation()

void l1tVertexFinder::VertexFinder::fastHistoLooseAssociation

(

)

High pT Vertex Algorithm.

Definition at line 608 of file VertexFinder.cc.

References funct::abs(), HLT_2024v13_cff::track, and bphysicsOniaDQM_cfi::vertex.

Referenced by VertexProducer::produce().

                                               {
    float vxPt = 0.;
    RecoVertex leading_vertex;

    for (float z = settings_->vx_histogram_min(); z < settings_->vx_histogram_max();
         z += settings_->vx_histogram_binwidth()) {
      RecoVertex vertex;
      for (const L1Track& track : fitTracks_) {
        if (std::abs(z - track.z0()) < settings_->vx_width()) {
          vertex.insert(&track);
        }
      }
      computeAndSetVertexParameters(vertex, {}, {});
      vertex.setZ0(z);
      if (vertex.pt() > vxPt) {
        leading_vertex = vertex;
        vxPt = vertex.pt();
      }
    }

    vertices_.emplace_back(leading_vertex);
    pv_index_ = 0;  // by default fastHistoLooseAssociation algorithm finds only hard PV
  }                 // end of fastHistoLooseAssociation

findPrimaryVertex()

void l1tVertexFinder::VertexFinder::findPrimaryVertex

(

)

Find the primary vertex.

Definition at line 445 of file VertexFinder.cc.

References HLT_2024v13_cff::distance, l1tVertexFinder::RecoVertex< T >::pt(), and l1t::VertexWord::pt().

Referenced by VertexProducer::produce().

                                       {
    if (settings_->vx_precision() == Precision::Emulation) {
      pv_index_ = std::distance(verticesEmulation_.begin(),
                                std::max_element(verticesEmulation_.begin(),
                                                 verticesEmulation_.end(),
                                                 [](const l1t::VertexWord& vertex0, const l1t::VertexWord& vertex1) {
                                                   return (vertex0.pt() < vertex1.pt());
                                                 }));
    } else {
      pv_index_ = std::distance(
          vertices_.begin(),
          std::max_element(
              vertices_.begin(), vertices_.end(), [](const RecoVertex<>& vertex0, const RecoVertex<>& vertex1) {
                return (vertex0.pt() < vertex1.pt());
              }));
    }
  }

fitTracks()

const FitTrackCollection& l1tVertexFinder::VertexFinder::fitTracks ( ) const

inline

Accessors.

Storage for tracks out of the L1 Track finder

Definition at line 50 of file VertexFinder.h.

References fitTracks_.

Referenced by VertexFinder().

{ return fitTracks_; }

GapClustering()

void l1tVertexFinder::VertexFinder::GapClustering

(

)

Gap Clustering Algorithm.

Definition at line 65 of file VertexFinder.cc.

References mps_fire::i, l1tVertexFinder::Vertex::insert(), l1tVertexFinder::Vertex::numTracks(), or, jetUpdater_cfi::sort, and HltBtagValidation_cff::Vertex.

Referenced by VertexProducer::produce().

                                   {
    sort(fitTracks_.begin(), fitTracks_.end(), SortTracksByZ0());
    iterations_ = 0;
    RecoVertex Vertex;
    for (unsigned int i = 0; i < fitTracks_.size(); ++i) {
      Vertex.insert(&fitTracks_[i]);
      iterations_++;
      if ((i + 1 < fitTracks_.size() and fitTracks_[i + 1].z0() - fitTracks_[i].z0() > settings_->vx_distance()) or
          i == fitTracks_.size() - 1) {
        if (Vertex.numTracks() >= settings_->vx_minTracks()) {
          computeAndSetVertexParameters(Vertex, {}, {});
          vertices_.push_back(Vertex);
        }
        Vertex.clear();
      }
    }
  }

HPV()

void l1tVertexFinder::VertexFinder::HPV

(

)

High pT Vertex Algorithm.

Definition at line 376 of file VertexFinder.cc.

References funct::abs(), dqmdumpme::first, jetUpdater_cfi::sort, HLT_2024v13_cff::track, and bphysicsOniaDQM_cfi::vertex.

Referenced by VertexProducer::produce().

                         {
    iterations_ = 0;
    sort(fitTracks_.begin(), fitTracks_.end(), SortTracksByPt());

    RecoVertex vertex;
    bool first = true;
    float z = 99.;
    for (const L1Track& track : fitTracks_) {
      if (track.pt() < 50.) {
        if (first) {
          first = false;
          z = track.z0();
          vertex.insert(&track);
        } else {
          if (std::abs(track.z0() - z) < settings_->vx_distance())
            vertex.insert(&track);
        }
      }
    }

    computeAndSetVertexParameters(vertex, {}, {});
    vertex.setZ0(z);
    vertices_.push_back(vertex);
  }

iterationsPerTrack()

unsigned int l1tVertexFinder::VertexFinder::iterationsPerTrack ( ) const

inline

Number of iterations.

Definition at line 52 of file VertexFinder.h.

References fitTracks_, and iterations_.

{ return double(iterations_) / double(fitTracks_.size()); }

Kmeans()

void l1tVertexFinder::VertexFinder::Kmeans

(

)

Kmeans Algorithm.

Definition at line 401 of file VertexFinder.cc.

References funct::abs(), HLT_2024v13_cff::distance, mps_fire::i, l1ctLayer2EG_cff::id, NA, and HLT_2024v13_cff::track.

Referenced by VertexProducer::produce().

                            {
    unsigned int NumberOfClusters = settings_->vx_kmeans_nclusters();

    vertices_.resize(NumberOfClusters);
    float ClusterSeparation = 30. / NumberOfClusters;

    for (unsigned int i = 0; i < NumberOfClusters; ++i) {
      float ClusterCentre = -15. + ClusterSeparation * (i + 0.5);
      vertices_[i].setZ0(ClusterCentre);
    }
    unsigned int iterations = 0;
    // Initialise Clusters
    while (iterations < settings_->vx_kmeans_iterations()) {
      for (unsigned int i = 0; i < NumberOfClusters; ++i) {
        vertices_[i].clear();
      }

      for (const L1Track& track : fitTracks_) {
        float distance = 9999;
        if (iterations == settings_->vx_kmeans_iterations() - 3)
          distance = settings_->vx_distance() * 2;
        if (iterations > settings_->vx_kmeans_iterations() - 3)
          distance = settings_->vx_distance();
        unsigned int ClusterId;
        bool NA = true;
        for (unsigned int id = 0; id < NumberOfClusters; ++id) {
          if (std::abs(track.z0() - vertices_[id].z0()) < distance) {
            distance = std::abs(track.z0() - vertices_[id].z0());
            ClusterId = id;
            NA = false;
          }
        }
        if (!NA) {
          vertices_[ClusterId].insert(&track);
        }
      }
      for (unsigned int i = 0; i < NumberOfClusters; ++i) {
        if (vertices_[i].numTracks() >= settings_->vx_minTracks())
          computeAndSetVertexParameters(vertices_[i], {}, {});
      }
      iterations++;
    }
  }

maxDistance()

float l1tVertexFinder::VertexFinder::maxDistance	(	RecoVertex<>	cluster0,
		RecoVertex<>	cluster1
	)

Find maximum distance in two clusters of tracks.

Definition at line 83 of file VertexFinder.cc.

References funct::abs(), HLT_2024v13_cff::distance, and l1tVertexFinder::RecoVertex< T >::tracks().

                                                                              {
    float distance = 0;
    for (const L1Track* track0 : cluster0.tracks()) {
      for (const L1Track* track1 : cluster1.tracks()) {
        if (std::abs(track0->z0() - track1->z0()) > distance) {
          distance = std::abs(track0->z0() - track1->z0());
        }
      }
    }

    return distance;
  }

meanDistance()

float l1tVertexFinder::VertexFinder::meanDistance	(	RecoVertex<>	cluster0,
		RecoVertex<>	cluster1
	)

Find average distance in two clusters of tracks.

Definition at line 109 of file VertexFinder.cc.

References funct::abs(), HLT_2024v13_cff::distance, l1tVertexFinder::RecoVertex< T >::numTracks(), and l1tVertexFinder::RecoVertex< T >::tracks().

                                                                               {
    float distanceSum = 0;

    for (const L1Track* track0 : cluster0.tracks()) {
      for (const L1Track* track1 : cluster1.tracks()) {
        distanceSum += std::abs(track0->z0() - track1->z0());
      }
    }

    float distance = distanceSum / (cluster0.numTracks() * cluster1.numTracks());
    return distance;
  }

minDistance()

float l1tVertexFinder::VertexFinder::minDistance	(	RecoVertex<>	cluster0,
		RecoVertex<>	cluster1
	)

Find minimum distance in two clusters of tracks.

Definition at line 96 of file VertexFinder.cc.

References funct::abs(), HLT_2024v13_cff::distance, and l1tVertexFinder::RecoVertex< T >::tracks().

                                                                              {
    float distance = 9999;
    for (const L1Track* track0 : cluster0.tracks()) {
      for (const L1Track* track1 : cluster1.tracks()) {
        if (std::abs(track0->z0() - track1->z0()) < distance) {
          distance = std::abs(track0->z0() - track1->z0());
        }
      }
    }

    return distance;
  }

NNVtxEmulation()

void l1tVertexFinder::VertexFinder::NNVtxEmulation	(	tensorflow::Session *	TrackWeightSesh = nullptr,
		tensorflow::Session *	PatternRecSesh = nullptr,
		tensorflow::Session *	AssociationSesh = nullptr
	)

NNVtx algorithm.

Definition at line 1081 of file VertexFinder.cc.

References funct::abs(), reco::ceil(), counter, nano_mu_digi_cff::float, mps_fire::i, hgcal_clustering::max_index(), tensorflow::run(), HLT_2024v13_cff::track, parallelization::uint, and AlignmentTracksFromVertexSelector_cfi::vertices.

Referenced by VertexProducer::produce().

                                                                        {
    // #### Weight Tracks: ####
    // Loop over tracks -> weight the network -> set track weights
    tensorflow::Tensor inputTrkWeight(tensorflow::DT_FLOAT, {1, 3});  //Single batch of 3 values
    uint counter = 0;

    for (auto& track : fitTracks_) {
      // Quantised Network: Use values from L1GTTInputProducer pT, MVA1, eta
      auto& gttTrack = fitTracks_.at(counter);

      TTTrack_TrackWord::tanl_t etaEmulationBits = gttTrack.getTTTrackPtr()->getTanlWord();
      ap_fixed<16, 3> etaEmulation;
      etaEmulation.V = (etaEmulationBits.range());

      ap_uint<14> ptEmulationBits = gttTrack.getTTTrackPtr()->getTrackWord()(
          TTTrack_TrackWord::TrackBitLocations::kRinvMSB - 1, TTTrack_TrackWord::TrackBitLocations::kRinvLSB);
      ap_ufixed<14, 9> ptEmulation;
      ptEmulation.V = (ptEmulationBits.range());

      ap_ufixed<22, 9> ptEmulation_rescale;
      ptEmulation_rescale = ptEmulation.to_double();

      ap_ufixed<22, 9> etaEmulation_rescale;
      etaEmulation_rescale = abs(etaEmulation.to_double());

      ap_ufixed<22, 9> MVAEmulation_rescale;
      MVAEmulation_rescale = gttTrack.getTTTrackPtr()->getMVAQualityBits();

      inputTrkWeight.tensor<float, 2>()(0, 0) = ptEmulation_rescale.to_double();
      inputTrkWeight.tensor<float, 2>()(0, 1) = MVAEmulation_rescale.to_double();
      inputTrkWeight.tensor<float, 2>()(0, 2) = etaEmulation_rescale.to_double();

      // CNN output: track weight
      std::vector<tensorflow::Tensor> outputTrkWeight;
      tensorflow::run(TrackWeightSesh, {{"weight:0", inputTrkWeight}}, {"Identity:0"}, &outputTrkWeight);
      // Set track weight pack into tracks:

      ap_ufixed<16, 5> NNOutput;
      NNOutput = (double)outputTrkWeight[0].tensor<float, 2>()(0, 0);

      //std::cout<<"NNOutput_weight_network = "<< NNOutput <<std::endl;

      track.setWeight(NNOutput.to_double());

      ++counter;
    }
    // #### Find Vertices: ####
    tensorflow::Tensor inputPV(tensorflow::DT_FLOAT,
                               {1, settings_->vx_histogram_numbins(), 1});  //Single batch with 256 bins and 1 weight
    std::vector<tensorflow::Tensor> outputPV;
    RecoVertexCollection vertices(settings_->vx_histogram_numbins());
    std::map<float, int> vertexMap;
    std::map<int, float> histogram;
    std::map<int, float> nnOutput;

    float binWidth = settings_->vx_histogram_binwidth();

    int track_z = 0;

    for (int z = 0; z < settings_->vx_histogram_numbins(); z += 1) {
      counter = 0;
      double vxWeight = 0;

      for (const L1Track& track : fitTracks_) {
        auto& gttTrack = fitTracks_.at(counter);
        double temp_z0 = gttTrack.getTTTrackPtr()->z0();

        track_z = std::floor((temp_z0 + settings_->vx_histogram_max()) / binWidth);

        if (track_z >= z && track_z < (z + 1)) {
          vertices.at(z).insert(&track);
          vxWeight += track.weight();
        }
        ++counter;
      }
      // Get centre of bin before setting z0
      vertices.at(z).setZ0(((z + 0.5) * binWidth) - settings_->vx_histogram_max());

      vertexMap[vxWeight] = z;
      inputPV.tensor<float, 3>()(0, z, 0) = vxWeight;
      //Fill histogram for 3 bin sliding window:
      histogram[z] = vxWeight;
    }

    // Run PV Network:
    tensorflow::run(PatternRecSesh, {{"hist:0", inputPV}}, {"Identity:0"}, &outputPV);
    // Threshold needed due to rounding differences in internal CNN layer emulation versus firmware
    const float histogrammingThreshold_ = 0.0;
    for (int i(0); i < settings_->vx_histogram_numbins(); ++i) {
      if (outputPV[0].tensor<float, 3>()(0, i, 0) >= histogrammingThreshold_) {
        nnOutput[i] = outputPV[0].tensor<float, 3>()(0, i, 0);
      }
    }

    //Find max then find all occurances of it in histogram and average their position -> python argmax layer
    //Performance is not optimised for multiple peaks in histogram or spread peaks these are edge cases, need to revisit
    int max_index = 0;
    int num_maxes = 0;
    float max_element = 0.0;
    for (int i(0); i < settings_->vx_histogram_numbins(); ++i) {
      if (nnOutput[i] > max_element) {
        max_element = nnOutput[i];
      }
    }

    for (int i(0); i < settings_->vx_histogram_numbins(); ++i) {
      if (nnOutput[i] == max_element) {
        num_maxes++;
        max_index += i;
      }
    }
    int PV_index = ceil((float)max_index / (float)num_maxes);

    edm::LogVerbatim("VertexFinder") << " NNVtxEmulation Chosen PV: prob: " << nnOutput[PV_index]
                                     << " bin = " << PV_index << " z0 = " << vertices.at(PV_index).z0() << '\n';

    verticesEmulation_.emplace_back(1, vertices.at(PV_index).z0(), 0, vertices.at(PV_index).pt(), 0, 0, 0);

  }  // end of NNVtx Algorithm

numInputTracks()

unsigned int l1tVertexFinder::VertexFinder::numInputTracks ( ) const

inline

Storage for tracks out of the L1 Track finder.

Definition at line 54 of file VertexFinder.h.

References fitTracks_.

{ return fitTracks_.size(); }

numIterations()

unsigned int l1tVertexFinder::VertexFinder::numIterations ( ) const

inline

Number of iterations.

Definition at line 56 of file VertexFinder.h.

References iterations_.

{ return iterations_; }

numVertices()

unsigned int l1tVertexFinder::VertexFinder::numVertices ( ) const

inline

Number of reconstructed vertices.

Definition at line 58 of file VertexFinder.h.

References vertices_.

{ return vertices_.size(); }

numVerticesEmulation()

unsigned int l1tVertexFinder::VertexFinder::numVerticesEmulation ( ) const

inline

Number of emulation vertices.

Definition at line 60 of file VertexFinder.h.

References verticesEmulation_.

{ return verticesEmulation_.size(); }

PrimaryVertex()

template<typename T >

T l1tVertexFinder::VertexFinder::PrimaryVertex ( ) const

inline

Reconstructed primary vertex.

Definition at line 63 of file VertexFinder.h.

References l1tVertexFinder::Emulation, pv_index_, settings_, l1tVertexFinder::Simulation, vertices_, verticesEmulation_, and l1tVertexFinder::AlgoSettings::vx_precision().

                            {
      if ((settings_->vx_precision() == Precision::Simulation) && (pv_index_ < vertices_.size()))
        return vertices_[pv_index_];
      else if ((settings_->vx_precision() == Precision::Emulation) && (pv_index_ < vertices_.size()))
        return verticesEmulation_[pv_index_];
      else {
        edm::LogWarning("VertexFinder") << "PrimaryVertex::No primary vertex has been found.";
        return RecoVertex<>();
      }
    }

primaryVertexId()

unsigned int l1tVertexFinder::VertexFinder::primaryVertexId ( ) const

inline

Reconstructed Primary Vertex Id.

Definition at line 74 of file VertexFinder.h.

References pv_index_.

{ return pv_index_; }

printHistogram()

template<class data_type , typename stream_type >

void l1tVertexFinder::VertexFinder::printHistogram	(	stream_type &	stream,
		std::vector< data_type >	data,
		int	width = 80,
		int	minimum = 0,
		int	maximum = -1,
		std::string	title = "",
		std::string	color = ""
	)

Print an ASCII histogram.

Definition at line 465 of file VertexFinder.cc.

References edmScanValgrind::buffer, data, nano_mu_digi_cff::float, mps_fire::i, writeEcalDQMStatus::interval, mps_splice::line, SiStripPI::max, EcalTangentSkim_cfg::o, MillePedeFileConverter_cfg::out, l1tEGammaCrystalsEmulatorProducer_cfi::scale, cms::cuda::stream, AlCaHLTBitMon_QueryRunRegistry::string, runGCPTkAlMap::title, contentValuesCheck::values, ApeEstimator_cff::width, x, and SiPixelPI::zero.

                                                     {
    int tableSize = data.size();

    if (maximum == -1) {
      maximum = float(*std::max_element(std::begin(data), std::end(data))) * 1.05;
    } else if (maximum <= minimum) {
      maximum = float(*std::max_element(std::begin(data), std::end(data))) * 1.05;
      minimum = float(*std::min_element(std::begin(data), std::end(data)));
    }

    if (minimum < 0) {
      minimum *= 1.05;
    } else {
      minimum = 0;
    }

    std::vector<std::string> intervals(tableSize, "");
    std::vector<std::string> values(tableSize, "");
    char buffer[128];
    int intervalswidth = 0, valueswidth = 0, tmpwidth = 0;
    for (int i = 0; i < tableSize; i++) {
      //Format the bin labels
      tmpwidth = sprintf(buffer, "[%-.5g, %-.5g)", float(i), float(i + 1));
      intervals[i] = buffer;
      if (i == (tableSize - 1)) {
        intervals[i][intervals[i].size() - 1] = ']';
      }
      if (tmpwidth > intervalswidth)
        intervalswidth = tmpwidth;

      //Format the values
      tmpwidth = sprintf(buffer, "%-.5g", float(data[i]));
      values[i] = buffer;
      if (tmpwidth > valueswidth)
        valueswidth = tmpwidth;
    }

    sprintf(buffer, "%-.5g", float(minimum));
    std::string minimumtext = buffer;
    sprintf(buffer, "%-.5g", float(maximum));
    std::string maximumtext = buffer;

    int plotwidth =
        std::max(int(minimumtext.size() + maximumtext.size()), width - (intervalswidth + 1 + valueswidth + 1 + 2));
    std::string scale =
        minimumtext + std::string(plotwidth + 2 - minimumtext.size() - maximumtext.size(), ' ') + maximumtext;

    float norm = float(plotwidth) / float(maximum - minimum);
    int zero = std::round((0.0 - minimum) * norm);
    std::vector<char> line(plotwidth, '-');

    if ((minimum != 0) && (0 <= zero) && (zero < plotwidth)) {
      line[zero] = '+';
    }
    std::string capstone =
        std::string(intervalswidth + 1 + valueswidth + 1, ' ') + "+" + std::string(line.begin(), line.end()) + "+";

    std::vector<std::string> out;
    if (!title.empty()) {
      out.push_back(title);
      out.push_back(std::string(title.size(), '='));
    }
    out.push_back(std::string(intervalswidth + valueswidth + 2, ' ') + scale);
    out.push_back(capstone);
    for (int i = 0; i < tableSize; i++) {
      std::string interval = intervals[i];
      std::string value = values[i];
      data_type x = data[i];
      std::fill_n(line.begin(), plotwidth, ' ');

      int pos = std::round((float(x) - minimum) * norm);
      if (x < 0) {
        std::fill_n(line.begin() + pos, zero - pos, '*');
      } else {
        std::fill_n(line.begin() + zero, pos - zero, '*');
      }

      if ((minimum != 0) && (0 <= zero) && (zero < plotwidth)) {
        line[zero] = '|';
      }

      sprintf(buffer,
              "%-*s %-*s |%s|",
              intervalswidth,
              interval.c_str(),
              valueswidth,
              value.c_str(),
              std::string(line.begin(), line.end()).c_str());
      out.push_back(buffer);
    }
    out.push_back(capstone);
    if (!color.empty())
      stream << color;
    for (const auto& o : out) {
      stream << o << "\n";
    }
    if (!color.empty())
      stream << "\e[0m";
    stream << "\n";
  }

PVR()

void l1tVertexFinder::VertexFinder::PVR

(

)

Find maximum distance in two clusters of tracks.

Definition at line 246 of file VertexFinder.cc.

References funct::abs(), ALPAKA_ACCELERATOR_NAMESPACE::vertexFinder::it, or, HLT_2024v13_cff::track, bphysicsOniaDQM_cfi::vertex, and l1tVertexFinder::L1Track::z0().

Referenced by VertexProducer::produce().

                         {
    bool start = true;
    FitTrackCollection discardedTracks, acceptedTracks;
    iterations_ = 0;
    for (const L1Track& track : fitTracks_) {
      acceptedTracks.push_back(track);
    }

    while (discardedTracks.size() >= settings_->vx_minTracks() or start == true) {
      start = false;
      bool removing = true;
      discardedTracks.clear();
      while (removing) {
        float oldDistance = 0.;

        if (settings_->debug() > 2)
          edm::LogInfo("VertexFinder") << "PVR::AcceptedTracks " << acceptedTracks.size();

        float z0start = 0;
        for (const L1Track& track : acceptedTracks) {
          z0start += track.z0();
          iterations_++;
        }

        z0start /= acceptedTracks.size();
        if (settings_->debug() > 2)
          edm::LogInfo("VertexFinder") << "PVR::z0 vertex " << z0start;
        FitTrackCollection::iterator badTrackIt = acceptedTracks.end();
        removing = false;

        for (FitTrackCollection::iterator it = acceptedTracks.begin(); it < acceptedTracks.end(); ++it) {
          const L1Track* track = &*it;
          iterations_++;
          if (std::abs(track->z0() - z0start) > settings_->vx_distance() and
              std::abs(track->z0() - z0start) > oldDistance) {
            badTrackIt = it;
            oldDistance = std::abs(track->z0() - z0start);
            removing = true;
          }
        }

        if (removing) {
          const L1Track badTrack = *badTrackIt;
          if (settings_->debug() > 2)
            edm::LogInfo("VertexFinder") << "PVR::Removing track " << badTrack.z0() << " at distance " << oldDistance;
          discardedTracks.push_back(badTrack);
          acceptedTracks.erase(badTrackIt);
        }
      }

      if (acceptedTracks.size() >= settings_->vx_minTracks()) {
        RecoVertex vertex;
        for (const L1Track& track : acceptedTracks) {
          vertex.insert(&track);
        }
        computeAndSetVertexParameters(vertex, {}, {});
        vertices_.push_back(vertex);
      }
      if (settings_->debug() > 2)
        edm::LogInfo("VertexFinder") << "PVR::DiscardedTracks size " << discardedTracks.size();
      acceptedTracks.clear();
      acceptedTracks = discardedTracks;
    }
  }

sortVerticesInPt()

void l1tVertexFinder::VertexFinder::sortVerticesInPt

(

)

Sort vertices in pT.

Definition at line 572 of file VertexFinder.cc.

References l1tVertexFinder::RecoVertex< T >::pt(), l1t::VertexWord::pt(), and jetUpdater_cfi::sort.

Referenced by VertexProducer::produce().

                                      {
    if (settings_->vx_precision() == Precision::Emulation) {
      std::sort(
          verticesEmulation_.begin(),
          verticesEmulation_.end(),
          [](const l1t::VertexWord& vertex0, const l1t::VertexWord& vertex1) { return (vertex0.pt() > vertex1.pt()); });
    } else {
      std::sort(vertices_.begin(), vertices_.end(), [](const RecoVertex<>& vertex0, const RecoVertex<>& vertex1) {
        return (vertex0.pt() > vertex1.pt());
      });
    }
  }

sortVerticesInZ0()

void l1tVertexFinder::VertexFinder::sortVerticesInZ0

(

)

Sort vertices in z.

Definition at line 585 of file VertexFinder.cc.

References jetUpdater_cfi::sort, l1tVertexFinder::RecoVertex< T >::z0(), and l1t::VertexWord::z0().

                                      {
    if (settings_->vx_precision() == Precision::Emulation) {
      std::sort(
          verticesEmulation_.begin(),
          verticesEmulation_.end(),
          [](const l1t::VertexWord& vertex0, const l1t::VertexWord& vertex1) { return (vertex0.z0() < vertex1.z0()); });
    } else {
      std::sort(vertices_.begin(), vertices_.end(), [](const RecoVertex<>& vertex0, const RecoVertex<>& vertex1) {
        return (vertex0.z0() < vertex1.z0());
      });
    }
  }

strided_iota()

template<typename ForwardIterator , typename T >

void l1tVertexFinder::VertexFinder::strided_iota ( ForwardIterator  first, ForwardIterator  last, T  value, T  stride  )

inline

Definition at line 125 of file VertexFinder.h.

References dqmdumpme::first, dqmdumpme::last, riemannFit::stride, and relativeConstraints::value.

                                                                                      {
      while (first != last) {
        *first++ = value;
        value += stride;
      }
    }

vertices()

const std::vector<RecoVertex<> >& l1tVertexFinder::VertexFinder::vertices ( ) const

inline

Returns the z positions of the reconstructed primary vertices.

Definition at line 76 of file VertexFinder.h.

References vertices_.

Referenced by VertexProducer::produce().

{ return vertices_; }

verticesEmulation()

const l1t::VertexWordCollection& l1tVertexFinder::VertexFinder::verticesEmulation ( ) const

inline

Returns the emulation primary vertices.

Definition at line 78 of file VertexFinder.h.

References verticesEmulation_.

Referenced by VertexProducer::produce().

{ return verticesEmulation_; }

Member Data Documentation

fitTracks_

FitTrackCollection l1tVertexFinder::VertexFinder::fitTracks_

private

Definition at line 152 of file VertexFinder.h.

Referenced by fitTracks(), iterationsPerTrack(), numInputTracks(), and VertexFinder().

iterations_

unsigned int l1tVertexFinder::VertexFinder::iterations_

private

Definition at line 154 of file VertexFinder.h.

Referenced by iterationsPerTrack(), and numIterations().

numMatchedVertices_

unsigned int l1tVertexFinder::VertexFinder::numMatchedVertices_

private

Definition at line 151 of file VertexFinder.h.

pv_index_

unsigned int l1tVertexFinder::VertexFinder::pv_index_

private

Definition at line 153 of file VertexFinder.h.

Referenced by PrimaryVertex(), and primaryVertexId().

settings_

const AlgoSettings* l1tVertexFinder::VertexFinder::settings_

private

Definition at line 148 of file VertexFinder.h.

Referenced by PrimaryVertex(), and VertexFinder().

vertices_

RecoVertexCollection l1tVertexFinder::VertexFinder::vertices_

private

Definition at line 149 of file VertexFinder.h.

Referenced by numVertices(), PrimaryVertex(), and vertices().

verticesEmulation_

l1t::VertexWordCollection l1tVertexFinder::VertexFinder::verticesEmulation_

private

Definition at line 150 of file VertexFinder.h.

Referenced by numVerticesEmulation(), PrimaryVertex(), and verticesEmulation().