Add Giovanna's PFA as an option for vertex simulation

linacre · linacre · commit a68f85714f33 · 2025-01-17T18:05:32.000Z
diff --git a/L1Trigger/VertexFinder/interface/AlgoSettings.h b/L1Trigger/VertexFinder/interface/AlgoSettings.h
@@ -11,6 +11,8 @@
 namespace l1tVertexFinder {
 
   enum class Algorithm {
+    PFA,
+    PFASingleVertex,
     fastHisto,
     fastHistoEmulation,
     fastHistoLooseAssociation,
@@ -49,6 +51,25 @@ namespace l1tVertexFinder {
     float vx_chi2cut() const { return vx_chi2cut_; }
     // Do track quality cuts in emulation algorithms
     bool vx_DoQualityCuts() const { return vx_DoQualityCuts_; }
+    // PFA scan parameters (min, max, width)
+    std::vector<double> vx_pfa_scanparameters() const { return vx_pfa_scanparameters_; }
+    double vx_pfa_min() const { return vx_pfa_scanparameters_.at(0); }
+    double vx_pfa_max() const { return vx_pfa_scanparameters_.at(1); }
+    double vx_pfa_binwidth() const { return vx_pfa_scanparameters_.at(2); }
+    // Include eta-dependence of the estimated track resolution used in PFA
+    bool vx_pfa_etadependentresolution() const { return vx_pfa_etadependentresolution_; }
+    // Scale factor for the PFA track resolution parameter (where the nominal values with and without eta-dependence are hard-coded using the fit results from Giovanna's thesis)
+    double vx_pfa_resolutionSF() const { return vx_pfa_resolutionSF_; }
+    // PFA Gaussian width cutoff
+    double vx_pfa_width() const { return vx_pfa_width_; }
+    // Enable 2-step process where the weighted pT sum is only calculated at positions where the weighted multiplicity is maximum ("local maxima"). In the second step, the local maximum with the largest weighted pT sum is chosen as the vertex.
+    bool vx_pfa_usemultiplicitymaxima() const { return vx_pfa_usemultiplicitymaxima_; }
+    // Weight function to use in PFA. 0: Gaussian, 1: Gaussian without width normalisation, 2: Complementary error function
+    unsigned int vx_pfa_weightfunction() const { return vx_pfa_weightfunction_; }
+    // Instead of taking the z0 value from the discrete PFA scan (0), calculate it from the Gaussian and pT-weighted sum of track z0 (1) or the optimal (1/variance) weighted mean of associated tracks, weighted also by pT and association probability (2)
+    unsigned int vx_pfa_weightedz0() const { return vx_pfa_weightedz0_; }
+    // Use vx_TrackMinPt cut specified below (otherwise no additional track selection is applied)
+    bool vx_pfa_doqualitycuts() const { return vx_pfa_doqualitycuts_; }
     // Window size of the sliding window
     unsigned int vx_windowSize() const { return vx_windowSize_; }
     // fastHisto histogram parameters (min, max, width)
@@ -112,6 +133,14 @@ namespace l1tVertexFinder {
     unsigned int vx_weightedmean_;
     float vx_chi2cut_;
     bool vx_DoQualityCuts_;
+    std::vector<double> vx_pfa_scanparameters_;
+    bool vx_pfa_etadependentresolution_;
+    double vx_pfa_resolutionSF_;
+    double vx_pfa_width_;
+    bool vx_pfa_usemultiplicitymaxima_;
+    unsigned int vx_pfa_weightfunction_;
+    unsigned int vx_pfa_weightedz0_;
+    bool vx_pfa_doqualitycuts_;
     bool vx_DoPtComp_;
     bool vx_DoTightChi2_;
     std::vector<double> vx_histogram_parameters_;
diff --git a/L1Trigger/VertexFinder/interface/VertexFinder.h b/L1Trigger/VertexFinder/interface/VertexFinder.h
@@ -135,6 +135,12 @@ namespace l1tVertexFinder {
     void computeAndSetVertexParameters(RecoVertex<>& vertex,
                                        const std::vector<float>& bin_centers,
                                        const std::vector<unsigned int>& counts);
+
+    double computeAndSetVertexParametersPFA(RecoVertex<>& vertex);
+    /// Peak finding algorithm
+    void PFA();
+    /// Peak finding algorithm, single vertex
+    void PFASingleVertex();
     /// DBSCAN algorithm
     void DBSCAN();
     /// High pT Vertex Algorithm
diff --git a/L1Trigger/VertexFinder/plugins/VertexProducer.cc b/L1Trigger/VertexFinder/plugins/VertexProducer.cc
@@ -23,6 +23,12 @@ VertexProducer::VertexProducer(const edm::ParameterSet& iConfig)
   // Get configuration parameters
 
   switch (settings_.vx_algo()) {
+    case Algorithm::PFA:
+      edm::LogInfo("VertexProducer") << "VertexProducer::Finding vertices using the PFA algorithm";
+      break;
+    case Algorithm::PFASingleVertex:
+      edm::LogInfo("VertexProducer") << "VertexProducer::Finding vertices using the PFASingleVertex algorithm";
+      break;
     case Algorithm::fastHisto:
       edm::LogInfo("VertexProducer") << "VertexProducer::Finding vertices using the fastHisto binning algorithm";
       break;
@@ -110,6 +116,12 @@ void VertexProducer::produce(edm::StreamID, edm::Event& iEvent, const edm::Event
   VertexFinder vf(l1Tracks, settings_);
 
   switch (settings_.vx_algo()) {
+    case Algorithm::PFA:
+      vf.PFA();
+      break;
+    case Algorithm::PFASingleVertex:
+      vf.PFASingleVertex();
+      break;
     case Algorithm::fastHisto: {
       const TrackerTopology& tTopo = iSetup.getData(tTopoToken);
       vf.fastHisto(&tTopo);
diff --git a/L1Trigger/VertexFinder/python/l1tVertexProducer_cfi.py b/L1Trigger/VertexFinder/python/l1tVertexProducer_cfi.py
@@ -8,7 +8,7 @@
   # === Vertex Reconstruction configuration
   VertexReconstruction = cms.PSet(
         # Vertex Reconstruction Algorithm
-        Algorithm = cms.string("fastHisto"),
+        Algorithm = cms.string("PFA"),
         # Vertex distance [cm]
         VertexDistance = cms.double(.15),
         # Assumed Vertex Resolution [cm]
@@ -21,15 +21,31 @@
         #   0 = unweighted
         #   1 = pT weighted
         #   2 = pT^2 weighted
-        WeightedMean = cms.uint32(1),
+        WeightedMean = cms.uint32(2),
         # 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),
+        FH_DoPtComp = cms.bool(False),
         # chi2dof < 5 for tracks with Pt > 10
         FH_DoTightChi2 = cms.bool(False),
+        # PFA algorithm scan parameters (min,max,width) [cm]
+        PFA_ScanParameters = cms.vdouble(-20.46912512, 20.46912512, 0.03997876),
+        # Include eta-dependence of the estimated track resolution used in PFA
+        PFA_EtaDependentResolution = cms.bool(True),
+        # Scale factor for the PFA track resolution parameter (where the nominal values with and without eta-dependence are hard-coded using the fit results from Giovanna's thesis)
+        PFA_ResolutionSF = cms.double(2.),
+        # PFA Gaussian width cutoff [cm]
+        PFA_VertexWidth = cms.double(1.31), # Giovanna's recommendation of 3*sigma(lowest-resolution tracks).
+        # Enable 2-step process where the weighted pT sum is only calculated at positions where the weighted multiplicity is maximum ("local maxima"). In the second step, the local maximum with the largest weighted pT sum is chosen as the vertex.
+        PFA_UseMultiplicityMaxima = cms.bool(False),
+        # Weight function to use in PFA. 0: Gaussian, 1: Gaussian without width normalisation, 2: Complementary error function
+        PFA_WeightFunction = cms.uint32(1),
+        # Instead of taking the z0 value from the discrete PFA scan (0), calculate it from the Gaussian and pT-weighted sum of track z0 (1) or the optimal (1/variance) weighted mean of associated tracks, weighted also by pT and association probability (2)
+        PFA_WeightedZ0 = cms.uint32(1),
+        # Use VxMinTrackPt cut specified below (otherwise no additional track selection is applied)
+        PFA_DoQualityCuts = cms.bool(False),
         # fastHisto algorithm histogram parameters (min,max,width) [cm]
         # TDR settings: [-14.95, 15.0, 0.1]
         # L1TkPrimaryVertexProducer: [-30.0, 30.0, 0.09983361065]
@@ -38,7 +54,7 @@
         # Track word limits (256 binns): [-20.46912512, 20.46912512, 0.15991504]
         FH_HistogramParameters = cms.vdouble(-20.46912512, 20.46912512, 0.15991504),
         # The number of vertixes to return (i.e. N windows with the highest combined pT)
-        FH_NVtx = cms.uint32(1),
+        FH_NVtx = cms.uint32(10),
         # fastHisto algorithm assumed vertex half-width [cm]
         FH_VertexWidth = cms.double(.15),
         # Window size of the sliding window
@@ -62,11 +78,11 @@
         # Option '0' was used for the TDR, but '1' is used for the firmware
         VxMaxTrackPtBehavior = cms.int32(1),
         # Maximum chi2 of tracks used to create vertex
-        VxMaxTrackChi2 = cms.double(100.),
+        VxMaxTrackChi2 = cms.double(99999999.),
         # Minimum number of stubs associated to a track
-        VxMinNStub = cms.uint32(4),
+        VxMinNStub = cms.uint32(0),
         # Minimum number of stubs in PS modules associated to a track
-        VxMinNStubPS = cms.uint32(3),
+        VxMinNStubPS = cms.uint32(0),
         # Track weight NN graph 
         TrackWeightGraph = cms.FileInPath("L1Trigger/VertexFinder/data/NNVtx_WeightModelGraph.pb"),
         # Pattern recognition NN graph
diff --git a/L1Trigger/VertexFinder/src/AlgoSettings.cc b/L1Trigger/VertexFinder/src/AlgoSettings.cc
@@ -13,6 +13,14 @@ namespace l1tVertexFinder {
         vx_weightedmean_(vertex_.getParameter<unsigned int>("WeightedMean")),
         vx_chi2cut_(vertex_.getParameter<double>("AVR_chi2cut")),
         vx_DoQualityCuts_(vertex_.getParameter<bool>("EM_DoQualityCuts")),
+        vx_pfa_scanparameters_(vertex_.getParameter<std::vector<double> >("PFA_ScanParameters")),
+        vx_pfa_etadependentresolution_(vertex_.getParameter<bool>("PFA_EtaDependentResolution")),
+        vx_pfa_resolutionSF_(vertex_.getParameter<double>("PFA_ResolutionSF")),
+        vx_pfa_width_(vertex_.getParameter<double>("PFA_VertexWidth")),
+        vx_pfa_usemultiplicitymaxima_(vertex_.getParameter<bool>("PFA_UseMultiplicityMaxima")),
+        vx_pfa_weightfunction_(vertex_.getParameter<unsigned int>("PFA_WeightFunction")),
+        vx_pfa_weightedz0_(vertex_.getParameter<unsigned int>("PFA_WeightedZ0")),
+        vx_pfa_doqualitycuts_(vertex_.getParameter<bool>("PFA_DoQualityCuts")),
         vx_DoPtComp_(vertex_.getParameter<bool>("FH_DoPtComp")),
         vx_DoTightChi2_(vertex_.getParameter<bool>("FH_DoTightChi2")),
         vx_histogram_parameters_(vertex_.getParameter<std::vector<double> >("FH_HistogramParameters")),
@@ -57,6 +65,8 @@ namespace l1tVertexFinder {
   }
 
   const std::map<std::string, Algorithm> AlgoSettings::algoNameMap = {
+      {"PFA", Algorithm::PFA},
+      {"PFASingleVertex", Algorithm::PFASingleVertex},
       {"fastHisto", Algorithm::fastHisto},
       {"fastHistoEmulation", Algorithm::fastHistoEmulation},
       {"fastHistoLooseAssociation", Algorithm::fastHistoLooseAssociation},
@@ -70,6 +80,8 @@ namespace l1tVertexFinder {
       {"NNEmulation", Algorithm::NNEmulation}};
 
   const std::map<Algorithm, Precision> AlgoSettings::algoPrecisionMap = {
+      {Algorithm::PFA, Precision::Simulation},
+      {Algorithm::PFASingleVertex, Precision::Simulation},
       {Algorithm::fastHisto, Precision::Simulation},
       {Algorithm::fastHistoEmulation, Precision::Emulation},
       {Algorithm::fastHistoLooseAssociation, Precision::Simulation},
diff --git a/L1Trigger/VertexFinder/src/VertexFinder.cc b/L1Trigger/VertexFinder/src/VertexFinder.cc
