File indexing completed on 2022-03-26 02:43:37
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014 #include <memory>
0015
0016
0017 #include "FWCore/Framework/interface/Frameworkfwd.h"
0018 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0019 #include "FWCore/Framework/interface/Event.h"
0020 #include "FWCore/Framework/interface/EventSetup.h"
0021 #include "FWCore/Framework/interface/MakerMacros.h"
0022 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0023 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0024 #include "FWCore/Framework/interface/ESHandle.h"
0025
0026
0027 #include "JetMETCorrections/JetCorrector/interface/JetCorrector.h"
0028 #include "CondFormats/EcalObjects/interface/EcalChannelStatus.h"
0029 #include "CondFormats/DataRecord/interface/EcalChannelStatusRcd.h"
0030
0031
0032 #include "DataFormats/JetReco/interface/CaloJetCollection.h"
0033 #include "DataFormats/JetReco/interface/PFJetCollection.h"
0034 #include "DataFormats/JetReco/interface/JetID.h"
0035 #include "DataFormats/METReco/interface/PFMETCollection.h"
0036 #include "DataFormats/MuonReco/interface/Muon.h"
0037 #include "DataFormats/MuonReco/interface/MuonFwd.h"
0038 #include "DataFormats/METReco/interface/PFMET.h"
0039 #include "DataFormats/METReco/interface/CaloMETCollection.h"
0040 #include "DataFormats/METReco/interface/CaloMET.h"
0041
0042
0043 #include "FWCore/ServiceRegistry/interface/Service.h"
0044 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0045 #include "TH1.h"
0046 #include "TTree.h"
0047 #include "TF1.h"
0048 #include <TVector2.h>
0049
0050
0051 #include "L1Trigger/L1TNtuples/interface/L1AnalysisRecoJetDataFormat.h"
0052 #include "L1Trigger/L1TNtuples/interface/L1AnalysisRecoMetDataFormat.h"
0053
0054
0055
0056
0057
0058 class L1JetRecoTreeProducer : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0059 public:
0060 explicit L1JetRecoTreeProducer(const edm::ParameterSet&);
0061 ~L1JetRecoTreeProducer() override;
0062
0063 private:
0064 void beginJob(void) override;
0065 void analyze(const edm::Event&, const edm::EventSetup&) override;
0066 void endJob() override;
0067
0068 void doPFJets(edm::Handle<reco::PFJetCollection> pfJets);
0069 void doPFJetCorr(edm::Handle<reco::PFJetCollection> pfJets, edm::Handle<reco::JetCorrector> pfJetCorr);
0070 void doCaloJets(edm::Handle<reco::CaloJetCollection> caloJets);
0071 void doCaloJetCorr(edm::Handle<reco::CaloJetCollection> caloJets, edm::Handle<reco::JetCorrector> caloJetCorr);
0072 void doCaloMet(edm::Handle<reco::CaloMETCollection> caloMet);
0073 void doCaloMetBE(edm::Handle<reco::CaloMETCollection> caloMetBE);
0074
0075 void doPFMet(edm::Handle<reco::PFMETCollection> pfMet);
0076 void doPFMetNoMu(edm::Handle<reco::PFMETCollection> pfMet, edm::Handle<reco::MuonCollection>);
0077
0078 bool pfJetID(const reco::PFJet& jet);
0079 bool caloJetID(const reco::CaloJet& jet);
0080
0081 public:
0082 L1Analysis::L1AnalysisRecoJetDataFormat* jet_data;
0083 L1Analysis::L1AnalysisRecoMetDataFormat* met_data;
0084
0085 private:
0086
0087 TTree* tree_;
0088
0089
0090 edm::EDGetTokenT<reco::PFJetCollection> pfJetToken_;
0091 edm::EDGetTokenT<reco::CaloJetCollection> caloJetToken_;
0092 edm::EDGetTokenT<edm::ValueMap<reco::JetID> > caloJetIDToken_;
0093 edm::EDGetTokenT<reco::JetCorrector> pfJECToken_;
0094 edm::EDGetTokenT<reco::JetCorrector> caloJECToken_;
0095
0096 edm::EDGetTokenT<reco::PFMETCollection> pfMetToken_;
0097 edm::EDGetTokenT<reco::CaloMETCollection> caloMetToken_;
0098 edm::EDGetTokenT<reco::CaloMETCollection> caloMetBEToken_;
0099
0100 edm::EDGetTokenT<reco::MuonCollection> muonToken_;
0101
0102
0103 bool pfJetsMissing_;
0104 double jetptThreshold_;
0105 double jetetaMax_;
0106 unsigned int maxCl_;
0107 unsigned int maxJet_;
0108 unsigned int maxVtx_;
0109 unsigned int maxTrk_;
0110
0111 bool pfJetCorrMissing_;
0112 bool caloJetCorrMissing_;
0113 bool caloJetsMissing_;
0114 bool caloJetIDMissing_;
0115 bool pfMetMissing_;
0116 bool caloMetMissing_;
0117 bool caloMetBEMissing_;
0118
0119 bool muonsMissing_;
0120 };
0121
0122 L1JetRecoTreeProducer::L1JetRecoTreeProducer(const edm::ParameterSet& iConfig)
0123 : pfJetsMissing_(false),
0124 pfJetCorrMissing_(false),
0125 caloJetCorrMissing_(false),
0126 caloJetsMissing_(false),
0127 caloJetIDMissing_(false),
0128 pfMetMissing_(false),
0129 caloMetMissing_(false),
0130 caloMetBEMissing_(false),
0131 muonsMissing_(false) {
0132 caloJetToken_ =
0133 consumes<reco::CaloJetCollection>(iConfig.getUntrackedParameter("caloJetToken", edm::InputTag("ak4CaloJets")));
0134 pfJetToken_ =
0135 consumes<reco::PFJetCollection>(iConfig.getUntrackedParameter("pfJetToken", edm::InputTag("ak4PFJetsCHS")));
0136 caloJetIDToken_ =
0137 consumes<edm::ValueMap<reco::JetID> >(iConfig.getUntrackedParameter("caloJetIDToken", edm::InputTag("ak4JetID")));
0138 pfJECToken_ = consumes<reco::JetCorrector>(
0139 iConfig.getUntrackedParameter<edm::InputTag>("pfJECToken", edm::InputTag("ak4PFCHSL1FastL2L3ResidualCorrector")));
0140 caloJECToken_ = consumes<reco::JetCorrector>(iConfig.getUntrackedParameter<edm::InputTag>(
0141 "caloJECToken", edm::InputTag("ak4CaloL1FastL2L3ResidualCorrector")));
0142
0143 pfMetToken_ = consumes<reco::PFMETCollection>(iConfig.getUntrackedParameter("pfMetToken", edm::InputTag("pfMetT1")));
0144 caloMetToken_ =
0145 consumes<reco::CaloMETCollection>(iConfig.getUntrackedParameter("caloMetToken", edm::InputTag("caloMet")));
0146 caloMetBEToken_ =
0147 consumes<reco::CaloMETCollection>(iConfig.getUntrackedParameter("caloMetBEToken", edm::InputTag("caloMetBE")));
0148
0149 muonToken_ = consumes<reco::MuonCollection>(iConfig.getUntrackedParameter("muonToken", edm::InputTag("muons")));
0150
0151 usesResource(TFileService::kSharedResource);
0152
0153 jetptThreshold_ = iConfig.getParameter<double>("jetptThreshold");
0154 jetetaMax_ = iConfig.getParameter<double>("jetetaMax");
0155 maxJet_ = iConfig.getParameter<unsigned int>("maxJet");
0156
0157 jet_data = new L1Analysis::L1AnalysisRecoJetDataFormat();
0158 met_data = new L1Analysis::L1AnalysisRecoMetDataFormat();
0159
0160
0161 edm::Service<TFileService> fs_;
0162 tree_ = fs_->make<TTree>("JetRecoTree", "JetRecoTree");
0163 tree_->Branch("Jet", "L1Analysis::L1AnalysisRecoJetDataFormat", &jet_data, 32000, 3);
0164 tree_->Branch("Sums", "L1Analysis::L1AnalysisRecoMetDataFormat", &met_data, 32000, 3);
0165 }
0166
0167 L1JetRecoTreeProducer::~L1JetRecoTreeProducer() {
0168
0169
0170 }
0171
0172
0173
0174
0175
0176
0177 void L1JetRecoTreeProducer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0178 jet_data->Reset();
0179 met_data->Reset();
0180
0181
0182 edm::Handle<reco::PFJetCollection> pfJets;
0183 iEvent.getByToken(pfJetToken_, pfJets);
0184
0185
0186 edm::Handle<reco::CaloJetCollection> caloJets;
0187 iEvent.getByToken(caloJetToken_, caloJets);
0188
0189
0190 edm::Handle<reco::PFMETCollection> pfMet;
0191 iEvent.getByToken(pfMetToken_, pfMet);
0192
0193
0194 edm::Handle<edm::ValueMap<reco::JetID> > jetsID;
0195 iEvent.getByToken(caloJetIDToken_, jetsID);
0196
0197 edm::Handle<reco::JetCorrector> pfJetCorr;
0198 iEvent.getByToken(pfJECToken_, pfJetCorr);
0199
0200 edm::Handle<reco::JetCorrector> caloJetCorr;
0201 iEvent.getByToken(caloJECToken_, caloJetCorr);
0202
0203 edm::Handle<reco::CaloMETCollection> caloMet;
0204 iEvent.getByToken(caloMetToken_, caloMet);
0205
0206 edm::Handle<reco::CaloMETCollection> caloMetBE;
0207 iEvent.getByToken(caloMetBEToken_, caloMetBE);
0208
0209
0210 edm::Handle<reco::MuonCollection> muons;
0211 iEvent.getByToken(muonToken_, muons);
0212
0213 if (pfJets.isValid()) {
0214 jet_data->nJets = 0;
0215
0216 doPFJets(pfJets);
0217
0218 } else {
0219 if (!pfJetsMissing_) {
0220 edm::LogWarning("MissingProduct") << "PFJets not found. Branch will not be filled" << std::endl;
0221 }
0222 pfJetsMissing_ = true;
0223 }
0224
0225 if (pfJetCorr.isValid()) {
0226 doPFJetCorr(pfJets, pfJetCorr);
0227
0228 } else {
0229 if (!pfJetCorrMissing_) {
0230 edm::LogWarning("MissingProduct") << "PF Jet Corrector not found. Branch will not be filled" << std::endl;
0231 }
0232 pfJetCorrMissing_ = true;
0233 }
0234
0235 if (caloJets.isValid()) {
0236 jet_data->nCaloJets = 0;
0237
0238 doCaloJets(caloJets);
0239
0240 } else {
0241 if (!caloJetsMissing_) {
0242 edm::LogWarning("MissingProduct") << "Calo Jets not found. Branch will not be filled" << std::endl;
0243 }
0244 caloJetsMissing_ = true;
0245 }
0246
0247 if (caloJetCorr.isValid()) {
0248 doCaloJetCorr(caloJets, caloJetCorr);
0249
0250 } else {
0251 if (!caloJetCorrMissing_) {
0252 edm::LogWarning("MissingProduct") << "Calo Jet Corrector not found. Branch will not be filled" << std::endl;
0253 }
0254 caloJetCorrMissing_ = true;
0255 }
0256
0257 if (!jetsID.isValid()) {
0258 if (!caloJetIDMissing_) {
0259 edm::LogWarning("MissingProduct") << "Calo Jet ID not found. Branch will not be filled" << std::endl;
0260 }
0261 caloJetIDMissing_ = true;
0262 }
0263
0264 if (pfMet.isValid()) {
0265 doPFMet(pfMet);
0266
0267 if (muons.isValid()) {
0268 doPFMetNoMu(pfMet, muons);
0269
0270 } else {
0271 if (!muonsMissing_) {
0272 edm::LogWarning("MissingProduct") << "Muons not found. PFMetNoMu branch will not be filled" << std::endl;
0273 }
0274 muonsMissing_ = true;
0275 }
0276 } else {
0277 if (!pfMetMissing_) {
0278 edm::LogWarning("MissingProduct") << "PFMet not found. Branch will not be filled" << std::endl;
0279 }
0280 pfMetMissing_ = true;
0281 }
0282
0283 if (caloMet.isValid()) {
0284 doCaloMet(caloMet);
0285
0286 } else {
0287 if (!caloMetMissing_) {
0288 edm::LogWarning("MissingProduct") << "CaloMet not found. Branch will not be filled" << std::endl;
0289 }
0290 caloMetMissing_ = true;
0291 }
0292
0293 if (caloMetBE.isValid()) {
0294 doCaloMetBE(caloMetBE);
0295
0296 } else {
0297 if (!caloMetBEMissing_) {
0298 edm::LogWarning("MissingProduct") << "CaloMetBE not found. Branch will not be filled" << std::endl;
0299 }
0300 caloMetBEMissing_ = true;
0301 }
0302
0303 tree_->Fill();
0304 }
0305
0306 void L1JetRecoTreeProducer::doCaloJets(edm::Handle<reco::CaloJetCollection> caloJets) {
0307 for (auto it = caloJets->begin(); it != caloJets->end() && jet_data->nCaloJets < maxJet_; ++it) {
0308 if (!caloJetIDMissing_)
0309 if (!caloJetID(*it))
0310 continue;
0311
0312 jet_data->caloEt.push_back(it->et());
0313 jet_data->caloEta.push_back(it->eta());
0314 jet_data->caloPhi.push_back(it->phi());
0315 jet_data->caloE.push_back(it->energy());
0316
0317 jet_data->eEMF.push_back(it->emEnergyFraction());
0318 jet_data->eEmEB.push_back(it->emEnergyInEB());
0319 jet_data->eEmEE.push_back(it->emEnergyInEE());
0320 jet_data->eEmHF.push_back(it->emEnergyInHF());
0321 jet_data->eHadHB.push_back(it->hadEnergyInHB());
0322 jet_data->eHadHE.push_back(it->hadEnergyInHE());
0323 jet_data->eHadHO.push_back(it->hadEnergyInHO());
0324 jet_data->eHadHF.push_back(it->hadEnergyInHF());
0325 jet_data->eMaxEcalTow.push_back(it->maxEInEmTowers());
0326 jet_data->eMaxHcalTow.push_back(it->maxEInHadTowers());
0327 jet_data->towerArea.push_back(it->towersArea());
0328 jet_data->n60.push_back(it->n60());
0329
0330 jet_data->nCaloJets++;
0331 }
0332 }
0333
0334 void L1JetRecoTreeProducer::doPFJets(edm::Handle<reco::PFJetCollection> pfJets) {
0335 for (auto it = pfJets->begin(); it != pfJets->end() && jet_data->nJets < maxJet_; ++it) {
0336 if (!pfJetID(*it))
0337 continue;
0338
0339 jet_data->et.push_back(it->et());
0340 jet_data->eta.push_back(it->eta());
0341 jet_data->phi.push_back(it->phi());
0342 jet_data->e.push_back(it->energy());
0343
0344 jet_data->chef.push_back(it->chargedHadronEnergyFraction());
0345 jet_data->nhef.push_back(it->neutralHadronEnergyFraction());
0346 jet_data->pef.push_back(it->photonEnergyFraction());
0347 jet_data->eef.push_back(it->electronEnergyFraction());
0348 jet_data->mef.push_back(it->muonEnergyFraction());
0349 jet_data->hfhef.push_back(it->HFHadronEnergyFraction());
0350 jet_data->hfemef.push_back(it->HFEMEnergyFraction());
0351 jet_data->chMult.push_back(it->chargedHadronMultiplicity());
0352 jet_data->nhMult.push_back(it->neutralHadronMultiplicity());
0353 jet_data->phMult.push_back(it->photonMultiplicity());
0354 jet_data->elMult.push_back(it->electronMultiplicity());
0355 jet_data->muMult.push_back(it->muonMultiplicity());
0356 jet_data->hfhMult.push_back(it->HFHadronMultiplicity());
0357 jet_data->hfemMult.push_back(it->HFEMMultiplicity());
0358
0359 jet_data->cemef.push_back(it->chargedEmEnergyFraction());
0360 jet_data->cmef.push_back(it->chargedMuEnergyFraction());
0361 jet_data->nemef.push_back(it->neutralEmEnergyFraction());
0362 jet_data->cMult.push_back(it->chargedMultiplicity());
0363 jet_data->nMult.push_back(it->neutralMultiplicity());
0364
0365 jet_data->nJets++;
0366 }
0367 }
0368
0369 void L1JetRecoTreeProducer::doPFJetCorr(edm::Handle<reco::PFJetCollection> pfJets,
0370 edm::Handle<reco::JetCorrector> pfJetCorr) {
0371 float corrFactor = 1.;
0372 unsigned int nJets = 0;
0373
0374 float mHx = 0;
0375 float mHy = 0;
0376
0377 met_data->Ht = 0;
0378 met_data->mHt = -999.;
0379 met_data->mHtPhi = -999.;
0380
0381 for (auto it = pfJets->begin(); it != pfJets->end() && nJets < maxJet_; ++it) {
0382 if (!pfJetID(*it))
0383 continue;
0384
0385 corrFactor = pfJetCorr.product()->correction(*it);
0386
0387 jet_data->etCorr.push_back(it->et() * corrFactor);
0388 jet_data->corrFactor.push_back(corrFactor);
0389
0390 nJets++;
0391
0392 if (it->pt() * corrFactor > jetptThreshold_ && std::abs(it->eta()) < jetetaMax_) {
0393 mHx += -1. * it->px() * corrFactor;
0394 mHy += -1. * it->py() * corrFactor;
0395 met_data->Ht += it->pt() * corrFactor;
0396 }
0397 }
0398
0399 TVector2 tv2 = TVector2(mHx, mHy);
0400 met_data->mHt = tv2.Mod();
0401 met_data->mHtPhi = tv2.Phi();
0402 }
0403
0404 void L1JetRecoTreeProducer::doCaloJetCorr(edm::Handle<reco::CaloJetCollection> caloJets,
0405 edm::Handle<reco::JetCorrector> caloJetCorr) {
0406 float caloCorrFactor = 1.;
0407 unsigned int nCaloJets = 0;
0408
0409 met_data->caloHt = 0;
0410
0411 for (auto it = caloJets->begin(); it != caloJets->end() && nCaloJets < maxJet_; ++it) {
0412 if (!caloJetIDMissing_)
0413 if (!caloJetID(*it))
0414 continue;
0415
0416 caloCorrFactor = caloJetCorr.product()->correction(*it);
0417
0418 jet_data->caloEtCorr.push_back(it->et() * caloCorrFactor);
0419 jet_data->caloCorrFactor.push_back(caloCorrFactor);
0420
0421 nCaloJets++;
0422
0423 if (it->pt() * caloCorrFactor > jetptThreshold_ && std::abs(it->eta()) < jetetaMax_) {
0424 met_data->caloHt += it->pt() * caloCorrFactor;
0425 }
0426 }
0427 }
0428
0429 void L1JetRecoTreeProducer::doPFMet(edm::Handle<reco::PFMETCollection> pfMet) {
0430 const reco::PFMETCollection* metCol = pfMet.product();
0431 const reco::PFMET theMet = metCol->front();
0432
0433 met_data->met = theMet.et();
0434 met_data->metPhi = theMet.phi();
0435 met_data->sumEt = theMet.sumEt();
0436 met_data->metPx = theMet.px();
0437 met_data->metPy = theMet.py();
0438 }
0439
0440 void L1JetRecoTreeProducer::doPFMetNoMu(edm::Handle<reco::PFMETCollection> pfMet,
0441 edm::Handle<reco::MuonCollection> muons) {
0442 const reco::PFMETCollection* metCol = pfMet.product();
0443 const reco::PFMET theMet = metCol->front();
0444 reco::PFMET thePFMetNoMu = metCol->front();
0445
0446 double pfMetNoMuPx = theMet.px();
0447 double pfMetNoMuPy = theMet.py();
0448
0449 double muPx(0.), muPy(0.);
0450
0451 for (auto it = muons->begin(); it != muons->end(); ++it) {
0452 if (it->isPFMuon()) {
0453 muPx += it->px();
0454 muPy += it->py();
0455 }
0456 }
0457
0458 pfMetNoMuPx += muPx;
0459 pfMetNoMuPy += muPy;
0460
0461 math::XYZTLorentzVector pfMetNoMuP4(pfMetNoMuPx, pfMetNoMuPy, 0, hypot(pfMetNoMuPx, pfMetNoMuPy));
0462
0463 thePFMetNoMu.setP4(pfMetNoMuP4);
0464
0465 met_data->pfMetNoMu = thePFMetNoMu.et();
0466 met_data->pfMetNoMuPhi = thePFMetNoMu.phi();
0467 met_data->pfMetNoMuPx = thePFMetNoMu.px();
0468 met_data->pfMetNoMuPy = thePFMetNoMu.py();
0469 }
0470
0471 void L1JetRecoTreeProducer::doCaloMet(edm::Handle<reco::CaloMETCollection> caloMet) {
0472 const reco::CaloMETCollection* metCol = caloMet.product();
0473 const reco::CaloMET theMet = metCol->front();
0474
0475 met_data->caloMet = theMet.et();
0476 met_data->caloMetPhi = theMet.phi();
0477 met_data->caloSumEt = theMet.sumEt();
0478 }
0479
0480 void L1JetRecoTreeProducer::doCaloMetBE(edm::Handle<reco::CaloMETCollection> caloMetBE) {
0481 const reco::CaloMETCollection* metCol = caloMetBE.product();
0482 const reco::CaloMET theMet = metCol->front();
0483
0484 met_data->caloMetBE = theMet.et();
0485 met_data->caloMetPhiBE = theMet.phi();
0486 met_data->caloSumEtBE = theMet.sumEt();
0487 }
0488
0489 bool L1JetRecoTreeProducer::pfJetID(const reco::PFJet& jet) {
0490 bool tmp = true;
0491 if (std::abs(jet.eta()) < 2.7) {
0492 tmp &= jet.neutralHadronEnergyFraction() < 0.9;
0493 tmp &= jet.neutralEmEnergyFraction() < 0.9;
0494 tmp &= (jet.chargedMultiplicity() + jet.neutralMultiplicity()) > 1;
0495 tmp &= jet.muonEnergyFraction() < 0.8;
0496 tmp &= jet.chargedHadronEnergyFraction() > 0.0;
0497 tmp &= jet.chargedMultiplicity() > 0;
0498 tmp &= jet.chargedEmEnergyFraction() < 0.9;
0499 }
0500 if (std::abs(jet.eta()) > 2.7 && std::abs(jet.eta()) < 3.0) {
0501 tmp &= jet.neutralEmEnergyFraction() > 0.01;
0502 tmp &= jet.neutralHadronEnergyFraction() < 0.98;
0503 tmp &= jet.neutralMultiplicity() > 2;
0504 }
0505 if (std::abs(jet.eta()) > 3.0) {
0506 tmp &= jet.neutralEmEnergyFraction() < 0.9;
0507 tmp &= jet.neutralMultiplicity() > 10;
0508 }
0509
0510
0511
0512
0513
0514 return tmp;
0515 }
0516
0517 bool L1JetRecoTreeProducer::caloJetID(const reco::CaloJet& jet) {
0518 bool tmp = true;
0519
0520 return tmp;
0521 }
0522
0523
0524 void L1JetRecoTreeProducer::beginJob(void) {}
0525
0526
0527 void L1JetRecoTreeProducer::endJob() {}
0528
0529
0530 DEFINE_FWK_MODULE(L1JetRecoTreeProducer);
