L1TTrackMatch/plugins/L1TkHTMissEmulatorProducer.cc

0001
0002 /**\class L1TrackerHTMissEmulatorProducer L1TrackerHTMissEmulatorProducer.cc
0003  L1Trigger/L1TTrackMatch/plugins/L1TrackerHTMissEmulatorProducer.cc
0004  Description: Takes L1TTkJets and performs a integer emulation of Track-based missing HT, outputting a collection of EtSum
0005 */
0006
0007 // Original Author:  Hardik Routray
0008 //         Created:  Mon, 11 Oct 2021
0009
0010 // system include files
0011 #include <memory>
0012 #include <numeric>
0013
0014 // user include files
0015 #include "FWCore/Framework/interface/Frameworkfwd.h"
0016 #include "FWCore/Framework/interface/stream/EDProducer.h"
0017 #include "FWCore/Framework/interface/Event.h"
0018 #include "FWCore/Framework/interface/MakerMacros.h"
0019 #include "FWCore/Framework/interface/EventSetup.h"
0020 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0021 #include "DataFormats/Math/interface/LorentzVector.h"
0022 #include "DataFormats/L1TCorrelator/interface/TkHTMiss.h"
0023 #include "DataFormats/L1TCorrelator/interface/TkHTMissFwd.h"
0024 #include "DataFormats/L1Trigger/interface/EtSum.h"
0025 #include "DataFormats/L1Trigger/interface/TkJetWord.h"
0026 #include "L1Trigger/L1TTrackMatch/interface/L1TkHTMissEmulatorProducer.h"
0027
0028 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0029
0030 using namespace l1t;
0031
0032 class L1TkHTMissEmulatorProducer : public edm::stream::EDProducer<> {
0033 public:
0034   explicit L1TkHTMissEmulatorProducer(const edm::ParameterSet&);
0035   ~L1TkHTMissEmulatorProducer() override = default;
0036
0037 private:
0038   virtual void beginJob();
0039   void produce(edm::Event&, const edm::EventSetup&) override;
0040   virtual void endJob();
0041
0042   // ----------member data ---------------------------
0043
0044   bool debug_ = false;
0045   bool displaced_;
0046
0047   int cosLUTbins;
0048   int aSteps = 8;
0049
0050   l1tmhtemu::pt_t jetMinPt_;
0051   l1tmhtemu::eta_t jetMaxEta_;
0052   l1tmhtemu::ntracks_t minNtracksHighPt_;
0053   l1tmhtemu::ntracks_t minNtracksLowPt_;
0054
0055   std::vector<l1tmhtemu::phi_t> cosLUT_;
0056   std::vector<l1tmhtemu::MHTphi_t> atanLUT_;
0057   std::vector<l1tmhtemu::Et_t> magNormalisationLUT_;
0058
0059   std::string L1MHTCollectionName_;
0060
0061   const edm::EDGetTokenT<TkJetWordCollection> jetToken_;
0062 };
0063
0064 L1TkHTMissEmulatorProducer::L1TkHTMissEmulatorProducer(const edm::ParameterSet& iConfig)
0065     : jetToken_(consumes<TkJetWordCollection>(iConfig.getParameter<edm::InputTag>("L1TkJetEmulationInputTag"))) {
0066   debug_ = iConfig.getParameter<bool>("debug");
0067   displaced_ = iConfig.getParameter<bool>("displaced");
0068
0069   jetMinPt_ = l1tmhtemu::digitizeSignedValue<l1tmhtemu::pt_t>(
0070       (float)iConfig.getParameter<double>("jet_minPt"), l1tmhtemu::kInternalPtWidth, l1tmhtemu::kStepPt);
0071   jetMaxEta_ = l1tmhtemu::digitizeSignedValue<l1tmhtemu::eta_t>(
0072       (float)iConfig.getParameter<double>("jet_maxEta"), l1tmhtemu::kInternalEtaWidth, l1tmhtemu::kStepEta);
0073   minNtracksHighPt_ = (l1tmhtemu::ntracks_t)iConfig.getParameter<int>("jet_minNtracksHighPt");
0074   minNtracksLowPt_ = (l1tmhtemu::ntracks_t)iConfig.getParameter<int>("jet_minNtracksLowPt");
0075
0076   cosLUTbins = floor(l1tmhtemu::kMaxCosLUTPhi / l1tmhtemu::kStepPhi);
0077   cosLUT_ = l1tmhtemu::generateCosLUT(cosLUTbins);
0078
0079   atanLUT_ = l1tmhtemu::generateaTanLUT(aSteps);
0080   magNormalisationLUT_ = l1tmhtemu::generatemagNormalisationLUT(aSteps);
0081
0082   // Name of output ED Product
0083   L1MHTCollectionName_ = (std::string)iConfig.getParameter<std::string>("L1MHTCollectionName");
0084
0085   produces<std::vector<EtSum>>(L1MHTCollectionName_);
0086
0087   if (debug_) {
0088     edm::LogVerbatim("L1TrackerHTMissEmulatorProducer")
0089         << "-------------------------------------------------------------------------\n"
0090         << "====BITWIDTHS====\n"
0091         << "pt: " << l1t::TkJetWord::TkJetBitWidths::kPtSize << " eta: " << l1t::TkJetWord::TkJetBitWidths::kGlbEtaSize
0092         << " phi:" << l1t::TkJetWord::TkJetBitWidths::kGlbPhiSize << "\n"
0093         << "====CUT AP_INTS====\n"
0094         << "minpt: " << jetMinPt_ << " maxeta: " << jetMaxEta_ << " minNtracksHighPt: " << minNtracksHighPt_
0095         << " minNtracksLowPt: " << minNtracksLowPt_ << "\n"
0096         << "====CUT AP_INTS TO FLOATS====\n"
0097         << "minpt: " << (float)jetMinPt_ * l1tmhtemu::kStepPt << " maxeta: " << (float)jetMaxEta_ * l1tmhtemu::kStepEta
0098         << " minNtracksHighPt: " << (int)minNtracksHighPt_ << " minNtracksLowPt: " << (int)minNtracksLowPt_ << "\n"
0099         << "-------------------------------------------------------------------------\n";
0100   }
0101 }
0102
0103 void L1TkHTMissEmulatorProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0104   using namespace edm;
0105   std::unique_ptr<std::vector<l1t::EtSum>> MHTCollection(new std::vector<l1t::EtSum>(0));
0106
0107   // L1 track-trigger jets
0108   edm::Handle<TkJetWordCollection> L1TkJetsHandle;
0109   iEvent.getByToken(jetToken_, L1TkJetsHandle);
0110   std::vector<TkJetWord>::const_iterator jetIter;
0111
0112   if (!L1TkJetsHandle.isValid() && !displaced_) {
0113     LogError("TkHTMissEmulatorProducer") << "\nWarning: TkJetCollection not found in the event. Exit\n";
0114     return;
0115   }
0116
0117   if (!L1TkJetsHandle.isValid() && displaced_) {
0118     LogError("TkHTMissEmulatorProducer") << "\nWarning: TkJetExtendedCollection not found in the event. Exit\n";
0119     return;
0120   }
0121
0122   // floats used for debugging
0123   float sumPx_ = 0;
0124   float sumPy_ = 0;
0125   float HT_ = 0;
0126
0127   l1tmhtemu::Et_t sumPx = 0;
0128   l1tmhtemu::Et_t sumPy = 0;
0129   l1tmhtemu::MHT_t HT = 0;
0130
0131   // loop over jets
0132   int jetn = 0;
0133
0134   for (jetIter = L1TkJetsHandle->begin(); jetIter != L1TkJetsHandle->end(); ++jetIter) {
0135     // floats used for debugging
0136     float tmp_jet_px_ = jetIter->pt() * cos(jetIter->glbphi());
0137     float tmp_jet_py_ = jetIter->pt() * sin(jetIter->glbphi());
0138     //float tmp_jet_et_ = jetIter->pt();  // FIXME Get Et from the emulated jets
0139     float tmp_jet_pt_ = jetIter->pt();
0140
0141     l1tmhtemu::pt_t tmp_jet_pt =
0142         l1tmhtemu::digitizeSignedValue<l1tmhtemu::pt_t>(jetIter->pt(), l1tmhtemu::kInternalPtWidth, l1tmhtemu::kStepPt);
0143     l1tmhtemu::eta_t tmp_jet_eta = l1tmhtemu::digitizeSignedValue<l1tmhtemu::eta_t>(
0144         jetIter->glbeta(), l1tmhtemu::kInternalEtaWidth, l1tmhtemu::kStepEta);
0145     l1tmhtemu::phi_t tmp_jet_phi = l1tmhtemu::digitizeSignedValue<l1tmhtemu::phi_t>(
0146         jetIter->glbphi(), l1tmhtemu::kInternalPhiWidth, l1tmhtemu::kStepPhi);
0147     l1tmhtemu::ntracks_t tmp_jet_nt = l1tmhtemu::ntracks_t(jetIter->nt());
0148
0149     l1tmhtemu::phi_t tmp_jet_cos_phi = l1tmhtemu::phi_t(-999);
0150     l1tmhtemu::phi_t tmp_jet_sin_phi = l1tmhtemu::phi_t(-999);
0151
0152     if (tmp_jet_phi >= 0) {
0153       tmp_jet_cos_phi = cosLUT_[tmp_jet_phi];
0154
0155       if (cosLUTbins / 2 + 1 - tmp_jet_phi >= 0)
0156         tmp_jet_sin_phi = cosLUT_[cosLUTbins / 2 + 1 - tmp_jet_phi];
0157       else
0158         tmp_jet_sin_phi = cosLUT_[-1 * (cosLUTbins / 2 + 1 - tmp_jet_phi)];
0159
0160     } else {
0161       tmp_jet_cos_phi = cosLUT_[-1 * tmp_jet_phi];
0162
0163       if (cosLUTbins / 2 + 1 - (-1 * tmp_jet_phi) >= 0)
0164         tmp_jet_sin_phi = -1 * cosLUT_[cosLUTbins / 2 + 1 - (-1 * tmp_jet_phi)];
0165       else
0166         tmp_jet_sin_phi = -1 * cosLUT_[-1 * (cosLUTbins / 2 + 1 - (-1 * tmp_jet_phi))];
0167     }
0168
0169     l1tmhtemu::Et_t tmp_jet_px = tmp_jet_pt * tmp_jet_cos_phi;
0170     l1tmhtemu::Et_t tmp_jet_py = tmp_jet_pt * tmp_jet_sin_phi;
0171
0172     jetn++;
0173
0174     if (debug_) {
0175       edm::LogVerbatim("L1TrackerHTMissEmulatorProducer")
0176           << "****JET EMULATION" << jetn << "****\n"
0177           << "FLOATS ORIGINAL\n"
0178           << "PT: " << jetIter->pt() << "| ETA: " << jetIter->glbeta() << "| PHI: " << jetIter->glbphi()
0179           << "| NTRACKS: " << jetIter->nt() << "| COS(PHI): " << cos(jetIter->glbphi())
0180           << "| SIN(PHI): " << sin(jetIter->glbphi()) << "| Px: " << jetIter->pt() * cos(jetIter->glbphi())
0181           << "| Py: " << jetIter->pt() * sin(jetIter->glbphi()) << "\n"
0182           << "AP_INTS RAW\n"
0183           << "PT: " << jetIter->ptWord() << "| ETA: " << jetIter->glbEtaWord() << "| PHI: " << jetIter->glbPhiWord()
0184           << "| NTRACKS: " << jetIter->ntWord() << "\n"
0185           << "AP_INTS NEW\n"
0186           << "PT: " << tmp_jet_pt << "| ETA: " << tmp_jet_eta << "| PHI: " << tmp_jet_phi << "| NTRACKS: " << tmp_jet_nt
0187           << "| COS(PHI): " << tmp_jet_cos_phi << "| SIN(PHI): " << tmp_jet_sin_phi << "| Px: " << tmp_jet_px
0188           << "| Py: " << tmp_jet_py << "\n"
0189           << "AP_INTS NEW TO FLOATS\n"
0190           << "PT: " << (float)tmp_jet_pt * l1tmhtemu::kStepPt << "| ETA: " << (float)tmp_jet_eta * l1tmhtemu::kStepEta
0191           << "| PHI: " << (float)tmp_jet_phi * l1tmhtemu::kStepPhi << "| NTRACKS: " << (int)tmp_jet_nt
0192           << "| COS(PHI): " << (float)tmp_jet_cos_phi * l1tmhtemu::kStepPhi
0193           << "| SIN(PHI): " << (float)tmp_jet_sin_phi * l1tmhtemu::kStepPhi
0194           << "| Px: " << (float)tmp_jet_px * l1tmhtemu::kStepPt * l1tmhtemu::kStepPhi
0195           << "| Py: " << (float)tmp_jet_py * l1tmhtemu::kStepPt * l1tmhtemu::kStepPhi << "\n"
0196           << "-------------------------------------------------------------------------\n";
0197     }
0198
0199     if (tmp_jet_pt < jetMinPt_)
0200       continue;
0201     if (tmp_jet_eta > jetMaxEta_ or tmp_jet_eta < -1 * jetMaxEta_)
0202       continue;
0203     if (tmp_jet_nt < minNtracksLowPt_ && tmp_jet_pt > 200)
0204       continue;
0205     if (tmp_jet_nt < minNtracksHighPt_ && tmp_jet_pt > 400)
0206       continue;
0207
0208     if (debug_) {
0209       sumPx_ += tmp_jet_px_;
0210       sumPy_ += tmp_jet_py_;
0211       HT_ += tmp_jet_pt_;
0212     }
0213
0214     sumPx += tmp_jet_pt * tmp_jet_cos_phi;
0215     sumPy += tmp_jet_pt * tmp_jet_sin_phi;
0216     HT += tmp_jet_pt;
0217
0218   }  // end jet loop
0219
0220   // define missing HT
0221
0222   // Perform cordic sqrt, take x,y and converts to polar coordinate r,phi where
0223   // r=sqrt(x**2+y**2) and phi = atan(y/x)
0224   l1tmhtemu::EtMiss EtMiss = l1tmhtemu::cordicSqrt(sumPx, sumPy, aSteps, atanLUT_, magNormalisationLUT_);
0225   math::XYZTLorentzVector missingEt(-sumPx, -sumPy, 0, EtMiss.Et);
0226
0227   l1tmhtemu::MHTphi_t phi = 0;
0228
0229   if ((sumPx < 0) && (sumPy < 0))
0230     phi = EtMiss.Phi - l1tmhtemu::kMHTPhiBins / 2;
0231   else if ((sumPx >= 0) && (sumPy >= 0))
0232     phi = (EtMiss.Phi) + l1tmhtemu::kMHTPhiBins / 2;
0233   else if ((sumPx >= 0) && (sumPy < 0))
0234     phi = EtMiss.Phi - l1tmhtemu::kMHTPhiBins / 2;
0235   else if ((sumPx < 0) && (sumPy >= 0))
0236     phi = EtMiss.Phi - 3 * l1tmhtemu::kMHTPhiBins / 2;
0237
0238   if (debug_) {
0239     edm::LogVerbatim("L1TrackerHTMissEmulatorProducer")
0240         << "-------------------------------------------------------------------------\n"
0241         << "====MHT FLOATS====\n"
0242         << "sumPx: " << sumPx_ << "| sumPy: " << sumPy_ << "| ET: " << sqrt(sumPx_ * sumPx_ + sumPy_ * sumPy_)
0243         << "| HT: " << HT_ << "| PHI: " << atan2(sumPy_, sumPx_) << "\n"
0244         << "====MHT AP_INTS====\n"
0245         << "sumPx: " << sumPx << "| sumPy: " << sumPy << "| ET: " << EtMiss.Et << "| HT: " << HT << "| PHI: " << phi
0246         << "\n"
0247         << "====MHT AP_INTS TO FLOATS====\n"
0248         << "sumPx: " << (float)sumPx * l1tmhtemu::kStepPt * l1tmhtemu::kStepPhi
0249         << "| sumPy: " << (float)sumPy * l1tmhtemu::kStepPt * l1tmhtemu::kStepPhi << "| ET: " << EtMiss.Et.to_double()
0250         << "| HT: " << (float)HT * l1tmhtemu::kStepPt << "| PHI: " << (float)phi * l1tmhtemu::kStepMHTPhi - M_PI << "\n"
0251         << "-------------------------------------------------------------------------\n";
0252   }
0253   //rescale HT to correct output range
0254   HT = HT / (int)(1 / l1tmhtemu::kStepPt);
0255
0256   EtSum L1HTSum(missingEt, EtSum::EtSumType::kMissingHt, (int)HT.range(), 0, (int)phi, (int)jetn);
0257
0258   MHTCollection->push_back(L1HTSum);
0259   iEvent.put(std::move(MHTCollection), L1MHTCollectionName_);
0260
0261 }  //end producer
0262
0263 void L1TkHTMissEmulatorProducer::beginJob() {}
0264
0265 void L1TkHTMissEmulatorProducer::endJob() {}
0266
0267 DEFINE_FWK_MODULE(L1TkHTMissEmulatorProducer);