Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-06-10 01:53:54

0001 // system include files
0002 #include <memory>
0003 #include <algorithm>
0004 
0005 // user include files
0006 #include "FWCore/Framework/interface/Frameworkfwd.h"
0007 #include "FWCore/Framework/interface/stream/EDProducer.h"
0008 #include "FWCore/Framework/interface/Event.h"
0009 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0010 #include "FWCore/ParameterSet/interface/FileInPath.h"
0011 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0012 
0013 #include "DataFormats/HcalDetId/interface/HcalTrigTowerDetId.h"
0014 #include "DataFormats/HcalDigi/interface/HcalDigiCollections.h"
0015 #include "CalibFormats/CaloTPG/interface/CaloTPGTranscoder.h"
0016 #include "CalibFormats/CaloTPG/interface/CaloTPGRecord.h"
0017 #include "L1Trigger/L1TCalorimeter/interface/CaloTools.h"
0018 
0019 #include "DataFormats/L1THGCal/interface/HGCalTower.h"
0020 
0021 #include "DataFormats/L1TCalorimeterPhase2/interface/CaloTower.h"
0022 
0023 #include "DataFormats/Math/interface/deltaPhi.h"
0024 
0025 #include "L1Trigger/Phase2L1ParticleFlow/src/corrector.h"
0026 #include "L1Trigger/Phase2L1ParticleFlow/interface/ParametricResolution.h"
0027 #include "L1Trigger/Phase2L1ParticleFlow/interface/CaloClusterer.h"
0028 
0029 //--------------------------------------------------------------------------------------------------
0030 class L1TPFCaloProducer : public edm::stream::EDProducer<> {
0031 public:
0032   explicit L1TPFCaloProducer(const edm::ParameterSet &);
0033 
0034 private:
0035   bool ecalOnly_, debug_;
0036   std::vector<edm::EDGetTokenT<reco::CandidateView>> ecalCands_;
0037   std::vector<edm::EDGetTokenT<reco::CandidateView>> hcalCands_;
0038 
0039   std::vector<edm::EDGetTokenT<HcalTrigPrimDigiCollection>> hcalDigis_;
0040   edm::ESGetToken<CaloTPGTranscoder, CaloTPGRecord> decoderTag_;
0041   bool hcalDigisBarrel_, hcalDigisHF_;
0042   std::vector<edm::EDGetTokenT<l1tp2::CaloTowerCollection>> phase2barrelTowers_;
0043   std::vector<edm::EDGetTokenT<l1t::HGCalTowerBxCollection>> hcalHGCTowers_;
0044   bool hcalHGCTowersHadOnly_;
0045 
0046   l1tpf::corrector emCorrector_;
0047   l1tpf::corrector hcCorrector_;
0048   l1tpf::corrector hadCorrector_;
0049 
0050   l1tpf_calo::SingleCaloClusterer ecalClusterer_, hcalClusterer_;
0051   std::unique_ptr<l1tpf_calo::SimpleCaloLinkerBase> caloLinker_;
0052 
0053   l1tpf::ParametricResolution resol_;
0054 
0055   void produce(edm::Event &, const edm::EventSetup &) override;
0056 
0057   void readHcalDigis_(edm::Event &event, const edm::EventSetup &);
0058   void readPhase2BarrelCaloTowers_(edm::Event &event, const edm::EventSetup &);
0059   void readHcalHGCTowers_(edm::Event &event, const edm::EventSetup &);
0060   struct SimpleHGCTC {
0061     float et, eta, phi;
0062     SimpleHGCTC(float aet, float aeta, float aphi) : et(aet), eta(aeta), phi(aphi) {}
0063   };
0064 };
0065 
0066 //
0067 // constructors and destructor
0068 //
0069 L1TPFCaloProducer::L1TPFCaloProducer(const edm::ParameterSet &iConfig)
0070     : ecalOnly_(iConfig.existsAs<bool>("ecalOnly") ? iConfig.getParameter<bool>("ecalOnly") : false),
0071       debug_(iConfig.getUntrackedParameter<int>("debug", 0)),
0072       decoderTag_(esConsumes<CaloTPGTranscoder, CaloTPGRecord>(edm::ESInputTag("", ""))),
0073       emCorrector_(iConfig.getParameter<std::string>("emCorrector"), -1, debug_),
0074       hcCorrector_(iConfig.getParameter<std::string>("hcCorrector"), -1, debug_),
0075       hadCorrector_(iConfig.getParameter<std::string>("hadCorrector"),
0076                     iConfig.getParameter<double>("hadCorrectorEmfMax"),
0077                     debug_),
0078       ecalClusterer_(iConfig.getParameter<edm::ParameterSet>("ecalClusterer")),
0079       hcalClusterer_(iConfig.getParameter<edm::ParameterSet>("hcalClusterer")),
0080       caloLinker_(l1tpf_calo::makeCaloLinker(
0081           iConfig.getParameter<edm::ParameterSet>("linker"), ecalClusterer_, hcalClusterer_)),
0082       resol_(iConfig.getParameter<edm::ParameterSet>("resol")) {
0083   produces<l1t::PFClusterCollection>("ecalCells");
0084 
0085   produces<l1t::PFClusterCollection>("emCalibrated");
0086   produces<l1t::PFClusterCollection>("emUncalibrated");
0087 
0088   for (auto &tag : iConfig.getParameter<std::vector<edm::InputTag>>("ecalCandidates")) {
0089     ecalCands_.push_back(consumes<reco::CandidateView>(tag));
0090   }
0091 
0092   if (ecalOnly_)
0093     return;
0094 
0095   produces<l1t::PFClusterCollection>("hcalCells");
0096 
0097   produces<l1t::PFClusterCollection>("hcalUnclustered");
0098   produces<l1t::PFClusterCollection>("hcalUncalibrated");
0099   produces<l1t::PFClusterCollection>("hcalCalibrated");
0100 
0101   produces<l1t::PFClusterCollection>("uncalibrated");
0102   produces<l1t::PFClusterCollection>("calibrated");
0103 
0104   for (auto &tag : iConfig.getParameter<std::vector<edm::InputTag>>("hcalCandidates")) {
0105     hcalCands_.push_back(consumes<reco::CandidateView>(tag));
0106   }
0107 
0108   for (auto &tag : iConfig.getParameter<std::vector<edm::InputTag>>("hcalDigis")) {
0109     hcalDigis_.push_back(consumes<HcalTrigPrimDigiCollection>(tag));
0110   }
0111   if (!hcalDigis_.empty()) {
0112     hcalDigisBarrel_ = iConfig.getParameter<bool>("hcalDigisBarrel");
0113     hcalDigisHF_ = iConfig.getParameter<bool>("hcalDigisHF");
0114   }
0115 
0116   for (auto &tag : iConfig.getParameter<std::vector<edm::InputTag>>("phase2barrelCaloTowers")) {
0117     phase2barrelTowers_.push_back(consumes<l1tp2::CaloTowerCollection>(tag));
0118   }
0119 
0120   for (auto &tag : iConfig.getParameter<std::vector<edm::InputTag>>("hcalHGCTowers")) {
0121     hcalHGCTowers_.push_back(consumes<l1t::HGCalTowerBxCollection>(tag));
0122   }
0123   if (!hcalHGCTowers_.empty())
0124     hcalHGCTowersHadOnly_ = iConfig.getParameter<bool>("hcalHGCTowersHadOnly");
0125 }
0126 
0127 // ------------ method called to produce the data  ------------
0128 void L1TPFCaloProducer::produce(edm::Event &iEvent, const edm::EventSetup &iSetup) {
0129   /// ----------------ECAL INFO-------------------
0130   edm::Handle<reco::CandidateView> ecals;
0131   for (const auto &token : ecalCands_) {
0132     iEvent.getByToken(token, ecals);
0133     for (const reco::Candidate &it : *ecals) {
0134       if (debug_)
0135         edm::LogWarning("L1TPFCaloProducer")
0136             << "adding ECal input pt " << it.pt() << ", eta " << it.eta() << ", phi " << it.phi() << "\n";
0137       ecalClusterer_.add(it);
0138     }
0139   }
0140 
0141   /// ----------------HCAL INFO-------------------
0142   if (!ecalOnly_) {
0143     edm::Handle<reco::CandidateView> hcals;
0144     for (const auto &token : hcalCands_) {
0145       iEvent.getByToken(token, hcals);
0146       for (const reco::Candidate &it : *hcals) {
0147         if (debug_)
0148           edm::LogWarning("L1TPFCaloProducer")
0149               << "adding HCal cand input pt " << it.pt() << ", eta " << it.eta() << ", phi " << it.phi() << "\n";
0150         hcalClusterer_.add(it);
0151       }
0152     }
0153     if (!hcalDigis_.empty()) {
0154       readHcalDigis_(iEvent, iSetup);
0155     }
0156     if (!phase2barrelTowers_.empty()) {
0157       readPhase2BarrelCaloTowers_(iEvent, iSetup);
0158     }
0159     if (!hcalHGCTowers_.empty()) {
0160       readHcalHGCTowers_(iEvent, iSetup);
0161     }
0162   }
0163 
0164   /// --------------- CLUSTERING ------------------
0165   ecalClusterer_.run();
0166 
0167   auto ecalCellsH = iEvent.put(ecalClusterer_.fetchCells(), "ecalCells");
0168 
0169   iEvent.put(ecalClusterer_.fetch(ecalCellsH), "emUncalibrated");
0170 
0171   if (emCorrector_.valid()) {
0172     ecalClusterer_.correct(
0173         [&](const l1tpf_calo::Cluster &c) -> float { return emCorrector_.correctedPt(0., c.et, std::abs(c.eta)); });
0174   }
0175 
0176   std::unique_ptr<l1t::PFClusterCollection> corrEcal = ecalClusterer_.fetch(ecalCellsH);
0177 
0178   if (debug_) {
0179     for (const l1t::PFCluster &it : *corrEcal) {
0180       edm::LogWarning("L1TPFCaloProducer")
0181           << "corrected ECal cluster pt " << it.pt() << ", eta " << it.eta() << ", phi " << it.phi() << "\n";
0182     }
0183   }
0184 
0185   auto ecalClustH = iEvent.put(std::move(corrEcal), "emCalibrated");
0186 
0187   if (ecalOnly_) {
0188     ecalClusterer_.clear();
0189     return;
0190   }
0191 
0192   hcalClusterer_.run();
0193 
0194   auto hcalCellsH = iEvent.put(hcalClusterer_.fetchCells(), "hcalCells");
0195 
0196   // this we put separately for debugging
0197   iEvent.put(hcalClusterer_.fetchCells(/*unclustered=*/true), "hcalUnclustered");
0198 
0199   iEvent.put(hcalClusterer_.fetch(hcalCellsH), "hcalUncalibrated");
0200 
0201   if (hcCorrector_.valid()) {
0202     hcalClusterer_.correct(
0203         [&](const l1tpf_calo::Cluster &c) -> float { return hcCorrector_.correctedPt(c.et, 0., std::abs(c.eta)); });
0204   }
0205 
0206   auto hcalClustH = iEvent.put(hcalClusterer_.fetch(hcalCellsH), "hcalCalibrated");
0207 
0208   // Calorimeter linking
0209   caloLinker_->run();
0210 
0211   iEvent.put(caloLinker_->fetch(ecalClustH, hcalClustH), "uncalibrated");
0212 
0213   if (hadCorrector_.valid()) {
0214     caloLinker_->correct([&](const l1tpf_calo::CombinedCluster &c) -> float {
0215       if (debug_)
0216         edm::LogWarning("L1TPFCaloProducer") << "raw linked cluster pt " << c.et << ", eta " << c.eta << ", phi "
0217                                              << c.phi << ", emPt " << c.ecal_et << "\n";
0218       return hadCorrector_.correctedPt(c.et, c.ecal_et, std::abs(c.eta));
0219     });
0220   }
0221 
0222   std::unique_ptr<l1t::PFClusterCollection> clusters = caloLinker_->fetch(ecalClustH, hcalClustH);
0223   for (l1t::PFCluster &c : *clusters) {
0224     c.setPtError(resol_(c.pt(), std::abs(c.eta())));
0225     if (debug_)
0226       edm::LogWarning("L1TPFCaloProducer") << "calibrated linked cluster pt " << c.pt() << ", eta " << c.eta()
0227                                            << ", phi " << c.phi() << ", emPt " << c.emEt() << "\n";
0228   }
0229   iEvent.put(std::move(clusters), "calibrated");
0230 
0231   ecalClusterer_.clear();
0232   hcalClusterer_.clear();
0233   caloLinker_->clear();
0234 }
0235 
0236 void L1TPFCaloProducer::readHcalDigis_(edm::Event &iEvent, const edm::EventSetup &iSetup) {
0237   const auto &decoder = iSetup.getData(decoderTag_);
0238   edm::Handle<HcalTrigPrimDigiCollection> hcalTPs;
0239   for (const auto &token : hcalDigis_) {
0240     iEvent.getByToken(token, hcalTPs);
0241     for (const auto &itr : *hcalTPs) {
0242       HcalTrigTowerDetId id = itr.id();
0243       double et = decoder.hcaletValue(itr.id(), itr.t0());
0244       if (et <= 0)
0245         continue;
0246       float towerEta = l1t::CaloTools::towerEta(id.ieta());
0247       float towerPhi = l1t::CaloTools::towerPhi(id.ieta(), id.iphi());
0248       if (!hcalDigisBarrel_ && std::abs(towerEta) < 2)  // |eta| < 2 => barrel (there's no HE in Phase2)
0249         continue;
0250       if (!hcalDigisHF_ && std::abs(towerEta) > 2)  // |eta| > 2 => HF
0251         continue;
0252       if (debug_)
0253         edm::LogWarning("L1TPFCaloProducer")
0254             << "adding HCal digi input pt " << et << ", eta " << towerEta << ", phi " << towerPhi << "\n";
0255       hcalClusterer_.add(et, towerEta, towerPhi);
0256     }
0257   }
0258 }
0259 
0260 void L1TPFCaloProducer::readPhase2BarrelCaloTowers_(edm::Event &event, const edm::EventSetup &) {
0261   edm::Handle<l1tp2::CaloTowerCollection> towers;
0262   for (const auto &token : phase2barrelTowers_) {
0263     event.getByToken(token, towers);
0264     for (const auto &t : *towers) {
0265       // sanity check from https://github.com/cms-l1t-offline/cmssw/blob/l1t-phase2-v3.0.2/L1Trigger/L1CaloTrigger/plugins/L1TowerCalibrator.cc#L259-L263
0266       if ((int)t.towerIEta() == -1016 && (int)t.towerIPhi() == -962)
0267         continue;
0268       if (debug_ && (t.hcalTowerEt() > 0 || t.ecalTowerEt() > 0)) {
0269         edm::LogWarning("L1TPFCaloProducer")
0270             << "adding phase2 L1 CaloTower eta " << t.towerEta() << "   phi " << t.towerPhi() << "   ieta "
0271             << t.towerIEta() << "   iphi " << t.towerIPhi() << "   ecal " << t.ecalTowerEt() << "    hcal "
0272             << t.hcalTowerEt() << "\n";
0273       }
0274       hcalClusterer_.add(t.hcalTowerEt(), t.towerEta(), t.towerPhi());
0275       ecalClusterer_.add(t.ecalTowerEt(), t.towerEta(), t.towerPhi());
0276     }
0277   }
0278 }
0279 
0280 void L1TPFCaloProducer::readHcalHGCTowers_(edm::Event &iEvent, const edm::EventSetup &iSetup) {
0281   edm::Handle<l1t::HGCalTowerBxCollection> hgcTowers;
0282 
0283   for (const auto &token : hcalHGCTowers_) {
0284     iEvent.getByToken(token, hgcTowers);
0285     for (auto it = hgcTowers->begin(0), ed = hgcTowers->end(0); it != ed; ++it) {
0286       if (debug_)
0287         edm::LogWarning("L1TPFCaloProducer")
0288             << "adding HGC Tower hadEt " << it->etHad() << ", emEt " << it->etEm() << ", pt " << it->pt() << ", eta "
0289             << it->eta() << ", phi " << it->phi() << "\n";
0290       hcalClusterer_.add(it->etHad(), it->eta(), it->phi());
0291       if (!hcalHGCTowersHadOnly_)
0292         ecalClusterer_.add(it->etEm(), it->eta(), it->phi());
0293     }
0294   }
0295 }
0296 
0297 //define this as a plug-in
0298 #include "FWCore/Framework/interface/MakerMacros.h"
0299 DEFINE_FWK_MODULE(L1TPFCaloProducer);