L1CaloTrigger/plugins/Phase1L1TJetSumsProducer.cc

0001 // -*- C++ -*-
0002 //
0003 // Package: L1CaloTrigger
0004 // Class: Phase1L1TJetSumsProducer
0005 //
0006 /**\class Phase1L1TJetSumsProducer Phase1L1TJetSumsProducer.cc L1Trigger/L1CaloTrigger/plugin/Phase1L1TJetSumsProducer.cc
0007
0008 Description: Computes HT and MHT from phase-1-like jets
0009
0010 *** INPUT PARAMETERS ***
0011   * sin/cosPhi: Value of sin/cos phi in the middle of each bin of the grid.
0012   * etaBinning: vdouble with eta binning (allows non-homogeneous binning in eta)
0013   * nBinsPhi: uint32, number of bins in phi
0014   * phiLow: double, min phi (typically -pi)
0015   * phiUp: double, max phi (typically +pi)
0016   * {m}htPtThreshold: Minimum jet pt for HT/MHT calculation
0017   * {m}htAbsEtaCut:
0018   * pt/eta/philsb : lsb of quantities used in firmware implementation
0019   * outputCollectionName: string, tag for the output collection
0020    * inputCollectionTag: tag for input jet collection
0021
0022 */
0023 //
0024 // Original Simone Bologna
0025 //
0026
0027 #include "FWCore/Framework/interface/Frameworkfwd.h"
0028 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0029 #include "FWCore/Framework/interface/one/EDProducer.h"
0030 #include "DataFormats/JetReco/interface/CaloJet.h"
0031 #include "DataFormats/L1Trigger/interface/EtSum.h"
0032 #include "DataFormats/Common/interface/View.h"
0033 #include "DataFormats/Candidate/interface/Candidate.h"
0034 #include "FWCore/Framework/interface/MakerMacros.h"
0035 #include "FWCore/Framework/interface/Event.h"
0036 #include "DataFormats/Math/interface/LorentzVector.h"
0037 #include "FWCore/ServiceRegistry/interface/Service.h"
0038 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0039
0040 #include <cmath>
0041
0042 class Phase1L1TJetSumsProducer : public edm::one::EDProducer<edm::one::SharedResources> {
0043 public:
0044   explicit Phase1L1TJetSumsProducer(const edm::ParameterSet&);
0045   ~Phase1L1TJetSumsProducer() override = default;
0046
0047   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0048
0049 private:
0050   void produce(edm::Event&, const edm::EventSetup&) override;
0051
0052   // computes ht, adds jet pt to ht only if the pt of the jet is above the ht calculation threshold
0053   l1t::EtSum computeHT(const edm::Handle<std::vector<reco::CaloJet> > inputJets) const;
0054
0055   // computes MHT
0056   // adds jet pt to mht only if the pt of the jet is above the mht calculation threshold
0057   // performs some calculations with digitised/integer quantities to ensure agreement with firmware
0058   l1t::EtSum computeMHT(const edm::Handle<std::vector<reco::CaloJet> > inputJets) const;
0059
0060   const edm::EDGetTokenT<std::vector<reco::CaloJet> > inputJetCollectionTag_;
0061
0062   // holds the sin and cos for HLs LUT emulation
0063   const std::vector<double> sinPhi_;
0064   const std::vector<double> cosPhi_;
0065   const unsigned int nBinsPhi_;
0066
0067   // lower phi value
0068   const double phiLow_;
0069   // higher phi value
0070   const double phiUp_;
0071   // size of a phi bin
0072   const double phiStep_;
0073   // threshold for ht calculation
0074   const double htPtThreshold_;
0075   // threshold for ht calculation
0076   const double mhtPtThreshold_;
0077   // jet eta cut for ht calculation
0078   const double htAbsEtaCut_;
0079   // jet eta cut for mht calculation
0080   const double mhtAbsEtaCut_;
0081   // LSB of pt quantity
0082   const double ptlsb_;
0083   // label of sums
0084   const std::string outputCollectionName_;
0085 };
0086
0087 // initialises plugin configuration and prepares ROOT file for saving the sums
0088 Phase1L1TJetSumsProducer::Phase1L1TJetSumsProducer(const edm::ParameterSet& iConfig)
0089     : inputJetCollectionTag_{consumes<std::vector<reco::CaloJet> >(
0090           iConfig.getParameter<edm::InputTag>("inputJetCollectionTag"))},
0091       sinPhi_(iConfig.getParameter<std::vector<double> >("sinPhi")),
0092       cosPhi_(iConfig.getParameter<std::vector<double> >("cosPhi")),
0093       nBinsPhi_(iConfig.getParameter<unsigned int>("nBinsPhi")),
0094       phiLow_(iConfig.getParameter<double>("phiLow")),
0095       phiUp_(iConfig.getParameter<double>("phiUp")),
0096       phiStep_((phiUp_ - phiLow_) / nBinsPhi_),
0097       htPtThreshold_(iConfig.getParameter<double>("htPtThreshold")),
0098       mhtPtThreshold_(iConfig.getParameter<double>("mhtPtThreshold")),
0099       htAbsEtaCut_(iConfig.getParameter<double>("htAbsEtaCut")),
0100       mhtAbsEtaCut_(iConfig.getParameter<double>("mhtAbsEtaCut")),
0101       ptlsb_(iConfig.getParameter<double>("ptlsb")),
0102       outputCollectionName_(iConfig.getParameter<std::string>("outputCollectionName")) {
0103   usesResource("TFileService");
0104   produces<std::vector<l1t::EtSum> >(outputCollectionName_).setBranchAlias(outputCollectionName_);
0105 }
0106
0107 void Phase1L1TJetSumsProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0108   const edm::Handle<std::vector<reco::CaloJet> >& jetCollectionHandle = iEvent.getHandle(inputJetCollectionTag_);
0109
0110   // computing sums and storing them in sum object
0111   l1t::EtSum lHT = computeHT(jetCollectionHandle);
0112   l1t::EtSum lMHT = computeMHT(jetCollectionHandle);
0113
0114   //packing sums in vector for event saving
0115   std::unique_ptr<std::vector<l1t::EtSum> > lSumVectorPtr(new std::vector<l1t::EtSum>(0));
0116   lSumVectorPtr->push_back(lHT);
0117   lSumVectorPtr->push_back(lMHT);
0118   iEvent.put(std::move(lSumVectorPtr), outputCollectionName_);
0119
0120   return;
0121 }
0122
0123 l1t::EtSum Phase1L1TJetSumsProducer::computeHT(const edm::Handle<std::vector<reco::CaloJet> > inputJets) const {
0124   double lHT = 0;
0125   for (const auto& jet : *inputJets) {
0126     double lJetPt = jet.pt();
0127     double lJetPhi = jet.phi();
0128     double lJetEta = jet.eta();
0129     if ((lJetPhi < phiLow_) || (lJetPhi >= phiUp_))
0130       continue;
0131
0132     lHT += (lJetPt >= htPtThreshold_ && std::fabs(lJetEta) < htAbsEtaCut_) ? lJetPt : 0;
0133   }
0134
0135   reco::Candidate::PolarLorentzVector lHTVector;
0136   lHTVector.SetPt(lHT);
0137   lHTVector.SetEta(0);
0138   lHTVector.SetPhi(0);
0139   l1t::EtSum lHTSum(lHTVector, l1t::EtSum::EtSumType::kTotalHt, 0, 0, 0, 0);
0140   return lHTSum;
0141 }
0142
0143 l1t::EtSum Phase1L1TJetSumsProducer::computeMHT(const edm::Handle<std::vector<reco::CaloJet> > inputJets) const {
0144   int lTotalJetPx = 0;
0145   int lTotalJetPy = 0;
0146
0147   std::vector<unsigned int> jetPtInPhiBins(nBinsPhi_, 0);
0148
0149   for (const auto& jet : *inputJets) {
0150     double lJetPhi = jet.phi();
0151
0152     if ((lJetPhi < phiLow_) || (lJetPhi >= phiUp_))
0153       continue;
0154
0155     unsigned int iPhi = (lJetPhi - phiLow_) / phiStep_;
0156
0157     if (jet.pt() >= mhtPtThreshold_ && std::fabs(jet.eta()) < mhtAbsEtaCut_) {
0158       unsigned int digiJetPt = floor(jet.pt() / ptlsb_);
0159       jetPtInPhiBins[iPhi] += digiJetPt;
0160     }
0161   }
0162
0163   for (unsigned int iPhi = 0; iPhi < jetPtInPhiBins.size(); ++iPhi) {
0164     unsigned int digiJetPtSum = jetPtInPhiBins[iPhi];
0165
0166     // retrieving sin cos from LUT emulator
0167     double lSinPhi = sinPhi_[iPhi];
0168     double lCosPhi = cosPhi_[iPhi];
0169
0170     // checking if above threshold
0171     lTotalJetPx += trunc(digiJetPtSum * lCosPhi);
0172     lTotalJetPy += trunc(digiJetPtSum * lSinPhi);
0173   }
0174
0175   double lMHT = floor(sqrt(lTotalJetPx * lTotalJetPx + lTotalJetPy * lTotalJetPy)) * ptlsb_;
0176   math::PtEtaPhiMLorentzVector lMHTVector(lMHT, 0, acos(lTotalJetPx / (lMHT / ptlsb_)), 0);
0177   l1t::EtSum lMHTSum(lMHTVector, l1t::EtSum::EtSumType::kMissingHt, 0, 0, 0, 0);
0178
0179   return lMHTSum;
0180 }
0181
0182 void Phase1L1TJetSumsProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0183   edm::ParameterSetDescription desc;
0184   desc.add<edm::InputTag>("inputJetCollectionTag",
0185                           edm::InputTag("l1tPhase1JetCalibrator", "Phase1L1TJetFromPfCandidates"));
0186   desc.add<std::vector<double> >("sinPhi");
0187   desc.add<std::vector<double> >("cosPhi");
0188   desc.add<unsigned int>("nBinsPhi", 72);
0189   desc.add<double>("phiLow", -M_PI);
0190   desc.add<double>("phiUp", M_PI);
0191   desc.add<double>("htPtThreshold", 30);
0192   desc.add<double>("mhtPtThreshold", 30);
0193   desc.add<double>("htAbsEtaCut", 3);
0194   desc.add<double>("mhtAbsEtaCut", 3);
0195   desc.add<double>("ptlsb", 0.25), desc.add<string>("outputCollectionName", "Sums");
0196   descriptions.add("Phase1L1TJetSumsProducer", desc);
0197 }
0198
0199 // creates the plugin for later use in python
0200 DEFINE_FWK_MODULE(Phase1L1TJetSumsProducer);