Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-05-14 02:38:03

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 
0011 #include "DataFormats/L1TParticleFlow/interface/PFCandidate.h"
0012 #include "DataFormats/L1Trigger/interface/EtSum.h"
0013 #include "DataFormats/Math/interface/LorentzVector.h"
0014 
0015 using namespace l1t;
0016 
0017 class L1METPFProducer : public edm::global::EDProducer<> {
0018 public:
0019   explicit L1METPFProducer(const edm::ParameterSet&);
0020   ~L1METPFProducer() override;
0021 
0022   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0023 
0024 private:
0025   void produce(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const override;
0026   edm::EDGetTokenT<vector<l1t::PFCandidate>> _l1PFToken;
0027 
0028   int maxCands_ = 128;
0029 
0030   // quantization controllers
0031   typedef ap_ufixed<14, 12, AP_RND, AP_WRAP> pt_t;  // LSB is 0.25 and max is 4 TeV
0032   typedef ap_int<12> phi_t;                         // LSB is pi/720 ~ 0.0044 and max is +/-8.9
0033   const float ptLSB_ = 0.25;                        // GeV
0034   const float phiLSB_ = M_PI / 720;                 // rad
0035 
0036   // derived, helper types
0037   typedef ap_fixed<pt_t::width + 1, pt_t::iwidth + 1, AP_RND, AP_SAT> pxy_t;
0038   typedef ap_fixed<2 * pt_t::width, 2 * pt_t::iwidth, AP_RND, AP_SAT> pt2_t;
0039   // derived, helper constants
0040   const float maxPt_ = ((1 << pt_t::width) - 1) * ptLSB_;
0041   const phi_t hwPi_ = round(M_PI / phiLSB_);
0042   const phi_t hwPiOverTwo_ = round(M_PI / (2 * phiLSB_));
0043 
0044   typedef ap_ufixed<pt_t::width, 0> inv_t;  // can't easily use the MAXPT/pt trick with ap_fixed
0045 
0046   // to make configurable...
0047   const int dropBits_ = 2;
0048   const int dropFactor_ = (1 << dropBits_);
0049   const int invTableBits_ = 10;
0050   const int invTableSize_ = (1 << invTableBits_);
0051 
0052   void Project(pt_t pt, phi_t phi, pxy_t& pxy, bool isX, bool debug = false) const;
0053   void PhiFromXY(pxy_t px, pxy_t py, phi_t& phi, bool debug = false) const;
0054 
0055   void CalcMetHLS(const std::vector<float>& pt,
0056                   const std::vector<float>& phi,
0057                   reco::Candidate::PolarLorentzVector& metVector) const;
0058 };
0059 
0060 L1METPFProducer::L1METPFProducer(const edm::ParameterSet& cfg)
0061     : _l1PFToken(consumes<std::vector<l1t::PFCandidate>>(cfg.getParameter<edm::InputTag>("L1PFObjects"))),
0062       maxCands_(cfg.getParameter<int>("maxCands")) {
0063   produces<std::vector<l1t::EtSum>>();
0064 }
0065 
0066 void L1METPFProducer::produce(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const {
0067   edm::Handle<l1t::PFCandidateCollection> l1PFCandidates;
0068   iEvent.getByToken(_l1PFToken, l1PFCandidates);
0069 
0070   std::vector<float> pt;
0071   std::vector<float> phi;
0072 
0073   for (int i = 0; i < int(l1PFCandidates->size()) && (i < maxCands_ || maxCands_ < 0); i++) {
0074     const auto& l1PFCand = l1PFCandidates->at(i);
0075     pt.push_back(l1PFCand.pt());
0076     phi.push_back(l1PFCand.phi());
0077   }
0078 
0079   reco::Candidate::PolarLorentzVector metVector;
0080 
0081   CalcMetHLS(pt, phi, metVector);
0082 
0083   l1t::EtSum theMET(metVector, l1t::EtSum::EtSumType::kTotalHt, 0, 0, 0, 0);
0084 
0085   std::unique_ptr<std::vector<l1t::EtSum>> metCollection(new std::vector<l1t::EtSum>(0));
0086   metCollection->push_back(theMET);
0087   iEvent.put(std::move(metCollection));
0088 }
0089 
0090 void L1METPFProducer::CalcMetHLS(const std::vector<float>& pt,
0091                                  const std::vector<float>& phi,
0092                                  reco::Candidate::PolarLorentzVector& metVector) const {
0093   pxy_t hw_px = 0;
0094   pxy_t hw_py = 0;
0095   pxy_t hw_sumx = 0;
0096   pxy_t hw_sumy = 0;
0097 
0098   for (uint i = 0; i < pt.size(); i++) {
0099     pt_t hw_pt = min(pt[i], maxPt_);
0100     phi_t hw_phi = float(TVector2::Phi_mpi_pi(phi[i]) / phiLSB_);
0101 
0102     Project(hw_pt, hw_phi, hw_px, true);
0103     Project(hw_pt, hw_phi, hw_py, false);
0104 
0105     hw_sumx = hw_sumx - hw_px;
0106     hw_sumy = hw_sumy - hw_py;
0107   }
0108 
0109   pt2_t hw_met = pt2_t(hw_sumx) * pt2_t(hw_sumx) + pt2_t(hw_sumy) * pt2_t(hw_sumy);
0110   hw_met = sqrt(int(hw_met));  // stand-in for HLS::sqrt
0111 
0112   phi_t hw_met_phi = 0;
0113   PhiFromXY(hw_sumx, hw_sumy, hw_met_phi);
0114 
0115   metVector.SetPt(hw_met.to_double());
0116   metVector.SetPhi(hw_met_phi.to_double() * phiLSB_);
0117   metVector.SetEta(0);
0118 }
0119 
0120 void L1METPFProducer::Project(pt_t pt, phi_t phi, pxy_t& pxy, bool isX, bool debug) const {
0121   /*
0122       Convert pt and phi to px (py)
0123       1) Map phi to the first quadrant to reduce LUT size
0124       2) Lookup sin(phiQ1), where the result is in [0,maxPt]
0125       which is used to encode [0,1].
0126       3) Multiply pt by sin(phiQ1) to get px. Result will be px*maxPt, but
0127       wrapping multiplication is 'mod maxPt' so the correct value is returned.
0128       4) Check px=-|px|.
0129     */
0130 
0131   // set phi to first quadrant
0132   phi_t phiQ1 = (phi > 0) ? phi : phi_t(-phi);  // Q1/Q4
0133   if (phiQ1 >= hwPiOverTwo_)
0134     phiQ1 = hwPi_ - phiQ1;
0135 
0136   if (phiQ1 > hwPiOverTwo_) {
0137     edm::LogWarning("L1METPFProducer") << "unexpected phi (high)";
0138     phiQ1 = hwPiOverTwo_;
0139   } else if (phiQ1 < 0) {
0140     edm::LogWarning("L1METPFProducer") << "unexpected phi (low)";
0141     phiQ1 = 0;
0142   }
0143   if (isX) {
0144     typedef ap_ufixed<14, 12, AP_RND, AP_WRAP> pt_t;  // LSB is 0.25 and max is 4 TeV
0145     ap_ufixed<pt_t::width, 0> cosPhi = cos(phiQ1.to_double() / hwPiOverTwo_.to_double() * M_PI / 2);
0146     pxy = pt * cosPhi;
0147     if (phi > hwPiOverTwo_ || phi < -hwPiOverTwo_)
0148       pxy = -pxy;
0149   } else {
0150     ap_ufixed<pt_t::width, 0> sinPhi = sin(phiQ1.to_double() / hwPiOverTwo_.to_double() * M_PI / 2);
0151     pxy = pt * sinPhi;
0152     if (phi < 0)
0153       pxy = -pxy;
0154   }
0155 }
0156 
0157 void L1METPFProducer::PhiFromXY(pxy_t px, pxy_t py, phi_t& phi, bool debug) const {
0158   if (px == 0 && py == 0) {
0159     phi = 0;
0160     return;
0161   }
0162   if (px == 0) {
0163     phi = py > 0 ? hwPiOverTwo_ : phi_t(-hwPiOverTwo_);
0164     return;
0165   }
0166   if (py == 0) {
0167     phi = px > 0 ? phi_t(0) : phi_t(-hwPi_);
0168     return;
0169   }
0170 
0171   // get q1 coordinates
0172   pt_t x = px > 0 ? pt_t(px) : pt_t(-px);  //px>=0 ? px : -px;
0173   pt_t y = py > 0 ? pt_t(py) : pt_t(-py);  //px>=0 ? px : -px;
0174   // transform so a<b
0175   pt_t a = x < y ? x : y;
0176   pt_t b = x < y ? y : x;
0177 
0178   if (b.to_double() > maxPt_ / dropFactor_)
0179     b = maxPt_ / dropFactor_;
0180   // map [0,max/4) to inv table size
0181   int index = round((b.to_double() / (maxPt_ / dropFactor_)) * invTableSize_);
0182   float bcheck = (float(index) / invTableSize_) * (maxPt_ / dropFactor_);
0183   inv_t inv_b = 1. / ((float(index) / invTableSize_) * (maxPt_ / dropFactor_));
0184 
0185   inv_t a_over_b = a * inv_b;
0186 
0187   if (debug) {
0188     LogDebug("L1METPFProducer") << "  a, b = \n  " << a.to_double() << " , " << b.to_double()
0189                                 << ";   index, inv = " << index << ", " << inv_b.to_double()
0190                                 << "; ratio= " << a_over_b.to_double() << " \n"
0191                                 << std::endl;
0192     LogDebug("L1METPFProducer") << "bcheck, 1/bc = " << bcheck << ", " << 1. / bcheck << " -- " << invTableSize_ << " "
0193                                 << maxPt_ << " " << dropFactor_ << " \n"
0194                                 << std::endl;
0195   }
0196 
0197   const int atanTableBits_ = 7;
0198   const int atanTableSize_ = (1 << atanTableBits_);
0199   index = round(a_over_b.to_double() * atanTableSize_);
0200   phi = atan(float(index) / atanTableSize_) / phiLSB_;
0201 
0202   if (debug) {
0203     LogDebug("L1METPFProducer") << "    atan index, phi = " << index << ", " << phi.to_double() << " ("
0204                                 << phi.to_double() * (M_PI / hwPi_.to_double())
0205                                 << " rad) real atan(a/b)= " << atan(a.to_double() / b.to_double()) << " \n"
0206                                 << std::endl;
0207   }
0208 
0209   // rotate from (0,pi/4) to full quad1
0210   if (y > x)
0211     phi = hwPiOverTwo_ - phi;  //phi = pi/2 - phi
0212   // other quadrants
0213   if (px < 0 && py > 0)
0214     phi = hwPi_ - phi;  // Q2 phi = pi - phi
0215   if (px > 0 && py < 0)
0216     phi = -phi;  // Q4 phi = -phi
0217   if (px < 0 && py < 0)
0218     phi = -(hwPi_ - phi);  // Q3 composition of both
0219 
0220   if (debug) {
0221     LogDebug("L1METPFProducer") << "    phi hw, float, real = " << phi.to_double() << ", "
0222                                 << phi.to_double() * (M_PI / hwPi_.to_double()) << "     ("
0223                                 << atan2(py.to_double(), px.to_double()) << " rad from x,y = " << px.to_double() << ", "
0224                                 << py.to_double() << ") \n"
0225                                 << std::endl;
0226   }
0227 }
0228 
0229 void L1METPFProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0230   edm::ParameterSetDescription desc;
0231   desc.add<int>("maxCandidates", 128);
0232   desc.add<edm::InputTag>("L1PFObjects", edm::InputTag("L1PFProducer", "l1pfCandidates"));
0233   descriptions.add("L1METPFProducer", desc);
0234 }
0235 
0236 L1METPFProducer::~L1METPFProducer() {}
0237 
0238 #include "FWCore/Framework/interface/MakerMacros.h"
0239 DEFINE_FWK_MODULE(L1METPFProducer);