Correlator PR corresponding to l1t-offline PR 1063

DataFormats/L1TParticleFlow/interface/PFJet.h

@@ -3,14 +3,15 @@
 
 #include <vector>
 #include "DataFormats/L1Trigger/interface/L1Candidate.h"
+#include "DataFormats/L1TParticleFlow/interface/PFCandidate.h"
 #include "DataFormats/Common/interface/Ptr.h"
 
 namespace l1t {
 
   class PFJet : public L1Candidate {
   public:
     /// constituent information. note that this is not going to be available in the hardware!
-    typedef std::vector<edm::Ptr<l1t::L1Candidate>> Constituents;
+    typedef std::vector<edm::Ptr<l1t::PFCandidate>> Constituents;
 
     PFJet() {}
     PFJet(float pt, float eta, float phi, float mass = 0, int hwpt = 0, int hweta = 0, int hwphi = 0)
@@ -28,13 +29,13 @@ namespace l1t {
     /// constituent information. note that this is not going to be available in the hardware!
     const Constituents& constituents() const { return constituents_; }
     /// adds a candidate to this cluster; note that this only records the information, it's up to you to also set the 4-vector appropriately
-    void addConstituent(const edm::Ptr<l1t::L1Candidate>& cand) { constituents_.emplace_back(cand); }
+    void addConstituent(const edm::Ptr<l1t::PFCandidate>& cand) { constituents_.emplace_back(cand); }
 
     // candidate interface
     size_t numberOfDaughters() const override { return constituents_.size(); }
     const reco::Candidate* daughter(size_type i) const override { return constituents_[i].get(); }
     using reco::LeafCandidate::daughter;  // avoid hiding the base
-    edm::Ptr<l1t::L1Candidate> daughterPtr(size_type i) const { return constituents_[i]; }
+    edm::Ptr<l1t::PFCandidate> daughterPtr(size_type i) const { return constituents_[i]; }
 
     // Get and set the encodedJet_ bits. The Jet is encoded in 128 bits as a 2-element array of uint64_t
     std::array<uint64_t, 2> encodedJet() { return encodedJet_; }

DataFormats/L1TParticleFlow/interface/layer1_emulator.h

@@ -175,18 +175,20 @@ namespace l1ct {
       hwIsoVars[1] = 0;
       hwIsoVars[2] = 0;
       hwIsoVars[3] = 0;
+      hwIsoVars[4] = 0;
+      hwIsoVars[5] = 0;
     }
 
     using EGIsoObj::floatIso;
 
-    enum IsoType { TkIso = 0, PfIso = 1, TkIsoPV = 2, PfIsoPV = 3 };
+    enum IsoType { TkIso = 0, PfIso = 1, TkIsoPV = 2, PfIsoPV = 3, PuppiIso = 4, PuppiIsoPV = 5 };
 
     float floatIso(IsoType type) const { return Scales::floatIso(hwIsoVars[type]); }
     float floatRelIso(IsoType type) const { return Scales::floatIso(hwIsoVars[type]) / floatPt(); }
     float hwIsoVar(IsoType type) const { return hwIsoVars[type]; }
     void setHwIso(IsoType type, iso_t value) { hwIsoVars[type] = value; }
 
-    iso_t hwIsoVars[4];
+    iso_t hwIsoVars[6];
   };
 
   struct EGIsoEleObjEmu : public EGIsoEleObj {
@@ -207,18 +209,19 @@ namespace l1ct {
     void clearIsoVars() {
       hwIsoVars[0] = 0;
       hwIsoVars[1] = 0;
+      hwIsoVars[2] = 0;
     }
 
     using EGIsoEleObj::floatIso;
 
-    enum IsoType { TkIso = 0, PfIso = 1 };
+    enum IsoType { TkIso = 0, PfIso = 1, PuppiIso = 2 };
 
     float floatIso(IsoType type) const { return Scales::floatIso(hwIsoVars[type]); }
     float floatRelIso(IsoType type) const { return Scales::floatIso(hwIsoVars[type]) / floatPt(); }
     float hwIsoVar(IsoType type) const { return hwIsoVars[type]; }
     void setHwIso(IsoType type, iso_t value) { hwIsoVars[type] = value; }
 
-    iso_t hwIsoVars[2];
+    iso_t hwIsoVars[3];
   };
 
   struct PVObjEmu : public PVObj {

DataFormats/L1TParticleFlow/src/classes_def.xml

@@ -33,7 +33,8 @@
   <class name="l1t::RegionalOutput<l1t::PFCandidateCollection>" />
   <class name="edm::Wrapper<l1t::RegionalOutput<l1t::PFCandidateCollection>>" />
 
-  <class name="l1t::PFJet"  ClassVersion="4">
+  <class name="l1t::PFJet"  ClassVersion="5">
+        <version ClassVersion="5" checksum="2270932343"/>
         <version ClassVersion="4" checksum="1424452548"/>
         <version ClassVersion="3" checksum="133342988"/>
   </class>

L1Trigger/Configuration/python/L1Trigger_EventContent_cff.py

@@ -199,6 +199,8 @@ def _appendPhase2Digis(obj):
         'keep *_l1tPFTracksFromL1TracksHGCal_*_*',
         'keep *_l1tSCPFL1PuppiCorrectedEmulator_*_*',
         'keep *_l1tSCPFL1PuppiCorrectedEmulatorMHT_*_*',
+        'keep *_l1tSCPFL1PuppiExtendedCorrectedEmulator_*_*',
+        'keep *_l1tSCPFL1PuppiExtendedCorrectedEmulatorMHT_*_*',
         'keep *_l1tPhase1JetProducer9x9_*_*', 
         'keep *_l1tPhase1JetCalibrator9x9_*_*',
         'keep *_l1tPhase1JetSumsProducer9x9_*_*',
@@ -210,13 +212,18 @@ def _appendPhase2Digis(obj):
         'keep *_l1tLayer1HGCalNoTK_*_*',
         'keep *_l1tLayer1HF_*_*',
         'keep *_l1tLayer1_*_*',
+        'keep *_l1tLayer1BarrelExtended_*_*',
+        'keep *_l1tLayer1HGCalExtended_*_*',
+        'keep *_l1tLayer1Extended_*_*',
         'keep *_l1tLayer1EG_*_*',
         'keep *_l1tLayer2EG_*_*',
         'keep *_l1tMETPFProducer_*_*',
         'keep *_l1tNNTauProducer_*_*',
         'keep *_l1tNNTauProducerPuppi_*_*',
         'keep *_l1tHPSPFTauProducerPF_*_*',
         'keep *_l1tHPSPFTauProducerPuppi_*_*',
+        'keep *_l1tBJetProducerPuppi_*_*',
+        'keep *_l1tBJetProducerPuppiCorrectedEmulator_*_*',
         'keep *_TTStubsFromPhase2TrackerDigis_*_*',
         'keep *_TTClustersFromPhase2TrackerDigis_*_*',
         'keep *_l1tTTTracksFromExtendedTrackletEmulation_*_*',

L1Trigger/Configuration/python/SimL1Emulator_cff.py

@@ -218,6 +218,13 @@
 from L1Trigger.Phase2L1ParticleFlow.L1NNTauProducer_cff import *
 _phase2_siml1emulator.add(l1tNNTauProducerPuppi)
 
+
+# BJets
+# ########################################################################
+from L1Trigger.Phase2L1ParticleFlow.L1BJetProducer_cff import *
+_phase2_siml1emulator.add(L1TBJetsTask)
+
+
 # --> add modules
 from Configuration.Eras.Modifier_phase2_trigger_cff import phase2_trigger
 phase2_trigger.toReplaceWith( SimL1EmulatorTask , _phase2_siml1emulator)

L1Trigger/Phase2L1ParticleFlow/interface/BJetId.h

@@ -0,0 +1,42 @@
+#ifndef L1TRIGGER_PHASE2L1PARTICLEFLOWS_BJETID_H
+#define L1TRIGGER_PHASE2L1PARTICLEFLOWS_BJETID_H
+
+#include <string>
+#include "PhysicsTools/TensorFlow/interface/TensorFlow.h"
+#include "DataFormats/L1TParticleFlow/interface/PFCandidate.h"
+#include "DataFormats/L1TParticleFlow/interface/PFJet.h"
+
+struct BJetTFCache {
+  BJetTFCache(const std::string &graphPath) : graphDef(tensorflow::loadGraphDef(graphPath)) {
+    session = tensorflow::createSession(graphDef.get());
+  }
+  ~BJetTFCache() { tensorflow::closeSession(session); }
+  std::unique_ptr<tensorflow::GraphDef> graphDef;
+  tensorflow::Session *session;
+};
+
+class BJetId {
+public:
+  BJetId(const std::string &iInput, const std::string &iOutput, const BJetTFCache *cache, int iNParticles);
+  ~BJetId();
+
+  void setNNVectorVar();
+  float EvaluateNN();
+  float compute(const l1t::PFJet &iJet, float vz, bool useRawPt);
+
+private:
+  std::vector<float> NNvectorVar_;
+  std::string fInput_;
+  std::string fOutput_;
+  int fNParticles_;
+  unique_ptr<float[]> fPt_;
+  unique_ptr<float[]> fEta_;
+  unique_ptr<float[]> fPhi_;
+  unique_ptr<float[]> fId_;
+  unique_ptr<int[]> fCharge_;
+  unique_ptr<float[]> fDZ_;
+  unique_ptr<float[]> fDX_;
+  unique_ptr<float[]> fDY_;
+  tensorflow::Session *sessionRef_;
+};
+#endif

L1Trigger/Phase2L1ParticleFlow/interface/CaloClusterer.h

@@ -101,6 +101,9 @@ namespace l1tpf_calo {
     T &operator()(float eta, float phi) { return data_[grid_->find_cell(eta, phi)]; }
     const T &operator()(float eta, float phi) const { return data_[grid_->find_cell(eta, phi)]; }
 
+    float eta(float eta, float phi) const { return grid().eta(grid_->find_cell(eta, phi)); }
+    float phi(float eta, float phi) const { return grid().phi(grid_->find_cell(eta, phi)); }
+
     const Grid &grid() const { return *grid_; }
 
     unsigned int size() const { return data_.size(); }
@@ -146,6 +149,7 @@ namespace l1tpf_calo {
     std::vector<T> data_;
     const T empty_;
   };
+  typedef GridData<float> EtaPhiCenterGrid;
   typedef GridData<float> EtGrid;
   typedef GridData<int> IndexGrid;
 
@@ -169,9 +173,11 @@ namespace l1tpf_calo {
 
   struct CombinedCluster : public Cluster {
     float ecal_et, hcal_et;
+    float ecal_eta, ecal_phi;
     void clear() {
       Cluster::clear();
       ecal_et = hcal_et = 0;
+      ecal_eta = ecal_phi = 0;
     }
   };
 
@@ -183,14 +189,25 @@ namespace l1tpf_calo {
     ~SingleCaloClusterer();
     void clear();
     void add(const reco::Candidate &c) { add(c.pt(), c.eta(), c.phi()); }
-    void add(float pt, float eta, float phi) { rawet_(eta, phi) += pt; }
+    void add(float pt, float eta, float phi) {
+      rawet_(eta, phi) += pt;
+      if (preciseEtaPhi_) {
+        float newet = rawet_(eta, phi);
+        float prevw = (newet - pt) / newet;
+        float nextw = pt / newet;
+        eta_center_(eta, phi) = eta_center_(eta, phi) * prevw + eta * nextw;
+        phi_center_(eta, phi) = phi_center_(eta, phi) * prevw + phi * nextw;
+      }
+    }
     void run();
 
     /// possibly grow clusters by adding unclustered energy on the sides
     //  note: there can be some double-counting as the same unclustered energy can go into more clusters
     void grow();
 
     const EtGrid &raw() const { return rawet_; }
+    const EtaPhiCenterGrid &etaCenter() const { return eta_center_; }
+    const EtaPhiCenterGrid &phiCenter() const { return phi_center_; }
     const IndexGrid &indexGrid() const { return clusterIndex_; }
     const std::vector<Cluster> &clusters() const { return clusters_; }
     const Cluster &cluster(int i) const {
@@ -199,6 +216,8 @@ namespace l1tpf_calo {
 
     /// non-const access to the energy: be careful to use it only before 'run()'
     EtGrid &raw() { return rawet_; }
+    EtaPhiCenterGrid &etaCenter() { return eta_center_; }
+    EtaPhiCenterGrid &phiCenter() { return phi_center_; }
 
     // for the moment, generic interface that takes a cluster and returns the corrected pt
     template <typename Corrector>
@@ -223,13 +242,20 @@ namespace l1tpf_calo {
     }; /* if there's more than one local maximum neighbour, they all take half of the value (no fp division) */
     const Grid *grid_;
     EtGrid rawet_, unclustered_;
+    EtaPhiCenterGrid eta_center_;
+    EtaPhiCenterGrid phi_center_;
     PreClusterGrid precluster_;
     IndexGrid clusterIndex_, cellKey_;
+    bool preciseEtaPhi_;
+    std::vector<double> etaBounds_;
+    std::vector<double> phiBounds_;
+    std::vector<unsigned int> maxClustersEtaPhi_;  //eta x phi
     std::vector<Cluster> clusters_;
     const Cluster nullCluster_;
     float zsEt_, seedEt_, minClusterEt_, minEtToGrow_;
     EnergyShareAlgo energyShareAlgo_;
     bool energyWeightedPosition_;  // do the energy-weighted cluster position instead of the cell center
+    std::vector<int> neighborCells_;
   };
 
   class SimpleCaloLinkerBase {
@@ -262,6 +288,9 @@ namespace l1tpf_calo {
     const SingleCaloClusterer &ecal_, &hcal_;
     IndexGrid clusterIndex_;
     std::vector<CombinedCluster> clusters_;
+    std::vector<double> etaBounds_;
+    std::vector<double> phiBounds_;
+    std::vector<unsigned int> maxClustersEtaPhi_;  //eta x phi
     float hoeCut_, minPhotonEt_, minHadronRawEt_, minHadronEt_;
     bool noEmInHGC_;
   };
@@ -287,6 +316,31 @@ namespace l1tpf_calo {
     SingleCaloClusterer combClusterer_;
   };
 
+  class CombinedCaloLinker : public SimpleCaloLinkerBase {
+  public:
+    CombinedCaloLinker(const edm::ParameterSet &pset, const SingleCaloClusterer &ecal, const SingleCaloClusterer &hcal);
+    ~CombinedCaloLinker() override;
+    void clear() override;
+    void run() override;
+
+  protected:
+    SingleCaloClusterer combClusterer_;
+  };
+
+  class GridSelector {
+  public:
+    GridSelector(std::vector<double> etaBounds, std::vector<double> phiBounds, std::vector<unsigned int> maxClusters);
+    ~GridSelector() {}
+    void fill(float pt, float eta, float phi, unsigned int index);
+    std::vector<unsigned int> returnSorted();
+
+  private:
+    const std::vector<double> etaBounds_;
+    const std::vector<double> phiBounds_;
+    const std::vector<unsigned int> maxClustersEtaPhi_;
+    std::vector<std::vector<std::pair<float, unsigned int>>> regionPtIndices_;  //pt and index pairs in each region
+  };
+
   // makes a calo linker (pointer will be owned by the callee)
   std::unique_ptr<SimpleCaloLinkerBase> makeCaloLinker(const edm::ParameterSet &pset,
                                                        const SingleCaloClusterer &ecal,