File indexing completed on 2023-03-17 11:13:13
0001 #include <vector>
0002 #include <ap_int.h>
0003 #include <ap_fixed.h>
0004 #include <TVector2.h>
0005
0006 #include "FWCore/Framework/interface/global/EDProducer.h"
0007 #include "FWCore/Framework/interface/Event.h"
0008 #include "FWCore/Framework/interface/MakerMacros.h"
0009 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0010 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0011
0012 #include "DataFormats/L1TParticleFlow/interface/PFCandidate.h"
0013 #include "DataFormats/L1Trigger/interface/EtSum.h"
0014 #include "DataFormats/Math/interface/LorentzVector.h"
0015
0016 using namespace l1t;
0017
0018 class L1MetPfProducer : public edm::global::EDProducer<> {
0019 public:
0020 explicit L1MetPfProducer(const edm::ParameterSet&);
0021 ~L1MetPfProducer() override;
0022
0023 private:
0024 void produce(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const override;
0025 edm::EDGetTokenT<vector<l1t::PFCandidate>> _l1PFToken;
0026
0027 int maxCands_ = 128;
0028
0029
0030 typedef ap_ufixed<14, 12, AP_RND, AP_WRAP> pt_t;
0031 typedef ap_int<12> phi_t;
0032 const float ptLSB_ = 0.25;
0033 const float phiLSB_ = M_PI / 720;
0034
0035
0036 typedef ap_fixed<pt_t::width + 1, pt_t::iwidth + 1, AP_RND, AP_SAT> pxy_t;
0037 typedef ap_fixed<2 * pt_t::width, 2 * pt_t::iwidth, AP_RND, AP_SAT> pt2_t;
0038
0039 const float maxPt_ = ((1 << pt_t::width) - 1) * ptLSB_;
0040 const phi_t hwPi_ = round(M_PI / phiLSB_);
0041 const phi_t hwPiOverTwo_ = round(M_PI / (2 * phiLSB_));
0042
0043 typedef ap_ufixed<pt_t::width, 0> inv_t;
0044
0045
0046 const int dropBits_ = 2;
0047 const int dropFactor_ = (1 << dropBits_);
0048 const int invTableBits_ = 10;
0049 const int invTableSize_ = (1 << invTableBits_);
0050
0051 void Project(pt_t pt, phi_t phi, pxy_t& pxy, bool isX, bool debug = false) const;
0052 void PhiFromXY(pxy_t px, pxy_t py, phi_t& phi, bool debug = false) const;
0053
0054 void CalcMetHLS(std::vector<float> pt, std::vector<float> phi, reco::Candidate::PolarLorentzVector& metVector) const;
0055 };
0056
0057 L1MetPfProducer::L1MetPfProducer(const edm::ParameterSet& cfg)
0058 : _l1PFToken(consumes<std::vector<l1t::PFCandidate>>(cfg.getParameter<edm::InputTag>("L1PFObjects"))),
0059 maxCands_(cfg.getParameter<int>("maxCands")) {
0060 produces<std::vector<l1t::EtSum>>();
0061 }
0062
0063 void L1MetPfProducer::produce(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const {
0064 edm::Handle<l1t::PFCandidateCollection> l1PFCandidates;
0065 iEvent.getByToken(_l1PFToken, l1PFCandidates);
0066
0067 std::vector<float> pt;
0068 std::vector<float> phi;
0069
0070 for (int i = 0; i < int(l1PFCandidates->size()) && (i < maxCands_ || maxCands_ < 0); i++) {
0071 const auto& l1PFCand = l1PFCandidates->at(i);
0072 pt.push_back(l1PFCand.pt());
0073 phi.push_back(l1PFCand.phi());
0074 }
0075
0076 reco::Candidate::PolarLorentzVector metVector;
0077
0078 CalcMetHLS(pt, phi, metVector);
0079
0080 l1t::EtSum theMET(metVector, l1t::EtSum::EtSumType::kTotalHt, 0, 0, 0, 0);
0081
0082 std::unique_ptr<std::vector<l1t::EtSum>> metCollection(new std::vector<l1t::EtSum>(0));
0083 metCollection->push_back(theMET);
0084 iEvent.put(std::move(metCollection));
0085 }
0086
0087 void L1MetPfProducer::CalcMetHLS(std::vector<float> pt,
0088 std::vector<float> phi,
0089 reco::Candidate::PolarLorentzVector& metVector) const {
0090 pxy_t hw_px = 0;
0091 pxy_t hw_py = 0;
0092 pxy_t hw_sumx = 0;
0093 pxy_t hw_sumy = 0;
0094
0095 for (uint i = 0; i < pt.size(); i++) {
0096 pt_t hw_pt = min(pt[i], maxPt_);
0097 phi_t hw_phi = float(TVector2::Phi_mpi_pi(phi[i]) / phiLSB_);
0098
0099 Project(hw_pt, hw_phi, hw_px, true);
0100 Project(hw_pt, hw_phi, hw_py, false);
0101
0102 hw_sumx = hw_sumx - hw_px;
0103 hw_sumy = hw_sumy - hw_py;
0104 }
0105
0106 pt2_t hw_met = pt2_t(hw_sumx) * pt2_t(hw_sumx) + pt2_t(hw_sumy) * pt2_t(hw_sumy);
0107 hw_met = sqrt(int(hw_met));
0108
0109 phi_t hw_met_phi = 0;
0110 PhiFromXY(hw_sumx, hw_sumy, hw_met_phi);
0111
0112 metVector.SetPt(hw_met.to_double());
0113 metVector.SetPhi(hw_met_phi.to_double() * phiLSB_);
0114 metVector.SetEta(0);
0115 }
0116
0117 void L1MetPfProducer::Project(pt_t pt, phi_t phi, pxy_t& pxy, bool isX, bool debug) const {
0118
0119
0120
0121
0122
0123
0124
0125
0126
0127
0128
0129 phi_t phiQ1 = (phi > 0) ? phi : phi_t(-phi);
0130 if (phiQ1 >= hwPiOverTwo_)
0131 phiQ1 = hwPi_ - phiQ1;
0132
0133 if (phiQ1 > hwPiOverTwo_) {
0134 edm::LogWarning("L1MetPfProducer") << "unexpected phi (high)";
0135 phiQ1 = hwPiOverTwo_;
0136 } else if (phiQ1 < 0) {
0137 edm::LogWarning("L1MetPfProducer") << "unexpected phi (low)";
0138 phiQ1 = 0;
0139 }
0140 if (isX) {
0141 typedef ap_ufixed<14, 12, AP_RND, AP_WRAP> pt_t;
0142 ap_ufixed<pt_t::width, 0> cosPhi = cos(phiQ1.to_double() / hwPiOverTwo_.to_double() * M_PI / 2);
0143 pxy = pt * cosPhi;
0144 if (phi > hwPiOverTwo_ || phi < -hwPiOverTwo_)
0145 pxy = -pxy;
0146 } else {
0147 ap_ufixed<pt_t::width, 0> sinPhi = sin(phiQ1.to_double() / hwPiOverTwo_.to_double() * M_PI / 2);
0148 pxy = pt * sinPhi;
0149 if (phi < 0)
0150 pxy = -pxy;
0151 }
0152 }
0153
0154 void L1MetPfProducer::PhiFromXY(pxy_t px, pxy_t py, phi_t& phi, bool debug) const {
0155 if (px == 0 && py == 0) {
0156 phi = 0;
0157 return;
0158 }
0159 if (px == 0) {
0160 phi = py > 0 ? hwPiOverTwo_ : phi_t(-hwPiOverTwo_);
0161 return;
0162 }
0163 if (py == 0) {
0164 phi = px > 0 ? phi_t(0) : phi_t(-hwPi_);
0165 return;
0166 }
0167
0168
0169 pt_t x = px > 0 ? pt_t(px) : pt_t(-px);
0170 pt_t y = py > 0 ? pt_t(py) : pt_t(-py);
0171
0172 pt_t a = x < y ? x : y;
0173 pt_t b = x < y ? y : x;
0174
0175 if (b.to_double() > maxPt_ / dropFactor_)
0176 b = maxPt_ / dropFactor_;
0177
0178 int index = round((b.to_double() / (maxPt_ / dropFactor_)) * invTableSize_);
0179 float bcheck = (float(index) / invTableSize_) * (maxPt_ / dropFactor_);
0180 inv_t inv_b = 1. / ((float(index) / invTableSize_) * (maxPt_ / dropFactor_));
0181
0182 inv_t a_over_b = a * inv_b;
0183
0184 if (debug) {
0185 printf(" a, b = %f, %f; index, inv = %d, %f; ratio = %f \n",
0186 a.to_double(),
0187 b.to_double(),
0188 index,
0189 inv_b.to_double(),
0190 a_over_b.to_double());
0191 printf(" bcheck, 1/bc = %f, %f -- %d %f %d \n", bcheck, 1. / bcheck, invTableSize_, maxPt_, dropFactor_);
0192 }
0193
0194 int atanTableBits_ = 7;
0195 int atanTableSize_ = (1 << atanTableBits_);
0196 index = round(a_over_b.to_double() * atanTableSize_);
0197 phi = atan(float(index) / atanTableSize_) / phiLSB_;
0198
0199 if (debug) {
0200 printf(" atan index, phi = %d, %f (%f rad) real atan(a/b)= %f \n",
0201 index,
0202 phi.to_double(),
0203 phi.to_double() * (M_PI / hwPi_.to_double()),
0204 atan(a.to_double() / b.to_double()));
0205 }
0206
0207
0208 if (y > x)
0209 phi = hwPiOverTwo_ - phi;
0210
0211 if (px < 0 && py > 0)
0212 phi = hwPi_ - phi;
0213 if (px > 0 && py < 0)
0214 phi = -phi;
0215 if (px < 0 && py < 0)
0216 phi = -(hwPi_ - phi);
0217
0218 if (debug) {
0219 printf(" phi hw, float, real = %f, %f (%f rad from x,y = %f, %f) \n",
0220 phi.to_double(),
0221 phi.to_double() * (M_PI / hwPi_.to_double()),
0222 atan2(py.to_double(), px.to_double()),
0223 px.to_double(),
0224 py.to_double());
0225 }
0226 }
0227
0228 L1MetPfProducer::~L1MetPfProducer() {}
0229
0230 #include "FWCore/Framework/interface/MakerMacros.h"
0231 DEFINE_FWK_MODULE(L1MetPfProducer);