From 85300fc43d1d90d5a12538ed54f5c5cad2db51a4 Mon Sep 17 00:00:00 2001 From: ccaillol Date: Mon, 28 Mar 2022 14:53:36 +0200 Subject: [PATCH 1/5] L1T vetex finder phase2 developments --- L1Trigger/VertexFinder/BuildFile.xml | 1 + .../VertexFinder/interface/AlgoSettings.h | 31 +- L1Trigger/VertexFinder/interface/RecoVertex.h | 12 + .../VertexFinder/interface/VertexFinder.h | 76 ++- .../VertexFinder/interface/VertexProducer.h | 5 +- .../plugins/TPStubValueMapProducer.cc | 1 - .../VertexFinder/plugins/VertexNTupler.cc | 103 ++-- .../VertexFinder/plugins/VertexProducer.cc | 75 ++- .../VertexFinder/python/VertexNTupler_cff.py | 4 +- .../VertexFinder/python/VertexProducer_cff.py | 8 +- L1Trigger/VertexFinder/src/AlgoSettings.cc | 29 +- L1Trigger/VertexFinder/src/InputData.cc | 4 +- L1Trigger/VertexFinder/src/VertexFinder.cc | 498 +++++++++++++++++- .../VertexFinder/test/vertexNTupler_cfg.py | 95 +++- 14 files changed, 797 insertions(+), 145 deletions(-) diff --git a/L1Trigger/VertexFinder/BuildFile.xml b/L1Trigger/VertexFinder/BuildFile.xml index 82b4a81c0ce38..2ee9f1b116ce8 100644 --- a/L1Trigger/VertexFinder/BuildFile.xml +++ b/L1Trigger/VertexFinder/BuildFile.xml @@ -1,3 +1,4 @@ + diff --git a/L1Trigger/VertexFinder/interface/AlgoSettings.h b/L1Trigger/VertexFinder/interface/AlgoSettings.h index 5d24dc2ad65dd..b3fb487c2765d 100644 --- a/L1Trigger/VertexFinder/interface/AlgoSettings.h +++ b/L1Trigger/VertexFinder/interface/AlgoSettings.h @@ -10,17 +10,20 @@ namespace l1tVertexFinder { enum class Algorithm { - FastHisto, - FastHistoLooseAssociation, + fastHisto, + fastHistoEmulation, + fastHistoLooseAssociation, GapClustering, - AgglomerativeHierarchical, + agglomerativeHierarchical, DBSCAN, PVR, - AdaptiveVertexReconstruction, + adaptiveVertexReconstruction, HPV, Kmeans }; + enum class Precision { Simulation, Emulation }; + class AlgoSettings { public: AlgoSettings(const edm::ParameterSet& iConfig); @@ -29,7 +32,8 @@ namespace l1tVertexFinder { //=== Vertex Reconstruction configuration // Vertex Reconstruction algo Algorithm vx_algo() const { return vx_algo_; } - /// For Agglomerative cluster algorithm, select a definition of distance between clusters + Precision vx_precision() const { return vx_precision_; } + // For agglomerative cluster algorithm, select a definition of distance between clusters unsigned int vx_distanceType() const { return vx_distanceType_; } // Assumed Vertex Distance float vx_distance() const { return vx_distance_; } @@ -39,21 +43,23 @@ namespace l1tVertexFinder { unsigned int vx_minTracks() const { return vx_minTracks_; } // Compute the z0 position of the vertex with a mean weighted with track momenta unsigned int vx_weightedmean() const { return vx_weightedmean_; } - /// Chi2 cut for the Adaptive Vertex Recostruction Algorithm + // Chi2 cut for the adaptive Vertex Recostruction Algorithm float vx_chi2cut() const { return vx_chi2cut_; } - /// Window size of the sliding window + // Do track quality cuts in emulation algorithms + bool vx_DoQualityCuts() const { return vx_DoQualityCuts_; } + // Window size of the sliding window unsigned int vx_windowSize() const { return vx_windowSize_; } - /// FastHisto histogram parameters (min, max, width) + // fastHisto histogram parameters (min, max, width) std::vector vx_histogram_parameters() const { return vx_histogram_parameters_; } double vx_histogram_min() const { return vx_histogram_parameters_.at(0); } double vx_histogram_max() const { return vx_histogram_parameters_.at(1); } double vx_histogram_binwidth() const { return vx_histogram_parameters_.at(2); } - /// FastHisto assumed vertex width + // fastHisto assumed vertex width float vx_width() const { return vx_width_; } - /// FastHisto track selection control + // fastHisto track selection control bool vx_DoPtComp() const { return vx_DoPtComp_; } bool vx_DoTightChi2() const { return vx_DoTightChi2_; } - /// Number of vertices to return for FastHisto + // Number of vertices to return for fastHisto unsigned int vx_nvtx() const { return vx_nvtx_; } float vx_dbscan_pt() const { return vx_dbscan_pt_; } unsigned int vx_dbscan_mintracks() const { return vx_dbscan_mintracks_; } @@ -82,18 +88,21 @@ namespace l1tVertexFinder { private: static const std::map algoNameMap; + static const std::map algoPrecisionMap; // Parameter sets for differents types of configuration parameter. edm::ParameterSet vertex_; // Vertex Reconstruction configuration Algorithm vx_algo_; + Precision vx_precision_; float vx_distance_; float vx_resolution_; unsigned int vx_distanceType_; unsigned int vx_minTracks_; unsigned int vx_weightedmean_; float vx_chi2cut_; + bool vx_DoQualityCuts_; bool vx_DoPtComp_; bool vx_DoTightChi2_; std::vector vx_histogram_parameters_; diff --git a/L1Trigger/VertexFinder/interface/RecoVertex.h b/L1Trigger/VertexFinder/interface/RecoVertex.h index ca27deb3d5292..801b45124aba0 100644 --- a/L1Trigger/VertexFinder/interface/RecoVertex.h +++ b/L1Trigger/VertexFinder/interface/RecoVertex.h @@ -69,6 +69,8 @@ namespace l1tVertexFinder { unsigned int nHighPt = -999, double highestPt = -999., bool pv = false); + /// Get the equivalent l1t::Vertex + l1t::Vertex vertex() const; /// Vertex z0 position [cm] double z0() const { return z0_; } /// Vertex z0 width [cm] @@ -151,6 +153,16 @@ namespace l1tVertexFinder { pv_ = pv; } + template + l1t::Vertex RecoVertex::vertex() const { + std::vector> tracks; + tracks.reserve(tracks_.size()); + for (const auto& t : tracks_) { + tracks.push_back(t->getTTTrackPtr()); + } + return l1t::Vertex(pT_, z0_, tracks); + } + // Template specializations template <> RecoVertexWithTP& RecoVertexWithTP::operator+=(const RecoVertexWithTP& rhs); diff --git a/L1Trigger/VertexFinder/interface/VertexFinder.h b/L1Trigger/VertexFinder/interface/VertexFinder.h index 7218c2c5ef3fd..e2caa52a68405 100644 --- a/L1Trigger/VertexFinder/interface/VertexFinder.h +++ b/L1Trigger/VertexFinder/interface/VertexFinder.h @@ -2,9 +2,9 @@ #define __L1Trigger_VertexFinder_VertexFinder_h__ #include "DataFormats/Common/interface/Ptr.h" +#include "DataFormats/L1Trigger/interface/VertexWord.h" #include "DataFormats/Math/interface/deltaPhi.h" #include "DataFormats/TrackerCommon/interface/TrackerTopology.h" -#include "FWCore/Framework/interface/ESHandle.h" #include "FWCore/MessageLogger/interface/MessageLogger.h" #include "L1Trigger/VertexFinder/interface/AlgoSettings.h" #include "L1Trigger/VertexFinder/interface/RecoVertex.h" @@ -15,6 +15,9 @@ namespace l1tVertexFinder { + // Returns the number of bits needed to represent and integer decimal + static constexpr unsigned BitsToRepresent(unsigned x) { return x < 2 ? 1 : 1 + BitsToRepresent(x >> 1); } + typedef std::vector FitTrackCollection; typedef std::vector> RecoVertexCollection; @@ -40,20 +43,27 @@ namespace l1tVertexFinder { /// Accessors + /// Storage for tracks out of the L1 Track finder + const FitTrackCollection& fitTracks() const { return fitTracks_; } /// Number of iterations - unsigned int IterationsPerTrack() const { return double(iterations_) / double(fitTracks_.size()); } + unsigned int iterationsPerTrack() const { return double(iterations_) / double(fitTracks_.size()); } /// Storage for tracks out of the L1 Track finder unsigned int numInputTracks() const { return fitTracks_.size(); } /// Number of iterations - unsigned int NumIterations() const { return iterations_; } + unsigned int numIterations() const { return iterations_; } /// Number of reconstructed vertices unsigned int numVertices() const { return vertices_.size(); } - /// Reconstructed Primary Vertex - RecoVertex<> primaryVertex() const { - if (pv_index_ < vertices_.size()) + /// Number of emulation vertices + unsigned int numVerticesEmulation() const { return verticesEmulation_.size(); } + /// Reconstructed primary vertex + template + T PrimaryVertex() const { + 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."; + edm::LogWarning("VertexFinder") << "PrimaryVertex::No primary vertex has been found."; return RecoVertex<>(); } } @@ -61,8 +71,8 @@ namespace l1tVertexFinder { unsigned int primaryVertexId() const { return pv_index_; } /// Returns the z positions of the reconstructed primary vertices const std::vector>& vertices() const { return vertices_; } - /// Storage for tracks out of the L1 Track finder - const FitTrackCollection& fitTracks() const { return fitTracks_; } + /// Returns the emulation primary vertices + const l1t::VertexWordCollection& verticesEmulation() const { return verticesEmulation_; } /// Find the primary vertex void findPrimaryVertex(); @@ -86,22 +96,23 @@ namespace l1tVertexFinder { void fastHistoLooseAssociation(); /// Histogramming algorithm void fastHisto(const TrackerTopology* tTopo); + /// Histogramming algorithm (emulation) + void fastHistoEmulation(); + /// Sort vertices in pT - void SortVerticesInPt() { - std::sort(vertices_.begin(), vertices_.end(), [](const RecoVertex<>& vertex0, const RecoVertex<>& vertex1) { - return (vertex0.pt() > vertex1.pt()); - }); - } + void sortVerticesInPt(); /// Sort vertices in z - void SortVerticesInZ0() { - std::sort(vertices_.begin(), vertices_.end(), [](const RecoVertex<>& vertex0, const RecoVertex<>& vertex1) { - return (vertex0.z0() < vertex1.z0()); - }); - } - /// Number of iterations - unsigned int numIterations() const { return iterations_; } - /// Number of iterations - unsigned int iterationsPerTrack() const { return double(iterations_) / double(fitTracks_.size()); } + void sortVerticesInZ0(); + + /// Print an ASCII histogram + template + void printHistogram(stream_type& stream, + std::vector data, + int width = 80, + int minimum = 0, + int maximum = -1, + std::string title = "", + std::string color = ""); template void strided_iota(ForwardIterator first, ForwardIterator last, T value, T stride) { @@ -113,38 +124,23 @@ namespace l1tVertexFinder { /// Vertexing algorithms - /// Adaptive Vertex Reconstruction algorithm - void AdaptiveVertexReconstruction(); - /// Simple Merge Algorithm - void AgglomerativeHierarchicalClustering(); - /// Find distance between centres of two clusters - float CentralDistance(RecoVertex<> cluster0, RecoVertex<> cluster1); /// Compute the vertex parameters void computeAndSetVertexParameters(RecoVertex<>& vertex, const std::vector& bin_centers, const std::vector& counts); /// DBSCAN algorithm void DBSCAN(); - /// TDR histogramming algorithmn - void FastHistoLooseAssociation(); - /// Histogramming algorithm - void FastHisto(const TrackerTopology* tTopo); /// High pT Vertex Algorithm void HPV(); /// Kmeans Algorithm void Kmeans(); /// Find maximum distance in two clusters of tracks - float MaxDistance(RecoVertex<> cluster0, RecoVertex<> cluster1); - /// Find minimum distance in two clusters of tracks - float MinDistance(RecoVertex<> cluster0, RecoVertex<> cluster1); - /// Find average distance in two clusters of tracks - float MeanDistance(RecoVertex<> cluster0, RecoVertex<> cluster1); - /// Principal Vertex Reconstructor algorithm void PVR(); private: const AlgoSettings* settings_; - std::vector> vertices_; + RecoVertexCollection vertices_; + l1t::VertexWordCollection verticesEmulation_; unsigned int numMatchedVertices_; FitTrackCollection fitTracks_; unsigned int pv_index_; diff --git a/L1Trigger/VertexFinder/interface/VertexProducer.h b/L1Trigger/VertexFinder/interface/VertexProducer.h index 8b9abac9a5fe6..f9d683ff0334a 100644 --- a/L1Trigger/VertexFinder/interface/VertexProducer.h +++ b/L1Trigger/VertexFinder/interface/VertexProducer.h @@ -3,6 +3,8 @@ #include "DataFormats/L1Trigger/interface/Vertex.h" #include "FWCore/Framework/interface/global/EDProducer.h" +#include "DataFormats/L1TrackTrigger/interface/TTTypes.h" +#include "DataFormats/L1Trigger/interface/VertexWord.h" #include "DataFormats/TrackerCommon/interface/TrackerTopology.h" #include "FWCore/Framework/interface/global/EDProducer.h" #include "FWCore/Framework/interface/Event.h" @@ -13,6 +15,7 @@ #include "Geometry/Records/interface/TrackerTopologyRcd.h" #include "L1Trigger/VertexFinder/interface/AlgoSettings.h" #include "L1Trigger/VertexFinder/interface/RecoVertex.h" +#include "L1Trigger/VertexFinder/interface/VertexFinder.h" #include #include @@ -36,7 +39,7 @@ class VertexProducer : public edm::global::EDProducer<> { private: const edm::EDGetTokenT l1TracksToken_; - const edm::ESGetToken trackerTopologyToken_; + const edm::ESGetToken tTopoToken; const std::string outputCollectionName_; l1tVertexFinder::AlgoSettings settings_; diff --git a/L1Trigger/VertexFinder/plugins/TPStubValueMapProducer.cc b/L1Trigger/VertexFinder/plugins/TPStubValueMapProducer.cc index e799bd734c16a..408564410affb 100644 --- a/L1Trigger/VertexFinder/plugins/TPStubValueMapProducer.cc +++ b/L1Trigger/VertexFinder/plugins/TPStubValueMapProducer.cc @@ -77,7 +77,6 @@ class TPStubValueMapProducer : public edm::global::EDProducer<> { const edm::EDGetTokenT stubToken_; const edm::EDGetTokenT stubTruthToken_; const edm::EDGetTokenT clusterTruthToken_; - edm::ESGetToken tTopoToken_; edm::ESGetToken tGeomToken_; }; diff --git a/L1Trigger/VertexFinder/plugins/VertexNTupler.cc b/L1Trigger/VertexFinder/plugins/VertexNTupler.cc index 822e2c269fcee..b8bb5edb9f584 100644 --- a/L1Trigger/VertexFinder/plugins/VertexNTupler.cc +++ b/L1Trigger/VertexFinder/plugins/VertexNTupler.cc @@ -5,12 +5,12 @@ #include "DataFormats/HepMCCandidate/interface/GenParticle.h" #include "DataFormats/L1TrackTrigger/interface/TTTypes.h" #include "DataFormats/L1Trigger/interface/Vertex.h" +#include "DataFormats/L1Trigger/interface/VertexWord.h" #include "DataFormats/JetReco/interface/GenJet.h" #include "DataFormats/Phase2TrackerDigi/interface/Phase2TrackerDigi.h" #include "DataFormats/TrackerCommon/interface/TrackerTopology.h" -#include "FWCore/Framework/interface/EDAnalyzer.h" +#include "FWCore/Framework/interface/one/EDAnalyzer.h" #include "FWCore/Framework/interface/MakerMacros.h" -#include "FWCore/Framework/interface/ESHandle.h" #include "FWCore/Framework/interface/Event.h" #include "FWCore/Framework/interface/EventSetup.h" #include "FWCore/MessageLogger/interface/MessageLogger.h" @@ -40,12 +40,24 @@ using namespace std; namespace l1tVertexFinder { - class VertexNTupler : public edm::EDAnalyzer { + class VertexNTupler : public edm::one::EDAnalyzer<> { public: explicit VertexNTupler(const edm::ParameterSet&); ~VertexNTupler() override; private: + struct EmulationVerticesBranchData { + std::vector numTracks; + std::vector z0; + std::vector sumPt; + + void clear() { + numTracks.clear(); + z0.clear(); + sumPt.clear(); + } + }; + struct GenJetsBranchData { std::vector energy; std::vector pt; @@ -78,11 +90,8 @@ namespace l1tVertexFinder { } }; - struct RecoVerticesBranchData { - std::vector numTracks; + struct RecoVerticesBranchData : public EmulationVerticesBranchData { std::vector> trackIdxs; - std::vector z0; - std::vector sumPt; void clear() { numTracks.clear(); @@ -166,7 +175,7 @@ namespace l1tVertexFinder { std::map> l1TracksTokenMap_; std::map>> l1TracksMapTokenMap_; std::map>> l1VerticesTokenMap_; - std::map>> l1VerticesExtraTokenMap_; + std::map>> l1VerticesEmulationTokenMap_; std::vector>> l1VerticesExtraTokens_; TTree* outputTree_; @@ -176,8 +185,7 @@ namespace l1tVertexFinder { // storage class for configuration parameters AnalysisSettings settings_; - //edm::Service fs_; - + edm::Service fs_; // Histograms for Vertex Reconstruction float numTrueInteractions_, hepMCVtxZ0_, genVtxZ0_; @@ -191,18 +199,12 @@ namespace l1tVertexFinder { std::map l1TracksBranchData_; std::map l1VerticesBranchData_; std::map l1VerticesInputMap_; - - std::unordered_map> l1Vertices_branchMap_numTracks_; - std::unordered_map> l1Vertices_branchMap_z0_; - std::unordered_map> l1Vertices_branchMap_z0_etaWeighted_; - std::unordered_map> l1Vertices_branchMap_sumPt_; + std::map l1VerticesEmulationBranchData_; std::vector> l1Vertices_extra_numTracks_; std::vector> l1Vertices_extra_z0_; std::vector> l1Vertices_extra_z0_etaWeighted_; std::vector> l1Vertices_extra_sumPt_; - - bool available_; // ROOT file for histograms is open. }; VertexNTupler::VertexNTupler(const edm::ParameterSet& iConfig) @@ -215,7 +217,7 @@ namespace l1tVertexFinder { consumes>(iConfig.getParameter("l1TracksTPInputTags"))), vTPsToken_(consumes>( iConfig.getParameter("l1TracksTPValueMapInputTags"))), - //outputTree_(fs_->make("l1VertexReco", "L1 vertex-related info")), + outputTree_(fs_->make("l1VertexReco", "L1 vertex-related info")), printResults_(iConfig.getParameter("printResults")), settings_(iConfig) { const std::vector trackBranchNames( @@ -225,13 +227,6 @@ namespace l1tVertexFinder { const std::vector trackMapInputTags( iConfig.getParameter>("l1TracksTruthMapInputTags")); - edm::Service fs_; - available_ = fs_.isAvailable(); - if (not available_) - return; // No ROOT file open. - - outputTree_ = fs_->make("l1VertexReco", "L1 vertex-related info"); - if (trackBranchNames.size() != trackInputTags.size()) throw cms::Exception("The number of track branch names (" + std::to_string(trackBranchNames.size()) + ") specified in the config does not match the number of input tags (" + @@ -258,6 +253,28 @@ namespace l1tVertexFinder { ") specified in the config does not match the number of associated input track collection names (" + std::to_string(vertexTrackNames.size()) + ")"); + const std::vector emulationVertexBranchNames( + iConfig.getParameter>("emulationVertexBranchNames")); + const std::vector emulationVertexInputTags( + iConfig.getParameter>("emulationVertexInputTags")); + + const std::vector extraVertexInputTags( + iConfig.getParameter>("extraL1VertexInputTags")); + const std::vector extraVertexDescriptions( + iConfig.getParameter>("extraL1VertexDescriptions")); + + std::vector::const_iterator branchNameIt = emulationVertexBranchNames.begin(); + std::vector::const_iterator inputTagIt = emulationVertexInputTags.begin(); + for (; branchNameIt != emulationVertexBranchNames.end(); branchNameIt++, inputTagIt++) { + l1VerticesEmulationTokenMap_[*branchNameIt] = consumes>(*inputTagIt); + l1VerticesEmulationBranchData_[*branchNameIt] = EmulationVerticesBranchData(); + EmulationVerticesBranchData& branchData = l1VerticesEmulationBranchData_.at(*branchNameIt); + + outputTree_->Branch(("emulationVertices_" + *branchNameIt + "_numTracks").c_str(), &branchData.numTracks); + outputTree_->Branch(("emulationVertices_" + *branchNameIt + "_z0").c_str(), &branchData.z0); + outputTree_->Branch(("emulationVertices_" + *branchNameIt + "_sumPt").c_str(), &branchData.sumPt); + } + outputTree_->Branch("genJets_energy", &genJetsBranchData_.energy); outputTree_->Branch("genJets_pt", &genJetsBranchData_.pt); outputTree_->Branch("genJets_eta", &genJetsBranchData_.eta); @@ -302,12 +319,11 @@ namespace l1tVertexFinder { &branchData.truthMapIsUnknown); } - std::vector::const_iterator branchNameIt = vertexBranchNames.begin(); - std::vector::const_iterator inputTagIt = vertexInputTags.begin(); + branchNameIt = vertexBranchNames.begin(); + inputTagIt = vertexInputTags.begin(); std::vector::const_iterator l1VertexTrackNameIt = vertexTrackNames.begin(); for (; branchNameIt != vertexBranchNames.end(); branchNameIt++, inputTagIt++, l1VertexTrackNameIt++) { l1VerticesTokenMap_[*branchNameIt] = consumes>(*inputTagIt); - l1VerticesBranchData_[*branchNameIt] = RecoVerticesBranchData(); RecoVerticesBranchData& branchData = l1VerticesBranchData_.at(*branchNameIt); l1VerticesInputMap_[*branchNameIt] = *l1VertexTrackNameIt; @@ -334,10 +350,6 @@ namespace l1tVertexFinder { outputTree_->Branch("trueTracks_useForAlgEff", &trueTracksBranchData_.useForAlgEff); outputTree_->Branch("trueTracks_useForVtxReco", &trueTracksBranchData_.useForVertexReco); - const std::vector extraVertexInputTags( - iConfig.getParameter>("extraL1VertexInputTags")); - const std::vector extraVertexDescriptions( - iConfig.getParameter>("extraL1VertexDescriptions")); for (const auto& inputTag : extraVertexInputTags) l1VerticesExtraTokens_.push_back(consumes>(inputTag)); TObjArray* descriptionArray = new TObjArray(); @@ -429,6 +441,8 @@ namespace l1tVertexFinder { } } + //const Vertex& TruePrimaryVertex = inputData.getPrimaryVertex(); + // create a map for associating fat reco tracks with their underlying // TTTrack pointers std::map>, const L1TrackTruthMatched*>> @@ -581,6 +595,31 @@ namespace l1tVertexFinder { } } + // Emulation vertex branches + for (const auto& tokenMapEntry : l1VerticesEmulationTokenMap_) { + EmulationVerticesBranchData& branchData = l1VerticesEmulationBranchData_.at(tokenMapEntry.first); + branchData.clear(); + + edm::Handle> handle; + iEvent.getByToken(tokenMapEntry.second, handle); + for (const auto& vtx : *handle) { + branchData.numTracks.push_back(vtx.multiplicity()); + branchData.z0.push_back(vtx.z0()); + branchData.sumPt.push_back(vtx.pt()); + } + + if (printResults_) { + edm::LogInfo("VertexNTupler") << "analyze::" << handle->size() << " '" << tokenMapEntry.first + << "' vertices were found ... "; + for (const auto& vtx : *handle) { + edm::LogInfo("VertexNTupler") << "analyze::" + << " * z0 = " << vtx.z0() << "; contains " << vtx.multiplicity() + << " tracks ..."; + } + } + } + + // Extra vertex collections l1Vertices_extra_numTracks_.resize(l1VerticesExtraTokens_.size()); l1Vertices_extra_z0_.resize(l1VerticesExtraTokens_.size()); l1Vertices_extra_z0_etaWeighted_.resize(l1VerticesExtraTokens_.size()); diff --git a/L1Trigger/VertexFinder/plugins/VertexProducer.cc b/L1Trigger/VertexFinder/plugins/VertexProducer.cc index 3499dca617841..3fe3288bf9dd2 100644 --- a/L1Trigger/VertexFinder/plugins/VertexProducer.cc +++ b/L1Trigger/VertexFinder/plugins/VertexProducer.cc @@ -18,23 +18,27 @@ using namespace std; VertexProducer::VertexProducer(const edm::ParameterSet& iConfig) : l1TracksToken_(consumes(iConfig.getParameter("l1TracksInputTag"))), - trackerTopologyToken_(esConsumes()), + tTopoToken(esConsumes()), outputCollectionName_(iConfig.getParameter("l1VertexCollectionName")), settings_(AlgoSettings(iConfig)) { // Get configuration parameters switch (settings_.vx_algo()) { - case Algorithm::FastHisto: - edm::LogInfo("VertexProducer") << "VertexProducer::Finding vertices using the FastHisto binning algorithm"; + case Algorithm::fastHisto: + edm::LogInfo("VertexProducer") << "VertexProducer::Finding vertices using the fastHisto binning algorithm"; break; - case Algorithm::FastHistoLooseAssociation: + case Algorithm::fastHistoEmulation: edm::LogInfo("VertexProducer") - << "VertexProducer::Finding vertices using the FastHistoLooseAssociation binning algorithm"; + << "VertexProducer::Finding vertices using the emulation version of the fastHisto binning algorithm"; + break; + case Algorithm::fastHistoLooseAssociation: + edm::LogInfo("VertexProducer") + << "VertexProducer::Finding vertices using the fastHistoLooseAssociation binning algorithm"; break; case Algorithm::GapClustering: edm::LogInfo("VertexProducer") << "VertexProducer::Finding vertices using a gap clustering algorithm"; break; - case Algorithm::AgglomerativeHierarchical: + case Algorithm::agglomerativeHierarchical: edm::LogInfo("VertexProducer") << "VertexProducer::Finding vertices using a Simple Merge Clustering algorithm"; break; case Algorithm::DBSCAN: @@ -43,7 +47,7 @@ VertexProducer::VertexProducer(const edm::ParameterSet& iConfig) case Algorithm::PVR: edm::LogInfo("VertexProducer") << "VertexProducer::Finding vertices using a PVR algorithm"; break; - case Algorithm::AdaptiveVertexReconstruction: + case Algorithm::adaptiveVertexReconstruction: edm::LogInfo("VertexProducer") << "VertexProducer::Finding vertices using an AdaptiveVertexReconstruction algorithm"; break; @@ -56,7 +60,11 @@ VertexProducer::VertexProducer(const edm::ParameterSet& iConfig) } //--- Define EDM output to be written to file (if required) - produces(outputCollectionName_); + if (settings_.vx_algo() == Algorithm::fastHistoEmulation) { + produces(outputCollectionName_ + "Emulation"); + } else { + produces(outputCollectionName_); + } } void VertexProducer::produce(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const { @@ -65,31 +73,46 @@ void VertexProducer::produce(edm::StreamID, edm::Event& iEvent, const edm::Event std::vector l1Tracks; l1Tracks.reserve(l1TracksHandle->size()); + if (settings_.debug() > 1) { + edm::LogInfo("VertexProducer") << "produce::Processing " << l1TracksHandle->size() << " tracks"; + } for (const auto& track : l1TracksHandle->ptrs()) { auto l1track = L1Track(track); // Check the minimum pT of the tracks // This is left here because it represents the smallest pT to be sent by the track finding boards // This has less to do with the algorithms than the constraints of what will be sent to the vertexing algorithm - if (l1track.pt() > settings_.vx_TrackMinPt()) { + if (l1track.pt() >= settings_.vx_TrackMinPt()) { l1Tracks.push_back(l1track); + } else { + if (settings_.debug() > 2) { + edm::LogInfo("VertexProducer") << "produce::Removing track with too low of a pt (" << l1track.pt() << ")\n" + << " word = " << l1track.getTTTrackPtr()->getTrackWord().to_string(2); + } } } + if (settings_.debug() > 1) { + edm::LogInfo("VertexProducer") << "produce::Processing " << l1Tracks.size() << " tracks after minimum pt cut of" + << settings_.vx_TrackMinPt() << " GeV"; + } VertexFinder vf(l1Tracks, settings_); switch (settings_.vx_algo()) { - case Algorithm::FastHisto: { - edm::ESHandle tTopoHandle = iSetup.getHandle(trackerTopologyToken_); - vf.fastHisto(tTopoHandle.product()); + case Algorithm::fastHisto: { + const TrackerTopology& tTopo = iSetup.getData(tTopoToken); + vf.fastHisto(&tTopo); break; } - case Algorithm::FastHistoLooseAssociation: + case Algorithm::fastHistoEmulation: + vf.fastHistoEmulation(); + break; + case Algorithm::fastHistoLooseAssociation: vf.fastHistoLooseAssociation(); break; case Algorithm::GapClustering: vf.GapClustering(); break; - case Algorithm::AgglomerativeHierarchical: + case Algorithm::agglomerativeHierarchical: vf.agglomerativeHierarchicalClustering(); break; case Algorithm::DBSCAN: @@ -98,7 +121,7 @@ void VertexProducer::produce(edm::StreamID, edm::Event& iEvent, const edm::Event case Algorithm::PVR: vf.PVR(); break; - case Algorithm::AdaptiveVertexReconstruction: + case Algorithm::adaptiveVertexReconstruction: vf.adaptiveVertexReconstruction(); break; case Algorithm::HPV: @@ -109,20 +132,20 @@ void VertexProducer::produce(edm::StreamID, edm::Event& iEvent, const edm::Event break; } - vf.SortVerticesInPt(); + vf.sortVerticesInPt(); vf.findPrimaryVertex(); // //=== Store output EDM track and hardware stub collections. - std::unique_ptr lProduct(new std::vector()); - - for (const auto& vtx : vf.vertices()) { - std::vector> lVtxTracks; - lVtxTracks.reserve(vtx.tracks().size()); - for (const auto& t : vtx.tracks()) - lVtxTracks.push_back(t->getTTTrackPtr()); - lProduct->emplace_back(l1t::Vertex(vtx.pt(), vtx.z0(), lVtxTracks)); + if (settings_.vx_algo() == Algorithm::fastHistoEmulation) { + std::unique_ptr product_emulation = + std::make_unique(vf.verticesEmulation().begin(), vf.verticesEmulation().end()); + iEvent.put(std::move(product_emulation), outputCollectionName_ + "Emulation"); + } else { + std::unique_ptr product(new std::vector()); + for (const auto& vtx : vf.vertices()) { + product->emplace_back(vtx.vertex()); + } + iEvent.put(std::move(product), outputCollectionName_); } - iEvent.put(std::move(lProduct), outputCollectionName_); } - DEFINE_FWK_MODULE(VertexProducer); diff --git a/L1Trigger/VertexFinder/python/VertexNTupler_cff.py b/L1Trigger/VertexFinder/python/VertexNTupler_cff.py index 3e2bd18f499f3..185fa67d775d4 100644 --- a/L1Trigger/VertexFinder/python/VertexNTupler_cff.py +++ b/L1Trigger/VertexFinder/python/VertexNTupler_cff.py @@ -11,7 +11,9 @@ l1TracksBranchNames = cms.vstring('hybrid'), l1VertexInputTags = cms.VInputTag( cms.InputTag("VertexProducer", VertexProducer.l1VertexCollectionName.value()) ), l1VertexTrackInputs = cms.vstring('hybrid'), - l1VertexBranchNames = cms.vstring('FastHisto'), + l1VertexBranchNames = cms.vstring('fastHisto'), + emulationVertexInputTags = cms.VInputTag(), + emulationVertexBranchNames = cms.vstring(), extraL1VertexInputTags = cms.VInputTag(), extraL1VertexDescriptions = cms.vstring(), diff --git a/L1Trigger/VertexFinder/python/VertexProducer_cff.py b/L1Trigger/VertexFinder/python/VertexProducer_cff.py index de410f32adb36..f946fbb0f1390 100644 --- a/L1Trigger/VertexFinder/python/VertexProducer_cff.py +++ b/L1Trigger/VertexFinder/python/VertexProducer_cff.py @@ -9,7 +9,7 @@ # === Vertex Reconstruction configuration VertexReconstruction = cms.PSet( # Vertex Reconstruction Algorithm - Algorithm = cms.string("FastHisto"), + Algorithm = cms.string("fastHisto"), # Vertex distance [cm] VertexDistance = cms.double(.15), # Assumed Vertex Resolution [cm] @@ -25,18 +25,20 @@ WeightedMean = cms.uint32(1), # Chi2 cut for the Adaptive Vertex Reconstruction Algorithm AVR_chi2cut = cms.double(5.), + # Do track quality cuts in emulation algorithms + EM_DoQualityCuts = cms.bool(False), # Track-stubs Pt compatibility cut FH_DoPtComp = cms.bool(True), # chi2dof < 5 for tracks with Pt > 10 FH_DoTightChi2 = cms.bool(False), - # FastHisto algorithm histogram parameters (min,max,width) [cm] + # fastHisto algorithm histogram parameters (min,max,width) [cm] # TDR settings: [-14.95, 15.0, 0.1] # L1TkPrimaryVertexProducer: [-30.0, 30.0, 0.09983361065] # Firmware: [-14.4, 14.4, 0.4] FH_HistogramParameters = cms.vdouble(-30.0, 30.0, 0.09983361065), # The number of vertixes to return (i.e. N windows with the highest combined pT) FH_NVtx = cms.uint32(10), - # FastHisto algorithm assumed vertex half-width [cm] + # fastHisto algorithm assumed vertex half-width [cm] FH_VertexWidth = cms.double(.15), # Window size of the sliding window FH_WindowSize = cms.uint32(3), diff --git a/L1Trigger/VertexFinder/src/AlgoSettings.cc b/L1Trigger/VertexFinder/src/AlgoSettings.cc index 52f5aa461ae79..f111a8bcbb63c 100644 --- a/L1Trigger/VertexFinder/src/AlgoSettings.cc +++ b/L1Trigger/VertexFinder/src/AlgoSettings.cc @@ -12,6 +12,7 @@ namespace l1tVertexFinder { vx_minTracks_(vertex_.getParameter("MinTracks")), vx_weightedmean_(vertex_.getParameter("WeightedMean")), vx_chi2cut_(vertex_.getParameter("AVR_chi2cut")), + vx_DoQualityCuts_(vertex_.getParameter("EM_DoQualityCuts")), vx_DoPtComp_(vertex_.getParameter("FH_DoPtComp")), vx_DoTightChi2_(vertex_.getParameter("FH_DoTightChi2")), vx_histogram_parameters_(vertex_.getParameter >("FH_HistogramParameters")), @@ -44,17 +45,37 @@ namespace l1tVertexFinder { throw cms::Exception("Invalid algo name '" + algoName + "' specified for L1T vertex producer. Valid algo names are: " + validAlgoNames.str()); } + + const auto algoPrecisionMapIt = algoPrecisionMap.find(vx_algo_); + if (algoPrecisionMapIt != algoPrecisionMap.end()) { + vx_precision_ = algoPrecisionMapIt->second; + } else { + throw cms::Exception("Unknown precision {Simulation, Emulation} for algo name " + algoName); + } } const std::map AlgoSettings::algoNameMap = { - {"FastHisto", Algorithm::FastHisto}, - {"FastHistoLooseAssociation", Algorithm::FastHistoLooseAssociation}, + {"fastHisto", Algorithm::fastHisto}, + {"fastHistoEmulation", Algorithm::fastHistoEmulation}, + {"fastHistoLooseAssociation", Algorithm::fastHistoLooseAssociation}, {"GapClustering", Algorithm::GapClustering}, - {"Agglomerative", Algorithm::AgglomerativeHierarchical}, + {"agglomerative", Algorithm::agglomerativeHierarchical}, {"DBSCAN", Algorithm::DBSCAN}, {"PVR", Algorithm::PVR}, - {"Adaptive", Algorithm::AdaptiveVertexReconstruction}, + {"adaptive", Algorithm::adaptiveVertexReconstruction}, {"HPV", Algorithm::HPV}, {"K-means", Algorithm::Kmeans}}; + const std::map AlgoSettings::algoPrecisionMap = { + {Algorithm::fastHisto, Precision::Simulation}, + {Algorithm::fastHistoEmulation, Precision::Emulation}, + {Algorithm::fastHistoLooseAssociation, Precision::Simulation}, + {Algorithm::GapClustering, Precision::Simulation}, + {Algorithm::agglomerativeHierarchical, Precision::Simulation}, + {Algorithm::DBSCAN, Precision::Simulation}, + {Algorithm::PVR, Precision::Simulation}, + {Algorithm::adaptiveVertexReconstruction, Precision::Simulation}, + {Algorithm::HPV, Precision::Simulation}, + {Algorithm::Kmeans, Precision::Simulation}}; + } // end namespace l1tVertexFinder diff --git a/L1Trigger/VertexFinder/src/InputData.cc b/L1Trigger/VertexFinder/src/InputData.cc index 46fd49c203a20..a0413b72b4908 100644 --- a/L1Trigger/VertexFinder/src/InputData.cc +++ b/L1Trigger/VertexFinder/src/InputData.cc @@ -85,7 +85,7 @@ namespace l1tVertexFinder { } else { genPt_PU_ += tp->pt(); } - if (settings.debug() > 0) { + if (settings.debug() > 2) { edm::LogInfo("InputData") << "InputData::genPt in the event " << genPt_; } @@ -113,7 +113,7 @@ namespace l1tVertexFinder { recoVertices_.push_back(vertex); } - if (settings.debug() > 0) + if (settings.debug() > 2) edm::LogInfo("InputData") << "InputData::" << vertices_.size() << " pileup vertices in the event, " << recoVertices_.size() << " reconstructable"; diff --git a/L1Trigger/VertexFinder/src/VertexFinder.cc b/L1Trigger/VertexFinder/src/VertexFinder.cc index 875d175d24219..3c6fe46c7bb24 100644 --- a/L1Trigger/VertexFinder/src/VertexFinder.cc +++ b/L1Trigger/VertexFinder/src/VertexFinder.cc @@ -8,7 +8,7 @@ namespace l1tVertexFinder { const std::vector& bin_centers, const std::vector& counts) { double pt = 0.; - double z0 = 0.; + double z0 = -999.; double z0width = 0.; bool highPt = false; double highestPt = 0.; @@ -25,6 +25,11 @@ namespace l1tVertexFinder { 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++; @@ -126,7 +131,7 @@ namespace l1tVertexFinder { sort(fitTracks_.begin(), fitTracks_.end(), SortTracksByZ0()); - std::vector> vClusters; + RecoVertexCollection vClusters; vClusters.resize(fitTracks_.size()); for (unsigned int i = 0; i < fitTracks_.size(); ++i) { @@ -438,14 +443,158 @@ namespace l1tVertexFinder { } void VertexFinder::findPrimaryVertex() { - double vertexPt = 0; - pv_index_ = 0; + 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()); + })); + } + } + + // Possible Formatting Codes: https://misc.flogisoft.com/bash/tip_colors_and_formatting + template + void VertexFinder::printHistogram(stream_type& stream, + std::vector data, + int width, + int minimum, + int maximum, + std::string title, + std::string color) { + int tableSize = data.size(); + + if (minimum < 0) { + minimum *= 1.05; + } else { + minimum = 0; + } + + if (maximum == -1) { + maximum = float(*std::max_element(std::begin(data), std::end(data))) * 1.05; + } + + if (maximum <= minimum) { + float average = (minimum + maximum) / 2.0; + minimum = average - 0.5; + maximum = average + 0.5; + } - for (unsigned int i = 0; i < vertices_.size(); ++i) { - if (vertices_[i].pt() > vertexPt) { - vertexPt = vertices_[i].pt(); - pv_index_ = i; + std::vector intervals(tableSize, ""); + std::vector 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 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 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"; + } + + void VertexFinder::sortVerticesInPt() { + 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()); + }); + } + } + + void VertexFinder::sortVerticesInZ0() { + 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()); + }); } } @@ -552,8 +701,13 @@ namespace l1tVertexFinder { } // Assign the track to the correct vertex - auto upper_bound = std::lower_bound(bounds.begin(), bounds.end(), track.z0()); + // 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 @@ -588,11 +742,335 @@ namespace l1tVertexFinder { 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 tmp; + std::transform(std::begin(sums), std::end(sums), std::back_inserter(tmp), [](const RecoVertex<>& v) -> double { + return v.pt(); + }); + printHistogram(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 + void VertexFinder::fastHistoEmulation() { + // Relevant constants for the track word + enum TrackBitWidths { + kZ0Size = 12, // Width of z-position (40cm / 0.1) + kZ0MagSize = 5, // Width of z-position magnitude (signed) + kPtSize = 14, // Width of pt + kPtMagSize = 9, // Width of pt magnitude (unsigned) + kReducedPrecisionPt = 7, // Width of the reduced precision, integer only, pt + }; + + enum HistogramBitWidths { + kBinSize = 10, // Width of a single bin in z + kBinFixedSize = 7, // Width of a single z0 bin in fixed point representation + kBinFixedMagSize = 4, // 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(HistogramBitWidths::kWindowSize * (1 << HistogramBitWidths::kSumPtLinkSize)), + }; + + static constexpr unsigned int kTableSize = + ((1 << HistogramBitWidths::kSumPtLinkSize) - 1) * HistogramBitWidths::kWindowSize; + static constexpr double kZ0Scale = + (TTTrack_TrackWord::stepZ0 * + (1 << (TrackBitWidths::kZ0Size - TrackBitWidths::kZ0MagSize))); // scale = 1.27932032 + + typedef ap_ufixed pt_t; + typedef ap_fixed 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 + track_pt_fixed_t; + // Histogram bin index + typedef ap_uint histbin_t; + // Histogram bin in fixed point representation, before truncation + typedef ap_ufixed + histbin_fixed_t; + // This value is slightly arbitrary, but small enough that the windows sums aren't too big. + typedef ap_ufixed + link_pt_sum_fixed_t; + // Enough bits to store HistogramBitWidths::kWindowSize * (2**HistogramBitWidths::kSumPtLinkSize) + typedef ap_ufixed + window_pt_sum_fixed_t; + // pt weighted sum of bins in window + typedef ap_fixed + zsliding_t; + // Sum of histogram bins in window + typedef ap_uint slidingsum_t; + // Inverse of sum of bins in a given window + typedef ap_ufixed + 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 bin = (z0 * histbin_fixed_t(1.0 / settings_->vx_histogram_binwidth())) + + histbin_t(std::floor( + nbins / 2.)); // Rounding down (std::floor) taken care of by implicitly casting to ap_uint + 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 { + std::vector table_out(kTableSize, 0.); + for (unsigned int ii = 0; ii < kTableSize; ii++) { + // First, convert from table index to X-value (unsigned 8-bit, range 0 to +1533) + float in_val = 1533.0 * (ii / float(kTableSize)); + // Next, compute lookup table function + table_out.at(ii) = (in_val > 0) ? (1.0 / in_val) : 0.0; + } + return table_out; + }; + + auto inversion = [&](slidingsum_t& data_den) -> inverse_t { + std::vector 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) -> z0_t { + zsliding_t z = iz - histbin_t(std::floor(nbins / 2.)); + std::unique_ptr log; + if (settings_->debug() >= 1) { + log = std::make_unique("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) + << z0_t(z * zsliding_t(settings_->vx_histogram_binwidth())); + } + return z0_t(z * zsliding_t(settings_->vx_histogram_binwidth())); + }; + + auto weighted_position = [&](histbin_t b_max, + const std::vector& 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 log; + if (settings_->debug() >= 1) { + log = std::make_unique("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 w = 0; w < HistogramBitWidths::kWindowSize; ++w) { + zvtx_sliding_sum += (binpt.at(w) * w); + if (settings_->debug() >= 1) { + *log << "(" << w << " * " << binpt.at(w) << ")"; + if (w < HistogramBitWidths::kWindowSize - 1) { + *log << " + "; + } + } + } + + if (settings_->debug() >= 1) { + *log << " = " << zvtx_sliding_sum << "\n"; + } + + if (maximums != 0) { + inv = inversion(maximums); + 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::ceil((settings_->vx_histogram_max() - settings_->vx_histogram_min()) / settings_->vx_histogram_binwidth()); + unsigned int nsums = nbins - settings_->vx_windowSize(); + std::vector hist(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 saturdation taken care of by the data type specification + pt_t tkpt = 0; + tkpt.V = track.getTTTrackPtr()->getTrackWord()(TTTrack_TrackWord::TrackBitLocations::kRinvMSB - 1, + TTTrack_TrackWord::TrackBitLocations::kRinvLSB); + track_pt_fixed_t pt_tmp = tkpt; + //z0_t tkZ0 = track.getTTTrackPtr()->getZ0(); + z0_t tkZ0 = 0; + tkZ0.V = track.getTTTrackPtr()->getTrackWord()(TTTrack_TrackWord::TrackBitLocations::kZ0MSB, + TTTrack_TrackWord::TrackBitLocations::kZ0LSB); + ap_ufixed<32, 1> kZ0Scale_fixed = kZ0Scale; + tkZ0 *= kZ0Scale_fixed; + + if ((settings_->vx_DoQualityCuts() && track_quality_check(tkpt)) || (!settings_->vx_DoQualityCuts())) { + // + // Check bin validity of bin found for the current track + // + std::pair 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) + << ")\tpt_tmp = " << pt_tmp << "\tbin = " << bin.first.to_int() << "\n" + << "pt sum in bin " << bin.first.to_int() + << " BEFORE adding track = " << hist.at(bin.first).to_double(); + } + if (bin.second) { + hist.at(bin.first) = hist.at(bin.first) + pt_tmp; + } + if (settings_->debug() > 2) { + edm::LogInfo("VertexProducer") << "fastHistoEmulation::\npt sum in bin " << bin.first.to_int() + << " AFTER adding track = " << hist.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) + << ")\tpt_tmp = " << pt_tmp; + } + } + } // end loop over tracks + + // 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 hist_window_sums(nsums, 0); + for (unsigned int b = 0; b < nsums; ++b) { + for (unsigned int w = 0; w < HistogramBitWidths::kWindowSize; ++w) { + unsigned int index = b + w; + hist_window_sums.at(b) += hist.at(index); + } + } + + // Find the top N vertices + std::vector 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 binpt_max(HistogramBitWidths::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 + HistogramBitWidths::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(log, hist, 80, 0, -1, "fastHistoEmulation::hist", "\e[92m"); + printHistogram( + log, hist_window_sums, 80, 0, -1, "fastHistoEmulation::hist_window_sums", "\e[92m"); + } + printHistogram( + 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 + } // namespace l1tVertexFinder diff --git a/L1Trigger/VertexFinder/test/vertexNTupler_cfg.py b/L1Trigger/VertexFinder/test/vertexNTupler_cfg.py index fca514676091d..d2e5bf974395a 100644 --- a/L1Trigger/VertexFinder/test/vertexNTupler_cfg.py +++ b/L1Trigger/VertexFinder/test/vertexNTupler_cfg.py @@ -1,23 +1,56 @@ import FWCore.ParameterSet.Config as cms import FWCore.Utilities.FileUtils as FileUtils import FWCore.ParameterSet.VarParsing as VarParsing +import pkgutil +import sys # PART 1 : PARSE ARGUMENTS options = VarParsing.VarParsing ('analysis') +options.register('redir', 'root://cms-xrd-global.cern.ch/', VarParsing.VarParsing.multiplicity.singleton, VarParsing.VarParsing.varType.string, "The XRootD redirector to use") +options.register('nstart', 0,VarParsing.VarParsing.multiplicity.singleton, VarParsing.VarParsing.varType.int, "File index to start on") +options.register('nfiles', -1,VarParsing.VarParsing.multiplicity.singleton, VarParsing.VarParsing.varType.int, "Number of files to process per job") options.register('storeTracks', False, VarParsing.VarParsing.multiplicity.singleton, VarParsing.VarParsing.varType.bool, "Store tracks in NTuple") options.register('l1Tracks','TTTracksFromTrackletEmulation:Level1TTTracks', VarParsing.VarParsing.multiplicity.singleton, VarParsing.VarParsing.varType.string, 'L1 track collection to use') options.register('runVariations', False, VarParsing.VarParsing.multiplicity.singleton, VarParsing.VarParsing.varType.bool, "Run some pre-defined algorithmic variations") -options.register('threads',1,VarParsing.VarParsing.multiplicity.singleton, VarParsing.VarParsing.varType.int, "Number of threads/streams to run") +options.register('threads', 1,VarParsing.VarParsing.multiplicity.singleton, VarParsing.VarParsing.varType.int, "Number of threads to run") +options.register('streams', 0,VarParsing.VarParsing.multiplicity.singleton, VarParsing.VarParsing.varType.int, "Number of streams to run") +options.register('memoryProfiler', False, VarParsing.VarParsing.multiplicity.singleton, VarParsing.VarParsing.varType.bool, "Run the memory profile") +options.register('tmi', False, VarParsing.VarParsing.multiplicity.singleton, VarParsing.VarParsing.varType.bool, "Run a simple profiler") +options.register('trace', False, VarParsing.VarParsing.multiplicity.singleton, VarParsing.VarParsing.varType.bool, "Dump the paths and consumes") +options.register('dump', False, VarParsing.VarParsing.multiplicity.singleton, VarParsing.VarParsing.varType.bool, "Dump the configuration and exit") options.parseArguments() -inputFiles = [] +# handle site name usage +if options.redir[0]=="T": + options.redir = "root://cms-xrd-global.cern.ch//store/test/xrootd/"+options.redir + +# Load input files +inputFiles = cms.untracked.vstring() + for filePath in options.inputFiles: - if filePath.endswith(".root"): - inputFiles.append(filePath) + if filePath.endswith(".root") : + inputFiles.append( filePath ) + elif filePath.endswith(".txt"): + inputFiles += FileUtils.loadListFromFile( filePath ) + elif filePath.endswith("_cff.py"): + inputFilesImport = getattr(__import__(filePath.strip(".py").strip("python/").replace('/','.'),fromlist=["readFiles"]),"readFiles") + if options.nfiles==-1: + inputFiles.extend( inputFilesImport ) + else: + inputFiles.extend( inputFilesImport[options.nstart:(options.nstart+options.nfiles)] ) + elif pkgutil.find_loader("L1Trigger.VertexFinder."+filePath+"_cff") is not None: + inputFilesImport = getattr(__import__("L1Trigger.VertexFinder."+filePath+"_cff",fromlist=["readFiles"]),"readFiles") + if options.nfiles==-1: + inputFiles.extend( inputFilesImport ) + else: + inputFiles.extend( inputFilesImport[options.nstart:(options.nstart+options.nfiles)] ) else: - inputFiles += FileUtils.loadListFromFile(filePath) + raise RuntimeError("Must specify a list of ROOT files, a list of txt files containing a list of ROOT files, or a list of python input files.") + +if options.redir != "": + inputFiles = [(options.redir if val.startswith("/") else "")+val for val in inputFiles] if options.l1Tracks.count(':') != 1: raise RuntimeError("Value for 'l1Tracks' command-line argument (= '{}') should contain one colon".format(options.l1Tracks)) @@ -36,14 +69,14 @@ process.load('Configuration.StandardSequences.FrontierConditions_GlobalTag_cff') from Configuration.AlCa.GlobalTag import GlobalTag process.GlobalTag = GlobalTag(process.GlobalTag, 'auto:upgradePLS3', '') -process.load("FWCore.MessageLogger.MessageLogger_cfi") +process.load("FWCore.MessageService.MessageLogger_cfi") process.source = cms.Source("PoolSource", fileNames = cms.untracked.vstring(inputFiles) ) process.TFileService = cms.Service("TFileService", fileName = cms.string(options.outputFile)) process.maxEvents = cms.untracked.PSet( input = cms.untracked.int32(options.maxEvents) ) process.options = cms.untracked.PSet( numberOfThreads = cms.untracked.uint32(options.threads), - numberOfStreams = cms.untracked.uint32(options.threads if options.threads>0 else 0) + numberOfStreams = cms.untracked.uint32(options.streams if options.streams>0 else 0) ) process.load('L1Trigger.VertexFinder.VertexProducer_cff') @@ -60,20 +93,30 @@ process.Timing = cms.Service("Timing", summaryOnly = cms.untracked.bool(True)) producerSum = process.VertexProducer -additionalProducerAlgorithms = ["FastHistoLooseAssociation", "DBSCAN"] +additionalProducerAlgorithms = ["fastHistoEmulation", "fastHistoLooseAssociation", "DBSCAN"] for algo in additionalProducerAlgorithms: producerName = 'VertexProducer{0}'.format(algo) producerName = producerName.replace(".","p") # legalize the name producer = process.VertexProducer.clone() producer.VertexReconstruction.Algorithm = cms.string(algo) + + if "Emulation" in algo: + if "L1GTTInputProducer" not in process.producerNames(): + process.load('L1Trigger.L1TTrackMatch.L1GTTInputProducer_cfi') + producer.l1TracksInputTag = cms.InputTag("L1GTTInputProducer","Level1TTTracksConverted") + producerSum = process.L1GTTInputProducer + producerSum + + process.L1TVertexNTupler.emulationVertexInputTags.append( cms.InputTag(producerName, 'l1verticesEmulation') ) + process.L1TVertexNTupler.emulationVertexBranchNames.append(algo) + else: + process.L1TVertexNTupler.l1VertexInputTags.append( cms.InputTag(producerName, 'l1vertices') ) + process.L1TVertexNTupler.l1VertexBranchNames.append(algo) + process.L1TVertexNTupler.l1VertexTrackInputs.append('hybrid') + setattr(process, producerName, producer) producerSum += producer - process.L1TVertexNTupler.l1VertexInputTags.append( cms.InputTag(producerName, 'l1vertices') ) - process.L1TVertexNTupler.l1VertexBranchNames.append(algo) - process.L1TVertexNTupler.l1VertexTrackInputs.append('hybrid') - # PART 3: PERFORM SCAN OVER ALGO PARAMETER SPACE if options.runVariations: for i in range(1, 9): @@ -108,7 +151,31 @@ print "Total number of producers =", len(additionalProducerAlgorithms)+1 print " Producers = [{0}]".format(producerSum.dumpSequenceConfig().replace('&',', ')) -print " Algorithms = [FastHisto, {0}]".format(', '.join(additionalProducerAlgorithms)) - +print " Algorithms = [fastHisto, {0}]".format(', '.join(additionalProducerAlgorithms)) + +# PART 4: UTILITIES + +# MEMORY PROFILING +if options.memoryProfiler: + process.IgProfService = cms.Service("IgProfService", + reportEventInterval = cms.untracked.int32(1), + reportFirstEvent = cms.untracked.int32(1), + reportToFileAtPostEndJob = cms.untracked.string('| gzip -c > '+options.outputFile+'___memory___%I_EndOfJob.gz'), + reportToFileAtPostEvent = cms.untracked.string('| gzip -c > '+options.outputFile+'___memory___%I.gz') + ) + +# SIMPLER PROFILING +if options.tmi: + from Validation.Performance.TimeMemoryInfo import customise + process = customise(process) + +if options.trace: + process.add_(cms.Service("Tracer", dumpPathsAndConsumes = cms.untracked.bool(True))) + +# SETUP THE PATH process.p = cms.Path(producerSum + process.TPStubValueMapProducer + process.InputDataProducer + process.L1TVertexNTupler) +# DUMP AND EXIT +if options.dump: + print process.dumpPython() + sys.exit(0) From 928641a309464f702aa3cf78d5d96312b3c2daf3 Mon Sep 17 00:00:00 2001 From: ccaillol Date: Mon, 28 Mar 2022 15:34:47 +0200 Subject: [PATCH 2/5] fix clang warning --- L1Trigger/VertexFinder/plugins/VertexNTupler.cc | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/L1Trigger/VertexFinder/plugins/VertexNTupler.cc b/L1Trigger/VertexFinder/plugins/VertexNTupler.cc index b8bb5edb9f584..38786407ce6f6 100644 --- a/L1Trigger/VertexFinder/plugins/VertexNTupler.cc +++ b/L1Trigger/VertexFinder/plugins/VertexNTupler.cc @@ -40,7 +40,7 @@ using namespace std; namespace l1tVertexFinder { - class VertexNTupler : public edm::one::EDAnalyzer<> { + class VertexNTupler : public edm::one::EDAnalyzer { public: explicit VertexNTupler(const edm::ParameterSet&); ~VertexNTupler() override; @@ -185,7 +185,6 @@ namespace l1tVertexFinder { // storage class for configuration parameters AnalysisSettings settings_; - edm::Service fs_; // Histograms for Vertex Reconstruction float numTrueInteractions_, hepMCVtxZ0_, genVtxZ0_; @@ -217,7 +216,7 @@ namespace l1tVertexFinder { consumes>(iConfig.getParameter("l1TracksTPInputTags"))), vTPsToken_(consumes>( iConfig.getParameter("l1TracksTPValueMapInputTags"))), - outputTree_(fs_->make("l1VertexReco", "L1 vertex-related info")), + //outputTree_(fs_->make("l1VertexReco", "L1 vertex-related info")), printResults_(iConfig.getParameter("printResults")), settings_(iConfig) { const std::vector trackBranchNames( @@ -263,6 +262,10 @@ namespace l1tVertexFinder { const std::vector extraVertexDescriptions( iConfig.getParameter>("extraL1VertexDescriptions")); + usesResource(TFileService::kSharedResource); + edm::Service fs; + outputTree_=fs->make("l1VertexReco", "L1 vertex-related info"); + std::vector::const_iterator branchNameIt = emulationVertexBranchNames.begin(); std::vector::const_iterator inputTagIt = emulationVertexInputTags.begin(); for (; branchNameIt != emulationVertexBranchNames.end(); branchNameIt++, inputTagIt++) { From daf0e11fe09c083c351db95440f10b243a80b67e Mon Sep 17 00:00:00 2001 From: ccaillol Date: Mon, 28 Mar 2022 15:39:18 +0200 Subject: [PATCH 3/5] fix clang warning --- L1Trigger/VertexFinder/plugins/VertexNTupler.cc | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/L1Trigger/VertexFinder/plugins/VertexNTupler.cc b/L1Trigger/VertexFinder/plugins/VertexNTupler.cc index 38786407ce6f6..13c535d0c9180 100644 --- a/L1Trigger/VertexFinder/plugins/VertexNTupler.cc +++ b/L1Trigger/VertexFinder/plugins/VertexNTupler.cc @@ -216,7 +216,6 @@ namespace l1tVertexFinder { consumes>(iConfig.getParameter("l1TracksTPInputTags"))), vTPsToken_(consumes>( iConfig.getParameter("l1TracksTPValueMapInputTags"))), - //outputTree_(fs_->make("l1VertexReco", "L1 vertex-related info")), printResults_(iConfig.getParameter("printResults")), settings_(iConfig) { const std::vector trackBranchNames( @@ -264,7 +263,7 @@ namespace l1tVertexFinder { usesResource(TFileService::kSharedResource); edm::Service fs; - outputTree_=fs->make("l1VertexReco", "L1 vertex-related info"); + outputTree_ = fs->make("l1VertexReco", "L1 vertex-related info"); std::vector::const_iterator branchNameIt = emulationVertexBranchNames.begin(); std::vector::const_iterator inputTagIt = emulationVertexInputTags.begin(); From 6b26f6b5dd932d41ece3916042986cd80bc54667 Mon Sep 17 00:00:00 2001 From: ccaillol Date: Fri, 1 Apr 2022 09:20:50 +0200 Subject: [PATCH 4/5] implement perrotta's comments --- L1Trigger/VertexFinder/plugins/VertexNTupler.cc | 2 -- L1Trigger/VertexFinder/src/VertexFinder.cc | 11 +++++------ 2 files changed, 5 insertions(+), 8 deletions(-) diff --git a/L1Trigger/VertexFinder/plugins/VertexNTupler.cc b/L1Trigger/VertexFinder/plugins/VertexNTupler.cc index 13c535d0c9180..aa9f6ba047b08 100644 --- a/L1Trigger/VertexFinder/plugins/VertexNTupler.cc +++ b/L1Trigger/VertexFinder/plugins/VertexNTupler.cc @@ -443,8 +443,6 @@ namespace l1tVertexFinder { } } - //const Vertex& TruePrimaryVertex = inputData.getPrimaryVertex(); - // create a map for associating fat reco tracks with their underlying // TTTrack pointers std::map>, const L1TrackTruthMatched*>> diff --git a/L1Trigger/VertexFinder/src/VertexFinder.cc b/L1Trigger/VertexFinder/src/VertexFinder.cc index 3c6fe46c7bb24..e011b17d65bfc 100644 --- a/L1Trigger/VertexFinder/src/VertexFinder.cc +++ b/L1Trigger/VertexFinder/src/VertexFinder.cc @@ -471,6 +471,11 @@ namespace l1tVertexFinder { std::string color) { int tableSize = data.size(); + 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 { @@ -481,12 +486,6 @@ namespace l1tVertexFinder { maximum = float(*std::max_element(std::begin(data), std::end(data))) * 1.05; } - if (maximum <= minimum) { - float average = (minimum + maximum) / 2.0; - minimum = average - 0.5; - maximum = average + 0.5; - } - std::vector intervals(tableSize, ""); std::vector values(tableSize, ""); char buffer[128]; From 2940b83e7ba48e8fcc8aab278ff79dbafb1c87f6 Mon Sep 17 00:00:00 2001 From: ccaillol Date: Fri, 1 Apr 2022 10:09:14 +0200 Subject: [PATCH 5/5] implement perrotta's comment --- L1Trigger/VertexFinder/src/VertexFinder.cc | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/L1Trigger/VertexFinder/src/VertexFinder.cc b/L1Trigger/VertexFinder/src/VertexFinder.cc index e011b17d65bfc..cedea48ea0bcd 100644 --- a/L1Trigger/VertexFinder/src/VertexFinder.cc +++ b/L1Trigger/VertexFinder/src/VertexFinder.cc @@ -471,7 +471,9 @@ namespace l1tVertexFinder { std::string color) { int tableSize = data.size(); - if (maximum <= minimum) { + 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))); } @@ -482,10 +484,6 @@ namespace l1tVertexFinder { minimum = 0; } - if (maximum == -1) { - maximum = float(*std::max_element(std::begin(data), std::end(data))) * 1.05; - } - std::vector intervals(tableSize, ""); std::vector values(tableSize, ""); char buffer[128];