Physics/src/TopDiLeptonOfflineDQM.cc

0001 //#include <algorithm>
0002 #include "JetMETCorrections/Objects/interface/JetCorrectionsRecord.h"
0003 #include "DataFormats/JetReco/interface/CaloJet.h"
0004 #include "DataFormats/BTauReco/interface/JetTag.h"
0005 #include "DataFormats/Common/interface/Handle.h"
0006 #include "DQM/Physics/src/TopDiLeptonOfflineDQM.h"
0007
0008 #include <memory>
0009
0010 #include "DQM/Physics/interface/TopDQMHelpers.h"
0011 #include "DataFormats/JetReco/interface/PFJet.h"
0012 #include "FWCore/Framework/interface/ConsumesCollector.h"
0013 #include "FWCore/Framework/interface/EDConsumerBase.h"
0014 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0015 #include "FWCore/Utilities/interface/EDGetToken.h"
0016 #include "FWCore/Utilities/interface/InputTag.h"
0017
0018 namespace TopDiLeptonOffline {
0019
0020   MonitorEnsemble::MonitorEnsemble(const char* label, const edm::ParameterSet& cfg, edm::ConsumesCollector&& iC)
0021       : label_(label),
0022         eidCutValue_(0.),
0023         elecIso_(nullptr),
0024         elecSelect_(nullptr),
0025         muonIso_(nullptr),
0026         muonSelect_(nullptr),
0027         jetIDSelect_(nullptr),
0028         lowerEdge_(-1.),
0029         upperEdge_(-1.),
0030         elecMuLogged_(0),
0031         diMuonLogged_(0),
0032         diElecLogged_(0) {
0033     // sources have to be given; this PSet is not optional
0034     edm::ParameterSet sources = cfg.getParameter<edm::ParameterSet>("sources");
0035     muons_ = iC.consumes<edm::View<reco::PFCandidate>>(sources.getParameter<edm::InputTag>("muons"));
0036     elecs_ = iC.consumes<edm::View<reco::PFCandidate>>(sources.getParameter<edm::InputTag>("elecs"));
0037     jets_ = iC.consumes<edm::View<reco::Jet>>(sources.getParameter<edm::InputTag>("jets"));
0038     for (edm::InputTag const& tag : sources.getParameter<std::vector<edm::InputTag>>("mets"))
0039       mets_.push_back(iC.consumes<edm::View<reco::MET>>(tag));
0040
0041     // elecExtras are optional; they may be omitted or empty
0042     if (cfg.existsAs<edm::ParameterSet>("elecExtras")) {
0043       edm::ParameterSet elecExtras = cfg.getParameter<edm::ParameterSet>("elecExtras");
0044       // select is optional; in case it's not found no
0045       // selection will be applied
0046       if (elecExtras.existsAs<std::string>("select")) {
0047         elecSelect_ = std::make_unique<StringCutObjectSelector<reco::PFCandidate>>(
0048             elecExtras.getParameter<std::string>("select"));
0049       }
0050       // isolation is optional; in case it's not found no
0051       // isolation will be applied
0052       if (elecExtras.existsAs<std::string>("isolation")) {
0053         elecIso_ = std::make_unique<StringCutObjectSelector<reco::PFCandidate>>(
0054             elecExtras.getParameter<std::string>("isolation"));
0055       }
0056       // electronId is optional; in case it's not found the
0057       // InputTag will remain empty
0058       if (elecExtras.existsAs<edm::ParameterSet>("electronId")) {
0059         edm::ParameterSet elecId = elecExtras.getParameter<edm::ParameterSet>("electronId");
0060         electronId_ = iC.consumes<edm::ValueMap<float>>(elecId.getParameter<edm::InputTag>("src"));
0061         eidCutValue_ = elecId.getParameter<double>("cutValue");
0062         //  eidPattern_= elecId.getParameter<int>("pattern");
0063       }
0064     }
0065     // muonExtras are optional; they may be omitted or empty
0066     if (cfg.existsAs<edm::ParameterSet>("muonExtras")) {
0067       edm::ParameterSet muonExtras = cfg.getParameter<edm::ParameterSet>("muonExtras");
0068       // select is optional; in case it's not found no
0069       // selection will be applied
0070       if (muonExtras.existsAs<std::string>("select")) {
0071         muonSelect_ = std::make_unique<StringCutObjectSelector<reco::PFCandidate, true>>(
0072             muonExtras.getParameter<std::string>("select"));
0073       }
0074       // isolation is optional; in case it's not found no
0075       // isolation will be applied
0076       if (muonExtras.existsAs<std::string>("isolation")) {
0077         muonIso_ = std::make_unique<StringCutObjectSelector<reco::PFCandidate, true>>(
0078             muonExtras.getParameter<std::string>("isolation"));
0079       }
0080     }
0081     // jetExtras are optional; they may be omitted or empty
0082     if (cfg.existsAs<edm::ParameterSet>("jetExtras")) {
0083       edm::ParameterSet jetExtras = cfg.getParameter<edm::ParameterSet>("jetExtras");
0084       // jetCorrector is optional; in case it's not found
0085       // the InputTag will remain empty
0086       if (jetExtras.existsAs<std::string>("jetCorrector")) {
0087         jetCorrector_ =
0088             iC.consumes<reco::JetCorrector>(edm::InputTag(jetExtras.getParameter<std::string>("jetCorrector")));
0089       }
0090       // read jetID information if it exists
0091       if (jetExtras.existsAs<edm::ParameterSet>("jetID")) {
0092         edm::ParameterSet jetID = jetExtras.getParameter<edm::ParameterSet>("jetID");
0093         jetIDLabel_ = iC.consumes<reco::JetIDValueMap>(jetID.getParameter<edm::InputTag>("label"));
0094         jetIDSelect_ =
0095             std::make_unique<StringCutObjectSelector<reco::JetID>>(jetID.getParameter<std::string>("select"));
0096       }
0097       // select is optional; in case it's not found no
0098       // selection will be applied (only implemented for
0099       // CaloJets at the moment)
0100       if (jetExtras.existsAs<std::string>("select")) {
0101         jetSelect_ = jetExtras.getParameter<std::string>("select");
0102       }
0103     }
0104     // triggerExtras are optional; they may be omitted or empty
0105     if (cfg.existsAs<edm::ParameterSet>("triggerExtras")) {
0106       edm::ParameterSet triggerExtras = cfg.getParameter<edm::ParameterSet>("triggerExtras");
0107       triggerTable_ = iC.consumes<edm::TriggerResults>(triggerExtras.getParameter<edm::InputTag>("src"));
0108       elecMuPaths_ = triggerExtras.getParameter<std::vector<std::string>>("pathsELECMU");
0109       diMuonPaths_ = triggerExtras.getParameter<std::vector<std::string>>("pathsDIMUON");
0110     }
0111     // massExtras is optional; in case it's not found no mass
0112     // window cuts are applied for the same flavor monitor
0113     // histograms
0114     if (cfg.existsAs<edm::ParameterSet>("massExtras")) {
0115       edm::ParameterSet massExtras = cfg.getParameter<edm::ParameterSet>("massExtras");
0116       lowerEdge_ = massExtras.getParameter<double>("lowerEdge");
0117       upperEdge_ = massExtras.getParameter<double>("upperEdge");
0118     }
0119
0120     // setup the verbosity level for booking histograms;
0121     // per default the verbosity level will be set to
0122     // STANDARD. This will also be the chosen level in
0123     // the case when the monitoring PSet is not found
0124     verbosity_ = STANDARD;
0125     if (cfg.existsAs<edm::ParameterSet>("monitoring")) {
0126       edm::ParameterSet monitoring = cfg.getParameter<edm::ParameterSet>("monitoring");
0127       if (monitoring.getParameter<std::string>("verbosity") == "DEBUG")
0128         verbosity_ = DEBUG;
0129       if (monitoring.getParameter<std::string>("verbosity") == "VERBOSE")
0130         verbosity_ = VERBOSE;
0131       if (monitoring.getParameter<std::string>("verbosity") == "STANDARD")
0132         verbosity_ = STANDARD;
0133     }
0134     // and don't forget to do the histogram booking
0135     directory_ = cfg.getParameter<std::string>("directory");
0136   }
0137
0138   void MonitorEnsemble::book(DQMStore::IBooker& ibooker) {
0139     // set up the current directory path
0140     std::string current(directory_);
0141     current += label_;
0142     ibooker.setCurrentFolder(current);
0143
0144     // determine number of bins for trigger monitoring
0145     unsigned int nElecMu = elecMuPaths_.size();
0146     unsigned int nDiMuon = diMuonPaths_.size();
0147
0148     // --- [STANDARD] --- //
0149     // Run Number
0150     hists_["RunNumb_"] = ibooker.book1D("RunNumber", "Run Nr.", 1.e4, 1.5e5, 3.e5);
0151     // invariant mass of opposite charge lepton pair (only filled for same flavor)
0152     hists_["invMass_"] = ibooker.book1D("InvMass", "M(lep1, lep2)", 80, 0., 320.);
0153     // invariant mass of opposite charge lepton pair (only filled for same flavor)
0154     hists_["invMassLog_"] = ibooker.book1D("InvMassLog", "log_{10}(M(lep1, lep2))", 80, .1, 2.5);
0155     // invariant mass of same charge lepton pair (log10 for low mass region, only
0156     // filled for same flavor)
0157     hists_["invMassWC_"] = ibooker.book1D("InvMassWC", "M_{WC}(L1, L2)", 80, 0., 320.);
0158     // invariant mass of same charge lepton pair (log10 for low mass region, only
0159     // filled for same flavor)
0160     hists_["invMassWCLog_"] = ibooker.book1D("InvMassLogWC", "log_{10}(M_{WC})", 80, .1, 2.5);
0161     // decay channel [1]: muon/muon, [2]:elec/elec, [3]:elec/muon
0162     hists_["decayChannel_"] = ibooker.book1D("DecayChannel", "Decay Channel", 3, 0, 3);
0163     // trigger efficiency estimates for the electron muon channel
0164     hists_["elecMuEff_"] = ibooker.book1D("ElecMuEff", "Eff(e/#mu paths)", nElecMu, 0., nElecMu);
0165     // monitored trigger occupancy for the electron muon channel
0166     hists_["elecMuMon_"] = ibooker.book1D("ElecMuMon", "Mon(e/#mu paths)", nElecMu, 0., nElecMu);
0167     // trigger efficiency estimates for the di muon channel
0168     hists_["diMuonEff_"] = ibooker.book1D("DiMuonEff", "Eff(#mu/#mu paths)", nDiMuon, 0., nDiMuon);
0169     // monitored trigger occupancy for the di muon channel
0170     hists_["diMuonMon_"] = ibooker.book1D("DiMuonMon", "Mon(#mu/#mu paths)", nDiMuon, 0., nDiMuon);
0171     // pt of the leading lepton
0172     hists_["lep1Pt_"] = ibooker.book1D("Lep1Pt", "pt(lep1)", 50, 0., 200.);
0173     // pt of the 2. leading lepton
0174     hists_["lep2Pt_"] = ibooker.book1D("Lep2Pt", "pt(lep2)", 50, 0., 200.);
0175     // multiplicity of jets with pt>30 (corrected to L2+L3)
0176     hists_["jetMult_"] = ibooker.book1D("JetMult", "N_{30}(jet)", 21, -0.5, 20.5);
0177     // MET (calo)
0178     hists_["metCalo_"] = ibooker.book1D("METCalo", "MET_{Calo}", 50, 0., 200.);
0179
0180     // set bin labels for trigger monitoring
0181     triggerBinLabels(std::string("elecMu"), elecMuPaths_);
0182     triggerBinLabels(std::string("diMuon"), diMuonPaths_);
0183     // set bin labels for decayChannel_
0184     hists_["decayChannel_"]->setBinLabel(1, "#mu e", 1);
0185     hists_["decayChannel_"]->setBinLabel(2, "#mu #mu", 1);
0186     hists_["decayChannel_"]->setBinLabel(3, "e e", 1);
0187
0188     if (verbosity_ == STANDARD)
0189       return;
0190
0191     // --- [VERBOSE] --- //
0192     // mean eta of the candidate leptons
0193     hists_["sumEtaL1L2_"] = ibooker.book1D("SumEtaL1L2", "<#eta>(lep1, lep2)", 100, -5., 5.);
0194     // deltaEta between the 2 candidate leptons
0195     hists_["dEtaL1L2_"] = ibooker.book1D("DEtaL1L2", "#Delta#eta(lep1,lep2)", 80, -4., 4.);
0196     // deltaPhi between the 2 candidate leptons
0197     hists_["dPhiL1L2_"] = ibooker.book1D("DPhiL1L2", "#Delta#phi(lep1,lep2)", 64, -3.2, 3.2);
0198     // pt of the candidate electron (depending on the decay channel)
0199     hists_["elecPt_"] = ibooker.book1D("ElecPt", "pt(e)", 50, 0., 200.);
0200     // relative isolation of the candidate electron (depending on the decay
0201     // channel)
0202     hists_["elecRelIso_"] = ibooker.book1D("ElecRelIso", "Iso_{Rel}(e)", 50, 0., 1.);
0203     // pt of the canddiate muon (depending on the decay channel)
0204     hists_["muonPt_"] = ibooker.book1D("MuonPt", "pt(#mu)", 50, 0., 200.);
0205     // relative isolation of the candidate muon (depending on the decay channel)
0206     hists_["muonRelIso_"] = ibooker.book1D("MuonRelIso", "Iso_{Rel}(#mu) (#Delta#beta Corrected)", 50, 0., 1.);
0207     // pt of the 1. leading jet (corrected to L2+L3)
0208     hists_["jet1Pt_"] = ibooker.book1D("Jet1Pt", "pt_{L2L3}(jet1)", 60, 0., 300.);
0209     // pt of the 2. leading jet (corrected to L2+L3)
0210     hists_["jet2Pt_"] = ibooker.book1D("Jet2Pt", "pt_{L2L3}(jet2)", 60, 0., 300.);
0211     // MET (PF)
0212     hists_["metPflow_"] = ibooker.book1D("METPflow", "MET_{Pflow}", 50, 0., 200.);
0213     // MET (TC)
0214     hists_["metTC_"] = ibooker.book1D("METTC", "MET_{TC}", 50, 0., 200.);
0215     // dz for muons (to suppress cosmis)
0216     hists_["muonDelZ_"] = ibooker.book1D("MuonDelZ", "d_{z}(#mu)", 50, -25., 25.);
0217     // dxy for muons (to suppress cosmics)
0218     hists_["muonDelXY_"] = ibooker.book2D("MuonDelXY", "d_{xy}(#mu)", 50, -1., 1., 50, -1., 1.);
0219     // lepton multiplicity after std isolation
0220     hists_["lepMultIso_"] = ibooker.book2D("LepMultIso", "N_{Iso}(e) vs N_{Iso}(#mu)", 5, 0., 5., 5, 0., 5.);
0221
0222     // set axes titles for dxy for muons
0223     hists_["muonDelXY_"]->setAxisTitle("x [cm]", 1);
0224     hists_["muonDelXY_"]->setAxisTitle("y [cm]", 2);
0225     // set axes titles for lepton multiplicity after std isolation
0226     hists_["lepMultIso_"]->setAxisTitle("N_{Iso}(#mu)", 1);
0227     hists_["lepMultIso_"]->setAxisTitle("N_{Iso}(elec)", 2);
0228
0229     if (verbosity_ == VERBOSE)
0230       return;
0231
0232     // --- [DEBUG] --- //
0233     // electron multiplicity after std isolation
0234     hists_["elecMultIso_"] = ibooker.book1D("ElecMultIso", "N_{Iso}(e)", 11, -0.5, 10.5);
0235     // muon multiplicity after std isolation
0236     hists_["muonMultIso_"] = ibooker.book1D("MuonMultIso", "N_{Iso}(#mu)", 11, -0.5, 10.5);
0237     // charged hadron isolation component of the candidate muon (depending on the
0238     // decay channel)
0239     hists_["muonChHadIso_"] = ibooker.book1D("MuonChHadIsoComp", "ChHad_{IsoComponent}(#mu)", 50, 0., 5.);
0240     // neutral hadron isolation component of the candidate muon (depending on the
0241     // decay channel)
0242     hists_["muonNeHadIso_"] = ibooker.book1D("MuonNeHadIsoComp", "NeHad_{IsoComponent}(#mu)", 50, 0., 5.);
0243     // photon isolation component of the candidate muon (depending on the decay
0244     // channel)
0245     hists_["muonPhIso_"] = ibooker.book1D("MuonPhIsoComp", "Photon_{IsoComponent}(#mu)", 50, 0., 5.);
0246     // charged hadron isolation component of the candidate electron (depending on
0247     // the decay channel)
0248     hists_["elecChHadIso_"] = ibooker.book1D("ElectronChHadIsoComp", "ChHad_{IsoComponent}(e)", 50, 0., 5.);
0249     // neutral hadron isolation component of the candidate electron (depending on
0250     // the decay channel)
0251     hists_["elecNeHadIso_"] = ibooker.book1D("ElectronNeHadIsoComp", "NeHad_{IsoComponent}(e)", 50, 0., 5.);
0252     // photon isolation component of the candidate electron (depending on the
0253     // decay channel)
0254     hists_["elecPhIso_"] = ibooker.book1D("ElectronPhIsoComp", "Photon_{IsoComponent}(e)", 50, 0., 5.);
0255     // eta of the leading jet
0256     hists_["jet1Eta_"] = ibooker.book1D("Jet1Eta", "#eta(jet1)", 30, -5., 5.);
0257     // eta of the 2. leading jet
0258     hists_["jet2Eta_"] = ibooker.book1D("Jet2Eta", "#eta(jet2)", 30, -5., 5.);
0259     // pt of the 1. leading jet (not corrected)
0260     hists_["jet1PtRaw_"] = ibooker.book1D("Jet1PtRaw", "pt_{Raw}(jet1)", 60, 0., 300.);
0261     // pt of the 2. leading jet (not corrected)
0262     hists_["jet2PtRaw_"] = ibooker.book1D("Jet2PtRaw", "pt_{Raw}(jet2)", 60, 0., 300.);
0263     // deltaEta between the 2 leading jets
0264     hists_["dEtaJet1Jet2_"] = ibooker.book1D("DEtaJet1Jet2", "#Delta#eta(jet1,jet2)", 80, -4., 4.);
0265     // deltaEta between the lepton and the leading jet
0266     hists_["dEtaJet1Lep1_"] = ibooker.book1D("DEtaJet1Lep1", "#Delta#eta(jet1,lep1)", 80, -4., 4.);
0267     // deltaEta between the lepton and MET
0268     hists_["dEtaLep1MET_"] = ibooker.book1D("DEtaLep1MET", "#Delta#eta(lep1,MET)", 80, -4., 4.);
0269     // deltaEta between leading jet and MET
0270     hists_["dEtaJet1MET_"] = ibooker.book1D("DEtaJet1MET", "#Delta#eta(jet1,MET)", 80, -4., 4.);
0271     // deltaPhi of 2 leading jets
0272     hists_["dPhiJet1Jet2_"] = ibooker.book1D("DPhiJet1Jet2", "#Delta#phi(jet1,jet2)", 64, -3.2, 3.2);
0273     // deltaPhi of 1. lepton and 1. jet
0274     hists_["dPhiJet1Lep1_"] = ibooker.book1D("DPhiJet1Lep1", "#Delta#phi(jet1,lep1)", 64, -3.2, 3.2);
0275     // deltaPhi of 1. lepton and MET
0276     hists_["dPhiLep1MET_"] = ibooker.book1D("DPhiLep1MET", "#Delta#phi(lep1,MET)", 64, -3.2, 3.2);
0277     // deltaPhi of 1. jet and MET
0278     hists_["dPhiJet1MET_"] = ibooker.book1D("DPhiJet1MET", "#Delta#phi(jet1,MET)", 64, -3.2, 3.2);
0279     // selected dimuon events
0280     hists_["diMuonLogger_"] = ibooker.book2D("DiMuonLogger", "Logged DiMuon Events", 8, 0., 8., 10, 0., 10.);
0281     // selected dielec events
0282     hists_["diElecLogger_"] = ibooker.book2D("DiElecLogger", "Logged DiElec Events", 8, 0., 8., 10, 0., 10.);
0283     // selected elemu events
0284     hists_["elecMuLogger_"] = ibooker.book2D("ElecMuLogger", "Logged ElecMu Events", 8, 0., 8., 10, 0., 10.);
0285
0286     // set bin labels for trigger monitoring
0287     loggerBinLabels(std::string("diMuonLogger_"));
0288     loggerBinLabels(std::string("diElecLogger_"));
0289     loggerBinLabels(std::string("elecMuLogger_"));
0290     return;
0291   }
0292
0293   void MonitorEnsemble::fill(const edm::Event& event, const edm::EventSetup& setup) {
0294     // fetch trigger event if configured such
0295     edm::Handle<edm::TriggerResults> triggerTable;
0296     if (!triggerTable_.isUninitialized()) {
0297       if (!event.getByToken(triggerTable_, triggerTable))
0298         return;
0299     }
0300     /*
0301   ------------------------------------------------------------
0302
0303   Run and Inst. Luminosity information (Inst. Lumi. filled now with a dummy
0304   value=5.0)
0305
0306   ------------------------------------------------------------
0307   */
0308
0309     if (!event.eventAuxiliary().run())
0310       return;
0311     fill("RunNumb_", event.eventAuxiliary().run());
0312
0313     double dummy = 5.;
0314     fill("InstLumi_", dummy);
0315
0316     /*
0317   ------------------------------------------------------------
0318
0319   Muon Selection
0320
0321   ------------------------------------------------------------
0322   */
0323
0324     std::vector<const reco::PFCandidate*> isoMuons;
0325
0326     edm::Handle<edm::View<reco::PFCandidate>> muons;
0327     edm::View<reco::PFCandidate>::const_iterator muonit;
0328
0329     if (!event.getByToken(muons_, muons))
0330       return;
0331
0332     for (edm::View<reco::PFCandidate>::const_iterator muonit = muons->begin(); muonit != muons->end(); ++muonit) {
0333       if (muonit->muonRef().isNull())
0334         continue;
0335       reco::MuonRef muon = muonit->muonRef();
0336
0337       if (muon->innerTrack().isNull())
0338         continue;
0339
0340       if (muon->isGlobalMuon()) {
0341         fill("muonDelZ_", muon->innerTrack()->vz());  // CB using inner track!
0342         fill("muonDelXY_", muon->innerTrack()->vx(), muon->innerTrack()->vy());
0343
0344         // apply selection
0345         if (!muonSelect_ || (*muonSelect_)(*muonit)) {
0346           double chHadPt = muon->pfIsolationR04().sumChargedHadronPt;
0347           double neHadEt = muon->pfIsolationR04().sumNeutralHadronEt;
0348           double phoEt = muon->pfIsolationR04().sumPhotonEt;
0349
0350           double pfRelIso = (chHadPt + std::max(0., neHadEt + phoEt - 0.5 * muon->pfIsolationR04().sumPUPt)) /
0351                             muon->pt();  // CB dBeta corrected iso!
0352
0353           fill("muonRelIso_", pfRelIso);
0354
0355           fill("muonChHadIso_", chHadPt);
0356           fill("muonNeHadIso_", neHadEt);
0357           fill("muonPhIso_", phoEt);
0358
0359           if (!muonIso_ || (*muonIso_)(*muonit))
0360             isoMuons.push_back(&(*muonit));
0361         }
0362       }
0363     }
0364
0365     fill("muonMultIso_", isoMuons.size());
0366
0367     /*
0368   ------------------------------------------------------------
0369
0370   Electron Selection
0371
0372   ------------------------------------------------------------
0373   */
0374
0375     // buffer isolated electronss
0376     std::vector<const reco::PFCandidate*> isoElecs;
0377     edm::Handle<edm::ValueMap<float>> electronId;
0378     if (!electronId_.isUninitialized()) {
0379       if (!event.getByToken(electronId_, electronId))
0380         return;
0381     }
0382     edm::Handle<edm::View<reco::PFCandidate>> elecs;
0383     if (!event.getByToken(elecs_, elecs))
0384       return;
0385
0386     for (edm::View<reco::PFCandidate>::const_iterator elec = elecs->begin(); elec != elecs->end(); ++elec) {
0387       if (elec->gsfElectronRef().isNull()) {
0388         continue;
0389       }
0390       reco::GsfElectronRef gsf_el = elec->gsfElectronRef();
0391       // restrict to electrons with good electronId
0392       if (electronId_.isUninitialized() ? true : ((double)(*electronId)[gsf_el] >= eidCutValue_)) {
0393         // apply preselection
0394         if (!elecSelect_ || (*elecSelect_)(*elec)) {
0395           double el_ChHadIso = gsf_el->pfIsolationVariables().sumChargedHadronPt;
0396           double el_NeHadIso = gsf_el->pfIsolationVariables().sumNeutralHadronEt;
0397           double el_PhIso = gsf_el->pfIsolationVariables().sumPhotonEt;
0398           double el_pfRelIso =
0399               (el_ChHadIso + std::max(0., el_NeHadIso + el_PhIso - 0.5 * gsf_el->pfIsolationVariables().sumPUPt)) /
0400               gsf_el->pt();
0401           fill("elecRelIso_", el_pfRelIso);
0402           fill("elecChHadIso_", el_ChHadIso);
0403           fill("elecNeHadIso_", el_NeHadIso);
0404           fill("elecPhIso_", el_PhIso);
0405           if (!elecIso_ || (*elecIso_)(*elec))
0406             isoElecs.push_back(&(*elec));
0407         }
0408       }
0409     }
0410     fill("elecMultIso_", isoElecs.size());
0411
0412     /*
0413   ------------------------------------------------------------
0414
0415   Jet Selection
0416
0417   ------------------------------------------------------------
0418   */
0419
0420     const reco::JetCorrector* corrector = nullptr;
0421     if (!jetCorrector_.isUninitialized()) {
0422       // check whether a jet corrector is in the event or not
0423       edm::Handle<reco::JetCorrector> correctorHandle = event.getHandle(jetCorrector_);
0424       if (correctorHandle.isValid()) {
0425         corrector = correctorHandle.product();
0426       } else {
0427         edm::LogVerbatim("TopDiLeptonOfflineDQM") << "\n"
0428                                                   << "-----------------------------------------------------------------"
0429                                                      "-------------------- \n"
0430                                                   << " No JetCorrector available from Event:\n"
0431                                                   << "  - Jets will not be corrected.                                  "
0432                                                      "                     \n"
0433                                                   << "-----------------------------------------------------------------"
0434                                                      "-------------------- \n";
0435       }
0436     }
0437
0438     unsigned int mult = 0;
0439     // buffer leadingJets
0440     std::vector<reco::Jet> leadingJets;
0441     edm::Handle<edm::View<reco::Jet>> jets;
0442     if (!event.getByToken(jets_, jets))
0443       return;
0444
0445     edm::Handle<reco::JetIDValueMap> jetID;
0446     if (jetIDSelect_) {
0447       if (!event.getByToken(jetIDLabel_, jetID))
0448         return;
0449     }
0450
0451     for (edm::View<reco::Jet>::const_iterator jet = jets->begin(); jet != jets->end(); ++jet) {
0452       unsigned int idx = jet - jets->begin();
0453       if (jetIDSelect_ && dynamic_cast<const reco::CaloJet*>(jets->refAt(idx).get())) {
0454         if (!(*jetIDSelect_)((*jetID)[jets->refAt(idx)]))
0455           continue;
0456       }
0457       // chekc additional jet selection for calo, pf and bare reco jets
0458       if (dynamic_cast<const reco::CaloJet*>(&*jet)) {
0459         reco::CaloJet sel = dynamic_cast<const reco::CaloJet&>(*jet);
0460         sel.scaleEnergy(corrector ? corrector->correction(*jet) : 1.);
0461         StringCutObjectSelector<reco::CaloJet> jetSelect(jetSelect_);
0462         if (!jetSelect(sel)) {
0463           continue;
0464         }
0465       } else if (dynamic_cast<const reco::PFJet*>(&*jet)) {
0466         reco::PFJet sel = dynamic_cast<const reco::PFJet&>(*jet);
0467         sel.scaleEnergy(corrector ? corrector->correction(*jet) : 1.);
0468         StringCutObjectSelector<reco::PFJet> jetSelect(jetSelect_);
0469         if (!jetSelect(sel))
0470           continue;
0471       } else {
0472         reco::Jet sel = *jet;
0473         sel.scaleEnergy(corrector ? corrector->correction(*jet) : 1.);
0474         StringCutObjectSelector<reco::Jet> jetSelect(jetSelect_);
0475         if (!jetSelect(sel))
0476           continue;
0477       }
0478       // check for overlaps
0479       bool overlap = false;
0480       for (std::vector<const reco::PFCandidate*>::const_iterator elec = isoElecs.begin(); elec != isoElecs.end();
0481            ++elec) {
0482         if (reco::deltaR((*elec)->eta(), (*elec)->phi(), jet->eta(), jet->phi()) < 0.4) {
0483           overlap = true;
0484           break;
0485         }
0486       }
0487       if (overlap) {
0488         continue;
0489       }
0490       // prepare jet to fill monitor histograms
0491       reco::Jet monitorJet = *jet;
0492       monitorJet.scaleEnergy(corrector ? corrector->correction(*jet) : 1.);
0493       ++mult;  // determine jet multiplicity
0494       if (idx == 0) {
0495         leadingJets.push_back(monitorJet);
0496         fill("jet1Pt_", monitorJet.pt());
0497         fill("jet1PtRaw_", jet->pt());
0498         fill("jet1Eta_", jet->eta());
0499       }
0500       if (idx == 1) {
0501         leadingJets.push_back(monitorJet);
0502         fill("jet2Pt_", monitorJet.pt());
0503         fill("jet2PtRaw_", jet->pt());
0504         fill("jet2Eta_", jet->eta());
0505       }
0506     }
0507     if (leadingJets.size() > 1) {
0508       fill("dEtaJet1Jet2_", leadingJets[0].eta() - leadingJets[1].eta());
0509       fill("dPhiJet1Jet2_", reco::deltaPhi(leadingJets[0].phi(), leadingJets[1].phi()));
0510       if (!isoMuons.empty()) {
0511         if (isoElecs.empty() || isoMuons[0]->pt() > isoElecs[0]->pt()) {
0512           fill("dEtaJet1Lep1_", isoMuons[0]->eta() - leadingJets[0].eta());
0513           fill("dPhiJet1Lep1_", reco::deltaPhi(isoMuons[0]->phi(), leadingJets[0].phi()));
0514         }
0515       }
0516       if (!isoElecs.empty()) {
0517         if (isoMuons.empty() || isoElecs[0]->pt() > isoMuons[0]->pt()) {
0518           fill("dEtaJet1Lep1_", isoElecs[0]->eta() - leadingJets[0].eta());
0519           fill("dPhiJet1Lep1_", reco::deltaPhi(isoElecs[0]->phi(), leadingJets[0].phi()));
0520         }
0521       }
0522     }
0523     fill("jetMult_", mult);
0524
0525     /*
0526   ------------------------------------------------------------
0527
0528   MET Selection
0529
0530   ------------------------------------------------------------
0531   */
0532
0533     // buffer for event logging
0534     reco::MET caloMET;
0535     for (std::vector<edm::EDGetTokenT<edm::View<reco::MET>>>::const_iterator met_ = mets_.begin(); met_ != mets_.end();
0536          ++met_) {
0537       edm::Handle<edm::View<reco::MET>> met;
0538       if (!event.getByToken(*met_, met))
0539         continue;
0540
0541       if (met->begin() != met->end()) {
0542         unsigned int idx = met_ - mets_.begin();
0543         if (idx == 0) {
0544           caloMET = *met->begin();
0545           fill("metCalo_", met->begin()->et());
0546           if (!leadingJets.empty()) {
0547             fill("dEtaJet1MET_", leadingJets[0].eta() - met->begin()->eta());
0548             fill("dPhiJet1MET_", reco::deltaPhi(leadingJets[0].phi(), met->begin()->phi()));
0549           }
0550           if (!isoMuons.empty()) {
0551             if (isoElecs.empty() || isoMuons[0]->pt() > isoElecs[0]->pt()) {
0552               fill("dEtaLep1MET_", isoMuons[0]->eta() - met->begin()->eta());
0553               fill("dPhiLep1MET_", reco::deltaPhi(isoMuons[0]->phi(), met->begin()->phi()));
0554             }
0555           }
0556           if (!isoElecs.empty()) {
0557             if (isoMuons.empty() || isoElecs[0]->pt() > isoMuons[0]->pt()) {
0558               fill("dEtaLep1MET_", isoElecs[0]->eta() - met->begin()->eta());
0559               fill("dPhiLep1MET_", reco::deltaPhi(isoElecs[0]->phi(), met->begin()->phi()));
0560             }
0561           }
0562         }
0563         if (idx == 1) {
0564           fill("metTC_", met->begin()->et());
0565         }
0566         if (idx == 2) {
0567           fill("metPflow_", met->begin()->et());
0568         }
0569       }
0570     }
0571
0572     /*
0573   ------------------------------------------------------------
0574
0575   Event Monitoring
0576
0577   ------------------------------------------------------------
0578   */
0579
0580     // check number of isolated leptons
0581     fill("lepMultIso_", isoMuons.size(), isoElecs.size());
0582     // ELECMU channel
0583     if (decayChannel(isoMuons, isoElecs) == ELECMU) {
0584       fill("decayChannel_", 0.5);
0585       double mass = (isoElecs[0]->p4() + isoMuons[0]->p4()).mass();
0586       if ((lowerEdge_ == -1. && upperEdge_ == -1.) || (lowerEdge_ < mass && mass < upperEdge_)) {
0587         fill("dEtaL1L2_", isoElecs[0]->eta() - isoMuons[0]->eta());
0588         fill("sumEtaL1L2_", (isoElecs[0]->eta() + isoMuons[0]->eta()) / 2);
0589         fill("dPhiL1L2_", reco::deltaPhi(isoElecs[0]->phi(), isoMuons[0]->eta()));
0590         fill("elecPt_", isoElecs[0]->pt());
0591         fill("muonPt_", isoMuons[0]->pt());
0592         fill("lep1Pt_", isoElecs[0]->pt() > isoMuons[0]->pt() ? isoElecs[0]->pt() : isoMuons[0]->pt());
0593         fill("lep2Pt_", isoElecs[0]->pt() > isoMuons[0]->pt() ? isoMuons[0]->pt() : isoElecs[0]->pt());
0594         // fill plots for trigger monitoring
0595         if (!triggerTable_.isUninitialized())
0596           fill(event, *triggerTable, "elecMu", elecMuPaths_);
0597         if (elecMuLogged_ <= hists_.find("elecMuLogger_")->second->getNbinsY()) {
0598           // log runnumber, lumi block, event number & some
0599           // more pysics infomation for interesting events
0600           fill("elecMuLogger_", 0.5, elecMuLogged_ + 0.5, event.eventAuxiliary().run());
0601           fill("elecMuLogger_", 1.5, elecMuLogged_ + 0.5, event.eventAuxiliary().luminosityBlock());
0602           fill("elecMuLogger_", 2.5, elecMuLogged_ + 0.5, event.eventAuxiliary().event());
0603           fill("elecMuLogger_", 3.5, elecMuLogged_ + 0.5, isoMuons[0]->pt());
0604           fill("elecMuLogger_", 4.5, elecMuLogged_ + 0.5, isoElecs[0]->pt());
0605           if (!leadingJets.empty())
0606             fill("elecMuLogger_", 5.5, elecMuLogged_ + 0.5, leadingJets[0].pt());
0607           if (leadingJets.size() > 1)
0608             fill("elecMuLogger_", 6.5, elecMuLogged_ + 0.5, leadingJets[1].pt());
0609           fill("elecMuLogger_", 7.5, elecMuLogged_ + 0.5, caloMET.et());
0610           ++elecMuLogged_;
0611         }
0612       }
0613     }
0614
0615     // DIMUON channel
0616     if (decayChannel(isoMuons, isoElecs) == DIMUON) {
0617       fill("decayChannel_", 1.5);
0618       int charge = isoMuons[0]->charge() * isoMuons[1]->charge();
0619       double mass = (isoMuons[0]->p4() + isoMuons[1]->p4()).mass();
0620
0621       fill(charge < 0 ? "invMass_" : "invMassWC_", mass);
0622       fill(charge < 0 ? "invMassLog_" : "invMassWCLog_", log10(mass));
0623       if ((lowerEdge_ == -1. && upperEdge_ == -1.) || (lowerEdge_ < mass && mass < upperEdge_)) {
0624         fill("dEtaL1L2_", isoMuons[0]->eta() - isoMuons[1]->eta());
0625         fill("sumEtaL1L2_", (isoMuons[0]->eta() + isoMuons[1]->eta()) / 2);
0626         fill("dPhiL1L2_", reco::deltaPhi(isoMuons[0]->phi(), isoMuons[1]->phi()));
0627         fill("muonPt_", isoMuons[0]->pt());
0628         fill("muonPt_", isoMuons[1]->pt());
0629         fill("lep1Pt_", isoMuons[0]->pt());
0630         fill("lep2Pt_", isoMuons[1]->pt());
0631         // fill plots for trigger monitoring
0632         if (!triggerTable_.isUninitialized())
0633           fill(event, *triggerTable, "diMuon", diMuonPaths_);
0634         if (diMuonLogged_ <= hists_.find("diMuonLogger_")->second->getNbinsY()) {
0635           // log runnumber, lumi block, event number & some
0636           // more pysics infomation for interesting events
0637           fill("diMuonLogger_", 0.5, diMuonLogged_ + 0.5, event.eventAuxiliary().run());
0638           fill("diMuonLogger_", 1.5, diMuonLogged_ + 0.5, event.eventAuxiliary().luminosityBlock());
0639           fill("diMuonLogger_", 2.5, diMuonLogged_ + 0.5, event.eventAuxiliary().event());
0640           fill("diMuonLogger_", 3.5, diMuonLogged_ + 0.5, isoMuons[0]->pt());
0641           fill("diMuonLogger_", 4.5, diMuonLogged_ + 0.5, isoMuons[1]->pt());
0642           if (!leadingJets.empty())
0643             fill("diMuonLogger_", 5.5, diMuonLogged_ + 0.5, leadingJets[0].pt());
0644           if (leadingJets.size() > 1)
0645             fill("diMuonLogger_", 6.5, diMuonLogged_ + 0.5, leadingJets[1].pt());
0646           fill("diMuonLogger_", 7.5, diMuonLogged_ + 0.5, caloMET.et());
0647           ++diMuonLogged_;
0648         }
0649       }
0650     }
0651
0652     // DIELEC channel
0653     if (decayChannel(isoMuons, isoElecs) == DIELEC) {
0654       fill("decayChannel_", 2.5);
0655       int charge = isoElecs[0]->charge() * isoElecs[1]->charge();
0656       double mass = (isoElecs[0]->p4() + isoElecs[1]->p4()).mass();
0657       fill(charge < 0 ? "invMass_" : "invMassWC_", mass);
0658       fill(charge < 0 ? "invMassLog_" : "invMassWCLog_", log10(mass));
0659       if ((lowerEdge_ == -1. && upperEdge_ == -1.) || (lowerEdge_ < mass && mass < upperEdge_)) {
0660         fill("dEtaL1L2_", isoElecs[0]->eta() - isoElecs[1]->eta());
0661         fill("sumEtaL1L2_", (isoElecs[0]->eta() + isoElecs[1]->eta()) / 2);
0662         fill("dPhiL1L2_", reco::deltaPhi(isoElecs[0]->phi(), isoElecs[1]->phi()));
0663         fill("elecPt_", isoElecs[0]->pt());
0664         fill("elecPt_", isoElecs[1]->pt());
0665         fill("lep1Pt_", isoElecs[0]->pt());
0666         fill("lep2Pt_", isoElecs[1]->pt());
0667         if (diElecLogged_ <= hists_.find("diElecLogger_")->second->getNbinsY()) {
0668           // log runnumber, lumi block, event number & some
0669           // more pysics infomation for interesting events
0670           fill("diElecLogger_", 0.5, diElecLogged_ + 0.5, event.eventAuxiliary().run());
0671           fill("diElecLogger_", 1.5, diElecLogged_ + 0.5, event.eventAuxiliary().luminosityBlock());
0672           fill("diElecLogger_", 2.5, diElecLogged_ + 0.5, event.eventAuxiliary().event());
0673           fill("diElecLogger_", 3.5, diElecLogged_ + 0.5, isoElecs[0]->pt());
0674           fill("diElecLogger_", 4.5, diElecLogged_ + 0.5, isoElecs[1]->pt());
0675           if (!leadingJets.empty())
0676             fill("diElecLogger_", 5.5, diElecLogged_ + 0.5, leadingJets[0].pt());
0677           if (leadingJets.size() > 1)
0678             fill("diElecLogger_", 6.5, diElecLogged_ + 0.5, leadingJets[1].pt());
0679           fill("diElecLogger_", 7.5, diElecLogged_ + 0.5, caloMET.et());
0680           ++diElecLogged_;
0681         }
0682       }
0683     }
0684   }
0685 }  // namespace TopDiLeptonOffline
0686
0687 TopDiLeptonOfflineDQM::TopDiLeptonOfflineDQM(const edm::ParameterSet& cfg)
0688     : vertexSelect_(nullptr),
0689       beamspotSelect_(nullptr),
0690       MuonStep(nullptr),
0691       ElectronStep(nullptr),
0692       PvStep(nullptr),
0693       METStep(nullptr) {
0694   JetSteps.clear();
0695   CaloJetSteps.clear();
0696   PFJetSteps.clear();
0697   // configure the preselection
0698   edm::ParameterSet presel = cfg.getParameter<edm::ParameterSet>("preselection");
0699   if (presel.existsAs<edm::ParameterSet>("trigger")) {
0700     edm::ParameterSet trigger = presel.getParameter<edm::ParameterSet>("trigger");
0701     //    triggerTable_=trigger.getParameter<edm::InputTag>("src");
0702     triggerTable_ = consumes<edm::TriggerResults>(trigger.getParameter<edm::InputTag>("src"));
0703     triggerPaths_ = trigger.getParameter<std::vector<std::string>>("select");
0704   }
0705   if (presel.existsAs<edm::ParameterSet>("vertex")) {
0706     edm::ParameterSet vertex = presel.getParameter<edm::ParameterSet>("vertex");
0707     vertex_ = consumes<std::vector<reco::Vertex>>(vertex.getParameter<edm::InputTag>("src"));
0708     vertexSelect_ = std::make_unique<StringCutObjectSelector<reco::Vertex>>(vertex.getParameter<std::string>("select"));
0709   }
0710   if (presel.existsAs<edm::ParameterSet>("beamspot")) {
0711     edm::ParameterSet beamspot = presel.getParameter<edm::ParameterSet>("beamspot");
0712     beamspot_ = consumes<reco::BeamSpot>(beamspot.getParameter<edm::InputTag>("src"));
0713     beamspotSelect_ =
0714         std::make_unique<StringCutObjectSelector<reco::BeamSpot>>(beamspot.getParameter<std::string>("select"));
0715   }
0716
0717   // conifgure the selection
0718   sel_ = cfg.getParameter<std::vector<edm::ParameterSet>>("selection");
0719   setup_ = cfg.getParameter<edm::ParameterSet>("setup");
0720   for (unsigned int i = 0; i < sel_.size(); ++i) {
0721     selectionOrder_.push_back(sel_.at(i).getParameter<std::string>("label"));
0722     selection_[selectionStep(selectionOrder_.back())] =
0723         std::make_pair(sel_.at(i),
0724                        std::make_unique<TopDiLeptonOffline::MonitorEnsemble>(
0725                            selectionStep(selectionOrder_.back()).c_str(), setup_, consumesCollector()));
0726   }
0727   for (std::vector<std::string>::const_iterator selIt = selectionOrder_.begin(); selIt != selectionOrder_.end();
0728        ++selIt) {
0729     std::string key = selectionStep(*selIt), type = objectType(*selIt);
0730     if (selection_.find(key) != selection_.end()) {
0731       if (type == "muons") {
0732         MuonStep = std::make_unique<SelectionStep<reco::PFCandidate>>(selection_[key].first, consumesCollector());
0733       }
0734       if (type == "elecs") {
0735         ElectronStep = std::make_unique<SelectionStep<reco::PFCandidate>>(selection_[key].first, consumesCollector());
0736       }
0737       if (type == "pvs") {
0738         PvStep = std::make_unique<SelectionStep<reco::Vertex>>(selection_[key].first, consumesCollector());
0739       }
0740       if (type == "jets") {
0741         JetSteps.push_back(std::make_unique<SelectionStep<reco::Jet>>(selection_[key].first, consumesCollector()));
0742       }
0743       if (type == "jets/pf") {
0744         PFJetSteps.push_back(std::make_unique<SelectionStep<reco::PFJet>>(selection_[key].first, consumesCollector()));
0745       }
0746       if (type == "jets/calo") {
0747         CaloJetSteps.push_back(
0748             std::make_unique<SelectionStep<reco::CaloJet>>(selection_[key].first, consumesCollector()));
0749       }
0750       if (type == "met") {
0751         METStep = std::make_unique<SelectionStep<reco::MET>>(selection_[key].first, consumesCollector());
0752       }
0753     }
0754   }
0755 }
0756
0757 void TopDiLeptonOfflineDQM::bookHistograms(DQMStore::IBooker& ibooker, edm::Run const&, edm::EventSetup const&) {
0758   for (auto selIt = selection_.begin(); selIt != selection_.end(); ++selIt) {
0759     selIt->second.second->book(ibooker);
0760   }
0761 }
0762 void TopDiLeptonOfflineDQM::analyze(const edm::Event& event, const edm::EventSetup& setup) {
0763   if (!triggerTable_.isUninitialized()) {
0764     edm::Handle<edm::TriggerResults> triggerTable;
0765     if (!event.getByToken(triggerTable_, triggerTable))
0766       return;
0767     if (!accept(event, *triggerTable, triggerPaths_))
0768       return;
0769   }
0770   if (!vertex_.isUninitialized()) {
0771     edm::Handle<std::vector<reco::Vertex>> vertex;
0772     if (!event.getByToken(vertex_, vertex))
0773       return;
0774     if (vertex->empty() || !(*vertexSelect_)(vertex->front()))
0775       return;
0776   }
0777   if (!beamspot_.isUninitialized()) {
0778     edm::Handle<reco::BeamSpot> beamspot;
0779     if (!event.getByToken(beamspot_, beamspot))
0780       return;
0781     if (!(*beamspotSelect_)(*beamspot))
0782       return;
0783   }
0784   unsigned int nJetSteps = -1;
0785
0786   unsigned int nPFJetSteps = -1;
0787
0788   unsigned int nCaloJetSteps = -1;
0789   // apply selection steps
0790   for (std::vector<std::string>::const_iterator selIt = selectionOrder_.begin(); selIt != selectionOrder_.end();
0791        ++selIt) {
0792     std::string key = selectionStep(*selIt), type = objectType(*selIt);
0793     if (selection_.find(key) != selection_.end()) {
0794       if (type == "empty") {
0795         selection_[key].second->fill(event, setup);
0796       }
0797       if (type == "muons" && MuonStep != nullptr) {
0798         if (MuonStep->select(event)) {
0799           selection_[key].second->fill(event, setup);
0800         } else
0801           break;
0802       }
0803       if (type == "elecs" && ElectronStep != nullptr) {
0804         if (ElectronStep->select(event, "electron")) {
0805           selection_[key].second->fill(event, setup);
0806         } else
0807           break;
0808       }
0809       if (type == "jets" && !JetSteps.empty()) {
0810         nJetSteps++;
0811         if (JetSteps[nJetSteps] != nullptr) {
0812           if (JetSteps[nJetSteps]->select(event, setup)) {
0813             selection_[key].second->fill(event, setup);
0814           } else
0815             break;
0816         }
0817       }
0818
0819       if (type == "jets/pf" && !PFJetSteps.empty()) {
0820         nPFJetSteps++;
0821         if (PFJetSteps[nPFJetSteps] != nullptr) {
0822           if (PFJetSteps[nPFJetSteps]->select(event, setup)) {
0823             selection_[key].second->fill(event, setup);
0824           } else
0825             break;
0826         }
0827       }
0828
0829       if (type == "jets/calo" && !CaloJetSteps.empty()) {
0830         nCaloJetSteps++;
0831         if (CaloJetSteps[nCaloJetSteps] != nullptr) {
0832           if (CaloJetSteps[nCaloJetSteps]->select(event, setup)) {
0833             selection_[key].second->fill(event, setup);
0834           } else
0835             break;
0836         }
0837       }
0838
0839       if (type == "met" && METStep != nullptr) {
0840         if (METStep->select(event)) {
0841           selection_[key].second->fill(event, setup);
0842         } else
0843           break;
0844       }
0845     }
0846   }
0847 }