L1CaloTrigger/plugins/L1EGammaEEProducer.cc

0001 #include "FWCore/Framework/interface/stream/EDProducer.h"
0002 #include "FWCore/Framework/interface/Event.h"
0003 #include "FWCore/Framework/interface/MakerMacros.h"
0004
0005 #include "DataFormats/L1THGCal/interface/HGCalMulticluster.h"
0006 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0007 #include "DataFormats/L1Trigger/interface/EGamma.h"
0008 #include "L1Trigger/L1CaloTrigger/interface/L1EGammaEECalibrator.h"
0009 #include "DataFormats/Math/interface/deltaPhi.h"
0010
0011 namespace l1tp2 {
0012   // we sort the clusters in pt
0013   bool compare_cluster_pt(const l1t::HGCalMulticluster *cl1, const l1t::HGCalMulticluster *cl2) {
0014     return cl1->pt() > cl2->pt();
0015   }
0016 };  // namespace l1tp2
0017
0018 int etaBin(const l1t::HGCalMulticluster *cl) {
0019   static float constexpr eta_min = 1.;
0020   static float constexpr eta_max = 4.;
0021   static unsigned constexpr n_eta_bins = 150;
0022   int eta_bin = floor((std::abs(cl->eta()) - eta_min) / ((eta_max - eta_min) / n_eta_bins));
0023   if (cl->eta() < 0)
0024     return -1 * eta_bin;  // bin 0 doesn't exist
0025   return eta_bin;
0026 }
0027
0028 int get_phi_bin(const l1t::HGCalMulticluster *cl) {
0029   static constexpr float phi_min = -M_PI;
0030   static constexpr float phi_max = M_PI;
0031   static constexpr unsigned n_phi_bins = 63;
0032   return floor(std::abs(reco::deltaPhi(cl->phi(), phi_min)) / ((phi_max - phi_min) / n_phi_bins));
0033 }
0034
0035 pair<int, int> get_eta_phi_bin(const l1t::HGCalMulticluster *cl) { return std::make_pair(etaBin(cl), get_phi_bin(cl)); }
0036
0037 class L1EGammaEEProducer : public edm::stream::EDProducer<> {
0038 public:
0039   explicit L1EGammaEEProducer(const edm::ParameterSet &);
0040
0041 private:
0042   void produce(edm::Event &, const edm::EventSetup &) override;
0043
0044   edm::EDGetToken multiclusters_token_;
0045   L1EGammaEECalibrator calibrator_;
0046 };
0047
0048 L1EGammaEEProducer::L1EGammaEEProducer(const edm::ParameterSet &iConfig)
0049     : multiclusters_token_(
0050           consumes<l1t::HGCalMulticlusterBxCollection>(iConfig.getParameter<edm::InputTag>("Multiclusters"))),
0051       calibrator_(iConfig.getParameter<edm::ParameterSet>("calibrationConfig")) {
0052   produces<BXVector<l1t::EGamma>>("L1EGammaCollectionBXVWithCuts");
0053 }
0054
0055 void L1EGammaEEProducer::produce(edm::Event &iEvent, const edm::EventSetup &iSetup) {
0056   float minEt_ = 0;
0057
0058   std::unique_ptr<BXVector<l1t::EGamma>> l1EgammaBxCollection(new l1t::EGammaBxCollection);
0059
0060   // retrieve clusters 3D
0061   edm::Handle<l1t::HGCalMulticlusterBxCollection> multiclusters_h;
0062   iEvent.getByToken(multiclusters_token_, multiclusters_h);
0063   const l1t::HGCalMulticlusterBxCollection &multiclusters = *multiclusters_h;
0064
0065   std::vector<const l1t::HGCalMulticluster *> selected_multiclusters;
0066   std::map<std::pair<int, int>, std::vector<const l1t::HGCalMulticluster *>> etaphi_bins;
0067
0068   // here we loop on the TPGs
0069   for (auto cl3d = multiclusters.begin(0); cl3d != multiclusters.end(0); cl3d++) {
0070     if (cl3d->hwQual()) {
0071       if (cl3d->et() > minEt_) {
0072         int hw_quality = 1;  // baseline EG ID passed
0073         if (std::abs(cl3d->eta()) >= 1.52) {
0074           hw_quality = 2;  // baseline EG ID passed + cleanup of transition region
0075         }
0076
0077         float calib_factor = calibrator_.calibrationFactor(cl3d->pt(), cl3d->eta());
0078         l1t::EGamma eg =
0079             l1t::EGamma(reco::Candidate::PolarLorentzVector(cl3d->pt() / calib_factor, cl3d->eta(), cl3d->phi(), 0.));
0080         eg.setHwQual(hw_quality);
0081         eg.setHwIso(1);
0082         eg.setIsoEt(-1);  // just temporarily as a dummy value
0083         l1EgammaBxCollection->push_back(0, eg);
0084         if (hw_quality == 2) {
0085           // we build the EM interpreted EG object
0086           l1t::EGamma eg_emint = l1t::EGamma(reco::Candidate::PolarLorentzVector(
0087               cl3d->iPt(l1t::HGCalMulticluster::EnergyInterpretation::EM), cl3d->eta(), cl3d->phi(), 0.));
0088           eg_emint.setHwQual(4);
0089           eg_emint.setHwIso(1);
0090           eg_emint.setIsoEt(-1);  // just temporarily as a dummy value
0091           l1EgammaBxCollection->push_back(0, eg_emint);
0092           // we also prepare for the brem recovery procedure
0093           selected_multiclusters.push_back(&(*cl3d));
0094           auto eta_phi_bin = get_eta_phi_bin(&(*cl3d));
0095           auto bucket = etaphi_bins.find(eta_phi_bin);
0096           if (bucket == etaphi_bins.end()) {
0097             std::vector<const l1t::HGCalMulticluster *> vec;
0098             vec.push_back(&(*cl3d));
0099             etaphi_bins[eta_phi_bin] = vec;
0100           } else {
0101             bucket->second.push_back(&(*cl3d));
0102           }
0103         }
0104       }
0105     }
0106   }
0107
0108   std::sort(selected_multiclusters.begin(), selected_multiclusters.end(), l1tp2::compare_cluster_pt);
0109   std::set<const l1t::HGCalMulticluster *> used_clusters;
0110   for (const auto &cl3d : selected_multiclusters) {
0111     if (used_clusters.find(cl3d) == used_clusters.end()) {
0112       float pt = cl3d->pt();
0113       // we drop the Had component of the energy
0114       if (cl3d->hOverE() != -1)
0115         pt = cl3d->pt() / (1 + cl3d->hOverE());
0116       reco::Candidate::PolarLorentzVector mom(pt, cl3d->eta(), cl3d->phi(), 0.);
0117       reco::Candidate::PolarLorentzVector mom_eint(
0118           cl3d->iPt(l1t::HGCalMulticluster::EnergyInterpretation::EM), cl3d->eta(), cl3d->phi(), 0.);
0119
0120       // this is not yet used
0121       used_clusters.insert(cl3d);
0122       auto eta_phi_bin = get_eta_phi_bin(cl3d);
0123
0124       for (int eta_bin : {eta_phi_bin.first - 1, eta_phi_bin.first, eta_phi_bin.first + 1}) {
0125         for (int phi_bin : {eta_phi_bin.second - 1, eta_phi_bin.second, eta_phi_bin.second + 1}) {
0126           auto bucket = etaphi_bins.find(std::make_pair(eta_bin, phi_bin));
0127           if (bucket != etaphi_bins.end()) {
0128             // this bucket is not empty
0129             for (const auto &other_cl_ptr : bucket->second) {
0130               if (used_clusters.find(other_cl_ptr) == used_clusters.end()) {
0131                 if (std::abs(other_cl_ptr->eta() - cl3d->eta()) < 0.02) {
0132                   if (std::abs(reco::deltaPhi(other_cl_ptr->phi(), cl3d->phi())) < 0.1) {
0133                     float pt_other = other_cl_ptr->pt();
0134                     if (other_cl_ptr->hOverE() != -1)
0135                       pt_other = other_cl_ptr->pt() / (1 + other_cl_ptr->hOverE());
0136                     mom += reco::Candidate::PolarLorentzVector(pt_other, other_cl_ptr->eta(), other_cl_ptr->phi(), 0.);
0137                     mom_eint += reco::Candidate::PolarLorentzVector(
0138                         other_cl_ptr->iPt(l1t::HGCalMulticluster::EnergyInterpretation::EM),
0139                         other_cl_ptr->eta(),
0140                         other_cl_ptr->phi(),
0141                         0.);
0142                     used_clusters.insert(other_cl_ptr);
0143                   }
0144                 }
0145               }
0146             }
0147           }
0148         }
0149       }
0150       float calib_factor = calibrator_.calibrationFactor(mom.pt(), mom.eta());
0151       l1t::EGamma eg =
0152           l1t::EGamma(reco::Candidate::PolarLorentzVector(mom.pt() / calib_factor, mom.eta(), mom.phi(), 0.));
0153       eg.setHwQual(3);
0154       eg.setHwIso(1);
0155       l1EgammaBxCollection->push_back(0, eg);
0156
0157       l1t::EGamma eg_emint_brec =
0158           l1t::EGamma(reco::Candidate::PolarLorentzVector(mom_eint.pt(), mom_eint.eta(), mom_eint.phi(), 0.));
0159       eg_emint_brec.setHwQual(5);
0160       eg_emint_brec.setHwIso(1);
0161       l1EgammaBxCollection->push_back(0, eg_emint_brec);
0162     }
0163   }
0164
0165   iEvent.put(std::move(l1EgammaBxCollection), "L1EGammaCollectionBXVWithCuts");
0166 }
0167
0168 DEFINE_FWK_MODULE(L1EGammaEEProducer);