L1TSC82ProngJetModel

DPGAnalysis/Phase2L1TNanoAOD/python/l1tPh2Nanotables_cff.py

@@ -394,12 +394,13 @@
 )
 
 sc8JetTable = pfJetTable.clone(
-    src = 'l1tSC8PFL1PuppiCorrectedEmulator',
+    src = cms.InputTag('l1tSC82ProngJetProducer', 'l1tSC82ProngJets'),
     name = "L1puppiJetSC8",
     doc = "SeededCone 0.8 Puppi jet,  origin: Correlator",
     variables = cms.PSet(
         pfJetTable.variables.clone(),
-        mass = Var("mass", float)
+        mass = Var("mass", float),
+        nprongTagScore = Var('getTagScore("nprong")', float)
     )
 )
 
@@ -580,7 +581,7 @@
     # ## jets
     sc4JetTable,
     sc8JetTable,
-    sc4ExtJetTable, 
+    sc4ExtJetTable,
     sc4NGJetTable,
     histoJetTable,
     caloJetTable,

DataFormats/L1TParticleFlow/interface/jets.h

@@ -15,7 +15,7 @@ namespace l1ct {
   // all possible tag categories (can be extended for new / separate taggers)
   class JetTagClass {
   public:
-    enum JetTagClassValue : uint8_t { b, c, uds, g, tau_p, tau_n, mu, e };
+    enum JetTagClassValue : uint8_t { b, c, uds, g, tau_p, tau_n, mu, e, nprong };
     JetTagClass() = default;
     JetTagClass(JetTagClassValue aJetTagClassValue) : value_(aJetTagClassValue) {}
     JetTagClass(std::string aJetTagClassValueString) {
@@ -47,7 +47,7 @@ namespace l1ct {
 
   // Define a separate class/struct for jet tag handling
   struct JetTagClassHandler {
-    static const unsigned NTagFields = 8;
+    static const unsigned NTagFields = 9;
     static const JetTagClass tagClassesDefault_[NTagFields];
 
     JetTagClass tagClassesArray[NTagFields];
@@ -195,7 +195,7 @@ namespace l1ct {
       j.v3.phi = CTtoGT_phi(hwPhi);
       j.v3.eta = CTtoGT_eta(hwEta);
       j.z0(l1ct::z0_t::width - 1, 0) = hwZ0(l1ct::z0_t::width - 1, 0);
-      j.hwNProngScore = 0;
+      j.hwNProngScore = hwTagScores[0];  // currently only one WideJet tagger
       j.hwMassSq = CTtoGT_massSq(hwMassSq);
 
       return j;

DataFormats/L1TParticleFlow/src/jets.cpp

@@ -8,7 +8,8 @@ const std::unordered_map<std::string, l1ct::JetTagClass::JetTagClassValue> l1ct:
     {"tau_p", l1ct::JetTagClass::JetTagClassValue::tau_p},
     {"tau_n", l1ct::JetTagClass::JetTagClassValue::tau_n},
     {"mu", l1ct::JetTagClass::JetTagClassValue::mu},
-    {"e", l1ct::JetTagClass::JetTagClassValue::e}};
+    {"e", l1ct::JetTagClass::JetTagClassValue::e},
+    {"nprong", l1ct::JetTagClass::JetTagClassValue::nprong}};
 
 const l1ct::JetTagClass l1ct::JetTagClassHandler::tagClassesDefault_[NTagFields] = {l1ct::JetTagClass("b"),
                                                                                     l1ct::JetTagClass("c"),
@@ -17,4 +18,5 @@ const l1ct::JetTagClass l1ct::JetTagClassHandler::tagClassesDefault_[NTagFields]
                                                                                     l1ct::JetTagClass("tau_p"),
                                                                                     l1ct::JetTagClass("tau_n"),
                                                                                     l1ct::JetTagClass("mu"),
-                                                                                    l1ct::JetTagClass("e")};
+                                                                                    l1ct::JetTagClass("e"),
+                                                                                    l1ct::JetTagClass("nprong")};

L1Trigger/Phase2L1ParticleFlow/interface/L1TSC82ProngJetID.h

@@ -0,0 +1,45 @@
+#ifndef L1TRIGGER_PHASE2L1PARTICLEFLOWS_L1TSC82ProngJetID_H
+#define L1TRIGGER_PHASE2L1PARTICLEFLOWS_L1TSC82ProngJetID_H
+
+#include "DataFormats/L1TParticleFlow/interface/PFJet.h"
+#include <memory>
+#include <vector>
+
+//HLS4ML compiled emulator modeling
+#include "ap_fixed.h"
+#include "hls4ml/emulator.h"
+
+class L1TSC82ProngJetID {
+public:
+  L1TSC82ProngJetID(const std::shared_ptr<hls4mlEmulator::Model> model, int iNParticles);
+
+  typedef ap_fixed<24, 12, AP_RND, AP_SAT, 0> inputtype;
+  typedef ap_ufixed<20, 10, AP_RND, AP_SAT, 0> prong_score;
+
+  void setNNVectorVar();
+  std::vector<float> EvaluateNNFixed();
+  std::vector<float> computeFixed(const l1t::PFJet &iJet);
+
+private:
+  std::vector<inputtype> NNvectorVar_;
+  int fNParticles_;
+  std::unique_ptr<float[]> fPt_;
+  std::unique_ptr<float[]> fPt_rel_;
+  std::unique_ptr<float[]> fDEta_;
+  std::unique_ptr<float[]> fDPhi_;
+  std::unique_ptr<float[]> fPt_log_;
+  std::unique_ptr<float[]> fMass_;
+  std::unique_ptr<float[]> fZ0_;
+  std::unique_ptr<float[]> fDxy_;
+  std::unique_ptr<int[]> fIs_filled_;
+  std::unique_ptr<float[]> fPuppi_weight_;
+  std::unique_ptr<int[]> fEmID_;
+  std::unique_ptr<float[]> fQuality_;
+
+  std::unique_ptr<int[]> fCharge_;
+  std::unique_ptr<int[]> fId_;
+  std::shared_ptr<hls4mlEmulator::Model> modelRef_;
+
+  //bool isDebugEnabled_;
+};
+#endif

L1Trigger/Phase2L1ParticleFlow/plugins/L1TSC82ProngJetModelProducer.cc

@@ -0,0 +1,108 @@
+#include "FWCore/Framework/interface/Frameworkfwd.h"
+#include "FWCore/Framework/interface/stream/EDProducer.h"
+#include "FWCore/Framework/interface/Event.h"
+#include "FWCore/Framework/interface/MakerMacros.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+#include "FWCore/Utilities/interface/InputTag.h"
+
+#include "DataFormats/L1TParticleFlow/interface/PFJet.h"
+#include "DataFormats/JetReco/interface/Jet.h"
+#include "DataFormats/L1TParticleFlow/interface/PFCandidate.h"
+#include "L1Trigger/Phase2L1ParticleFlow/interface/L1TSC82ProngJetID.h"
+#include "DataFormats/Common/interface/ValueMap.h"
+
+#include "DataFormats/L1Trigger/interface/VertexWord.h"
+
+#include <cmath>
+#include <vector>
+
+#include <string>
+#include "ap_fixed.h"
+#include "hls4ml/emulator.h"
+
+using namespace l1t;
+
+class L1TSC82ProngJetProducer : public edm::stream::EDProducer<> {
+public:
+  explicit L1TSC82ProngJetProducer(const edm::ParameterSet&);
+  ~L1TSC82ProngJetProducer() override = default;
+
+  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
+
+private:
+  std::unique_ptr<L1TSC82ProngJetID> fJetId_;
+  void produce(edm::Event& iEvent, const edm::EventSetup& iSetup) override;
+
+  edm::EDGetTokenT<edm::View<l1t::PFJet>> const jets_;
+  double const fMinPt_;
+  double const fMaxEta_;
+  unsigned int const fMaxJets_;
+  int const fNParticles_;
+
+  hls4mlEmulator::ModelLoader loader;
+  std::shared_ptr<hls4mlEmulator::Model> model;
+};
+
+L1TSC82ProngJetProducer::L1TSC82ProngJetProducer(const edm::ParameterSet& cfg)
+    : jets_(consumes<edm::View<l1t::PFJet>>(cfg.getParameter<edm::InputTag>("jets"))),
+      fMinPt_(cfg.getParameter<double>("minPt")),
+      fMaxEta_(cfg.getParameter<double>("maxEta")),
+      fMaxJets_(cfg.getParameter<int>("maxJets")),
+      fNParticles_(cfg.getParameter<int>("nParticles")),
+      loader(hls4mlEmulator::ModelLoader(cfg.getParameter<std::string>("l1tSC82ProngJetModelPath"))) {
+  model = loader.load_model();
+  fJetId_ = std::make_unique<L1TSC82ProngJetID>(model, fNParticles_);
+  produces<l1t::PFJetCollection>("l1tSC82ProngJets");
+}
+
+void L1TSC82ProngJetProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
+  edm::Handle<edm::View<l1t::PFJet>> jets;
+  iEvent.getByToken(jets_, jets);
+  std::vector<l1t::PFJet> taggedJets;
+
+  for (const auto& srcjet : *jets) {
+    l1ct::Jet ctHWJet = l1ct::Jet::unpack(srcjet.encodedJet(l1t::PFJet::HWEncoding::CT));
+
+    if (srcjet.pt() < fMinPt_ || std::abs(srcjet.eta()) > fMaxEta_ || taggedJets.size() >= fMaxJets_) {
+      ctHWJet.clear();
+      continue;
+    }
+    std::vector<float> JetProngScore_float = fJetId_->computeFixed(srcjet);
+    l1gt::WideJet gtwHWJet = l1gt::WideJet::unpack(srcjet.getHWJetGTWide());
+
+    l1t::PFJet edmJet(
+        srcjet.pt(), srcjet.eta(), srcjet.phi(), srcjet.mass(), gtwHWJet.v3.pt.V, gtwHWJet.v3.eta.V, gtwHWJet.v3.phi.V);
+
+    std::vector<l1ct::JetTagClass> classes{l1ct::JetTagClass(l1ct::JetTagClass::JetTagClassValue::nprong)};
+
+    edmJet.addTagScores(JetProngScore_float, classes, 1.);
+    edmJet.setEncodedJet(l1t::PFJet::HWEncoding::CT, ctHWJet.pack());
+    edmJet.setEncodedJet(l1t::PFJet::HWEncoding::GTWide, gtwHWJet.pack());
+
+    std::vector<edm::Ptr<l1t::PFCandidate>> constituents;
+    std::for_each(srcjet.constituents().begin(), srcjet.constituents().end(), [&](auto constituent) {
+      edmJet.addConstituent(constituent);
+    });
+
+    taggedJets.push_back(edmJet);
+  }
+  std::sort(taggedJets.begin(), taggedJets.end(), [](l1t::PFJet a, l1t::PFJet b) { return (a.pt() > b.pt()); });
+
+  std::unique_ptr<l1t::PFJetCollection> taggedJetsCollection(new l1t::PFJetCollection);
+  taggedJetsCollection->swap(taggedJets);
+  iEvent.put(std::move(taggedJetsCollection), "l1tSC82ProngJets");
+}
+
+void L1TSC82ProngJetProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
+  edm::ParameterSetDescription desc;
+  desc.add<edm::InputTag>("jets", edm::InputTag("l1tSC8PFL1PuppiEmulator"));
+  desc.add<std::string>("l1tSC82ProngJetModelPath", std::string("L1TSC82ProngJetModel_v0"));
+  desc.add<double>("minPt", 0.);
+  desc.add<double>("maxEta", 5.);
+  desc.add<int>("maxJets", 16);
+  desc.add<int>("nParticles", 8);
+  descriptions.add("l1tSC82ProngJetProducer", desc);
+}
+
+#include "FWCore/Framework/interface/MakerMacros.h"
+DEFINE_FWK_MODULE(L1TSC82ProngJetProducer);

L1Trigger/Phase2L1ParticleFlow/python/l1pfJetMet_cff.py

@@ -60,9 +60,12 @@
                                                                      correctorFile = cms.string("L1Trigger/Phase2L1ParticleFlow/data/jecs/jecs_20220308.root"),
                                                                      correctorDir = cms.string('L1PuppiSC4EmuJets'))
 
+
+l1tSC82ProngJetProducer = l1tSC4NGJetProducer.clone()
+
 L1TPFJetsTask = cms.Task(
     l1tLayer2Deregionizer, l1tSC4PFL1PF, l1tSC4PFL1Puppi, l1tSC4PFL1PuppiCorrected, l1tSC8PFL1Puppi, l1tSC8PFL1PuppiCorrected, l1tSC4PFL1PuppiEmulator, l1tSC4PFL1PuppiCorrectedEmulator, l1tSC4PFL1PuppiCorrectedEmulatorMHT,
-    l1tSC8PFL1PuppiEmulator, l1tSC8PFL1PuppiCorrectedEmulator
+    l1tSC8PFL1PuppiEmulator, l1tSC8PFL1PuppiCorrectedEmulator, l1tSC82ProngJetProducer
 )
 
 L1TPFJetsExtendedTask = cms.Task(
@@ -71,5 +74,5 @@
 
 L1TPFJetsEmulationTask = cms.Task(
     l1tLayer2Deregionizer, l1tSC4PFL1PuppiEmulator, l1tSC4PFL1PuppiCorrectedEmulator, l1tSC4PFL1PuppiCorrectedEmulatorMHT,
-    l1tSC8PFL1PuppiEmulator, l1tSC8PFL1PuppiCorrectedEmulator
+    l1tSC8PFL1PuppiEmulator, l1tSC8PFL1PuppiCorrectedEmulator, l1tSC82ProngJetProducer
 )

L1Trigger/Phase2L1ParticleFlow/src/L1TSC82ProngJetID.cc

@@ -0,0 +1,154 @@
+#include "L1Trigger/Phase2L1ParticleFlow/interface/L1TSC82ProngJetID.h"
+#include "DataFormats/Math/interface/deltaPhi.h"
+#include <algorithm>
+#include <cmath>
+
+L1TSC82ProngJetID::L1TSC82ProngJetID(const std::shared_ptr<hls4mlEmulator::Model> model, int iNParticles)
+    : modelRef_(model) {
+  NNvectorVar_.clear();
+  fNParticles_ = iNParticles;
+
+  fPt_ = std::make_unique<float[]>(fNParticles_);
+  fPt_rel_ = std::make_unique<float[]>(fNParticles_);
+  fDEta_ = std::make_unique<float[]>(fNParticles_);
+  fDPhi_ = std::make_unique<float[]>(fNParticles_);
+  fPt_log_ = std::make_unique<float[]>(fNParticles_);
+  fMass_ = std::make_unique<float[]>(fNParticles_);
+  fZ0_ = std::make_unique<float[]>(fNParticles_);
+  fDxy_ = std::make_unique<float[]>(fNParticles_);
+  fIs_filled_ = std::make_unique<int[]>(fNParticles_);
+  fPuppi_weight_ = std::make_unique<float[]>(fNParticles_);
+  fEmID_ = std::make_unique<int[]>(fNParticles_);
+  fQuality_ = std::make_unique<float[]>(fNParticles_);
+
+  fId_ = std::make_unique<int[]>(fNParticles_);
+  fCharge_ = std::make_unique<int[]>(fNParticles_);
+}
+
+void L1TSC82ProngJetID::setNNVectorVar() {
+  NNvectorVar_.clear();
+  for (int i0 = 0; i0 < fNParticles_; i0++) {
+    NNvectorVar_.push_back(fPt_.get()[i0]);                              // pt
+    NNvectorVar_.push_back(fPt_rel_.get()[i0]);                          //pT as a fraction of jet pT
+    NNvectorVar_.push_back(fPt_log_.get()[i0]);                          // pt log
+    NNvectorVar_.push_back(fDEta_.get()[i0]);                            //dEta from jet axis
+    NNvectorVar_.push_back(fDPhi_.get()[i0]);                            //dPhi from jet axis
+    NNvectorVar_.push_back(fMass_.get()[i0]);                            // Mass
+    NNvectorVar_.push_back(fId_.get()[i0] == l1t::PFCandidate::Photon);  // Photon
+    NNvectorVar_.push_back(fId_.get()[i0] == l1t::PFCandidate::Electron && fCharge_.get()[i0] > 0);       // Positron
+    NNvectorVar_.push_back(fId_.get()[i0] == l1t::PFCandidate::Electron && fCharge_.get()[i0] < 0);       // Electron
+    NNvectorVar_.push_back(fId_.get()[i0] == l1t::PFCandidate::Muon && fCharge_.get()[i0] > 0);           // Anti-muon
+    NNvectorVar_.push_back(fId_.get()[i0] == l1t::PFCandidate::Muon && fCharge_.get()[i0] < 0);           // Muon
+    NNvectorVar_.push_back(fId_.get()[i0] == l1t::PFCandidate::NeutralHadron);                            // Neutral Had
+    NNvectorVar_.push_back(fId_.get()[i0] == l1t::PFCandidate::ChargedHadron && fCharge_.get()[i0] > 0);  // Anti-Pion
+    NNvectorVar_.push_back(fId_.get()[i0] == l1t::PFCandidate::ChargedHadron && fCharge_.get()[i0] < 0);  // Pion
+    NNvectorVar_.push_back(fZ0_.get()[i0]);                                                               // z0
+    NNvectorVar_.push_back(fDxy_.get()[i0]);                                                              // dxy
+    NNvectorVar_.push_back(fIs_filled_.get()[i0]);                                                        // isfilled
+    NNvectorVar_.push_back(fPuppi_weight_.get()[i0]);  // puppi weight
+    NNvectorVar_.push_back(fEmID_.get()[i0]);          // emID
+    NNvectorVar_.push_back(fQuality_.get()[i0]);       // quality
+  }
+}
+
+std::vector<float> L1TSC82ProngJetID::EvaluateNNFixed() {
+  const unsigned int NInputs = 160;
+  const unsigned int FeaturesPerParticle = 20;
+
+  if (fNParticles_ * FeaturesPerParticle != NInputs) {
+    throw std::runtime_error(
+        "L1TSC82ProngJetID::EvaluateNNFixed: "
+        "Model expects exactly 8 particles (160 inputs), but fNParticles_ = " +
+        std::to_string(fNParticles_));
+  }
+
+  prong_score prong_scores;
+  inputtype fillzero = 0.0;
+
+  inputtype modelInput[NInputs];
+  std::fill(std::begin(modelInput), std::end(modelInput), fillzero);
+
+  for (unsigned int i = 0; i < NInputs; i++) {
+    modelInput[i] = NNvectorVar_[i];
+  }
+
+  modelRef_->prepare_input(modelInput);
+  modelRef_->predict();
+  modelRef_->read_result(&prong_scores);
+
+  std::vector<float> prong_score_;
+  prong_score_.push_back(prong_scores.to_float());
+
+  return prong_score_;
+}  //end EvaluateNNFixed
+
+std::vector<float> L1TSC82ProngJetID::computeFixed(const l1t::PFJet &iJet) {
+  for (int i0 = 0; i0 < fNParticles_; i0++) {
+    fPt_rel_.get()[i0] = 0;
+    fPt_.get()[i0] = 0;
+    fDEta_.get()[i0] = 0;
+    fDPhi_.get()[i0] = 0;
+    fPt_log_.get()[i0] = 0;
+    fMass_.get()[i0] = 0;
+    fZ0_.get()[i0] = 0;
+    fDxy_.get()[i0] = 0;
+    fIs_filled_.get()[i0] = 0;
+    fPuppi_weight_.get()[i0] = 0;
+    fEmID_.get()[i0] = 0;
+    fQuality_.get()[i0] = 0;
+
+    fId_.get()[i0] = 0;
+    fCharge_.get()[i0] = 0;
+  }
+  auto iParts = iJet.constituents();
+  std::sort(iParts.begin(), iParts.end(), [](edm::Ptr<l1t::PFCandidate> i, edm::Ptr<l1t::PFCandidate> j) {
+    return (i->pt() > j->pt());
+  });
+
+  l1ct::Jet ctJet = l1ct::Jet::unpack(iJet.getHWJetCT());
+  float jet_pt_ = ctJet.floatPt();
+  float jet_eta_ = ctJet.floatEta();
+  float jet_phi_ = ctJet.floatPhi();
+
+  for (unsigned int i0 = 0; i0 < iParts.size(); i0++) {
+    if (i0 >= (unsigned int)fNParticles_)
+      break;
+    fPt_.get()[i0] = iParts[i0]->pt();
+    fZ0_.get()[i0] = iParts[i0]->z0();
+    fPt_rel_.get()[i0] = iParts[i0]->pt() / jet_pt_;
+    fPt_log_.get()[i0] = std::log(iParts[i0]->pt());
+    fDEta_.get()[i0] = jet_eta_ - float(iParts[i0]->eta());
+
+    if (1 < iParts[i0]->phi() - jet_phi_) {
+      fDPhi_.get()[i0] = 1;
+    } else if (-1 > iParts[i0]->phi() - jet_phi_) {
+      fDPhi_.get()[i0] = -1;
+    } else {
+      fDPhi_.get()[i0] = iParts[i0]->phi() - jet_phi_;
+    }
+
+    fIs_filled_.get()[i0] = 1;
+
+    float massCand = 0.13f;
+    if (abs(iParts[i0]->charge())) {
+      if ((iParts[i0]->id() == l1t::PFCandidate::Muon)) {
+        massCand = 0.105;
+      } else if ((iParts[i0]->id() == l1t::PFCandidate::Electron)) {
+        massCand = 0.005;
+      }
+    } else {
+      massCand = iParts[i0]->id() == l1t::PFCandidate::Photon ? 0.0 : 0.5;
+    }
+
+    fMass_.get()[i0] = massCand;
+    fPuppi_weight_.get()[i0] = iParts[i0]->puppiWeight();
+    fId_.get()[i0] = iParts[i0]->id();
+    fCharge_.get()[i0] = iParts[i0]->charge();
+
+    fDxy_.get()[i0] = iParts[i0]->hwDxy();
+    fEmID_.get()[i0] = iParts[i0]->hwEmID();
+    fQuality_.get()[i0] = iParts[i0]->hwTkQuality();
+  }
+  setNNVectorVar();
+  return EvaluateNNFixed();
+}