[L1T] Update L1T MET and Add JUMP algorithm (Phase 2)

L1Trigger/Phase2L1ParticleFlow/interface/jetmet/L1PFJUMPEmulator.h

@@ -0,0 +1,123 @@
+#ifndef L1Trigger_Phase2L1ParticleFlow_JUMP_h
+#define L1Trigger_Phase2L1ParticleFlow_JUMP_h
+
+#include "DataFormats/L1TParticleFlow/interface/jets.h"
+#include "DataFormats/L1TParticleFlow/interface/sums.h"
+#include "DataFormats/L1TParticleFlow/interface/layer1_emulator.h"
+
+#include "L1Trigger/Phase2L1ParticleFlow/interface/dbgPrintf.h"
+#include "L1Trigger/Phase2L1ParticleFlow/interface/jetmet/L1PFMetEmulator.h"
+
+#include "FWCore/Utilities/interface/FileInPath.h"
+#include "FWCore/Utilities/interface/Exception.h"
+#include "nlohmann/json.hpp"
+
+#include <vector>
+#include <numeric>
+#include <array>
+#include <fstream>
+#include <algorithm>
+#include <cmath>
+#include "ap_int.h"
+#include "ap_fixed.h"
+
+namespace L1JUMPEmu {
+  /*
+    Emulator for the JUMP Algorithm
+    DPS Note publicly available on CDS: CMS DP-2025/023
+    JUMP: Jet Uncertainty-aware MET Prediction
+    - Approximate L1 Jet energy resolution by pT, eta value
+    - Apply the estimated resolution to MET
+  */
+  struct JER_param {
+    std::array<ap_fixed<11, 1>, 5> par0;   // eta.par0 (slope)
+    std::array<ap_fixed<8, 5>, 5> par1;    // eta.par1 (offset)
+    std::array<L1METEmu::eta_t, 4> edges;  // |eta| boundaries in HW scale
+  };
+
+  inline const JER_param& Get_jer_param() {
+    static JER_param P = []() {
+      JER_param t{};
+      edm::FileInPath fip("L1Trigger/Phase2L1ParticleFlow/data/met/l1jump_jer_v1.json");
+      std::ifstream in(fip.fullPath());
+      if (!in)
+        throw cms::Exception("FileNotFound") << fip.fullPath();
+      nlohmann::json j;
+      in >> j;
+
+      for (int i = 0; i < 5; ++i) {
+        t.par0[i] = ap_fixed<11, 1>(j["eta"]["par0"][i].get<double>());
+        t.par1[i] = ap_fixed<8, 5>(j["eta"]["par1"][i].get<double>());
+      }
+      for (int i = 0; i < 4; ++i) {
+        t.edges[i] = l1ct::Scales::makeGlbEta(j["eta_edges"][i].get<double>());
+      }
+      return t;
+    }();
+    return P;
+  }
+
+  inline void Get_dPt(const l1ct::Jet jet, L1METEmu::proj2_t& dPx_2, L1METEmu::proj2_t& dPy_2) {
+    /*
+      L1 Jet Energy Resolution parameterization
+      - Fitted σ(pT)/pT as a function of jet pT in each η region (detector boundary at η≈1.3, 1.7, 2.5, 3.0)
+      - Derived from simulated QCD multijet samples to calculate detector‐dependent resolution
+      - σ(pT) ≈ eta_par0[i] * pT + eta_par1[i]
+    */
+
+    const auto& J = Get_jer_param();
+
+    L1METEmu::eta_t abseta = abs(jet.hwEta.to_float());
+    int etabin = 0;
+    for (uint i = 0; i < 4;) {
+      if (abseta < J.edges[i]) {
+        etabin = i + 1;
+        break;
+      }
+    }
+
+    dPx_2 = 0;
+    dPy_2 = 0;
+    l1ct::Sum jet_resolution;
+    jet_resolution.hwPt = J.par0[etabin] * jet.hwPt + J.par1[etabin];
+    jet_resolution.hwPhi = jet.hwPhi;
+    L1METEmu::Particle_xy dpt_xy = L1METEmu::Get_xy(jet_resolution.hwPt, jet_resolution.hwPhi);
+
+    dPx_2 = dpt_xy.hwPx * dpt_xy.hwPx;
+    dPy_2 = dpt_xy.hwPy * dpt_xy.hwPy;
+  }
+
+  inline void Met_dPt(const std::vector<l1ct::Jet> jets, L1METEmu::proj2_t& dPx_2, L1METEmu::proj2_t& dPy_2) {
+    L1METEmu::proj2_t each_dPx2 = 0;
+    L1METEmu::proj2_t each_dPy2 = 0;
+
+    L1METEmu::proj2_t sum_dPx2 = 0;
+    L1METEmu::proj2_t sum_dPy2 = 0;
+
+    for (uint i = 0; i < jets.size(); i++) {
+      Get_dPt(jets[i], each_dPx2, each_dPy2);
+      sum_dPx2 += each_dPx2;
+      sum_dPy2 += each_dPy2;
+    }
+
+    dPx_2 = sum_dPx2;
+    dPy_2 = sum_dPy2;
+  }
+}  // namespace L1JUMPEmu
+
+inline void JUMP_emu(const l1ct::Sum inMet, const std::vector<l1ct::Jet> jets, l1ct::Sum& outMet) {
+  L1METEmu::Particle_xy inMet_xy = L1METEmu::Get_xy(inMet.hwPt, inMet.hwPhi);
+
+  L1METEmu::proj2_t dPx_2;
+  L1METEmu::proj2_t dPy_2;
+  L1JUMPEmu::Met_dPt(jets, dPx_2, dPy_2);
+
+  L1METEmu::Particle_xy outMet_xy;
+  outMet_xy.hwPx = (inMet_xy.hwPx > 0) ? inMet_xy.hwPx + L1METEmu::proj2_t(sqrt(dPx_2.to_float()))
+                                       : inMet_xy.hwPx - L1METEmu::proj2_t(sqrt(dPx_2.to_float()));
+  outMet_xy.hwPy = (inMet_xy.hwPy > 0) ? inMet_xy.hwPy + L1METEmu::proj2_t(sqrt(dPy_2.to_float()))
+                                       : inMet_xy.hwPy - L1METEmu::proj2_t(sqrt(dPy_2.to_float()));
+  L1METEmu::pxpy_to_ptphi(outMet_xy, outMet);
+}
+
+#endif

L1Trigger/Phase2L1ParticleFlow/interface/jetmet/L1PFMetEmulator.h

@@ -0,0 +1,149 @@
+#ifndef L1Trigger_Phase2L1ParticleFlow_MET_h
+#define L1Trigger_Phase2L1ParticleFlow_MET_h
+
+#include "DataFormats/L1TParticleFlow/interface/jets.h"
+#include "DataFormats/L1TParticleFlow/interface/sums.h"
+#include "DataFormats/L1TParticleFlow/interface/layer1_emulator.h"
+
+#include "L1Trigger/Phase2L1ParticleFlow/interface/dbgPrintf.h"
+
+#include "FWCore/Utilities/interface/FileInPath.h"
+#include "FWCore/Utilities/interface/Exception.h"
+#include "nlohmann/json.hpp"
+
+#ifndef CMSSW_GIT_HASH
+#include "hls_math.h"
+#endif
+
+#include <vector>
+#include <array>
+#include <fstream>
+#include <numeric>
+#include <algorithm>
+#include "ap_int.h"
+#include "ap_fixed.h"
+
+namespace L1METEmu {
+  // Define Data types for P2 L1 MET Emulator
+  typedef l1ct::pt_t pt_t;
+  typedef l1ct::glbphi_t phi_t;
+  typedef l1ct::glbeta_t eta_t;
+  typedef l1ct::dpt_t pxy_t;
+  typedef l1ct::Sum Met;
+
+  typedef ap_fixed<22, 12> proj_t;
+  typedef ap_fixed<32, 22> proj2_t;
+  typedef ap_fixed<32, 2> poly_t;
+
+  struct Particle_xy {
+    // 44bits
+    proj_t hwPx;
+    proj_t hwPy;
+  };
+
+  struct Poly2_param {
+    std::array<poly_t, 16> c0, c1, c2, s0, s1, s2;
+    std::array<phi_t, 17> phi_edges;
+  };
+
+  inline const Poly2_param& Get_Poly2_param() {
+    static Poly2_param P = [] {
+      Poly2_param t{};
+      edm::FileInPath f("L1Trigger/Phase2L1ParticleFlow/data/met/l1met_ptphi2pxpy_poly2_v1.json");
+      std::ifstream in(f.fullPath());
+      if (!in)
+        throw cms::Exception("FileNotFound") << f.fullPath();
+      nlohmann::json j;
+      in >> j;
+
+      for (int i = 0; i < 16; ++i) {
+        t.c0[i] = poly_t(j["cos"]["par0"][i].get<double>());
+        t.c1[i] = poly_t(j["cos"]["par1"][i].get<double>());
+        t.c2[i] = poly_t(j["cos"]["par2"][i].get<double>());
+        t.s0[i] = poly_t(j["sin"]["par0"][i].get<double>());
+        t.s1[i] = poly_t(j["sin"]["par1"][i].get<double>());
+        t.s2[i] = poly_t(j["sin"]["par2"][i].get<double>());
+      }
+      for (int i = 0; i < 17; ++i) {
+        t.phi_edges[i] = l1ct::Scales::makeGlbPhi(j["phi_edges"][i].get<double>() * M_PI);
+      }
+      return t;
+    }();
+    return P;
+  }
+
+  inline Particle_xy Get_xy(const l1ct::pt_t hwPt, const l1ct::glbphi_t hwPhi) {
+    /*
+      Convert pt, phi to px, py
+      - Use 2nd order Polynomial interpolation for cos, sin
+      - Divide the sine and cosine value from -2 pi to 2 pi into 16 parts
+      - Fitting the value with 2nd order function
+    */
+
+    const auto& P = L1METEmu::Get_Poly2_param();
+    int phibin = 0;
+    for (int i = 0; i < 16; i++) {
+      if (hwPhi >= P.phi_edges[i] && hwPhi < P.phi_edges[i + 1]) {
+        phibin = i;
+        break;
+      }
+    }
+
+    Particle_xy proj_xy;
+
+    poly_t cos_var = P.c0[phibin] + P.c1[phibin] * (hwPhi - P.phi_edges[phibin]) +
+                     P.c2[phibin] * (hwPhi - P.phi_edges[phibin]) * (hwPhi - P.phi_edges[phibin]);
+    poly_t sin_var = P.s0[phibin] + P.s1[phibin] * (hwPhi - P.phi_edges[phibin]) +
+                     P.s2[phibin] * (hwPhi - P.phi_edges[phibin]) * (hwPhi - P.phi_edges[phibin]);
+
+    proj_xy.hwPx = hwPt * cos_var;
+    proj_xy.hwPy = hwPt * sin_var;
+
+    return proj_xy;
+  }
+
+  inline void Sum_Particles(const std::vector<Particle_xy> particles_xy, Particle_xy& met_xy) {
+    met_xy.hwPx = 0;
+    met_xy.hwPy = 0;
+
+    for (uint i = 0; i < particles_xy.size(); ++i) {
+      met_xy.hwPx -= particles_xy[i].hwPx;
+      met_xy.hwPy -= particles_xy[i].hwPy;
+    }
+  }
+
+  inline void pxpy_to_ptphi(const Particle_xy met_xy, l1ct::Sum& hls_met) {
+    // convert x, y coordinate to pt, phi coordinate using math library
+    hls_met.clear();
+
+#ifdef CMSSW_GIT_HASH
+    hls_met.hwPt = hypot(met_xy.hwPx.to_float(), met_xy.hwPy.to_float());
+    hls_met.hwPhi = phi_t(
+        ap_fixed<26, 11>(l1ct::Scales::makeGlbPhi(atan2(met_xy.hwPy.to_float(), met_xy.hwPx.to_float()))));  // 720/pi
+#else
+    hls_met.hwPt = hls::hypot(met_xy.hwPx, met_xy.hwPy);
+    hls_met.hwPhi = phi_t(ap_fixed<26, 11>(l1ct::Scales::makeGlbPhi(hls::atan2(met_xy.hwPy, met_xy.hwPx))));
+#endif
+  }
+
+  inline void met_format(const l1ct::Sum d, ap_uint<l1gt::Sum::BITWIDTH>& q) {
+    // Change output formats to GT formats
+    q = d.toGT().pack();
+  }
+
+}  // namespace L1METEmu
+
+inline void puppimet_emu(const std::vector<l1ct::PuppiObjEmu> particles, l1ct::Sum& out_met) {
+  std::vector<L1METEmu::Particle_xy> particles_xy;
+
+  for (uint i = 0; i < particles.size(); i++) {
+    L1METEmu::Particle_xy each_particle_xy = L1METEmu::Get_xy(particles[i].hwPt, particles[i].hwPhi);
+    particles_xy.push_back(each_particle_xy);
+  }
+
+  L1METEmu::Particle_xy met_xy;
+  L1METEmu::Sum_Particles(particles_xy, met_xy);
+  L1METEmu::pxpy_to_ptphi(met_xy, out_met);
+}
+
+#endif

L1Trigger/Phase2L1ParticleFlow/plugins/L1JUMPProducer.cc

@@ -0,0 +1,122 @@
+#include <vector>
+#include <string>
+#include <ap_int.h>
+#include <ap_fixed.h>
+#include <TVector2.h>
+#include <iostream>
+
+#include "FWCore/Framework/interface/global/EDProducer.h"
+#include "FWCore/Framework/interface/Event.h"
+#include "FWCore/Utilities/interface/InputTag.h"
+#include "FWCore/Framework/interface/MakerMacros.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+#include "FWCore/MessageLogger/interface/MessageLogger.h"
+
+#include "DataFormats/L1TParticleFlow/interface/PFCandidate.h"
+#include "DataFormats/L1TParticleFlow/interface/PFJet.h"
+
+#include "DataFormats/L1TParticleFlow/interface/puppi.h"
+#include "DataFormats/L1TParticleFlow/interface/sums.h"
+#include "DataFormats/L1TParticleFlow/interface/jets.h"
+#include "DataFormats/L1TParticleFlow/interface/layer1_emulator.h"
+
+#include "DataFormats/L1Trigger/interface/EtSum.h"
+#include "DataFormats/Math/interface/LorentzVector.h"
+
+#include "L1Trigger/Phase2L1ParticleFlow/interface/jetmet/L1PFJUMPEmulator.h"
+
+class L1JUMPProducer : public edm::global::EDProducer<> {
+  /*
+    Producer for the JUMP Algorithm
+    JUMP: Jet Uncertainty-aware MET Prediction
+    - Approximate L1 Jet energy resolution by pT, eta value
+    - Apply the estimated resolution to MET
+  */
+public:
+  explicit L1JUMPProducer(const edm::ParameterSet&);
+  ~L1JUMPProducer() override;
+
+private:
+  void produce(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const override;
+  edm::EDGetTokenT<std::vector<l1t::EtSum>> metToken;
+  edm::EDGetTokenT<std::vector<l1t::PFJet>> jetsToken;
+
+  typedef l1ct::pt_t pt_t;
+  typedef l1ct::Jet Jet;
+
+  static constexpr float ptLSB_ = 0.25;
+  static constexpr float phiLSB_ = M_PI / 720;
+  static constexpr float maxPt_ = ((1 << pt_t::width) - 1) * ptLSB_;
+
+  void CalcJUMP_HLS(const l1t::EtSum& metVector,
+                    const std::vector<l1ct::Jet>& jets,
+                    reco::Candidate::PolarLorentzVector& JUMPVector) const;
+
+  std::vector<l1ct::Jet> convertEDMToHW(const std::vector<l1t::PFJet> edmJets) const;
+
+  double minJetPt;
+  double maxJetEta;
+};
+
+L1JUMPProducer::L1JUMPProducer(const edm::ParameterSet& cfg)
+    : metToken(consumes<std::vector<l1t::EtSum>>(cfg.getParameter<edm::InputTag>("RawMET"))),
+      jetsToken(consumes<std::vector<l1t::PFJet>>(cfg.getParameter<edm::InputTag>("L1PFJets"))),
+      minJetPt(cfg.getParameter<double>("MinJetpT")),
+      maxJetEta(cfg.getParameter<double>("MaxJetEta")) {
+  produces<std::vector<l1t::EtSum>>();
+}
+
+void L1JUMPProducer::produce(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const {
+  // Load Met
+  l1t::EtSum rawMET = iEvent.get(metToken)[0];
+  // Load Jets
+  l1t::PFJetCollection edmJets = iEvent.get(jetsToken);
+
+  std::vector<l1ct::Jet> hwJets = convertEDMToHW(edmJets);  // convert to the emulator format
+  // Apply pT and eta selections
+  std::vector<l1ct::Jet> hwJetsFiltered;
+  std::copy_if(hwJets.begin(), hwJets.end(), std::back_inserter(hwJetsFiltered), [&](auto jet) {
+    return jet.hwPt > l1ct::Scales::makePtFromFloat(float(minJetPt)) &&
+           std::abs(jet.hwEta) < l1ct::Scales::makeGlbEta(maxJetEta);
+  });
+
+  // JUMP Algorithm
+  reco::Candidate::PolarLorentzVector JUMPVector;
+  CalcJUMP_HLS(rawMET, hwJetsFiltered, JUMPVector);
+
+  l1t::EtSum theJUMP(JUMPVector, l1t::EtSum::EtSumType::kMissingEt, 0, 0, 0, 0);
+  auto JUMPCollection = std::make_unique<std::vector<l1t::EtSum>>(0);
+  JUMPCollection->push_back(theJUMP);
+  iEvent.put(std::move(JUMPCollection));
+}
+
+void L1JUMPProducer::CalcJUMP_HLS(const l1t::EtSum& metVector,
+                                  const std::vector<l1ct::Jet>& jets,
+                                  reco::Candidate::PolarLorentzVector& outMet_Vector) const {
+  // JUMP Calculate
+  l1ct::Sum inMet;
+  inMet.hwPt = metVector.pt();
+  inMet.hwPhi = l1ct::Scales::makeGlbPhi(metVector.phi());
+
+  l1ct::Sum outMet;
+
+  JUMP_emu(inMet, jets, outMet);
+
+  outMet_Vector.SetPt(outMet.hwPt.to_double());
+  outMet_Vector.SetPhi(outMet.hwPhi.to_double() * phiLSB_);
+  outMet_Vector.SetEta(0);
+}
+
+std::vector<l1ct::Jet> L1JUMPProducer::convertEDMToHW(const std::vector<l1t::PFJet> edmJets) const {
+  std::vector<l1ct::Jet> hwJets;
+  std::for_each(edmJets.begin(), edmJets.end(), [&](const l1t::PFJet jet) {
+    l1ct::Jet hwJet = l1ct::Jet::unpack(jet.getHWJetCT());
+    hwJets.push_back(hwJet);
+  });
+  return hwJets;
+}
+
+L1JUMPProducer::~L1JUMPProducer() {}
+
+#include "FWCore/Framework/interface/MakerMacros.h"
+DEFINE_FWK_MODULE(L1JUMPProducer);