Integrate the FastHisto algorithm from L1TkPrimaryVertexProducer into…

… the VertexFinder. This commit also changes the std::cout statements to edm::MessageLogger statements among other cleanup actions.

L1Trigger/VertexFinder/interface/AlgoSettings.h

@@ -11,6 +11,7 @@ namespace l1tVertexFinder {
 
   enum class Algorithm {
     FastHisto,
+    FastHistoLooseAssociation,
     GapClustering,
     AgglomerativeHierarchical,
     DBSCAN,
@@ -37,22 +38,34 @@ namespace l1tVertexFinder {
     // Minimum number of tracks to accept vertex
     unsigned int vx_minTracks() const { return vx_minTracks_; }
     // Compute the z0 position of the vertex with a mean weighted with track momenta
-    bool vx_weightedmean() const { return vx_weightedmean_; }
+    unsigned int vx_weightedmean() const { return vx_weightedmean_; }
     /// Chi2 cut for the Adaptive Vertex Recostruction Algorithm
     float vx_chi2cut() const { return vx_chi2cut_; }
+    /// Window size of the sliding window
+    unsigned int vx_windowSize() const { return vx_windowSize_; }
     /// FastHisto histogram parameters (min, max, width)
     std::vector<double> 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
     float vx_width() const { return vx_width_; }
+    /// 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
+    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_; }
 
     unsigned int vx_kmeans_iterations() const { return vx_kmeans_iterations_; }
     unsigned int vx_kmeans_nclusters() const { return vx_kmeans_nclusters_; }
     float vx_TrackMinPt() const { return vx_TrackMinPt_; }
+    float vx_TrackMaxPt() const { return vx_TrackMaxPt_; }
+    float vx_TrackMaxPtBehavior() const { return vx_TrackMaxPtBehavior_; }
+    float vx_TrackMaxChi2() const { return vx_TrackMaxChi2_; }
+    unsigned int vx_NStubMin() const { return vx_NStubMin_; }
+    unsigned int vx_NStubPSMin() const { return vx_NStubPSMin_; }
 
     //=== Debug printout
     unsigned int debug() const { return debug_; }
@@ -79,11 +92,20 @@ namespace l1tVertexFinder {
     float vx_resolution_;
     unsigned int vx_distanceType_;
     unsigned int vx_minTracks_;
-    bool vx_weightedmean_;
+    unsigned int vx_weightedmean_;
     float vx_chi2cut_;
+    bool vx_DoPtComp_;
+    bool vx_DoTightChi2_;
     std::vector<double> vx_histogram_parameters_;
+    unsigned int vx_nvtx_;
     float vx_width_;
+    unsigned int vx_windowSize_;
     float vx_TrackMinPt_;
+    float vx_TrackMaxPt_;
+    int vx_TrackMaxPtBehavior_;
+    float vx_TrackMaxChi2_;
+    unsigned int vx_NStubMin_;
+    unsigned int vx_NStubPSMin_;
     float vx_dbscan_pt_;
     float vx_dbscan_mintracks_;
     unsigned int vx_kmeans_iterations_;

L1Trigger/VertexFinder/interface/InputData.h

@@ -3,6 +3,7 @@
 
 #include "FWCore/Framework/interface/EDAnalyzer.h"
 #include "FWCore/Framework/interface/Frameworkfwd.h"
+#include "FWCore/MessageLogger/interface/MessageLogger.h"
 #include "FWCore/Utilities/interface/InputTag.h"
 #include "L1Trigger/VertexFinder/interface/Stub.h"
 #include "L1Trigger/VertexFinder/interface/TP.h"
@@ -23,10 +24,12 @@ namespace l1tVertexFinder {
     InputData(const edm::Event& iEvent,
               const edm::EventSetup& iSetup,
               const AnalysisSettings& settings,
-              const edm::EDGetTokenT<TrackingParticleCollection> tpInputTag,
-              const edm::EDGetTokenT<DetSetVec> stubInputTag,
-              const edm::EDGetTokenT<TTStubAssMap> stubTruthInputTag,
-              const edm::EDGetTokenT<TTClusterAssMap> clusterTruthInputTag);
+              const edm::EDGetTokenT<edm::HepMCProduct> hepMCToken,
+              const edm::EDGetTokenT<edm::View<reco::GenParticle>> genParticlesToken,
+              const edm::EDGetTokenT<TrackingParticleCollection> tpToken,
+              const edm::EDGetTokenT<DetSetVec> stubToken,
+              const edm::EDGetTokenT<TTStubAssMap> stubTruthToken,
+              const edm::EDGetTokenT<TTClusterAssMap> clusterTruthToken);
 
     // Sort Tracking Particles by vertex z position
     struct SortVertexByZ0 {
@@ -35,9 +38,13 @@ namespace l1tVertexFinder {
 
     // Get tracking particles
     const std::vector<TP>& getTPs() const { return vTPs_; }
-    /// Get Primary vertex information
+    /// Get primary vertex information (vertex from HepMCProduct)
+    const Vertex& getHepMCVertex() const { return hepMCVertex_; }
+    /// Get primary vertex information (vertex from gen particles)
+    const Vertex& getGenVertex() const { return genVertex_; } 
+    /// Get primary vertex information (vertex from tracking particles)
     const Vertex& getPrimaryVertex() const { return vertex_; }
-    /// Get PileUp Vertices information
+    /// Get pile-up vertices information
     const std::vector<Vertex>& getPileUpVertices() const { return vertices_; }
     /// Get reconstructable pile-up vertices information
     const std::vector<Vertex>& getRecoPileUpVertices() const { return recoVertices_; }
@@ -51,24 +58,13 @@ namespace l1tVertexFinder {
     const float genPt() const { return genPt_; }
     const float genPt_PU() const { return genPt_PU_; }
 
-  private:
-    // const edm::EDGetTokenT<TrackingParticleCollection> inputTag;
-
-    //  // Can optionally be used to sort stubs by bend.
-    //  struct SortStubsInBend {
-    //     inline bool operator() (const Stub* stub1, const Stub* stub2) {
-    //        return(std::abs(stub1->bend()) < std::abs(stub2->bend()));
-    //     }
-    //  };
-
   private:
     std::vector<TP> vTPs_;  // tracking particles
-    Vertex vertex_;
-    std::vector<Vertex> vertices_;
-    std::vector<Vertex> recoVertices_;
-
-    //--- of minor importance ...
-
+    Vertex hepMCVertex_;  // primary vertex from the edm::HepMCProduct
+    Vertex genVertex_;  // primary vertex formed from the gen particles collection
+    Vertex vertex_;  // primary vertex formed from the tracking particles
+    std::vector<Vertex> vertices_;  // pile-up vertices
+    std::vector<Vertex> recoVertices_;  // reconstructable pile-up vertices
     std::vector<Stub>
         vAllStubs_;  // all stubs, even those that would fail any tightened front-end readout electronic cuts specified in section StubCuts of Analyze_Defaults_cfi.py. (Only used to measure
                      // the efficiency of these cuts).

L1Trigger/VertexFinder/interface/RecoVertex.h

@@ -11,66 +11,60 @@ namespace l1tVertexFinder {
 
   class RecoVertex {
   public:
-    // Fill useful info about tracking particle.
+    /// Constructor and destructor
     RecoVertex() {
       z0_ = -999.;
-      pT_ = -9999.;
-      met_ = -999.;
+      z0width_ = -999.;
+      pT_ = -999.;
+      highestPt_ = -999.;
+      pv_ = false;
+      highPt_ = false;
+      numHighPtTracks_ = 0;
+      tracks_.clear();
     }
     ~RecoVertex() {}
 
-    /// Tracking Particles in vertex
-    const std::vector<const L1Track*>& tracks() const { return tracks_; }
-    /// Tracking Particles in vertex
-    const std::set<const TP*>& trueTracks() const { return trueTracks_; }
-    /// Number of tracks originating from this vertex
-    unsigned int numTracks() const { return tracks_.size(); }
-    /// Number of true particles assigned to this vertex
-    unsigned int numTrueTracks() const { return trueTracks_.size(); }
-    /// Assign fitted track to this vertex
-    void insert(const L1Track* fitTrack) {
-      tracks_.push_back(
-          fitTrack); /* if(fitTrack->getMatchedTP()!= nullptr and fitTrack->getMatchedTP()->physicsCollision()) trueTracks_.insert(fitTrack->getMatchedTP()); */
+    /// Operators
+    RecoVertex& operator+=(const RecoVertex& rhs){
+      this->tracks_.insert(std::end(this->tracks_), std::begin(rhs.tracks()), std::end(rhs.tracks()));
+      return *this;
     }
-    /// Compute vertex parameters
-    void computeParameters(bool weightedmean = false);
-    /// Set z0 position
-    void setZ(double z) { z0_ = z; }
-    /// Sum ot fitted tracks transverse momentum
-    double pT() const { return pT_; }
-    /// Vertex z0 position
-    double z0() const { return z0_; }
-    /// Vertex z0 width
-    double z0width() const { return z0width_; }
+
     /// Clear track vector
-    void clear() {
-      tracks_.clear();
-      trueTracks_.clear();
-    }
-    /// True if primary vertex
-    bool primaryVertex() const { return pv_; }
-    /// Set primary vertex tag
-    void isPrimary(bool is) { pv_ = is; }
+    void clear() { tracks_.clear(); }
+    /// Compute vertex parameters
+    void computeParameters(unsigned int weightedmean = false, double highPtThreshold = 50., int highPtBehavior = -1);
     /// Contain high-pT track?
     bool hasHighPt() const { return highPt_; }
     /// highest track pT in the vertex
     double highestPt() const { return highestPt_; }
+    /// Assign fitted track to this vertex
+    void insert(const L1Track* fitTrack) { tracks_.push_back(fitTrack); }
+    /// Set primary vertex tag
+    void isPrimary(bool is) { pv_ = is; }
     /// Number of high-pT tracks (pT > 10 GeV)
     unsigned int numHighPtTracks() const { return numHighPtTracks_; }
-    /// Vertec MET
-    double met() const { return met_; }
+    /// Number of tracks originating from this vertex
+    unsigned int numTracks() const { return tracks_.size(); }
+    /// True if primary vertex
+    bool primaryVertex() const { return pv_; }
+    /// Sum of fitted tracks transverse momentum [GeV]
+    double pT() const { return pT_; }
+    /// Tracking Particles in vertex
+    const std::vector<const L1Track*>& tracks() const { return tracks_; }
+    /// Set z0 position [cm]
+    void setZ(double z) { z0_ = z; }
+    /// Vertex z0 position [cm]
+    double z0() const { return z0_; }
+    /// Vertex z0 width [cm]
+    double z0width() const { return z0width_; }
 
   private:
     double z0_;
     double z0width_;
     double pT_;
-    double met_;
-    double metX_;
-    double metY_;
     double highestPt_;
-
     std::vector<const L1Track*> tracks_;
-    std::set<const TP*> trueTracks_;
     bool pv_;
     bool highPt_;
     unsigned int numHighPtTracks_;

L1Trigger/VertexFinder/interface/Stub.h

@@ -1,19 +1,17 @@
 #ifndef __L1Trigger_VertexFinder_Stub_h__
 #define __L1Trigger_VertexFinder_Stub_h__
 
-#include <array>
-#include <map>
-#include <set>
-
 #include "DataFormats/L1TrackTrigger/interface/TTStub.h"
 #include "DataFormats/L1TrackTrigger/interface/TTTypes.h"
 #include "DataFormats/Phase2TrackerDigi/interface/Phase2TrackerDigi.h"
-
-// TTStubAssociationMap.h forgets to two needed files, so must include them here ...
 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
 #include "SimTracker/TrackTriggerAssociation/interface/TTClusterAssociationMap.h"
 #include "SimTracker/TrackTriggerAssociation/interface/TTStubAssociationMap.h"
 
+#include <array>
+#include <map>
+#include <set>
+
 class TrackerGeometry;
 
 namespace l1tVertexFinder {

L1Trigger/VertexFinder/interface/Vertex.h

@@ -10,9 +10,15 @@ namespace l1tVertexFinder {
   class Vertex {
   public:
     // Fill useful info about tracking particle.
-    Vertex() {}
+    Vertex() {
+        Vertex(-999.);
+    }
 
-    Vertex(double vz) : vz_(vz) {}
+    Vertex(double vz) : vz_(vz) {
+        z0_ = -999.;
+        z0width_ = -999.;
+        pT_ = -999.;
+    }
 
     ~Vertex() {}
 
@@ -24,25 +30,22 @@ namespace l1tVertexFinder {
     void insert(TP tp) { tracks_.push_back(tp); }
     /// Compute vertex parameters
     void computeParameters();
-    /// Sum ot fitted tracks transverse momentum
+    /// Sum of fitted tracks transverse momentum [GeV]
     double pT() const { return pT_; }
-    /// Vertex z0 position
+    /// Vertex z0 position [cm]
     double z0() const { return z0_; }
-    /// Vertex z0 width
+    /// Vertex z0 width [cm]
     double z0width() const { return z0width_; }
-    /// Vertex z position
+    /// Vertex z position [cm]
     double vz() const { return vz_; }
-    /// Vertec MET
-    double met() const { return met_; }
+    /// Reset/initialize all of the member data
+    void reset();
 
   private:
     double vz_;
     double z0_;
     double z0width_;
     double pT_;
-    double met_;
-    double metX_;
-    double metY_;
 
     std::vector<TP> tracks_;
   };

L1Trigger/VertexFinder/interface/VertexFinder.h

@@ -3,9 +3,14 @@
 
 #include "DataFormats/Common/interface/Ptr.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"
 
+#include <algorithm>
+#include <iterator>
 #include <vector>
 
 namespace l1tVertexFinder {
@@ -47,7 +52,7 @@ namespace l1tVertexFinder {
       if (pv_index_ < vertices_.size())
         return vertices_[pv_index_];
       else {
-        std::cout << "No Primary Vertex has been found." << std::endl;
+        edm::LogWarning("VertexFinder")<<"PrimaryVertex::No Primary Vertex has been found.";
         return RecoVertex();
       }
     }
@@ -84,15 +89,25 @@ namespace l1tVertexFinder {
     void HPV();
     /// Kmeans Algorithm
     void Kmeans();
-    /// Histogramming algorithmn
-    void FastHisto();
+    /// TDR histogramming algorithmn
+    void FastHistoLooseAssociation();
+    /// Histogramming algorithm
+    void FastHisto(const TrackerTopology* tTopo);
     /// Sort Vertices in z
     void SortVerticesInZ0() { std::sort(vertices_.begin(), vertices_.end(), SortVertexByZ0()); }
     /// Number of iterations
     unsigned int NumIterations() const { return iterations_; }
     /// Number of iterations
     unsigned int IterationsPerTrack() const { return double(iterations_) / double(fitTracks_.size()); }
 
+    template<typename ForwardIterator, typename T>
+    void strided_iota(ForwardIterator first, ForwardIterator last, T value, T stride) {
+        while(first != last) {
+            *first++ = value;
+            value += stride;
+        }
+    }
+
   private:
     const AlgoSettings* settings_;
     std::vector<RecoVertex> vertices_;