From db7e90b9d176d025196bf902907fbd0e031073c7 Mon Sep 17 00:00:00 2001 From: JunwonTomOh Date: Fri, 13 Jun 2025 15:25:08 +0900 Subject: [PATCH 1/6] Update L1T MET and Add JUMP algorithm --- .../interface/jetmet/L1PFJUMPEmulator.h | 106 ++++++++ .../interface/jetmet/L1PFMetEmulator.h | 221 +++++++++++++++ .../plugins/L1JUMPProducer.cc | 124 +++++++++ .../plugins/L1MetPfProducer.cc | 256 ++++-------------- .../python/l1tJUMPProducer_cfi.py | 8 + 5 files changed, 516 insertions(+), 199 deletions(-) create mode 100644 L1Trigger/Phase2L1ParticleFlow/interface/jetmet/L1PFJUMPEmulator.h create mode 100644 L1Trigger/Phase2L1ParticleFlow/interface/jetmet/L1PFMetEmulator.h create mode 100644 L1Trigger/Phase2L1ParticleFlow/plugins/L1JUMPProducer.cc create mode 100644 L1Trigger/Phase2L1ParticleFlow/python/l1tJUMPProducer_cfi.py diff --git a/L1Trigger/Phase2L1ParticleFlow/interface/jetmet/L1PFJUMPEmulator.h b/L1Trigger/Phase2L1ParticleFlow/interface/jetmet/L1PFJUMPEmulator.h new file mode 100644 index 0000000000000..08a5969170355 --- /dev/null +++ b/L1Trigger/Phase2L1ParticleFlow/interface/jetmet/L1PFJUMPEmulator.h @@ -0,0 +1,106 @@ +#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 +#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 + */ + + inline void Get_dPt(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_par1[i] * pT + eta_par2[i] + */ + + ap_fixed<11, 1> eta_par1[5] = {0, 0.073, 0.247, 0.128, 0.091}; + ap_fixed<8, 5> eta_par2[5] = {0, 12.322, 6.061, 10.944, 12.660}; + + L1METEmu::eta_t eta_edges[4]; + float eta_boundaries[4] = {1.3, 1.7, 2.5, 3.0}; + for (uint i=0; i < 4; i++){ + eta_edges[i] = l1ct::Scales::makeGlbEta(eta_boundaries[i]); + } + + L1METEmu::eta_t abseta = abs(jet.hwEta.to_float()); + int etabin = 0; + if (abseta == 0.00) + etabin = 1; + else if (abseta < eta_edges[0]) + etabin = 1; + else if (abseta < eta_edges[1]) + etabin = 2; + else if (abseta < eta_edges[2]) + etabin = 3; + else if (abseta < eta_edges[3]) + etabin = 4; + else + etabin = 0; + + dPx_2 = 0; + dPy_2 = 0; + l1ct::Sum jet_resolution; + jet_resolution.hwPt = eta_par1[etabin] * jet.hwPt + eta_par2[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(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 (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(l1ct::Sum inMet, 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; + 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); + + return; +} + +#endif diff --git a/L1Trigger/Phase2L1ParticleFlow/interface/jetmet/L1PFMetEmulator.h b/L1Trigger/Phase2L1ParticleFlow/interface/jetmet/L1PFMetEmulator.h new file mode 100644 index 0000000000000..65f53b094b5ee --- /dev/null +++ b/L1Trigger/Phase2L1ParticleFlow/interface/jetmet/L1PFMetEmulator.h @@ -0,0 +1,221 @@ +#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" + +#ifndef CMSSW_GIT_HASH +#include "hls_math.h" +#endif + +#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; + }; + + inline Particle_xy Get_xy(l1ct::pt_t hwPt, 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 + */ + + poly_t cos2_par0[16] = {-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}; + poly_t cos2_par1[16] = {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}; + poly_t cos2_par2[16] = {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}; + + poly_t sin2_par0[16] = {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}; + poly_t sin2_par1[16] = {-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}; + poly_t sin2_par2[16] = {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}; + + phi_t phi2_edges[17]; + float phi2_points[17] = {-1.0*M_PI, -0.875*M_PI, -0.75*M_PI, -0.625*M_PI, -0.5*M_PI, -0.375*M_PI, -0.25*M_PI, -0.125*M_PI, 0.0, 0.125*M_PI, 0.25*M_PI, 0.375*M_PI, 0.5*M_PI, 0.625*M_PI, 0.75*M_PI, 0.875*M_PI, 1.0*M_PI}; + for (uint i=0; i < 17; i++){ + phi2_edges[i] = l1ct::Scales::makeGlbPhi(phi2_points[i]); + } + + int phibin = 0; + for (uint i = 0; i < 16; i++){ + if (hwPhi >= phi2_edges[i] && hwPhi < phi2_edges[i + 1]) { + phibin = i; + break; + } + } + + Particle_xy proj_xy; + + poly_t cos_var = + cos2_par0[phibin] + cos2_par1[phibin] * (hwPhi - phi2_edges[phibin]) + + cos2_par2[phibin] * (hwPhi - phi2_edges[phibin]) * (hwPhi - phi2_edges[phibin]); + poly_t sin_var = + sin2_par0[phibin] + sin2_par1[phibin] * (hwPhi - phi2_edges[phibin]) + + sin2_par2[phibin] * (hwPhi - phi2_edges[phibin]) * (hwPhi - phi2_edges[phibin]); + + proj_xy.hwPx = hwPt * cos_var; + proj_xy.hwPy = hwPt * sin_var; + + return proj_xy; + } + + inline void Sum_Particles(std::vector 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(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 + + return; + } + + inline void met_format(l1ct::Sum d, ap_uint& q) { + // Change output formats to GT formats + q = d.toGT().pack(); + } + +} // namespace L1METEmu + +inline void puppimet_emu(std::vector particles, l1ct::Sum& out_met) { + std::vector 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); + + return; +} + +#endif diff --git a/L1Trigger/Phase2L1ParticleFlow/plugins/L1JUMPProducer.cc b/L1Trigger/Phase2L1ParticleFlow/plugins/L1JUMPProducer.cc new file mode 100644 index 0000000000000..70af860bf7752 --- /dev/null +++ b/L1Trigger/Phase2L1ParticleFlow/plugins/L1JUMPProducer.cc @@ -0,0 +1,124 @@ +#include +#include +#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" + +using namespace l1t; + +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> metToken; + edm::EDGetTokenT> 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(l1t::EtSum& metVector, + std::vector& jets, + reco::Candidate::PolarLorentzVector& JUMPVector) const; + + std::vector convertEDMToHW(std::vector edmJets) const; + + double minJetPt; + double maxJetEta; +}; + +L1JUMPProducer::L1JUMPProducer(const edm::ParameterSet& cfg) + : metToken(consumes>(cfg.getParameter("RawMET"))), + jetsToken(consumes>(cfg.getParameter("L1PFJets"))), + minJetPt(cfg.getParameter("MinJetpT")), + maxJetEta(cfg.getParameter("MaxJetEta")) { + produces>(); +} + +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(l1t::EtSum& metVector, + 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(std::vector edmJets) const { + std::vector hwJets; + std::for_each(edmJets.begin(), edmJets.end(), [&](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); \ No newline at end of file diff --git a/L1Trigger/Phase2L1ParticleFlow/plugins/L1MetPfProducer.cc b/L1Trigger/Phase2L1ParticleFlow/plugins/L1MetPfProducer.cc index efcb7e78d02fc..6a1811dbee3a4 100644 --- a/L1Trigger/Phase2L1ParticleFlow/plugins/L1MetPfProducer.cc +++ b/L1Trigger/Phase2L1ParticleFlow/plugins/L1MetPfProducer.cc @@ -11,9 +11,16 @@ #include "FWCore/MessageLogger/interface/MessageLogger.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,32 +33,15 @@ 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 + typedef l1ct::pt_t pt_t; + typedef l1ct::glbphi_t phi_t; 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_); - // hls4ml emulator objects bool useMlModel_; std::shared_ptr model; @@ -63,20 +53,12 @@ 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, + void CalcMlMet(const std::vector& pfcands, reco::Candidate::PolarLorentzVector& metVector) const; }; @@ -85,7 +67,7 @@ L1MetPfProducer::L1MetPfProducer(const edm::ParameterSet& cfg) maxCands_(cfg.getParameter("maxCands")), modelVersion_(cfg.getParameter("modelVersion")) { produces>(); - useMlModel_ = (!modelVersion_.empty()); + useMlModel_ = (modelVersion_.length() > 0); if (useMlModel_) { hls4mlEmulator::ModelLoader loader(modelVersion_); model = loader.load_model(); @@ -104,38 +86,42 @@ 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 +139,26 @@ 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]; @@ -196,147 +195,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" diff --git a/L1Trigger/Phase2L1ParticleFlow/python/l1tJUMPProducer_cfi.py b/L1Trigger/Phase2L1ParticleFlow/python/l1tJUMPProducer_cfi.py new file mode 100644 index 0000000000000..d6b02f8b9aec7 --- /dev/null +++ b/L1Trigger/Phase2L1ParticleFlow/python/l1tJUMPProducer_cfi.py @@ -0,0 +1,8 @@ +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), +) \ No newline at end of file From 070951de521a48b6f556f761ab0de0162d670261 Mon Sep 17 00:00:00 2001 From: JunwonTomOh Date: Fri, 13 Jun 2025 09:00:17 +0200 Subject: [PATCH 2/6] Apply code checks and format --- .../interface/jetmet/L1PFJUMPEmulator.h | 6 +- .../interface/jetmet/L1PFMetEmulator.h | 232 ++++++++++-------- .../plugins/L1JUMPProducer.cc | 2 +- .../plugins/L1MetPfProducer.cc | 22 +- 4 files changed, 135 insertions(+), 127 deletions(-) diff --git a/L1Trigger/Phase2L1ParticleFlow/interface/jetmet/L1PFJUMPEmulator.h b/L1Trigger/Phase2L1ParticleFlow/interface/jetmet/L1PFJUMPEmulator.h index 08a5969170355..7ed55d6028227 100644 --- a/L1Trigger/Phase2L1ParticleFlow/interface/jetmet/L1PFJUMPEmulator.h +++ b/L1Trigger/Phase2L1ParticleFlow/interface/jetmet/L1PFJUMPEmulator.h @@ -36,7 +36,7 @@ namespace L1JUMPEmu { L1METEmu::eta_t eta_edges[4]; float eta_boundaries[4] = {1.3, 1.7, 2.5, 3.0}; - for (uint i=0; i < 4; i++){ + for (uint i = 0; i < 4; i++) { eta_edges[i] = l1ct::Scales::makeGlbEta(eta_boundaries[i]); } @@ -64,7 +64,6 @@ namespace L1JUMPEmu { dPx_2 = dpt_xy.hwPx * dpt_xy.hwPx; dPy_2 = dpt_xy.hwPy * dpt_xy.hwPy; - } inline void Met_dPt(std::vector jets, L1METEmu::proj2_t& dPx_2, L1METEmu::proj2_t& dPy_2) { @@ -74,7 +73,7 @@ namespace L1JUMPEmu { L1METEmu::proj2_t sum_dPx2 = 0; L1METEmu::proj2_t sum_dPy2 = 0; - for (uint i = 0; i < jets.size(); i++) { + for (uint i = 0; i < jets.size(); i++) { Get_dPt(jets[i], each_dPx2, each_dPy2); sum_dPx2 += each_dPx2; sum_dPy2 += each_dPy2; @@ -86,7 +85,6 @@ namespace L1JUMPEmu { } // namespace L1JUMPEmu inline void JUMP_emu(l1ct::Sum inMet, std::vector jets, l1ct::Sum& outMet) { - L1METEmu::Particle_xy inMet_xy = L1METEmu::Get_xy(inMet.hwPt, inMet.hwPhi); L1METEmu::proj2_t dPx_2; diff --git a/L1Trigger/Phase2L1ParticleFlow/interface/jetmet/L1PFMetEmulator.h b/L1Trigger/Phase2L1ParticleFlow/interface/jetmet/L1PFMetEmulator.h index 65f53b094b5ee..cc4fddbfc15c5 100644 --- a/L1Trigger/Phase2L1ParticleFlow/interface/jetmet/L1PFMetEmulator.h +++ b/L1Trigger/Phase2L1ParticleFlow/interface/jetmet/L1PFMetEmulator.h @@ -44,125 +44,139 @@ namespace L1METEmu { */ poly_t cos2_par0[16] = {-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}; + -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}; poly_t cos2_par1[16] = {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}; + 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}; poly_t cos2_par2[16] = {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}; + 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}; poly_t sin2_par0[16] = {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}; + -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}; poly_t sin2_par1[16] = {-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}; + -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}; poly_t sin2_par2[16] = {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}; + 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}; phi_t phi2_edges[17]; - float phi2_points[17] = {-1.0*M_PI, -0.875*M_PI, -0.75*M_PI, -0.625*M_PI, -0.5*M_PI, -0.375*M_PI, -0.25*M_PI, -0.125*M_PI, 0.0, 0.125*M_PI, 0.25*M_PI, 0.375*M_PI, 0.5*M_PI, 0.625*M_PI, 0.75*M_PI, 0.875*M_PI, 1.0*M_PI}; - for (uint i=0; i < 17; i++){ + float phi2_points[17] = {-1.0 * M_PI, + -0.875 * M_PI, + -0.75 * M_PI, + -0.625 * M_PI, + -0.5 * M_PI, + -0.375 * M_PI, + -0.25 * M_PI, + -0.125 * M_PI, + 0.0, + 0.125 * M_PI, + 0.25 * M_PI, + 0.375 * M_PI, + 0.5 * M_PI, + 0.625 * M_PI, + 0.75 * M_PI, + 0.875 * M_PI, + 1.0 * M_PI}; + for (uint i = 0; i < 17; i++) { phi2_edges[i] = l1ct::Scales::makeGlbPhi(phi2_points[i]); } int phibin = 0; - for (uint i = 0; i < 16; i++){ + for (uint i = 0; i < 16; i++) { if (hwPhi >= phi2_edges[i] && hwPhi < phi2_edges[i + 1]) { phibin = i; break; - } } + } Particle_xy proj_xy; - poly_t cos_var = - cos2_par0[phibin] + cos2_par1[phibin] * (hwPhi - phi2_edges[phibin]) + - cos2_par2[phibin] * (hwPhi - phi2_edges[phibin]) * (hwPhi - phi2_edges[phibin]); - poly_t sin_var = - sin2_par0[phibin] + sin2_par1[phibin] * (hwPhi - phi2_edges[phibin]) + - sin2_par2[phibin] * (hwPhi - phi2_edges[phibin]) * (hwPhi - phi2_edges[phibin]); + poly_t cos_var = cos2_par0[phibin] + cos2_par1[phibin] * (hwPhi - phi2_edges[phibin]) + + cos2_par2[phibin] * (hwPhi - phi2_edges[phibin]) * (hwPhi - phi2_edges[phibin]); + poly_t sin_var = sin2_par0[phibin] + sin2_par1[phibin] * (hwPhi - phi2_edges[phibin]) + + sin2_par2[phibin] * (hwPhi - phi2_edges[phibin]) * (hwPhi - phi2_edges[phibin]); proj_xy.hwPx = hwPt * cos_var; proj_xy.hwPy = hwPt * sin_var; @@ -173,26 +187,26 @@ namespace L1METEmu { inline void Sum_Particles(std::vector 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(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 + +#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.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 - +#endif + return; } @@ -210,7 +224,7 @@ inline void puppimet_emu(std::vector particles, l1ct::Sum& ou 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); diff --git a/L1Trigger/Phase2L1ParticleFlow/plugins/L1JUMPProducer.cc b/L1Trigger/Phase2L1ParticleFlow/plugins/L1JUMPProducer.cc index 70af860bf7752..4bf871bf30d43 100644 --- a/L1Trigger/Phase2L1ParticleFlow/plugins/L1JUMPProducer.cc +++ b/L1Trigger/Phase2L1ParticleFlow/plugins/L1JUMPProducer.cc @@ -28,7 +28,7 @@ using namespace l1t; 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 diff --git a/L1Trigger/Phase2L1ParticleFlow/plugins/L1MetPfProducer.cc b/L1Trigger/Phase2L1ParticleFlow/plugins/L1MetPfProducer.cc index 6a1811dbee3a4..22e815f6c4865 100644 --- a/L1Trigger/Phase2L1ParticleFlow/plugins/L1MetPfProducer.cc +++ b/L1Trigger/Phase2L1ParticleFlow/plugins/L1MetPfProducer.cc @@ -39,8 +39,8 @@ class L1MetPfProducer : public edm::global::EDProducer<> { // quantization controllers typedef l1ct::pt_t pt_t; - typedef l1ct::glbphi_t phi_t; - static constexpr float phiLSB_ = M_PI / 720; // rad + typedef l1ct::glbphi_t phi_t; + static constexpr float phiLSB_ = M_PI / 720; // rad // hls4ml emulator objects bool useMlModel_; @@ -53,13 +53,11 @@ class L1MetPfProducer : public edm::global::EDProducer<> { static constexpr int numCatInputs_ = 2; static constexpr int numInputs_ = numContInputs_ + numPxPyInputs_ + numCatInputs_; - void CalcMetHLS(const std::vector& pfcands, - reco::Candidate::PolarLorentzVector& metVector) const; + void CalcMetHLS(const std::vector& pfcands, reco::Candidate::PolarLorentzVector& metVector) const; int EncodePdgId(int pdgId) const; - void CalcMlMet(const std::vector& pfcands, - reco::Candidate::PolarLorentzVector& metVector) const; + void CalcMlMet(const std::vector& pfcands, reco::Candidate::PolarLorentzVector& metVector) const; }; L1MetPfProducer::L1MetPfProducer(const edm::ParameterSet& cfg) @@ -67,7 +65,7 @@ L1MetPfProducer::L1MetPfProducer(const edm::ParameterSet& cfg) maxCands_(cfg.getParameter("maxCands")), modelVersion_(cfg.getParameter("modelVersion")) { produces>(); - useMlModel_ = (modelVersion_.length() > 0); + useMlModel_ = (!modelVersion_.empty()); if (useMlModel_) { hls4mlEmulator::ModelLoader loader(modelVersion_); model = loader.load_model(); @@ -102,12 +100,12 @@ void L1MetPfProducer::produce(edm::StreamID, edm::Event& iEvent, const edm::Even iEvent.put(std::move(metCollection)); } - void L1MetPfProducer::CalcMetHLS(const std::vector& pfcands, - reco::Candidate::PolarLorentzVector& metVector) const { +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++){ + 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()); @@ -119,8 +117,7 @@ void L1MetPfProducer::produce(edm::StreamID, edm::Event& iEvent, const edm::Even 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)) { @@ -141,7 +138,6 @@ int L1MetPfProducer::EncodePdgId(int pdgId) const { void L1MetPfProducer::CalcMlMet(const std::vector& pfcands, reco::Candidate::PolarLorentzVector& metVector) const { - std::vector pt; std::vector eta; std::vector phi; From 104769c67b15f252baaa6864962d051a686c6ccf Mon Sep 17 00:00:00 2001 From: JunwonTomOh Date: Thu, 21 Aug 2025 04:18:45 +0200 Subject: [PATCH 3/6] Fix the code to read the factors from cms-data --- L1Trigger/Phase2L1ParticleFlow/BuildFile.xml | 1 + .../interface/jetmet/L1PFJUMPEmulator.h | 63 ++++-- .../interface/jetmet/L1PFMetEmulator.h | 184 +++++------------- .../plugins/L1JUMPProducer.cc | 16 +- 4 files changed, 101 insertions(+), 163 deletions(-) diff --git a/L1Trigger/Phase2L1ParticleFlow/BuildFile.xml b/L1Trigger/Phase2L1ParticleFlow/BuildFile.xml index e45b710a713f9..b104906011606 100644 --- a/L1Trigger/Phase2L1ParticleFlow/BuildFile.xml +++ b/L1Trigger/Phase2L1ParticleFlow/BuildFile.xml @@ -16,6 +16,7 @@ + diff --git a/L1Trigger/Phase2L1ParticleFlow/interface/jetmet/L1PFJUMPEmulator.h b/L1Trigger/Phase2L1ParticleFlow/interface/jetmet/L1PFJUMPEmulator.h index 7ed55d6028227..6a10d6f148761 100644 --- a/L1Trigger/Phase2L1ParticleFlow/interface/jetmet/L1PFJUMPEmulator.h +++ b/L1Trigger/Phase2L1ParticleFlow/interface/jetmet/L1PFJUMPEmulator.h @@ -8,9 +8,16 @@ #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" @@ -22,35 +29,55 @@ namespace L1JUMPEmu { - Approximate L1 Jet energy resolution by pT, eta value - Apply the estimated resolution to MET */ + struct JER_param { + std::array, 5> par0; // eta.par0 (slope) + std::array, 5> par1; // eta.par1 (offset) + std::array 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()); + t.par1[i] = ap_fixed<8, 5>(j["eta"]["par1"][i].get()); + } + for (int i = 0; i < 4; ++i) { + t.edges[i] = l1ct::Scales::makeGlbEta(j["eta_edges"][i].get()); + } + return t; + }(); + return P; + } - inline void Get_dPt(l1ct::Jet jet, L1METEmu::proj2_t& dPx_2, L1METEmu::proj2_t& dPy_2) { + 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_par1[i] * pT + eta_par2[i] + - σ(pT) ≈ eta_par0[i] * pT + eta_par1[i] */ - ap_fixed<11, 1> eta_par1[5] = {0, 0.073, 0.247, 0.128, 0.091}; - ap_fixed<8, 5> eta_par2[5] = {0, 12.322, 6.061, 10.944, 12.660}; - - L1METEmu::eta_t eta_edges[4]; - float eta_boundaries[4] = {1.3, 1.7, 2.5, 3.0}; - for (uint i = 0; i < 4; i++) { - eta_edges[i] = l1ct::Scales::makeGlbEta(eta_boundaries[i]); - } + const auto& J = Get_jer_param(); L1METEmu::eta_t abseta = abs(jet.hwEta.to_float()); int etabin = 0; if (abseta == 0.00) etabin = 1; - else if (abseta < eta_edges[0]) + else if (abseta < J.edges[0]) etabin = 1; - else if (abseta < eta_edges[1]) + else if (abseta < J.edges[1]) etabin = 2; - else if (abseta < eta_edges[2]) + else if (abseta < J.edges[2]) etabin = 3; - else if (abseta < eta_edges[3]) + else if (abseta < J.edges[3]) etabin = 4; else etabin = 0; @@ -58,7 +85,7 @@ namespace L1JUMPEmu { dPx_2 = 0; dPy_2 = 0; l1ct::Sum jet_resolution; - jet_resolution.hwPt = eta_par1[etabin] * jet.hwPt + eta_par2[etabin]; + 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); @@ -66,7 +93,7 @@ namespace L1JUMPEmu { dPy_2 = dpt_xy.hwPy * dpt_xy.hwPy; } - inline void Met_dPt(std::vector jets, L1METEmu::proj2_t& dPx_2, L1METEmu::proj2_t& dPy_2) { + 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; @@ -84,7 +111,7 @@ namespace L1JUMPEmu { } } // namespace L1JUMPEmu -inline void JUMP_emu(l1ct::Sum inMet, std::vector jets, l1ct::Sum& outMet) { +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; @@ -97,8 +124,6 @@ inline void JUMP_emu(l1ct::Sum inMet, std::vector jets, l1ct::Sum& ou 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); - - return; } #endif diff --git a/L1Trigger/Phase2L1ParticleFlow/interface/jetmet/L1PFMetEmulator.h b/L1Trigger/Phase2L1ParticleFlow/interface/jetmet/L1PFMetEmulator.h index cc4fddbfc15c5..127a9b2dbc76a 100644 --- a/L1Trigger/Phase2L1ParticleFlow/interface/jetmet/L1PFMetEmulator.h +++ b/L1Trigger/Phase2L1ParticleFlow/interface/jetmet/L1PFMetEmulator.h @@ -7,11 +7,17 @@ #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 +#include +#include #include #include #include "ap_int.h" @@ -35,7 +41,38 @@ namespace L1METEmu { proj_t hwPy; }; - inline Particle_xy Get_xy(l1ct::pt_t hwPt, l1ct::glbphi_t hwPhi) { + struct Poly2_param { + std::array c0, c1, c2, s0, s1, s2; + std::array 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()); + 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 (int i = 0; i < 17; ++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 @@ -43,129 +80,10 @@ namespace L1METEmu { - Fitting the value with 2nd order function */ - poly_t cos2_par0[16] = {-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}; - poly_t cos2_par1[16] = {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}; - poly_t cos2_par2[16] = {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}; - - poly_t sin2_par0[16] = {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}; - poly_t sin2_par1[16] = {-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}; - poly_t sin2_par2[16] = {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}; - - phi_t phi2_edges[17]; - float phi2_points[17] = {-1.0 * M_PI, - -0.875 * M_PI, - -0.75 * M_PI, - -0.625 * M_PI, - -0.5 * M_PI, - -0.375 * M_PI, - -0.25 * M_PI, - -0.125 * M_PI, - 0.0, - 0.125 * M_PI, - 0.25 * M_PI, - 0.375 * M_PI, - 0.5 * M_PI, - 0.625 * M_PI, - 0.75 * M_PI, - 0.875 * M_PI, - 1.0 * M_PI}; - for (uint i = 0; i < 17; i++) { - phi2_edges[i] = l1ct::Scales::makeGlbPhi(phi2_points[i]); - } - + const auto& P = L1METEmu::Get_Poly2_param(); int phibin = 0; - for (uint i = 0; i < 16; i++) { - if (hwPhi >= phi2_edges[i] && hwPhi < phi2_edges[i + 1]) { + for (int i = 0; i < 16; i++) { + if (hwPhi >= P.phi_edges[i] && hwPhi < P.phi_edges[i + 1]) { phibin = i; break; } @@ -173,10 +91,10 @@ namespace L1METEmu { Particle_xy proj_xy; - poly_t cos_var = cos2_par0[phibin] + cos2_par1[phibin] * (hwPhi - phi2_edges[phibin]) + - cos2_par2[phibin] * (hwPhi - phi2_edges[phibin]) * (hwPhi - phi2_edges[phibin]); - poly_t sin_var = sin2_par0[phibin] + sin2_par1[phibin] * (hwPhi - phi2_edges[phibin]) + - sin2_par2[phibin] * (hwPhi - phi2_edges[phibin]) * (hwPhi - phi2_edges[phibin]); + 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; @@ -184,7 +102,7 @@ namespace L1METEmu { return proj_xy; } - inline void Sum_Particles(std::vector particles_xy, Particle_xy& met_xy) { + inline void Sum_Particles(const std::vector particles_xy, Particle_xy& met_xy) { met_xy.hwPx = 0; met_xy.hwPy = 0; @@ -194,7 +112,7 @@ namespace L1METEmu { } } - inline void pxpy_to_ptphi(Particle_xy met_xy, l1ct::Sum& hls_met) { + 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(); @@ -206,18 +124,16 @@ namespace L1METEmu { 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 - - return; } - inline void met_format(l1ct::Sum d, ap_uint& q) { + 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(std::vector particles, l1ct::Sum& out_met) { +inline void puppimet_emu(const std::vector particles, l1ct::Sum& out_met) { std::vector particles_xy; for (uint i = 0; i < particles.size(); i++) { @@ -228,8 +144,6 @@ inline void puppimet_emu(std::vector particles, l1ct::Sum& ou L1METEmu::Particle_xy met_xy; L1METEmu::Sum_Particles(particles_xy, met_xy); L1METEmu::pxpy_to_ptphi(met_xy, out_met); - - return; } #endif diff --git a/L1Trigger/Phase2L1ParticleFlow/plugins/L1JUMPProducer.cc b/L1Trigger/Phase2L1ParticleFlow/plugins/L1JUMPProducer.cc index 4bf871bf30d43..49c5a3a62d104 100644 --- a/L1Trigger/Phase2L1ParticleFlow/plugins/L1JUMPProducer.cc +++ b/L1Trigger/Phase2L1ParticleFlow/plugins/L1JUMPProducer.cc @@ -25,8 +25,6 @@ #include "L1Trigger/Phase2L1ParticleFlow/interface/jetmet/L1PFJUMPEmulator.h" -using namespace l1t; - class L1JUMPProducer : public edm::global::EDProducer<> { /* Producer for the JUMP Algorithm @@ -50,11 +48,11 @@ class L1JUMPProducer : public edm::global::EDProducer<> { static constexpr float phiLSB_ = M_PI / 720; static constexpr float maxPt_ = ((1 << pt_t::width) - 1) * ptLSB_; - void CalcJUMP_HLS(l1t::EtSum& metVector, - std::vector& jets, + void CalcJUMP_HLS(const l1t::EtSum& metVector, + const std::vector& jets, reco::Candidate::PolarLorentzVector& JUMPVector) const; - std::vector convertEDMToHW(std::vector edmJets) const; + std::vector convertEDMToHW(const std::vector edmJets) const; double minJetPt; double maxJetEta; @@ -92,8 +90,8 @@ void L1JUMPProducer::produce(edm::StreamID, edm::Event& iEvent, const edm::Event iEvent.put(std::move(JUMPCollection)); } -void L1JUMPProducer::CalcJUMP_HLS(l1t::EtSum& metVector, - std::vector& jets, +void L1JUMPProducer::CalcJUMP_HLS(const l1t::EtSum& metVector, + const std::vector& jets, reco::Candidate::PolarLorentzVector& outMet_Vector) const { // JUMP Calculate l1ct::Sum inMet; @@ -109,9 +107,9 @@ void L1JUMPProducer::CalcJUMP_HLS(l1t::EtSum& metVector, outMet_Vector.SetEta(0); } -std::vector L1JUMPProducer::convertEDMToHW(std::vector edmJets) const { +std::vector L1JUMPProducer::convertEDMToHW(const std::vector edmJets) const { std::vector hwJets; - std::for_each(edmJets.begin(), edmJets.end(), [&](l1t::PFJet jet) { + std::for_each(edmJets.begin(), edmJets.end(), [&](const l1t::PFJet jet) { l1ct::Jet hwJet = l1ct::Jet::unpack(jet.getHWJetCT()); hwJets.push_back(hwJet); }); From 9a8ce1869417edf98b69dfaec35d627145f7102a Mon Sep 17 00:00:00 2001 From: JunwonTomOh Date: Mon, 1 Sep 2025 15:20:13 +0200 Subject: [PATCH 4/6] Modified the JUMP emulator --- L1Trigger/Phase2L1ParticleFlow/BuildFile.xml | 1 - .../interface/jetmet/L1PFJUMPEmulator.h | 18 ++++++------------ 2 files changed, 6 insertions(+), 13 deletions(-) diff --git a/L1Trigger/Phase2L1ParticleFlow/BuildFile.xml b/L1Trigger/Phase2L1ParticleFlow/BuildFile.xml index b104906011606..e45b710a713f9 100644 --- a/L1Trigger/Phase2L1ParticleFlow/BuildFile.xml +++ b/L1Trigger/Phase2L1ParticleFlow/BuildFile.xml @@ -16,7 +16,6 @@ - diff --git a/L1Trigger/Phase2L1ParticleFlow/interface/jetmet/L1PFJUMPEmulator.h b/L1Trigger/Phase2L1ParticleFlow/interface/jetmet/L1PFJUMPEmulator.h index 6a10d6f148761..0e3e2ec4f7506 100644 --- a/L1Trigger/Phase2L1ParticleFlow/interface/jetmet/L1PFJUMPEmulator.h +++ b/L1Trigger/Phase2L1ParticleFlow/interface/jetmet/L1PFJUMPEmulator.h @@ -69,18 +69,12 @@ namespace L1JUMPEmu { L1METEmu::eta_t abseta = abs(jet.hwEta.to_float()); int etabin = 0; - if (abseta == 0.00) - etabin = 1; - else if (abseta < J.edges[0]) - etabin = 1; - else if (abseta < J.edges[1]) - etabin = 2; - else if (abseta < J.edges[2]) - etabin = 3; - else if (abseta < J.edges[3]) - etabin = 4; - else - etabin = 0; + for (uint i = 0; i < 4;) { + if (abseta < J.edges[i]) { + etabin = i+1; + break; + } + } dPx_2 = 0; dPy_2 = 0; From 88d6bf5ea35b98871e09d3f72abf0ac0acba3a41 Mon Sep 17 00:00:00 2001 From: JunwonTomOh Date: Mon, 1 Sep 2025 15:30:47 +0200 Subject: [PATCH 5/6] Apply code format --- .../Phase2L1ParticleFlow/interface/jetmet/L1PFJUMPEmulator.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/L1Trigger/Phase2L1ParticleFlow/interface/jetmet/L1PFJUMPEmulator.h b/L1Trigger/Phase2L1ParticleFlow/interface/jetmet/L1PFJUMPEmulator.h index 0e3e2ec4f7506..21bbee23c9e40 100644 --- a/L1Trigger/Phase2L1ParticleFlow/interface/jetmet/L1PFJUMPEmulator.h +++ b/L1Trigger/Phase2L1ParticleFlow/interface/jetmet/L1PFJUMPEmulator.h @@ -71,7 +71,7 @@ namespace L1JUMPEmu { int etabin = 0; for (uint i = 0; i < 4;) { if (abseta < J.edges[i]) { - etabin = i+1; + etabin = i + 1; break; } } From af82bf2ca97c7bb94e56f486bafb654a450f9bb9 Mon Sep 17 00:00:00 2001 From: JunwonTomOh Date: Tue, 9 Sep 2025 17:43:04 +0200 Subject: [PATCH 6/6] Change emulators to define the json file directory --- .../interface/jetmet/L1PFJUMPEmulator.h | 62 ++++++++++++++----- .../interface/jetmet/L1PFMetEmulator.h | 61 ++++++++++++++---- .../plugins/L1JUMPProducer.cc | 6 +- .../plugins/L1MetPfProducer.cc | 11 +++- .../python/l1tJUMPProducer_cfi.py | 1 + .../python/l1tMETPFProducer_cfi.py | 1 + 6 files changed, 109 insertions(+), 33 deletions(-) diff --git a/L1Trigger/Phase2L1ParticleFlow/interface/jetmet/L1PFJUMPEmulator.h b/L1Trigger/Phase2L1ParticleFlow/interface/jetmet/L1PFJUMPEmulator.h index 21bbee23c9e40..2d66892906e40 100644 --- a/L1Trigger/Phase2L1ParticleFlow/interface/jetmet/L1PFJUMPEmulator.h +++ b/L1Trigger/Phase2L1ParticleFlow/interface/jetmet/L1PFJUMPEmulator.h @@ -1,5 +1,5 @@ -#ifndef L1Trigger_Phase2L1ParticleFlow_JUMP_h -#define L1Trigger_Phase2L1ParticleFlow_JUMP_h +#ifndef L1Trigger_Phase2L1ParticleFlow_JUMPEmulator_h +#define L1Trigger_Phase2L1ParticleFlow_JUMPEmulator_h #include "DataFormats/L1TParticleFlow/interface/jets.h" #include "DataFormats/L1TParticleFlow/interface/sums.h" @@ -30,26 +30,56 @@ namespace L1JUMPEmu { - Apply the estimated resolution to MET */ struct JER_param { - std::array, 5> par0; // eta.par0 (slope) - std::array, 5> par1; // eta.par1 (offset) - std::array edges; // |eta| boundaries in HW scale + 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{}; - edm::FileInPath fip("L1Trigger/Phase2L1ParticleFlow/data/met/l1jump_jer_v1.json"); - std::ifstream in(fip.fullPath()); - if (!in) - throw cms::Exception("FileNotFound") << fip.fullPath(); + 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 + std::ifstream in(path); + if (!in) { + throw std::runtime_error(std::string("File not found: ") + path); + } +#endif + nlohmann::json j; in >> j; - for (int i = 0; i < 5; ++i) { + 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 (int i = 0; i < 4; ++i) { + for (unsigned int i = 0; i < N; ++i) { t.edges[i] = l1ct::Scales::makeGlbEta(j["eta_edges"][i].get()); } return t; @@ -68,8 +98,8 @@ namespace L1JUMPEmu { const auto& J = Get_jer_param(); L1METEmu::eta_t abseta = abs(jet.hwEta.to_float()); - int etabin = 0; - for (uint i = 0; i < 4;) { + unsigned int etabin = 0; + for (unsigned int i = 0; i < J.eta_bins; i++) { if (abseta < J.edges[i]) { etabin = i + 1; break; @@ -87,14 +117,14 @@ namespace L1JUMPEmu { 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) { + 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 (uint i = 0; i < jets.size(); i++) { + 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; @@ -105,7 +135,7 @@ namespace L1JUMPEmu { } } // namespace L1JUMPEmu -inline void JUMP_emu(const l1ct::Sum inMet, const std::vector jets, l1ct::Sum& outMet) { +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; diff --git a/L1Trigger/Phase2L1ParticleFlow/interface/jetmet/L1PFMetEmulator.h b/L1Trigger/Phase2L1ParticleFlow/interface/jetmet/L1PFMetEmulator.h index 127a9b2dbc76a..7b295160e0d5a 100644 --- a/L1Trigger/Phase2L1ParticleFlow/interface/jetmet/L1PFMetEmulator.h +++ b/L1Trigger/Phase2L1ParticleFlow/interface/jetmet/L1PFMetEmulator.h @@ -1,5 +1,5 @@ -#ifndef L1Trigger_Phase2L1ParticleFlow_MET_h -#define L1Trigger_Phase2L1ParticleFlow_MET_h +#ifndef L1Trigger_Phase2L1ParticleFlow_METEmulator_h +#define L1Trigger_Phase2L1ParticleFlow_METEmulator_h #include "DataFormats/L1TParticleFlow/interface/jets.h" #include "DataFormats/L1TParticleFlow/interface/sums.h" @@ -42,21 +42,55 @@ namespace L1METEmu { }; struct Poly2_param { - std::array c0, c1, c2, s0, s1, s2; - std::array phi_edges; + 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{}; - edm::FileInPath f("L1Trigger/Phase2L1ParticleFlow/data/met/l1met_ptphi2pxpy_poly2_v1.json"); + std::string path = poly2_path_config().path; + +#ifdef CMSSW_GIT_HASH + edm::FileInPath f(path); std::ifstream in(f.fullPath()); - if (!in) + if (!in) { throw cms::Exception("FileNotFound") << f.fullPath(); + } +#else + std::ifstream in(path); + if (!in) { + throw std::runtime_error(std::string("File not found: ") + path); + } +#endif + nlohmann::json j; in >> j; - for (int i = 0; i < 16; ++i) { + 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()); @@ -64,7 +98,7 @@ namespace L1METEmu { t.s1[i] = poly_t(j["sin"]["par1"][i].get()); t.s2[i] = poly_t(j["sin"]["par2"][i].get()); } - for (int i = 0; i < 17; ++i) { + 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; @@ -82,7 +116,8 @@ namespace L1METEmu { const auto& P = L1METEmu::Get_Poly2_param(); int phibin = 0; - for (int i = 0; i < 16; i++) { + + 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; @@ -102,11 +137,11 @@ namespace L1METEmu { return proj_xy; } - inline void Sum_Particles(const std::vector particles_xy, Particle_xy& met_xy) { + inline void Sum_Particles(const std::vector& particles_xy, Particle_xy& met_xy) { met_xy.hwPx = 0; met_xy.hwPy = 0; - for (uint i = 0; i < particles_xy.size(); ++i) { + for (unsigned int i = 0; i < particles_xy.size(); ++i) { met_xy.hwPx -= particles_xy[i].hwPx; met_xy.hwPy -= particles_xy[i].hwPy; } @@ -133,10 +168,10 @@ namespace L1METEmu { } // namespace L1METEmu -inline void puppimet_emu(const std::vector particles, l1ct::Sum& out_met) { +inline void puppimet_emu(const std::vector& particles, l1ct::Sum& out_met) { std::vector particles_xy; - for (uint i = 0; i < particles.size(); i++) { + 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); } diff --git a/L1Trigger/Phase2L1ParticleFlow/plugins/L1JUMPProducer.cc b/L1Trigger/Phase2L1ParticleFlow/plugins/L1JUMPProducer.cc index 49c5a3a62d104..c5430083a366c 100644 --- a/L1Trigger/Phase2L1ParticleFlow/plugins/L1JUMPProducer.cc +++ b/L1Trigger/Phase2L1ParticleFlow/plugins/L1JUMPProducer.cc @@ -56,14 +56,17 @@ class L1JUMPProducer : public edm::global::EDProducer<> { 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")) { + maxJetEta(cfg.getParameter("MaxJetEta")), + jerFilePath_(cfg.getParameter("JERFile")) { produces>(); + L1JUMPEmu::SetJERFile(jerFilePath_); } void L1JUMPProducer::produce(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const { @@ -118,5 +121,4 @@ std::vector L1JUMPProducer::convertEDMToHW(const std::vector { 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; @@ -69,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()); } } @@ -77,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); } @@ -165,7 +173,7 @@ void L1MetPfProducer::CalcMlMet(const std::vector& pfcands, 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]; @@ -193,5 +201,4 @@ void L1MetPfProducer::CalcMlMet(const std::vector& pfcands, 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 index d6b02f8b9aec7..df8b34a13b5aa 100644 --- a/L1Trigger/Phase2L1ParticleFlow/python/l1tJUMPProducer_cfi.py +++ b/L1Trigger/Phase2L1ParticleFlow/python/l1tJUMPProducer_cfi.py @@ -5,4 +5,5 @@ 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",