File indexing completed on 2024-04-06 12:08:07
0001
0002
0003
0004
0005
0006
0007 #include "DQM/Physics/src/QcdPhotonsDQM.h"
0008
0009 #include "FWCore/Framework/interface/Event.h"
0010 #include "FWCore/Framework/interface/EventSetup.h"
0011 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0012
0013 #include "DQMServices/Core/interface/DQMStore.h"
0014 #include "FWCore/ServiceRegistry/interface/Service.h"
0015 #include "FWCore/Common/interface/TriggerNames.h"
0016
0017 #include "DataFormats/Common/interface/Handle.h"
0018
0019 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0020
0021
0022 #include "DataFormats/EgammaCandidates/interface/Photon.h"
0023 #include "DataFormats/JetReco/interface/Jet.h"
0024
0025
0026 #include "DataFormats/VertexReco/interface/Vertex.h"
0027
0028
0029 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgo.h"
0030 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgoRcd.h"
0031 #include "FWCore/ServiceRegistry/interface/Service.h"
0032
0033 #include "FWCore/Framework/interface/ESHandle.h"
0034 #include "CondFormats/EcalObjects/interface/EcalCondObjectContainer.h"
0035
0036
0037 #include "DataFormats/Math/interface/deltaR.h"
0038 #include "DataFormats/Math/interface/deltaPhi.h"
0039
0040 #include <vector>
0041
0042 #include <string>
0043 #include <cmath>
0044 using namespace std;
0045 using namespace edm;
0046 using namespace reco;
0047
0048 QcdPhotonsDQM::QcdPhotonsDQM(const ParameterSet& parameters) : ecalClusterToolsESGetTokens_{consumesCollector()} {
0049
0050 theTriggerPathToPass_ = parameters.getParameter<string>("triggerPathToPass");
0051 thePlotTheseTriggersToo_ = parameters.getParameter<vector<string> >("plotTheseTriggersToo");
0052 theJetCollectionLabel_ = parameters.getParameter<InputTag>("jetCollection");
0053 trigTagToken_ = consumes<edm::TriggerResults>(parameters.getUntrackedParameter<edm::InputTag>("trigTag"));
0054 thePhotonCollectionToken_ = consumes<reco::PhotonCollection>(parameters.getParameter<InputTag>("photonCollection"));
0055 theJetCollectionToken_ = consumes<edm::View<reco::Jet> >(parameters.getParameter<InputTag>("jetCollection"));
0056 theVertexCollectionToken_ = consumes<reco::VertexCollection>(parameters.getParameter<InputTag>("vertexCollection"));
0057 theMinJetPt_ = parameters.getParameter<double>("minJetPt");
0058 theMinPhotonEt_ = parameters.getParameter<double>("minPhotonEt");
0059 theRequirePhotonFound_ = parameters.getParameter<bool>("requirePhotonFound");
0060 thePlotPhotonMaxEt_ = parameters.getParameter<double>("plotPhotonMaxEt");
0061 thePlotPhotonMaxEta_ = parameters.getParameter<double>("plotPhotonMaxEta");
0062 thePlotJetMaxEta_ = parameters.getParameter<double>("plotJetMaxEta");
0063 theBarrelRecHitTag_ = parameters.getParameter<InputTag>("barrelRecHitTag");
0064 theEndcapRecHitTag_ = parameters.getParameter<InputTag>("endcapRecHitTag");
0065 theBarrelRecHitToken_ = consumes<EcalRecHitCollection>(parameters.getParameter<InputTag>("barrelRecHitTag"));
0066 theEndcapRecHitToken_ = consumes<EcalRecHitCollection>(parameters.getParameter<InputTag>("endcapRecHitTag"));
0067
0068
0069 h_deltaEt_photon_jet = nullptr;
0070 h_deltaPhi_jet_jet2 = nullptr;
0071 h_deltaPhi_photon_jet = nullptr;
0072 h_deltaPhi_photon_jet2 = nullptr;
0073 h_deltaR_jet_jet2 = nullptr;
0074 h_deltaR_photon_jet2 = nullptr;
0075 h_jet2_eta = nullptr;
0076 h_jet2_pt = nullptr;
0077 h_jet2_ptOverPhotonEt = nullptr;
0078 h_jet_count = nullptr;
0079 h_jet_eta = nullptr;
0080 h_jet_pt = nullptr;
0081 h_photon_count_bar = nullptr;
0082 h_photon_count_end = nullptr;
0083 h_photon_et = nullptr;
0084 h_photon_et_beforeCuts = nullptr;
0085 h_photon_et_jetco = nullptr;
0086 h_photon_et_jetcs = nullptr;
0087 h_photon_et_jetfo = nullptr;
0088 h_photon_et_jetfs = nullptr;
0089 h_photon_eta = nullptr;
0090 h_triggers_passed = nullptr;
0091 }
0092
0093 QcdPhotonsDQM::~QcdPhotonsDQM() {}
0094
0095 void QcdPhotonsDQM::bookHistograms(DQMStore::IBooker& ibooker, edm::Run const&, edm::EventSetup const&) {
0096 logTraceName = "QcdPhotonAnalyzer";
0097
0098 LogTrace(logTraceName) << "Parameters initialization";
0099
0100 ibooker.setCurrentFolder("Physics/QcdPhotons");
0101
0102 std::stringstream aStringStream;
0103 std::string aString;
0104 aStringStream << theMinJetPt_;
0105 aString = aStringStream.str();
0106
0107
0108 int numOfTriggersToMonitor = thePlotTheseTriggersToo_.size();
0109 h_triggers_passed = ibooker.book1D(
0110 "triggers_passed", "Events passing these trigger paths", numOfTriggersToMonitor, 0, numOfTriggersToMonitor);
0111 for (int i = 0; i < numOfTriggersToMonitor; i++) {
0112 h_triggers_passed->setBinLabel(i + 1, thePlotTheseTriggersToo_[i]);
0113 }
0114
0115
0116 h_photon_et_beforeCuts = ibooker.book1D(
0117 "photon_et_beforeCuts", "#gamma with highest E_{T};E_{T}(#gamma) (GeV)", 20, 0., thePlotPhotonMaxEt_);
0118 h_photon_et =
0119 ibooker.book1D("photon_et", "#gamma with highest E_{T};E_{T}(#gamma) (GeV)", 20, 0., thePlotPhotonMaxEt_);
0120 h_photon_eta = ibooker.book1D(
0121 "photon_eta", "#gamma with highest E_{T};#eta(#gamma)", 40, -thePlotPhotonMaxEta_, thePlotPhotonMaxEta_);
0122 h_photon_count_bar = ibooker.book1D(
0123 "photon_count_bar", "Number of #gamma's passing selection (Barrel);Number of #gamma's", 8, -0.5, 7.5);
0124 h_photon_count_end = ibooker.book1D(
0125 "photon_count_end", "Number of #gamma's passing selection (Endcap);Number of #gamma's", 8, -0.5, 7.5);
0126 h_jet_pt =
0127 ibooker.book1D("jet_pt",
0128 "Jet with highest p_{T} (from " + theJetCollectionLabel_.label() + ");p_{T}(1^{st} jet) (GeV)",
0129 20,
0130 0.,
0131 thePlotPhotonMaxEt_);
0132 h_jet_eta = ibooker.book1D("jet_eta",
0133 "Jet with highest p_{T} (from " + theJetCollectionLabel_.label() + ");#eta(1^{st} jet)",
0134 20,
0135 -thePlotJetMaxEta_,
0136 thePlotJetMaxEta_);
0137 h_deltaPhi_photon_jet =
0138 ibooker.book1D("deltaPhi_photon_jet",
0139 "#Delta#phi between Highest E_{T} #gamma and jet;#Delta#phi(#gamma,1^{st} jet)",
0140 20,
0141 0,
0142 3.1415926);
0143 h_deltaPhi_jet_jet2 = ibooker.book1D("deltaPhi_jet_jet2",
0144 "#Delta#phi between Highest E_{T} jet and 2^{nd} "
0145 "jet;#Delta#phi(1^{st} jet,2^{nd} jet)",
0146 20,
0147 0,
0148 3.1415926);
0149 h_deltaEt_photon_jet = ibooker.book1D("deltaEt_photon_jet",
0150 "(E_{T}(#gamma)-p_{T}(jet))/E_{T}(#gamma) when #Delta#phi(#gamma,1^{st} "
0151 "jet) > 2.8;#DeltaE_{T}(#gamma,1^{st} jet)/E_{T}(#gamma)",
0152 20,
0153 -1.0,
0154 1.0);
0155 h_jet_count =
0156 ibooker.book1D("jet_count",
0157 "Number of " + theJetCollectionLabel_.label() + " (p_{T} > " + aString + " GeV);Number of Jets",
0158 8,
0159 -0.5,
0160 7.5);
0161 h_jet2_pt = ibooker.book1D(
0162 "jet2_pt",
0163 "Jet with 2^{nd} highest p_{T} (from " + theJetCollectionLabel_.label() + ");p_{T}(2^{nd} jet) (GeV)",
0164 20,
0165 0.,
0166 thePlotPhotonMaxEt_);
0167 h_jet2_eta =
0168 ibooker.book1D("jet2_eta",
0169 "Jet with 2^{nd} highest p_{T} (from " + theJetCollectionLabel_.label() + ");#eta(2^{nd} jet)",
0170 20,
0171 -thePlotJetMaxEta_,
0172 thePlotJetMaxEta_);
0173 h_jet2_ptOverPhotonEt = ibooker.book1D(
0174 "jet2_ptOverPhotonEt", "p_{T}(2^{nd} highest jet) / E_{T}(#gamma);p_{T}(2^{nd} Jet)/E_{T}(#gamma)", 20, 0.0, 4.0);
0175 h_deltaPhi_photon_jet2 = ibooker.book1D("deltaPhi_photon_jet2",
0176 "#Delta#phi between Highest E_{T} #gamma and 2^{nd} "
0177 "highest jet;#Delta#phi(#gamma,2^{nd} jet)",
0178 20,
0179 0,
0180 3.1415926);
0181 h_deltaR_jet_jet2 = ibooker.book1D(
0182 "deltaR_jet_jet2", "#DeltaR between Highest Jet and 2^{nd} Highest;#DeltaR(1^{st} jet,2^{nd} jet)", 30, 0, 6.0);
0183 h_deltaR_photon_jet2 = ibooker.book1D("deltaR_photon_jet2",
0184 "#DeltaR between Highest E_{T} #gamma and 2^{nd} "
0185 "jet;#DeltaR(#gamma, 2^{nd} jet)",
0186 30,
0187 0,
0188 6.0);
0189
0190
0191 Float_t bins_et[] = {15, 20, 30, 50, 80};
0192 int num_bins_et = 4;
0193 h_photon_et_jetcs = ibooker.book1D("photon_et_jetcs",
0194 "#gamma with highest E_{T} (#eta(jet)<1.45, "
0195 "#eta(#gamma)#eta(jet)>0);E_{T}(#gamma) (GeV)",
0196 num_bins_et,
0197 bins_et);
0198 h_photon_et_jetco = ibooker.book1D("photon_et_jetco",
0199 "#gamma with highest E_{T} (#eta(jet)<1.45, "
0200 "#eta(#gamma)#eta(jet)<0);E_{T}(#gamma) (GeV)",
0201 num_bins_et,
0202 bins_et);
0203 h_photon_et_jetfs = ibooker.book1D("photon_et_jetfs",
0204 "#gamma with highest E_{T} (1.55<#eta(jet)<2.5, "
0205 "#eta(#gamma)#eta(jet)>0);E_{T}(#gamma) (GeV)",
0206 num_bins_et,
0207 bins_et);
0208 h_photon_et_jetfo = ibooker.book1D("photon_et_jetfo",
0209 "#gamma with highest E_{T} (1.55<#eta(jet)<2.5, "
0210 "#eta(#gamma)#eta(jet)<0);E_{T}(#gamma) (GeV)",
0211 num_bins_et,
0212 bins_et);
0213
0214 auto setSumw2 = [](MonitorElement* me) {
0215 if (me->getTH1F()->GetSumw2N() == 0) {
0216 me->enableSumw2();
0217 }
0218 };
0219
0220 setSumw2(h_photon_et_jetcs);
0221 setSumw2(h_photon_et_jetco);
0222 setSumw2(h_photon_et_jetfs);
0223 setSumw2(h_photon_et_jetfo);
0224 }
0225
0226 void QcdPhotonsDQM::analyze(const Event& iEvent, const EventSetup& iSetup) {
0227 LogTrace(logTraceName) << "Analysis of event # ";
0228
0229
0230
0231 Handle<TriggerResults> HLTresults;
0232 iEvent.getByToken(trigTagToken_, HLTresults);
0233 if (!HLTresults.isValid()) {
0234
0235 return;
0236 }
0237 const edm::TriggerNames& trigNames = iEvent.triggerNames(*HLTresults);
0238
0239 bool passed_HLT = false;
0240
0241
0242
0243 for (unsigned int i = 0; i < thePlotTheseTriggersToo_.size(); i++) {
0244 passed_HLT = false;
0245 for (unsigned int ti = 0; (ti < trigNames.size()) && !passed_HLT; ++ti) {
0246 size_t pos = trigNames.triggerName(ti).find(thePlotTheseTriggersToo_[i]);
0247 if (pos == 0)
0248 passed_HLT = HLTresults->accept(ti);
0249 }
0250 if (passed_HLT)
0251 h_triggers_passed->Fill(i);
0252 }
0253
0254
0255 Handle<PhotonCollection> photonCollection;
0256 iEvent.getByToken(thePhotonCollectionToken_, photonCollection);
0257
0258
0259 if (!photonCollection.isValid())
0260 return;
0261
0262
0263 passed_HLT = false;
0264 {
0265
0266 for (unsigned int ti = 0; ti < trigNames.size(); ++ti) {
0267 size_t pos = trigNames.triggerName(ti).find(theTriggerPathToPass_);
0268 if (pos == 0) {
0269 passed_HLT = HLTresults->accept(ti);
0270
0271 break;
0272 }
0273 }
0274
0275
0276 for (PhotonCollection::const_iterator recoPhoton = photonCollection->begin(); recoPhoton != photonCollection->end();
0277 recoPhoton++) {
0278
0279 if (recoPhoton->et() < theMinPhotonEt_)
0280 break;
0281
0282 h_photon_et_beforeCuts->Fill(recoPhoton->et());
0283 break;
0284 }
0285
0286 if (!passed_HLT) {
0287 return;
0288 }
0289 }
0290
0291
0292
0293
0294
0295
0296
0297
0298
0299 Handle<VertexCollection> vertexHandle;
0300 iEvent.getByToken(theVertexCollectionToken_, vertexHandle);
0301 VertexCollection vertexCollection = *(vertexHandle.product());
0302
0303
0304
0305
0306
0307
0308
0309
0310
0311
0312 int nvvertex = 0;
0313 for (unsigned int i = 0; i < vertexCollection.size(); ++i) {
0314 if (vertexCollection[i].isValid())
0315 nvvertex++;
0316 }
0317 if (nvvertex == 0)
0318 return;
0319
0320
0321
0322
0323
0324
0325
0326
0327
0328
0329
0330
0331 Handle<EcalRecHitCollection> EBReducedRecHits;
0332 iEvent.getByToken(theBarrelRecHitToken_, EBReducedRecHits);
0333 Handle<EcalRecHitCollection> EEReducedRecHits;
0334 iEvent.getByToken(theEndcapRecHitToken_, EEReducedRecHits);
0335 EcalClusterLazyTools lazyTool(
0336 iEvent, ecalClusterToolsESGetTokens_.get(iSetup), theBarrelRecHitToken_, theEndcapRecHitToken_);
0337
0338
0339 float photon_et = -9.0;
0340 float photon_eta = -9.0;
0341 float photon_phi = -9.0;
0342 bool photon_passPhotonID = false;
0343 bool found_lead_pho = false;
0344 int photon_count_bar = 0;
0345 int photon_count_end = 0;
0346
0347
0348 auto pho_maxet = std::max_element(
0349 photonCollection->begin(),
0350 photonCollection->end(),
0351 [](const PhotonCollection::value_type& a, const PhotonCollection::value_type& b) { return a.et() < b.et(); });
0352 if (pho_maxet != photonCollection->end() && pho_maxet->et() >= theMinPhotonEt_) {
0353
0354
0355
0356
0357
0358
0359
0360
0361
0362
0363
0364
0365
0366
0367
0368
0369
0370
0371
0372
0373
0374
0375
0376
0377
0378 bool pho_current_passPhotonID = false;
0379 bool pho_current_isEB = pho_maxet->isEB();
0380 bool pho_current_isEE = pho_maxet->isEE();
0381
0382 if (pho_current_isEB && (pho_maxet->sigmaIetaIeta() < 0.01 || pho_maxet->hadronicOverEm() < 0.05)) {
0383
0384 pho_current_passPhotonID = true;
0385 photon_count_bar++;
0386 } else if (pho_current_isEE && (pho_maxet->hadronicOverEm() < 0.05)) {
0387
0388 pho_current_passPhotonID = true;
0389 photon_count_end++;
0390 }
0391
0392 if (!found_lead_pho) {
0393 found_lead_pho = true;
0394 photon_passPhotonID = pho_current_passPhotonID;
0395 photon_et = pho_maxet->et();
0396 photon_eta = pho_maxet->eta();
0397 photon_phi = pho_maxet->phi();
0398 }
0399 }
0400
0401
0402
0403
0404 if (theRequirePhotonFound_ && (!photon_passPhotonID || photon_et < theMinPhotonEt_))
0405 return;
0406
0407
0408
0409 Handle<View<Jet> > jetCollection;
0410 iEvent.getByToken(theJetCollectionToken_, jetCollection);
0411 if (!jetCollection.isValid())
0412 return;
0413
0414 float jet_pt = -8.0;
0415 float jet_eta = -8.0;
0416 float jet_phi = -8.0;
0417 int jet_count = 0;
0418 float jet2_pt = -9.0;
0419 float jet2_eta = -9.0;
0420 float jet2_phi = -9.0;
0421
0422 for (unsigned int i_jet = 0; i_jet < jetCollection->size(); i_jet++) {
0423 const Jet* jet = &jetCollection->at(i_jet);
0424
0425 float jet_current_pt = jet->pt();
0426
0427
0428 if (deltaR(jet->eta(), jet->phi(), photon_eta, photon_phi) < 0.5)
0429 continue;
0430
0431 if (jet_current_pt < theMinJetPt_)
0432 break;
0433
0434 jet_count++;
0435 if (jet_current_pt > jet_pt) {
0436 jet2_pt = jet_pt;
0437 jet2_eta = jet_eta;
0438 jet2_phi = jet_phi;
0439 jet_pt = jet_current_pt;
0440 jet_eta = jet->eta();
0441 jet_phi = jet->phi();
0442 } else if (jet_current_pt > jet2_pt) {
0443 jet2_pt = jet_current_pt;
0444 jet2_eta = jet->eta();
0445 jet2_phi = jet->phi();
0446 }
0447 }
0448
0449
0450
0451
0452
0453
0454 if (jet_pt > 0.0) {
0455
0456 h_photon_et->Fill(photon_et);
0457 h_photon_eta->Fill(photon_eta);
0458 h_photon_count_bar->Fill(photon_count_bar);
0459 h_photon_count_end->Fill(photon_count_end);
0460
0461
0462 if (fabs(photon_eta) < 1.45 && photon_passPhotonID) {
0463 if (fabs(jet_eta) < 1.45) {
0464 if (photon_eta * jet_eta > 0) {
0465 h_photon_et_jetcs->Fill(photon_et);
0466 } else {
0467 h_photon_et_jetco->Fill(photon_et);
0468 }
0469 } else if (jet_eta > 1.55 && jet_eta < 2.5) {
0470 if (photon_eta * jet_eta > 0) {
0471 h_photon_et_jetfs->Fill(photon_et);
0472 } else {
0473 h_photon_et_jetfo->Fill(photon_et);
0474 }
0475 }
0476 }
0477
0478
0479 h_jet_pt->Fill(jet_pt);
0480 h_jet_eta->Fill(jet_eta);
0481 h_jet_count->Fill(jet_count);
0482 h_deltaPhi_photon_jet->Fill(abs(deltaPhi(photon_phi, jet_phi)));
0483 if (abs(deltaPhi(photon_phi, jet_phi)) > 2.8)
0484 h_deltaEt_photon_jet->Fill((photon_et - jet_pt) / photon_et);
0485
0486
0487 if (jet2_pt > 0.0) {
0488 h_jet2_pt->Fill(jet2_pt);
0489 h_jet2_eta->Fill(jet2_eta);
0490 h_jet2_ptOverPhotonEt->Fill(jet2_pt / photon_et);
0491 h_deltaPhi_photon_jet2->Fill(abs(deltaPhi(photon_phi, jet2_phi)));
0492 h_deltaPhi_jet_jet2->Fill(abs(deltaPhi(jet_phi, jet2_phi)));
0493 h_deltaR_jet_jet2->Fill(deltaR(jet_eta, jet_phi, jet2_eta, jet2_phi));
0494 h_deltaR_photon_jet2->Fill(deltaR(photon_eta, photon_phi, jet2_eta, jet2_phi));
0495 }
0496 }
0497
0498
0499 }
0500
0501
0502
0503
0504