diff --git a/L1Trigger/Phase2L1ParticleFlow/data/met/l1jump_jer_v1.json b/L1Trigger/Phase2L1ParticleFlow/data/met/l1jump_jer_v1.json new file mode 100644 index 0000000000000..4c85b333ed9ea --- /dev/null +++ b/L1Trigger/Phase2L1ParticleFlow/data/met/l1jump_jer_v1.json @@ -0,0 +1,27 @@ +{ + "version": 1, + "description": "JUMP L1 jet resolution: σ(pT)[GeV] = a(ηbin)·pT + b(ηbin); |η| bins at [1.3, 1.7, 2.5, 3.0]; used to build dPx,dPy for MET update.", + "eta_bins": 4, + "eta_edges": [ + 1.3, + 1.7, + 2.5, + 3.0 + ], + "eta": { + "par0": [ + 0, + 0.073, + 0.247, + 0.128, + 0.091 + ], + "par1": [ + 0, + 12.322, + 6.061, + 10.944, + 12.660 + ] + } +} \ No newline at end of file diff --git a/L1Trigger/Phase2L1ParticleFlow/data/met/l1met_ptphi2pxpy_poly2_v1.json b/L1Trigger/Phase2L1ParticleFlow/data/met/l1met_ptphi2pxpy_poly2_v1.json new file mode 100644 index 0000000000000..c8bcd1a9770b2 --- /dev/null +++ b/L1Trigger/Phase2L1ParticleFlow/data/met/l1met_ptphi2pxpy_poly2_v1.json @@ -0,0 +1,136 @@ +{ + "version": 1, + "description": "L1METEmu 2nd-order polynomial coefficients for cos/sin; 16 phi bins over [-pi, pi]; phi_edges in units of pi.", + "phi_bins": 16, + "phi_edges": [ + -1.0, + -0.875, + -0.75, + -0.625, + -0.5, + -0.375, + -0.25, + -0.125, + 0.0, + 0.125, + 0.25, + 0.375, + 0.5, + 0.625, + 0.75, + 0.875, + 1.0 + ], + "cos": { + "par0": [ + -1.00007, + -0.924181, + -0.707596, + -0.382902, + -0.000618262, + 0.382137, + 0.707056, + 0.923708, + 1.00007, + 0.924181, + 0.707594, + 0.383285, + 0.000188727, + -0.382139, + -0.706719, + -0.923708 + ], + "par1": [ + 9.164680268990924e-06, + 0.0017064607695524156, + 0.0031441321076514446, + 0.004079929656016374, + 0.004437063290882583, + 0.004095969231842202, + 0.0031107221424451436, + 0.001689531075808071, + -9.161756842493832e-06, + -0.001706456406229286, + -0.003143961938049376, + -0.004103015998697129, + -0.004411145151490469, + -0.0040958165155326525, + -0.0031310072316764474, + -0.001689531075808071 + ], + "par2": [ + 9.319674765430664e-06, + 7.871694899063284e-06, + 5.222989318251642e-06, + 2.0256106486379287e-06, + -1.9299417402361656e-06, + -5.35167113952279e-06, + -7.740062096537953e-06, + -9.348822844786505e-06, + -9.319674765430664e-06, + -7.871694899063284e-06, + -5.225331064666252e-06, + -1.780776301343235e-06, + 1.6556927733433181e-06, + 5.3495197789955455e-06, + 7.954684107366423e-06, + 9.348822844786505e-06 + ] + }, + "sin": { + "par0": [ + 0.000524872, + -0.382229, + -0.706791, + -0.923959, + -1.00008, + -0.924156, + -0.707264, + -0.383199, + -0.000525527, + 0.382228, + 0.706792, + 0.923752, + 1.00013, + 0.924155, + 0.707535, + 0.3832 + ], + "par1": [ + -0.004431478237276202, + -0.00409041472149773, + -0.0031267268116859314, + -0.00167440343451641, + 9.741773386162849e-06, + 0.0017049641497188307, + 0.00312406082125351, + 0.0040978672774037465, + 0.004431478237276202, + 0.00409041472149773, + 0.0031266351819002015, + 0.0016868781753450394, + -1.249302315254411e-05, + -0.001704846339994321, + -0.003140405829698437, + -0.0040978672774037465 + ], + "par2": [ + 1.870674613498914e-06, + 5.292404012785538e-06, + 7.909829192302831e-06, + 9.188746390688592e-06, + 9.313525301268721e-06, + 7.887020962996302e-06, + 5.435897856093815e-06, + 1.8358587462761668e-06, + -1.870668901922293e-06, + -5.292404012785538e-06, + -7.908420336736317e-06, + -9.320836119343602e-06, + -9.284396260501616e-06, + -7.88869635880513e-06, + -5.262894200243701e-06, + -1.835864457852788e-06 + ] + } +} \ No newline at end of file diff --git a/L1Trigger/Phase2L1ParticleFlow/interface/jetmet/L1PFJUMPEmulator.h b/L1Trigger/Phase2L1ParticleFlow/interface/jetmet/L1PFJUMPEmulator.h new file mode 100644 index 0000000000000..25d6e68e6cd80 --- /dev/null +++ b/L1Trigger/Phase2L1ParticleFlow/interface/jetmet/L1PFJUMPEmulator.h @@ -0,0 +1,158 @@ +#ifndef L1Trigger_Phase2L1ParticleFlow_JUMPEmulator_h +#define L1Trigger_Phase2L1ParticleFlow_JUMPEmulator_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 +#include +#include +#include +#include +#include +#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::vector> par0; // eta.par0 (slope) + std::vector> par1; // eta.par1 (offset) + std::vector edges; // |eta| boundaries in HW scale + unsigned int eta_bins = 0; + }; + + struct JER_Path { + std::string path = "L1Trigger/Phase2L1ParticleFlow/data/met/l1jump_jer_v1.json"; + }; + + inline JER_Path& jer_path_config() { + static JER_Path jump_p; + return jump_p; + } + + inline void SetJERFile(std::string jump_p) { jer_path_config().path = std::move(jump_p); } + + inline const JER_param& Get_jer_param() { + static JER_param P = []() { + JER_param t{}; + std::string path = jer_path_config().path; + +#ifdef CMSSW_GIT_HASH + edm::FileInPath f(path); + std::ifstream in(f.fullPath()); + if (!in) { + throw cms::Exception("FileNotFound") << f.fullPath(); + } +#else + path = "l1jump_jer_v1.json"; // For HLS Emulator + std::ifstream in(path); + if (!in) { + throw std::runtime_error(std::string("File not found: ") + path); + } +#endif + + nlohmann::json j; + in >> j; + + unsigned int N = j["eta_bins"].get(); + t.eta_bins = N; + + t.par0.resize(N + 1); + t.par1.resize(N + 1); + t.edges.resize(N); + + for (unsigned int i = 0; i < N + 1; ++i) { + t.par0[i] = ap_fixed<11, 1>(j["eta"]["par0"][i].get()); + t.par1[i] = ap_fixed<8, 5>(j["eta"]["par1"][i].get()); + } + for (unsigned int i = 0; i < N; ++i) { + t.edges[i] = l1ct::Scales::makeGlbEta(j["eta_edges"][i].get()); + } + 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()); + unsigned int etabin = 0; + for (unsigned int i = 0; i < J.eta_bins; i++) { + if (abseta < J.edges[i]) { + etabin = i + 1; + break; + } + } + // If abseta >= all bin edges, etabin remains 0, which is reserved for overflow (jets beyond max eta range). + // This is intentional: par0[0] and par1[0] should be set for this overflow bin. + + 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& 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 (unsigned int 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& 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; + float sqrt_dPx_2 = sqrt(dPx_2.to_float()); + float sqrt_dPy_2 = sqrt(dPy_2.to_float()); + outMet_xy.hwPx = (inMet_xy.hwPx > 0) ? inMet_xy.hwPx + L1METEmu::proj2_t(sqrt_dPx_2) + : inMet_xy.hwPx - L1METEmu::proj2_t(sqrt_dPx_2); + outMet_xy.hwPy = (inMet_xy.hwPy > 0) ? inMet_xy.hwPy + L1METEmu::proj2_t(sqrt_dPy_2) + : inMet_xy.hwPy - L1METEmu::proj2_t(sqrt_dPy_2); + L1METEmu::pxpy_to_ptphi(outMet_xy, outMet); +} + +#endif diff --git a/L1Trigger/Phase2L1ParticleFlow/interface/jetmet/L1PFMetEmulator.h b/L1Trigger/Phase2L1ParticleFlow/interface/jetmet/L1PFMetEmulator.h new file mode 100644 index 0000000000000..a0fcd0c736768 --- /dev/null +++ b/L1Trigger/Phase2L1ParticleFlow/interface/jetmet/L1PFMetEmulator.h @@ -0,0 +1,192 @@ +#ifndef L1Trigger_Phase2L1ParticleFlow_METEmulator_h +#define L1Trigger_Phase2L1ParticleFlow_METEmulator_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 "nlohmann/json.hpp" + +#ifdef CMSSW_GIT_HASH +#include "FWCore/Utilities/interface/FileInPath.h" +#include "FWCore/Utilities/interface/Exception.h" +#endif + +#ifndef CMSSW_GIT_HASH +#include "hls_math.h" +#endif + +#include +#include +#include +#include +#include +#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::vector c0, c1, c2, s0, s1, s2; + std::vector phi_edges; + unsigned int phi_bins = 0; + }; + + struct Poly2_Path { + std::string path = "L1Trigger/Phase2L1ParticleFlow/data/met/l1met_ptphi2pxpy_poly2_v1.json"; + }; + + inline Poly2_Path& poly2_path_config() { + static Poly2_Path met_p; + return met_p; + } + + inline void SetPoly2File(std::string met_p) { poly2_path_config().path = std::move(met_p); } + + inline const Poly2_param& Get_Poly2_param() { + static Poly2_param P = [] { + Poly2_param t{}; + std::string path = poly2_path_config().path; + +#ifdef CMSSW_GIT_HASH + edm::FileInPath f(path); + std::ifstream in(f.fullPath()); + if (!in) { + throw cms::Exception("FileNotFound") << f.fullPath(); + } +#else + path = "l1met_ptphi2pxpy_poly2_v1.json"; // For HLS Emulator + std::ifstream in(path); + if (!in) { + throw std::runtime_error(std::string("File not found: ") + path); + } +#endif + + nlohmann::json j; + in >> j; + + unsigned int N = j["phi_bins"].get(); + t.phi_bins = N; + + t.c0.resize(N); + t.c1.resize(N); + t.c2.resize(N); + t.s0.resize(N); + t.s1.resize(N); + t.s2.resize(N); + t.phi_edges.resize(N + 1); + + for (unsigned int i = 0; i < N; ++i) { + t.c0[i] = poly_t(j["cos"]["par0"][i].get()); + t.c1[i] = poly_t(j["cos"]["par1"][i].get()); + t.c2[i] = poly_t(j["cos"]["par2"][i].get()); + t.s0[i] = poly_t(j["sin"]["par0"][i].get()); + t.s1[i] = poly_t(j["sin"]["par1"][i].get()); + t.s2[i] = poly_t(j["sin"]["par2"][i].get()); + } + for (unsigned int i = 0; i < N + 1; ++i) { + t.phi_edges[i] = l1ct::Scales::makeGlbPhi(j["phi_edges"][i].get() * 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 -pi to pi into 16 parts + - Fitting the value with 2nd order function + */ + + const auto& P = L1METEmu::Get_Poly2_param(); + int phibin = 0; + + for (unsigned int i = 0; i < P.phi_bins; i++) { + if (hwPhi >= P.phi_edges[i] && hwPhi < P.phi_edges[i + 1]) { + phibin = i; + break; + } + } + // Handle the edge case where hwPhi is exactly equal to the last bin edge + if (hwPhi == P.phi_edges[P.phi_bins]) { + phibin = P.phi_bins - 1; + } + + 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& particles_xy, Particle_xy& met_xy) { + met_xy.hwPx = 0; + met_xy.hwPy = 0; + + for (unsigned int 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>(atan2(met_xy.hwPy.to_float(), met_xy.hwPx.to_float())) * + ap_fixed<26, 11>(229.29936)); // Scale for L1 phi value (720 / M_PI) +#else + hls_met.hwPt = hls::hypot(met_xy.hwPx, met_xy.hwPy); + hls_met.hwPhi = phi_t(ap_fixed<26, 11>(hls::atan2(met_xy.hwPy, met_xy.hwPx)) * ap_fixed<26, 11>(229.29936)); +#endif + } + + inline void met_format(const l1ct::Sum d, ap_uint& q) { + // Change output formats to GT formats + q = d.toGT().pack(); + } + +} // namespace L1METEmu + +inline void puppimet_emu(const std::vector& particles, l1ct::Sum& out_met) { + std::vector particles_xy; + + for (unsigned int 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 diff --git a/L1Trigger/Phase2L1ParticleFlow/plugins/L1JUMPProducer.cc b/L1Trigger/Phase2L1ParticleFlow/plugins/L1JUMPProducer.cc new file mode 100644 index 0000000000000..5235b6be399b8 --- /dev/null +++ b/L1Trigger/Phase2L1ParticleFlow/plugins/L1JUMPProducer.cc @@ -0,0 +1,131 @@ +#include +#include +#include +#include + +#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; + static void fillDescriptions(edm::ConfigurationDescriptions& descriptions); + +private: + void produce(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const override; + edm::EDGetTokenT> metToken; + edm::EDGetTokenT> jetsToken; + + typedef l1ct::pt_t pt_t; + typedef l1ct::Jet Jet; + + static constexpr float phiLSB_ = M_PI / 720; + + void CalcJUMP_HLS(const l1t::EtSum& metVector, + const std::vector& jets, + reco::Candidate::PolarLorentzVector& JUMPVector) const; + + std::vector convertEDMToHW(const std::vector& edmJets) const; + + double minJetPt; + double maxJetEta; + std::string jerFilePath_; +}; + +L1JUMPProducer::L1JUMPProducer(const edm::ParameterSet& cfg) + : metToken(consumes>(cfg.getParameter("RawMET"))), + jetsToken(consumes>(cfg.getParameter("L1PFJets"))), + minJetPt(cfg.getParameter("MinJetpT")), + maxJetEta(cfg.getParameter("MaxJetEta")), + jerFilePath_(cfg.getParameter("JERFile")) { + produces>(); + L1JUMPEmu::SetJERFile(jerFilePath_); +} + +void L1JUMPProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { + edm::ParameterSetDescription desc; + desc.add("RawMET", edm::InputTag("l1tMETPFProducer")); + desc.add("L1PFJets", edm::InputTag("l1tSC4PFL1PuppiCorrectedEmulator")); + desc.add("MinJetpT", 30); + desc.add("MaxJetEta", 3.0); + desc.add("JERFile", "L1Trigger/Phase2L1ParticleFlow/data/met/l1jump_jer_v1.json"); + descriptions.add("L1JUMPProducer", desc); +} + +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 hwJets = convertEDMToHW(edmJets); // convert to the emulator format + // Apply pT and eta selections + std::vector 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>(0); + JUMPCollection->push_back(theJUMP); + iEvent.put(std::move(JUMPCollection)); +} + +void L1JUMPProducer::CalcJUMP_HLS(const l1t::EtSum& metVector, + const std::vector& 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 L1JUMPProducer::convertEDMToHW(const std::vector& edmJets) const { + std::vector 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() {} + +DEFINE_FWK_MODULE(L1JUMPProducer); \ No newline at end of file diff --git a/L1Trigger/Phase2L1ParticleFlow/plugins/L1MetPfProducer.cc b/L1Trigger/Phase2L1ParticleFlow/plugins/L1MetPfProducer.cc index efcb7e78d02fc..97903712056a8 100644 --- a/L1Trigger/Phase2L1ParticleFlow/plugins/L1MetPfProducer.cc +++ b/L1Trigger/Phase2L1ParticleFlow/plugins/L1MetPfProducer.cc @@ -9,11 +9,20 @@ #include "FWCore/Framework/interface/MakerMacros.h" #include "FWCore/ParameterSet/interface/ParameterSet.h" #include "FWCore/MessageLogger/interface/MessageLogger.h" +#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h" +#include "FWCore/ParameterSet/interface/ParameterSetDescription.h" #include "DataFormats/L1TParticleFlow/interface/PFCandidate.h" +// For HLS MET Data Formats +#include "DataFormats/L1TParticleFlow/interface/puppi.h" +#include "DataFormats/L1TParticleFlow/interface/sums.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/L1PFMetEmulator.h" + #include "hls4ml/emulator.h" using namespace l1t; @@ -26,36 +35,20 @@ class L1MetPfProducer : public edm::global::EDProducer<> { private: void produce(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const override; - edm::EDGetTokenT> _l1PFToken; + edm::EDGetTokenT> _l1PFToken; int maxCands_ = 128; // quantization controllers - typedef ap_ufixed<14, 12, AP_RND, AP_WRAP> pt_t; // LSB is 0.25 and max is 4 TeV - typedef ap_int<12> phi_t; // LSB is pi/720 ~ 0.0044 and max is +/-8.9 - static constexpr float ptLSB_ = 0.25; // GeV - static constexpr float phiLSB_ = M_PI / 720; // rad - - // derived, helper types - typedef ap_fixed pxy_t; - typedef ap_fixed<2 * pt_t::width, 2 * pt_t::iwidth, AP_RND, AP_SAT> pt2_t; - // derived, helper constants - static constexpr float maxPt_ = ((1 << pt_t::width) - 1) * ptLSB_; - const phi_t hwPi_ = round(M_PI / phiLSB_); - const phi_t hwPiOverTwo_ = round(M_PI / (2 * phiLSB_)); - - typedef ap_ufixed inv_t; // can't easily use the MAXPT/pt trick with ap_fixed - - // to make configurable... - static constexpr int dropBits_ = 2; - static constexpr int dropFactor_ = (1 << dropBits_); - static constexpr int invTableBits_ = 10; - static constexpr int invTableSize_ = (1 << invTableBits_); + typedef l1ct::pt_t pt_t; + typedef l1ct::glbphi_t phi_t; + static constexpr float phiLSB_ = M_PI / 720; // rad // hls4ml emulator objects bool useMlModel_; std::shared_ptr model; std::string modelVersion_; + typedef ap_fixed<32, 16> input_t; typedef ap_fixed<32, 16> result_t; static constexpr int numContInputs_ = 4; @@ -63,21 +56,11 @@ class L1MetPfProducer : public edm::global::EDProducer<> { static constexpr int numCatInputs_ = 2; static constexpr int numInputs_ = numContInputs_ + numPxPyInputs_ + numCatInputs_; - void Project(pt_t pt, phi_t phi, pxy_t& pxy, bool isX, bool debug = false) const; - void PhiFromXY(pxy_t px, pxy_t py, phi_t& phi, bool debug = false) const; + void CalcMetHLS(const std::vector& pfcands, reco::Candidate::PolarLorentzVector& metVector) const; int EncodePdgId(int pdgId) const; - void CalcMetHLS(const std::vector& pt, - const std::vector& phi, - reco::Candidate::PolarLorentzVector& metVector) const; - void CalcMlMet(const std::vector& pt, - const std::vector& eta, - const std::vector& phi, - const std::vector& puppiWeight, - const std::vector& pdgId, - const std::vector& charge, - reco::Candidate::PolarLorentzVector& metVector) const; + void CalcMlMet(const std::vector& pfcands, reco::Candidate::PolarLorentzVector& metVector) const; }; L1MetPfProducer::L1MetPfProducer(const edm::ParameterSet& cfg) @@ -89,6 +72,9 @@ L1MetPfProducer::L1MetPfProducer(const edm::ParameterSet& cfg) if (useMlModel_) { hls4mlEmulator::ModelLoader loader(modelVersion_); model = loader.load_model(); + } else { + edm::FileInPath f = cfg.getParameter("Poly2File"); + L1METEmu::SetPoly2File(f.fullPath()); } } @@ -97,6 +83,8 @@ void L1MetPfProducer::fillDescriptions(edm::ConfigurationDescriptions& descripti desc.add("L1PFObjects", edm::InputTag("L1PFProducer", "l1pfCandidates")); desc.add("maxCands", 128); desc.add("modelVersion", ""); + desc.add("Poly2File", + edm::FileInPath("L1Trigger/Phase2L1ParticleFlow/data/met/l1met_ptphi2pxpy_poly2_v1.json")); descriptions.add("L1MetPfProducer", desc); } @@ -104,38 +92,41 @@ void L1MetPfProducer::produce(edm::StreamID, edm::Event& iEvent, const edm::Even edm::Handle l1PFCandidates; iEvent.getByToken(_l1PFToken, l1PFCandidates); - std::vector pt; - std::vector eta; - std::vector phi; - std::vector puppiWeight; - std::vector pdgId; - std::vector charge; - - for (int i = 0; i < int(l1PFCandidates->size()) && (i < maxCands_ || maxCands_ < 0); i++) { - const auto& l1PFCand = l1PFCandidates->at(i); - pt.push_back(l1PFCand.pt()); - eta.push_back(l1PFCand.eta()); - phi.push_back(l1PFCand.phi()); - puppiWeight.push_back(l1PFCand.puppiWeight()); - pdgId.push_back(l1PFCand.pdgId()); - charge.push_back(l1PFCand.charge()); - } - + const std::vector& pfcands = *l1PFCandidates; reco::Candidate::PolarLorentzVector metVector; if (useMlModel_) { - CalcMlMet(pt, eta, phi, puppiWeight, pdgId, charge, metVector); + CalcMlMet(pfcands, metVector); } else { - CalcMetHLS(pt, phi, metVector); + CalcMetHLS(pfcands, metVector); } - l1t::EtSum theMET(metVector, l1t::EtSum::EtSumType::kTotalHt, 0, 0, 0, 0); + l1t::EtSum theMET(metVector, l1t::EtSum::EtSumType::kMissingEt, 0, 0, 0, 0); auto metCollection = std::make_unique>(0); metCollection->push_back(theMET); iEvent.put(std::move(metCollection)); } +void L1MetPfProducer::CalcMetHLS(const std::vector& pfcands, + reco::Candidate::PolarLorentzVector& metVector) const { + std::vector particles; + l1ct::Sum hw_met; + + for (int i = 0; i < int(pfcands.size()) && (i < maxCands_ || maxCands_ < 0); i++) { + const auto& cand = pfcands[i]; + l1ct::PuppiObjEmu each_particle; + each_particle.initFromBits(cand.encodedPuppi64()); + particles.push_back(each_particle); + } + + puppimet_emu(particles, hw_met); + + metVector.SetPt(hw_met.hwPt.to_double()); + metVector.SetPhi(hw_met.hwPhi.to_double() * phiLSB_); + metVector.SetEta(0); +} + int L1MetPfProducer::EncodePdgId(int pdgId) const { switch (abs(pdgId)) { case 211: // charged hadron (pion) @@ -153,13 +144,25 @@ int L1MetPfProducer::EncodePdgId(int pdgId) const { } } -void L1MetPfProducer::CalcMlMet(const std::vector& pt, - const std::vector& eta, - const std::vector& phi, - const std::vector& puppiWeight, - const std::vector& pdgId, - const std::vector& charge, +void L1MetPfProducer::CalcMlMet(const std::vector& pfcands, reco::Candidate::PolarLorentzVector& metVector) const { + std::vector pt; + std::vector eta; + std::vector phi; + std::vector puppiWeight; + std::vector pdgId; + std::vector charge; + + for (int i = 0; i < int(pfcands.size()) && (i < maxCands_ || maxCands_ < 0); i++) { + const auto& l1PFCand = pfcands[i]; + pt.push_back(l1PFCand.pt()); + eta.push_back(l1PFCand.eta()); + phi.push_back(l1PFCand.phi()); + puppiWeight.push_back(l1PFCand.puppiWeight()); + pdgId.push_back(l1PFCand.pdgId()); + charge.push_back(l1PFCand.charge()); + } + const int inputSize = maxCands_ * numInputs_; input_t input[800]; @@ -170,7 +173,7 @@ void L1MetPfProducer::CalcMlMet(const std::vector& pt, input[i] = 0; } - for (uint i = 0; i < pt.size(); i++) { + for (unsigned int i = 0; i < pt.size(); i++) { // input_cont input[i * numContInputs_] = pt[i]; input[i * numContInputs_ + 1] = eta[i]; @@ -196,148 +199,6 @@ void L1MetPfProducer::CalcMlMet(const std::vector& pt, metVector.SetEta(0); } -void L1MetPfProducer::CalcMetHLS(const std::vector& pt, - const std::vector& phi, - reco::Candidate::PolarLorentzVector& metVector) const { - pxy_t hw_px = 0; - pxy_t hw_py = 0; - pxy_t hw_sumx = 0; - pxy_t hw_sumy = 0; - - for (uint i = 0; i < pt.size(); i++) { - pt_t hw_pt = min(pt[i], maxPt_); - phi_t hw_phi = float(TVector2::Phi_mpi_pi(phi[i]) / phiLSB_); - - Project(hw_pt, hw_phi, hw_px, true); - Project(hw_pt, hw_phi, hw_py, false); - - hw_sumx = hw_sumx - hw_px; - hw_sumy = hw_sumy - hw_py; - } - - pt2_t hw_met = pt2_t(hw_sumx) * pt2_t(hw_sumx) + pt2_t(hw_sumy) * pt2_t(hw_sumy); - hw_met = sqrt(int(hw_met)); // stand-in for HLS::sqrt - - phi_t hw_met_phi = 0; - PhiFromXY(hw_sumx, hw_sumy, hw_met_phi); - - metVector.SetPt(hw_met.to_double()); - metVector.SetPhi(hw_met_phi.to_double() * phiLSB_); - metVector.SetEta(0); -} - -void L1MetPfProducer::Project(pt_t pt, phi_t phi, pxy_t& pxy, bool isX, bool debug) const { - /* - Convert pt and phi to px (py) - 1) Map phi to the first quadrant to reduce LUT size - 2) Lookup sin(phiQ1), where the result is in [0,maxPt] - which is used to encode [0,1]. - 3) Multiply pt by sin(phiQ1) to get px. Result will be px*maxPt, but - wrapping multiplication is 'mod maxPt' so the correct value is returned. - 4) Check px=-|px|. - */ - - // set phi to first quadrant - phi_t phiQ1 = (phi > 0) ? phi : phi_t(-phi); // Q1/Q4 - if (phiQ1 >= hwPiOverTwo_) - phiQ1 = hwPi_ - phiQ1; - - if (phiQ1 > hwPiOverTwo_) { - edm::LogWarning("L1MetPfProducer") << "unexpected phi (high)"; - phiQ1 = hwPiOverTwo_; - } else if (phiQ1 < 0) { - edm::LogWarning("L1MetPfProducer") << "unexpected phi (low)"; - phiQ1 = 0; - } - if (isX) { - typedef ap_ufixed<14, 12, AP_RND, AP_WRAP> pt_t; // LSB is 0.25 and max is 4 TeV - ap_ufixed cosPhi = cos(phiQ1.to_double() / hwPiOverTwo_.to_double() * M_PI / 2); - pxy = pt * cosPhi; - if (phi > hwPiOverTwo_ || phi < -hwPiOverTwo_) - pxy = -pxy; - } else { - ap_ufixed sinPhi = sin(phiQ1.to_double() / hwPiOverTwo_.to_double() * M_PI / 2); - pxy = pt * sinPhi; - if (phi < 0) - pxy = -pxy; - } -} - -void L1MetPfProducer::PhiFromXY(pxy_t px, pxy_t py, phi_t& phi, bool debug) const { - if (px == 0 && py == 0) { - phi = 0; - return; - } - if (px == 0) { - phi = py > 0 ? hwPiOverTwo_ : phi_t(-hwPiOverTwo_); - return; - } - if (py == 0) { - phi = px > 0 ? phi_t(0) : phi_t(-hwPi_); - return; - } - - // get q1 coordinates - pt_t x = px > 0 ? pt_t(px) : pt_t(-px); //px>=0 ? px : -px; - pt_t y = py > 0 ? pt_t(py) : pt_t(-py); //px>=0 ? px : -px; - // transform so a maxPt_ / dropFactor_) - b = maxPt_ / dropFactor_; - // map [0,max/4) to inv table size - int index = round((b.to_double() / (maxPt_ / dropFactor_)) * invTableSize_); - float bcheck = (float(index) / invTableSize_) * (maxPt_ / dropFactor_); - inv_t inv_b = 1. / ((float(index) / invTableSize_) * (maxPt_ / dropFactor_)); - - inv_t a_over_b = a * inv_b; - - if (debug) { - printf(" a, b = %f, %f; index, inv = %d, %f; ratio = %f \n", - a.to_double(), - b.to_double(), - index, - inv_b.to_double(), - a_over_b.to_double()); - printf(" bcheck, 1/bc = %f, %f -- %d %f %d \n", bcheck, 1. / bcheck, invTableSize_, maxPt_, dropFactor_); - } - - int atanTableBits_ = 7; - int atanTableSize_ = (1 << atanTableBits_); - index = round(a_over_b.to_double() * atanTableSize_); - phi = atan(float(index) / atanTableSize_) / phiLSB_; - - if (debug) { - printf(" atan index, phi = %d, %f (%f rad) real atan(a/b)= %f \n", - index, - phi.to_double(), - phi.to_double() * (M_PI / hwPi_.to_double()), - atan(a.to_double() / b.to_double())); - } - - // rotate from (0,pi/4) to full quad1 - if (y > x) - phi = hwPiOverTwo_ - phi; //phi = pi/2 - phi - // other quadrants - if (px < 0 && py > 0) - phi = hwPi_ - phi; // Q2 phi = pi - phi - if (px > 0 && py < 0) - phi = -phi; // Q4 phi = -phi - if (px < 0 && py < 0) - phi = -(hwPi_ - phi); // Q3 composition of both - - if (debug) { - printf(" phi hw, float, real = %f, %f (%f rad from x,y = %f, %f) \n", - phi.to_double(), - phi.to_double() * (M_PI / hwPi_.to_double()), - atan2(py.to_double(), px.to_double()), - px.to_double(), - py.to_double()); - } -} - L1MetPfProducer::~L1MetPfProducer() {} -#include "FWCore/Framework/interface/MakerMacros.h" DEFINE_FWK_MODULE(L1MetPfProducer); diff --git a/L1Trigger/Phase2L1ParticleFlow/python/l1tJUMPProducer_cfi.py b/L1Trigger/Phase2L1ParticleFlow/python/l1tJUMPProducer_cfi.py new file mode 100644 index 0000000000000..7ba1375f4b6d7 --- /dev/null +++ b/L1Trigger/Phase2L1ParticleFlow/python/l1tJUMPProducer_cfi.py @@ -0,0 +1,9 @@ +import FWCore.ParameterSet.Config as cms + +l1tJUMPProducer = cms.EDProducer("L1JUMPProducer", + RawMET = cms.InputTag("l1tMETPFProducer"), + L1PFJets = cms.InputTag("l1tSC4PFL1PuppiCorrectedEmulator"), + MinJetpT = cms.double(30), + MaxJetEta = cms.double(3.0), + JERFile = cms.string("L1Trigger/Phase2L1ParticleFlow/data/met/l1jump_jer_v1.json") +) \ No newline at end of file diff --git a/L1Trigger/Phase2L1ParticleFlow/python/l1tMETPFProducer_cfi.py b/L1Trigger/Phase2L1ParticleFlow/python/l1tMETPFProducer_cfi.py index 7f047d3dfb675..de3da0e7b7c50 100644 --- a/L1Trigger/Phase2L1ParticleFlow/python/l1tMETPFProducer_cfi.py +++ b/L1Trigger/Phase2L1ParticleFlow/python/l1tMETPFProducer_cfi.py @@ -4,6 +4,7 @@ L1PFObjects = cms.InputTag("l1tLayer1","Puppi"), maxCands = cms.int32(128), modelVersion = cms.string(""), + Poly2File = cms.FileInPath("L1Trigger/Phase2L1ParticleFlow/data/met/l1met_ptphi2pxpy_poly2_v1.json"), ) l1tMETMLProducer = cms.EDProducer("L1MetPfProducer",