Physics/src/TopSingleLeptonDQM.cc

0001 #include "DQM/Physics/src/TopSingleLeptonDQM.h"
0002 #include "DataFormats/BTauReco/interface/JetTag.h"
0003 #include "DataFormats/Common/interface/Handle.h"
0004 #include "DataFormats/JetReco/interface/CaloJet.h"
0005 #include "DataFormats/JetReco/interface/PFJet.h"
0006 #include "DataFormats/Math/interface/deltaR.h"
0007 #include "JetMETCorrections/Objects/interface/JetCorrectionsRecord.h"
0008 #include <iostream>
0009 #include <memory>
0010
0011 #include "FWCore/Framework/interface/ConsumesCollector.h"
0012 #include "FWCore/Framework/interface/EDConsumerBase.h"
0013 #include "FWCore/Utilities/interface/EDGetToken.h"
0014 #include "FWCore/Utilities/interface/InputTag.h"
0015
0016 using namespace std;
0017 namespace TopSingleLepton {
0018
0019   // maximal number of leading jets
0020   // to be used for top mass estimate
0021   static const unsigned int MAXJETS = 4;
0022   // nominal mass of the W boson to
0023   // be used for the top mass estimate
0024   static const double WMASS = 80.4;
0025
0026   MonitorEnsemble::MonitorEnsemble(const char* label, const edm::ParameterSet& cfg, edm::ConsumesCollector&& iC)
0027       : label_(label),
0028         elecSelect_(nullptr),
0029         pvSelect_(nullptr),
0030         muonIso_(nullptr),
0031         muonSelect_(nullptr),
0032         jetIDSelect_(nullptr),
0033         jetlooseSelection_(nullptr),
0034         jetSelection_(nullptr),
0035         includeBTag_(false),
0036         lowerEdge_(-1.),
0037         upperEdge_(-1.),
0038         logged_(0) {
0039     // sources have to be given; this PSet is not optional
0040     edm::ParameterSet sources = cfg.getParameter<edm::ParameterSet>("sources");
0041     muons_ = iC.consumes<edm::View<reco::PFCandidate>>(sources.getParameter<edm::InputTag>("muons"));
0042     elecs_ = iC.consumes<edm::View<reco::PFCandidate>>(sources.getParameter<edm::InputTag>("elecs"));
0043     pvs_ = iC.consumes<edm::View<reco::Vertex>>(sources.getParameter<edm::InputTag>("pvs"));
0044     jets_ = iC.consumes<edm::View<reco::Jet>>(sources.getParameter<edm::InputTag>("jets"));
0045     for (edm::InputTag const& tag : sources.getParameter<std::vector<edm::InputTag>>("mets"))
0046       mets_.push_back(iC.consumes<edm::View<reco::MET>>(tag));
0047     // electronExtras are optional; they may be omitted or
0048     // empty
0049     if (cfg.existsAs<edm::ParameterSet>("elecExtras")) {
0050       // rho for PF isolation with EA corrections
0051       // eventrhoToken_ =
0052       // iC.consumes<double>(edm::InputTag("fixedGridRhoFastjetAll"));
0053
0054       edm::ParameterSet elecExtras = cfg.getParameter<edm::ParameterSet>("elecExtras");
0055       // select is optional; in case it's not found no
0056       // selection will be applied
0057       if (elecExtras.existsAs<std::string>("select")) {
0058         elecSelect_ = std::make_unique<StringCutObjectSelector<reco::PFCandidate>>(
0059             elecExtras.getParameter<std::string>("select"));
0060       }
0061
0062       if (elecExtras.existsAs<std::string>("rho")) {
0063         rhoTag = elecExtras.getParameter<edm::InputTag>("rho");
0064       }
0065       // electronId is optional; in case it's not found the
0066       // InputTag will remain empty
0067       if (elecExtras.existsAs<edm::ParameterSet>("electronId")) {
0068         edm::ParameterSet elecId = elecExtras.getParameter<edm::ParameterSet>("electronId");
0069         electronId_ = iC.consumes<edm::ValueMap<float>>(elecId.getParameter<edm::InputTag>("src"));
0070         eidCutValue_ = elecId.getParameter<double>("cutValue");
0071       }
0072     }
0073     // pvExtras are optional; they may be omitted or empty
0074     if (cfg.existsAs<edm::ParameterSet>("pvExtras")) {
0075       edm::ParameterSet pvExtras = cfg.getParameter<edm::ParameterSet>("pvExtras");
0076       // select is optional; in case it's not found no
0077       // selection will be applied
0078       if (pvExtras.existsAs<std::string>("select")) {
0079         pvSelect_ =
0080             std::make_unique<StringCutObjectSelector<reco::Vertex>>(pvExtras.getParameter<std::string>("select"));
0081       }
0082     }
0083     // muonExtras are optional; they may be omitted or empty
0084     if (cfg.existsAs<edm::ParameterSet>("muonExtras")) {
0085       edm::ParameterSet muonExtras = cfg.getParameter<edm::ParameterSet>("muonExtras");
0086       // select is optional; in case it's not found no
0087       // selection will be applied
0088       if (muonExtras.existsAs<std::string>("select")) {
0089         muonSelect_ = std::make_unique<StringCutObjectSelector<reco::PFCandidate>>(
0090             muonExtras.getParameter<std::string>("select"));
0091       }
0092       // isolation is optional; in case it's not found no
0093       // isolation will be applied
0094       if (muonExtras.existsAs<std::string>("isolation")) {
0095         muonIso_ = std::make_unique<StringCutObjectSelector<reco::PFCandidate>>(
0096             muonExtras.getParameter<std::string>("isolation"));
0097       }
0098     }
0099
0100     // jetExtras are optional; they may be omitted or
0101     // empty
0102     if (cfg.existsAs<edm::ParameterSet>("jetExtras")) {
0103       edm::ParameterSet jetExtras = cfg.getParameter<edm::ParameterSet>("jetExtras");
0104       // jetCorrector is optional; in case it's not found
0105       // the InputTag will remain empty
0106       if (jetExtras.existsAs<std::string>("jetCorrector")) {
0107         jetCorrector_ =
0108             iC.consumes<reco::JetCorrector>(edm::InputTag(jetExtras.getParameter<std::string>("jetCorrector")));
0109       }
0110       // read jetID information if it exists
0111       if (jetExtras.existsAs<edm::ParameterSet>("jetID")) {
0112         edm::ParameterSet jetID = jetExtras.getParameter<edm::ParameterSet>("jetID");
0113         jetIDLabel_ = iC.consumes<reco::JetIDValueMap>(jetID.getParameter<edm::InputTag>("label"));
0114         jetIDSelect_ =
0115             std::make_unique<StringCutObjectSelector<reco::JetID>>(jetID.getParameter<std::string>("select"));
0116       }
0117       // select is optional; in case it's not found no
0118       // selection will be applied (only implemented for
0119       // CaloJets at the moment)
0120       if (jetExtras.existsAs<std::string>("select")) {
0121         jetSelect_ = jetExtras.getParameter<std::string>("select");
0122         jetSelection_ = std::make_unique<StringCutObjectSelector<reco::PFJet>>(jetSelect_);
0123         jetlooseSelection_ = std::make_unique<StringCutObjectSelector<reco::PFJet>>(
0124             "chargedHadronEnergyFraction()>0 && chargedMultiplicity()>0 && chargedEmEnergyFraction()<0.99 && "
0125             "neutralHadronEnergyFraction()<0.99 && neutralEmEnergyFraction()<0.99 && "
0126             "(chargedMultiplicity()+neutralMultiplicity())>1");
0127       }
0128       // jetBDiscriminators are optional; in case they are
0129       // not found the InputTag will remain empty; they
0130       // consist of pairs of edm::JetFlavorAssociation's &
0131       // corresponding working points
0132       includeBTag_ = jetExtras.existsAs<edm::ParameterSet>("jetBTaggers");
0133       if (includeBTag_) {
0134         edm::ParameterSet btagCSV =
0135             jetExtras.getParameter<edm::ParameterSet>("jetBTaggers").getParameter<edm::ParameterSet>("cvsVertex");
0136         btagCSV_ = iC.consumes<reco::JetTagCollection>(btagCSV.getParameter<edm::InputTag>("label"));
0137         btagCSVWP_ = btagCSV.getParameter<double>("workingPoint");
0138       }
0139     }
0140
0141     // triggerExtras are optional; they may be omitted or empty
0142     if (cfg.existsAs<edm::ParameterSet>("triggerExtras")) {
0143       edm::ParameterSet triggerExtras = cfg.getParameter<edm::ParameterSet>("triggerExtras");
0144       triggerTable_ = iC.consumes<edm::TriggerResults>(triggerExtras.getParameter<edm::InputTag>("src"));
0145       triggerPaths_ = triggerExtras.getParameter<std::vector<std::string>>("paths");
0146     }
0147
0148     // massExtras is optional; in case it's not found no mass
0149     // window cuts are applied for the same flavor monitor
0150     // histograms
0151     if (cfg.existsAs<edm::ParameterSet>("massExtras")) {
0152       edm::ParameterSet massExtras = cfg.getParameter<edm::ParameterSet>("massExtras");
0153       lowerEdge_ = massExtras.getParameter<double>("lowerEdge");
0154       upperEdge_ = massExtras.getParameter<double>("upperEdge");
0155     }
0156
0157     // setup the verbosity level for booking histograms;
0158     // per default the verbosity level will be set to
0159     // STANDARD. This will also be the chosen level in
0160     // the case when the monitoring PSet is not found
0161     verbosity_ = STANDARD;
0162     if (cfg.existsAs<edm::ParameterSet>("monitoring")) {
0163       edm::ParameterSet monitoring = cfg.getParameter<edm::ParameterSet>("monitoring");
0164       if (monitoring.getParameter<std::string>("verbosity") == "DEBUG")
0165         verbosity_ = DEBUG;
0166       if (monitoring.getParameter<std::string>("verbosity") == "VERBOSE")
0167         verbosity_ = VERBOSE;
0168       if (monitoring.getParameter<std::string>("verbosity") == "STANDARD")
0169         verbosity_ = STANDARD;
0170     }
0171     // and don't forget to do the histogram booking
0172     directory_ = cfg.getParameter<std::string>("directory");
0173   }
0174
0175   void MonitorEnsemble::book(DQMStore::IBooker& ibooker) {
0176     // set up the current directory path
0177     std::string current(directory_);
0178     current += label_;
0179     ibooker.setCurrentFolder(current);
0180
0181     // determine number of bins for trigger monitoring
0182     // unsigned int nPaths = triggerPaths_.size();
0183
0184     // --- [STANDARD] --- //
0185     // Run Number
0186     // hists_["RunNumb_"] = ibooker.book1D("RunNumber", "Run Nr.", 10, 0, 10);
0187     // instantaneous luminosity
0188     // hists_["InstLumi_"] = ibooker.book1D("InstLumi", "Inst. Lumi.", 100, 0., 1.e3);
0189     // number of selected primary vertices
0190     hists_["pvMult_"] = ibooker.book1D("PvMult", "N_{good pvs}", 50, 0., 50.);
0191     // pt of the leading muon
0192     hists_["muonPt_"] = ibooker.book1D("MuonPt", "pt(#mu TightId, TightIso)", 40, 0., 200.);
0193     // muon multiplicity after  tight Id
0194     hists_["muonMult_"] = ibooker.book1D("MuonMult", "N_{loose}(#mu)", 10, 0., 10.);
0195     // muon multiplicity after  isolation
0196     //hists_["muonMultIso_"] = ibooker.book1D("MuonMultIso",
0197     //    "N_{TightIso}(#mu)", 10, 0., 10.);
0198     // muon multiplicity after  tight Id and isolation
0199     hists_["muonMultTight_"] = ibooker.book1D("MuonMultTight", "N_{TightIso,TightId}(#mu)", 10, 0., 10.);
0200     // pt of the leading electron
0201     hists_["elecPt_"] = ibooker.book1D("ElecPt", "pt(e TightId, TightIso)", 40, 0., 200.);
0202     // electron multiplicity before std isolation
0203     //hists_["elecMult_"] = ibooker.book1D("ElecMult", "N_{All}(e, tightId)", 10, 0., 10.);
0204     // electron multiplicity after  std isolation
0205     //hists_["elecMultIso_"] = ibooker.book1D("ElecMultIso", "N_{Iso}(e, tightId, tightIso)", 10, 0., 10.);
0206     // multiplicity of jets with pt>30
0207     hists_["jetMult_"] = ibooker.book1D("JetMult", "N_{30}(jet)", 10, 0., 10.);
0208     // multiplicity of loose jets with pt>30
0209     hists_["jetMultLoose_"] = ibooker.book1D("JetMultLoose", "N_{30, loose}(jet)", 10, 0., 10.);
0210
0211     // trigger efficiency estimates for single lepton triggers
0212     // hists_["triggerEff_"] = ibooker.book1D("TriggerEff",
0213     //    "Eff(trigger)", nPaths, 0., nPaths);
0214     // monitored trigger occupancy for single lepton triggers
0215     // hists_["triggerMon_"] = ibooker.book1D("TriggerMon",
0216     //    "Mon(trigger)", nPaths, 0., nPaths);
0217     // MET (pflow)
0218     hists_["metPflow_"] = ibooker.book1D("METPflow", "MET_{Pflow}", 50, 0., 200.);
0219     // W mass estimate
0220     hists_["massW_"] = ibooker.book1D("MassW", "M(W)", 60, 0., 300.);
0221     // Top mass estimate
0222     hists_["massTop_"] = ibooker.book1D("MassTop", "M(Top)", 50, 0., 500.);
0223     // b-tagged Top mass
0224     hists_["massBTop_"] = ibooker.book1D("MassBTop", "M(Top, 1 b-tag)", 50, 0., 500.);
0225     // set bin labels for trigger monitoring
0226     triggerBinLabels(std::string("trigger"), triggerPaths_);
0227
0228     if (verbosity_ == STANDARD)
0229       return;
0230
0231     // --- [VERBOSE] --- //
0232     // eta of the leading muon
0233     hists_["muonEta_"] = ibooker.book1D("MuonEta", "#eta(#mu TightId, TightIso)", 30, -3., 3.);
0234     // relative isolation of the candidate muon (depending on the decay channel)
0235     hists_["muonRelIso_"] = ibooker.book1D("MuonRelIso", "Iso_{Rel}(#mu TightId) (#Delta#beta Corrected)", 50, 0., 1.);
0236     // phi of the leading muon
0237     hists_["muonPhi_"] = ibooker.book1D("MuonPhi", "#phi(#mu TightId, TightIso)", 40, -4., 4.);
0238     // eta of the leading electron
0239     hists_["elecEta_"] = ibooker.book1D("ElecEta", "#eta(e tightId)", 30, -3., 3.);
0240     // std isolation variable of the leading electron
0241     hists_["elecRelIso_"] = ibooker.book1D("ElecRelIso", "Iso_{Rel}(e TightId, TightIso)", 50, 0., 1.);
0242     // phi of the leading electron
0243     hists_["elecPhi_"] = ibooker.book1D("ElecPhi", "#phi(e tightId)", 40, -4., 4.);
0244     // multiplicity of tight Id, tight Iso electorns
0245     hists_["elecMultTight_"] = ibooker.book1D("ElecMultTight", "N_{TightIso,TightId}(e)", 10, 0., 10.);
0246     // multiplicity of btagged jets (for track counting high efficiency) with
0247     // pt(L2L3)>30
0248     // hists_["jetMultBEff_"] = ibooker.book1D("JetMultBEff",
0249     //    "N_{30}(TCHE)", 10, 0., 10.);
0250     // btag discriminator for track counting high efficiency for jets with
0251     // pt(L2L3)>30
0252     //hists_["jetBDiscEff_"] = ibooker.book1D("JetBDiscEff",
0253     //    "Disc_{TCHE}(jet)", 100, 0., 10.);
0254     // eta of the 1. leading jet (corrected to L2+L3)
0255     hists_["jet1Eta_"] = ibooker.book1D("Jet1Eta", "#eta_{30,loose}(jet1)", 60, -3., 3.);
0256     // pt of the 1. leading jet (corrected to L2+L3)
0257     hists_["jet1Pt_"] = ibooker.book1D("Jet1Pt", "pt_{30,loose}(jet1)", 60, 0., 300.);
0258     // eta of the 2. leading jet (corrected to L2+L3)
0259     hists_["jet2Eta_"] = ibooker.book1D("Jet2Eta", "#eta_{30,loose}(jet2)", 60, -3., 3.);
0260     // pt of the 2. leading jet (corrected to L2+L3)
0261     hists_["jet2Pt_"] = ibooker.book1D("Jet2Pt", "pt_{30,loose}(jet2)", 60, 0., 300.);
0262     // eta of the 3. leading jet (corrected to L2+L3)
0263     hists_["jet3Eta_"] = ibooker.book1D("Jet3Eta", "#eta_{30,loose}(jet3)", 60, -3., 3.);
0264     // pt of the 3. leading jet (corrected to L2+L3)
0265     hists_["jet3Pt_"] = ibooker.book1D("Jet3Pt", "pt_{30,loose}(jet3)", 60, 0., 300.);
0266     // eta of the 4. leading jet (corrected to L2+L3)
0267     hists_["jet4Eta_"] = ibooker.book1D("Jet4Eta", "#eta_{30,loose}(jet4)", 60, -3., 3.);
0268     // pt of the 4. leading jet (corrected to L2+L3)
0269     hists_["jet4Pt_"] = ibooker.book1D("Jet4Pt", "pt_{30,loose}(jet4)", 60, 0., 300.);
0270     // dz for muons (to suppress cosmis)
0271     hists_["muonDelZ_"] = ibooker.book1D("MuonDelZ", "d_{z}(#mu)", 50, -25., 25.);
0272     // dxy for muons (to suppress cosmics)
0273     hists_["muonDelXY_"] = ibooker.book2D("MuonDelXY", "d_{xy}(#mu)", 50, -0.1, 0.1, 50, -0.1, 0.1);
0274
0275     // set axes titles for dxy for muons
0276     hists_["muonDelXY_"]->setAxisTitle("x [cm]", 1);
0277     hists_["muonDelXY_"]->setAxisTitle("y [cm]", 2);
0278
0279     if (verbosity_ == VERBOSE)
0280       return;
0281
0282     // --- [DEBUG] --- //
0283     // charged hadron isolation component of the candidate muon (depending on the
0284     // decay channel)
0285     hists_["muonChHadIso_"] = ibooker.book1D("MuonChHadIsoComp", "ChHad_{IsoComponent}(#mu TightId)", 50, 0., 5.);
0286     // neutral hadron isolation component of the candidate muon (depending on the
0287     // decay channel)
0288     hists_["muonNeHadIso_"] = ibooker.book1D("MuonNeHadIsoComp", "NeHad_{IsoComponent}(#mu TightId)", 50, 0., 5.);
0289     // photon isolation component of the candidate muon (depending on the decay
0290     // channel)
0291     hists_["muonPhIso_"] = ibooker.book1D("MuonPhIsoComp", "Photon_{IsoComponent}(#mu )", 50, 0., 5.);
0292     // charged hadron isolation component of the candidate electron (depending on
0293     // the decay channel)
0294     hists_["elecChHadIso_"] = ibooker.book1D("ElectronChHadIsoComp", "ChHad_{IsoComponent}(e tightId)", 50, 0., 5.);
0295     // neutral hadron isolation component of the candidate electron (depending on
0296     // the decay channel)
0297     hists_["elecNeHadIso_"] = ibooker.book1D("ElectronNeHadIsoComp", "NeHad_{IsoComponent}(e tightId)", 50, 0., 5.);
0298     // photon isolation component of the candidate electron (depending on the
0299     // decay channel)
0300     hists_["elecPhIso_"] = ibooker.book1D("ElectronPhIsoComp", "Photon_{IsoComponent}(e tightId)", 50, 0., 5.);
0301
0302     // multiplicity for combined secondary vertex
0303     hists_["jetMultBCSVM_"] = ibooker.book1D("JetMultBCSVM", "N_{30}(CSVM)", 10, 0., 10.);
0304     // btag discriminator for combined secondary vertex
0305     hists_["jetBCSV_"] = ibooker.book1D("JetDiscCSV", "BJet Disc_{CSV}(JET)", 100, -1., 2.);
0306     // pt of the 1. leading jet (uncorrected)
0307     // hists_["jet1PtRaw_"] = ibooker.book1D("Jet1PtRaw", "pt_{Raw}(jet1)", 60, 0., 300.);
0308     // pt of the 2. leading jet (uncorrected)
0309     // hists_["jet2PtRaw_"] = ibooker.book1D("Jet2PtRaw", "pt_{Raw}(jet2)", 60, 0., 300.);
0310     // pt of the 3. leading jet (uncorrected)
0311     // hists_["jet3PtRaw_"] = ibooker.book1D("Jet3PtRaw", "pt_{Raw}(jet3)", 60, 0., 300.);
0312     // pt of the 4. leading jet (uncorrected)
0313     // hists_["jet4PtRaw_"] = ibooker.book1D("Jet4PtRaw", "pt_{Raw}(jet4)", 60, 0., 300.);
0314     // selected events
0315     hists_["eventLogger_"] = ibooker.book2D("EventLogger", "Logged Events", 9, 0., 9., 10, 0., 10.);
0316
0317     // set axes titles for selected events
0318     hists_["eventLogger_"]->getTH1()->SetOption("TEXT");
0319     hists_["eventLogger_"]->setBinLabel(1, "Run", 1);
0320     hists_["eventLogger_"]->setBinLabel(2, "Block", 1);
0321     hists_["eventLogger_"]->setBinLabel(3, "Event", 1);
0322     hists_["eventLogger_"]->setBinLabel(4, "pt_{30,loose}(jet1)", 1);
0323     hists_["eventLogger_"]->setBinLabel(5, "pt_{30,loose}(jet2)", 1);
0324     hists_["eventLogger_"]->setBinLabel(6, "pt_{30,loose}(jet3)", 1);
0325     hists_["eventLogger_"]->setBinLabel(7, "pt_{30,loose}(jet4)", 1);
0326     hists_["eventLogger_"]->setBinLabel(8, "M_{W}", 1);
0327     hists_["eventLogger_"]->setBinLabel(9, "M_{Top}", 1);
0328     hists_["eventLogger_"]->setAxisTitle("logged evts", 2);
0329     return;
0330   }
0331
0332   void MonitorEnsemble::fill(const edm::Event& event, const edm::EventSetup& setup) {
0333     // fetch trigger event if configured such
0334     edm::Handle<edm::TriggerResults> triggerTable;
0335
0336     if (!triggerTable_.isUninitialized()) {
0337       if (!event.getByToken(triggerTable_, triggerTable))
0338         return;
0339     }
0340
0341     /*
0342   ------------------------------------------------------------
0343
0344   Primary Vertex Monitoring
0345
0346   ------------------------------------------------------------
0347   */
0348     // fill monitoring plots for primary verices
0349     edm::Handle<edm::View<reco::Vertex>> pvs;
0350
0351     if (!event.getByToken(pvs_, pvs))
0352       return;
0353     const reco::Vertex& Pvertex = pvs->front();
0354     unsigned int pvMult = 0;
0355     for (edm::View<reco::Vertex>::const_iterator pv = pvs->begin(); pv != pvs->end(); ++pv) {
0356       if (!pvSelect_ || (*pvSelect_)(*pv))
0357         pvMult++;
0358     }
0359     fill("pvMult_", pvMult);
0360
0361     /*
0362   ------------------------------------------------------------
0363
0364   Electron Monitoring
0365
0366   ------------------------------------------------------------
0367   */
0368
0369     // fill monitoring plots for electrons
0370     edm::Handle<edm::View<reco::PFCandidate>> elecs;
0371     edm::Handle<double> _rhoHandle;
0372     event.getByLabel(rhoTag, _rhoHandle);
0373     if (!event.getByToken(elecs_, elecs))
0374       return;
0375
0376     // check availability of electron id
0377     edm::Handle<edm::ValueMap<float>> electronId;
0378     if (!electronId_.isUninitialized()) {
0379       if (!event.getByToken(electronId_, electronId))
0380         return;
0381     }
0382
0383     // loop electron collection
0384     unsigned int eMult = 0, eMultIso = 0;
0385     std::vector<const reco::PFCandidate*> isoElecs;
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()
0393               ? true
0394               : ((double)(*electronId)[gsf_el] >=
0395                  eidCutValue_)) {  //This Electron Id is not currently used, but we can keep this for future needs
0396         if (!elecSelect_ || (*elecSelect_)(*elec)) {
0397           double el_ChHadIso = gsf_el->pfIsolationVariables().sumChargedHadronPt;
0398           double el_NeHadIso = gsf_el->pfIsolationVariables().sumNeutralHadronEt;
0399           double el_PhIso = gsf_el->pfIsolationVariables().sumPhotonEt;
0400           double absEta = std::fabs(gsf_el->superCluster()->eta());
0401
0402           //Effective Area computation
0403           double eA = 0;
0404           if (absEta < 1.000)
0405             eA = 0.1703;
0406           else if (absEta < 1.479)
0407             eA = 0.1715;
0408           else if (absEta < 2.000)
0409             eA = 0.1213;
0410           else if (absEta < 2.200)
0411             eA = 0.1230;
0412           else if (absEta < 2.300)
0413             eA = 0.1635;
0414           else if (absEta < 2.400)
0415             eA = 0.1937;
0416           else if (absEta < 5.000)
0417             eA = 0.2393;
0418
0419           double rho = _rhoHandle.isValid() ? (float)(*_rhoHandle) : 0;
0420           double el_pfRelIso = (el_ChHadIso + max(0., el_NeHadIso + el_PhIso - rho * eA)) / gsf_el->pt();
0421
0422           //Only TightId
0423           if (eMult == 0) {  // Restricted to the leading tight electron
0424             fill("elecRelIso_", el_pfRelIso);
0425             fill("elecChHadIso_", el_ChHadIso);
0426             fill("elecNeHadIso_", el_NeHadIso);
0427             fill("elecPhIso_", el_PhIso);
0428           }
0429           ++eMult;
0430
0431           if (!((el_pfRelIso < 0.0588 && absEta < 1.479) || (el_pfRelIso < 0.0571 && absEta > 1.479)))
0432             continue;  // PF Isolation with Effective Area Corrections according to https://twiki.cern.ch/twiki/bin/viewauth/CMS/CutBasedElectronIdentificationRun2
0433
0434           // TightId and TightIso
0435           if (eMultIso == 0) {  //Only leading
0436             fill("elecPt_", elec->pt());
0437             fill("elecEta_", elec->eta());
0438             fill("elecPhi_", elec->phi());
0439           }
0440           ++eMultIso;
0441         }
0442       }
0443     }
0444     //fill("elecMult_", eMult);
0445     fill("elecMultTight_", eMultIso);
0446
0447     /*
0448   ------------------------------------------------------------
0449
0450   Muon Monitoring
0451
0452   ------------------------------------------------------------
0453   */
0454
0455     // fill monitoring plots for muons
0456     unsigned int mMult = 0, mTight = 0, mTightId = 0;
0457
0458     edm::Handle<edm::View<reco::PFCandidate>> muons;
0459     edm::View<reco::PFCandidate>::const_iterator muonit;
0460
0461     if (!event.getByToken(muons_, muons))
0462       return;
0463
0464     for (edm::View<reco::PFCandidate>::const_iterator muonit = muons->begin(); muonit != muons->end(); ++muonit) {
0465       if (muonit->muonRef().isNull())
0466         continue;
0467       reco::MuonRef muon = muonit->muonRef();
0468
0469       // restrict to globalMuons
0470       if (muon->isGlobalMuon()) {
0471         fill("muonDelZ_", muon->innerTrack()->vz());  // CB using inner track!
0472         fill("muonDelXY_", muon->innerTrack()->vx(), muon->innerTrack()->vy());
0473
0474         // apply preselection
0475         if ((!muonSelect_ || (*muonSelect_)(*muonit))) {
0476           mMult++;
0477           double chHadPt = muon->pfIsolationR04().sumChargedHadronPt;
0478           double neHadEt = muon->pfIsolationR04().sumNeutralHadronEt;
0479           double phoEt = muon->pfIsolationR04().sumPhotonEt;
0480           double pfRelIso = (chHadPt + max(0., neHadEt + phoEt - 0.5 * muon->pfIsolationR04().sumPUPt)) /
0481                             muon->pt();  // CB dBeta corrected iso!
0482
0483           if (!(muon->isGlobalMuon() && muon->isPFMuon() && muon->globalTrack()->normalizedChi2() < 10. &&
0484                 muon->globalTrack()->hitPattern().numberOfValidMuonHits() > 0 && muon->numberOfMatchedStations() > 1 &&
0485                 muon->innerTrack()->hitPattern().numberOfValidPixelHits() > 0 &&
0486                 muon->innerTrack()->hitPattern().trackerLayersWithMeasurement() > 5 &&
0487                 fabs(muon->muonBestTrack()->dxy(Pvertex.position())) < 0.2 &&
0488                 fabs(muon->muonBestTrack()->dz(Pvertex.position())) < 0.5))
0489             continue;  //Only tight muons
0490
0491           if (mTightId == 0) {
0492             fill("muonRelIso_", pfRelIso);
0493             fill("muonChHadIso_", chHadPt);
0494             fill("muonNeHadIso_", neHadEt);
0495             fill("muonPhIso_", phoEt);
0496           }
0497           mTightId++;
0498
0499           if (!(pfRelIso < 0.15))
0500             continue;  //Tight Iso
0501
0502           if (mTight == 0) {  //Leading tightId tightIso muon
0503             fill("muonPt_", muon->pt());
0504             fill("muonEta_", muon->eta());
0505             fill("muonPhi_", muon->phi());
0506           }
0507           mTight++;
0508         }
0509       }
0510     }
0511     fill("muonMult_", mMult);
0512     fill("muonMultTight_", mTight);
0513
0514     /*
0515   ------------------------------------------------------------
0516
0517   Jet Monitoring
0518
0519   ------------------------------------------------------------
0520   */
0521
0522     const reco::JetCorrector* corrector = nullptr;
0523     if (!jetCorrector_.isUninitialized()) {
0524       // check whether a jet corrector is in the event or not
0525       edm::Handle<reco::JetCorrector> correctorHandle = event.getHandle(jetCorrector_);
0526       if (correctorHandle.isValid()) {
0527         corrector = correctorHandle.product();
0528       } else {
0529         edm::LogVerbatim("TopDiLeptonOfflineDQM") << "\n"
0530                                                   << "-----------------------------------------------------------------"
0531                                                      "-------------------- \n"
0532                                                   << " No JetCorrector available from Event:              "
0533                                                      "                     \n"
0534                                                   << "  - Jets will not be corrected.                                  "
0535                                                      "                     \n"
0536                                                   << "-----------------------------------------------------------------"
0537                                                      "-------------------- \n";
0538       }
0539     }
0540
0541     // check availability of the btaggers
0542     edm::Handle<reco::JetTagCollection> btagEff, btagPur, btagVtx, btagCSV;
0543     if (includeBTag_) {
0544       if (!event.getByToken(btagCSV_, btagCSV))
0545         return;
0546     }
0547
0548     // loop jet collection
0549     std::vector<reco::Jet> correctedJets;
0550     std::vector<double> JetTagValues;
0551     unsigned int mult = 0, multLoose = 0, multCSV = 0;
0552
0553     edm::Handle<edm::View<reco::Jet>> jets;
0554     if (!event.getByToken(jets_, jets)) {
0555       return;
0556     }
0557
0558     for (edm::View<reco::Jet>::const_iterator jet = jets->begin(); jet != jets->end(); ++jet) {
0559       bool isLoose = false;
0560       // check additional jet selection for pf jets
0561       if (dynamic_cast<const reco::PFJet*>(&*jet)) {
0562         reco::PFJet sel = dynamic_cast<const reco::PFJet&>(*jet);
0563         if ((*jetlooseSelection_)(sel))
0564           isLoose = true;
0565         sel.scaleEnergy(corrector ? corrector->correction(*jet) : 1.);
0566         if (!(*jetSelection_)(sel))
0567           continue;
0568       }
0569
0570       // prepare jet to fill monitor histograms
0571       reco::Jet monitorJet = *jet;
0572
0573       ++mult;  // determine jet (no Id) multiplicity
0574       monitorJet.scaleEnergy(corrector ? corrector->correction(*jet) : 1.);
0575
0576       if (isLoose) {  //Loose Id
0577         unsigned int idx = jet - jets->begin();
0578         correctedJets.push_back(monitorJet);
0579         if (includeBTag_) {
0580           // fill b-discriminators
0581           edm::RefToBase<reco::Jet> jetRef = jets->refAt(idx);
0582           fill("jetBCSV_", (*btagCSV)[jetRef]);
0583           if ((*btagCSV)[jetRef] > btagCSVWP_)
0584             ++multCSV;
0585           // Fill a vector with Jet b-tag WP for later M3+1tag calculation: CSV
0586           // tagger
0587           JetTagValues.push_back((*btagCSV)[jetRef]);
0588         }
0589
0590         // fill pt/eta for the leading four jets
0591         if (multLoose == 0) {
0592           fill("jet1Pt_", monitorJet.pt());
0593           fill("jet1Eta_", monitorJet.eta());
0594         };
0595         if (multLoose == 1) {
0596           fill("jet2Pt_", monitorJet.pt());
0597           fill("jet2Eta_", monitorJet.eta());
0598         }
0599         if (multLoose == 2) {
0600           fill("jet3Pt_", monitorJet.pt());
0601           fill("jet3Eta_", monitorJet.eta());
0602         }
0603         if (multLoose == 3) {
0604           fill("jet4Pt_", monitorJet.pt());
0605           fill("jet4Eta_", monitorJet.eta());
0606         }
0607         multLoose++;
0608       }
0609     }
0610     fill("jetMult_", mult);
0611     fill("jetMultLoose_", multLoose);
0612     fill("jetMultBCSVM_", multCSV);
0613
0614     /*
0615   ------------------------------------------------------------
0616
0617   MET Monitoring
0618
0619   ------------------------------------------------------------
0620   */
0621
0622     // fill monitoring histograms for met
0623     for (std::vector<edm::EDGetTokenT<edm::View<reco::MET>>>::const_iterator met_ = mets_.begin(); met_ != mets_.end();
0624          ++met_) {
0625       edm::Handle<edm::View<reco::MET>> met;
0626       if (!event.getByToken(*met_, met))
0627         continue;
0628       if (met->begin() != met->end()) {  //If we ever have to use more than one type of met again
0629         unsigned int idx = met_ - mets_.begin();
0630         if (idx == 0)
0631           fill("metPflow_", met->begin()->et());
0632       }
0633     }
0634
0635     /*
0636   ------------------------------------------------------------
0637
0638   Event Monitoring
0639
0640   ------------------------------------------------------------
0641   */
0642
0643     // fill W boson and top mass estimates
0644
0645     Calculate eventKinematics(MAXJETS, WMASS);
0646     double wMass = eventKinematics.massWBoson(correctedJets);
0647     double topMass = eventKinematics.massTopQuark(correctedJets);
0648     if (wMass >= 0 && topMass >= 0) {
0649       fill("massW_", wMass);
0650       fill("massTop_", topMass);
0651     }
0652
0653     // Fill M3 with Btag (CSV Tight) requirement
0654
0655     if (!includeBTag_)
0656       return;
0657     if (correctedJets.size() != JetTagValues.size())
0658       return;
0659     double btopMass = eventKinematics.massBTopQuark(correctedJets, JetTagValues, btagCSVWP_);
0660     if (btopMass >= 0)
0661       fill("massBTop_", btopMass);
0662
0663     // fill plots for trigger monitoring
0664     if ((lowerEdge_ == -1. && upperEdge_ == -1.) || (lowerEdge_ < wMass && wMass < upperEdge_)) {
0665       if (!triggerTable_.isUninitialized())
0666         fill(event, *triggerTable, "trigger", triggerPaths_);
0667       if (logged_ <= hists_.find("eventLogger_")->second->getNbinsY()) {
0668         // log runnumber, lumi block, event number & some
0669         // more pysics infomation for interesting events
0670         fill("eventLogger_", 0.5, logged_ + 0.5, event.eventAuxiliary().run());
0671         fill("eventLogger_", 1.5, logged_ + 0.5, event.eventAuxiliary().luminosityBlock());
0672         fill("eventLogger_", 2.5, logged_ + 0.5, event.eventAuxiliary().event());
0673         if (!correctedJets.empty())
0674           fill("eventLogger_", 3.5, logged_ + 0.5, correctedJets[0].pt());
0675         if (correctedJets.size() > 1)
0676           fill("eventLogger_", 4.5, logged_ + 0.5, correctedJets[1].pt());
0677         if (correctedJets.size() > 2)
0678           fill("eventLogger_", 5.5, logged_ + 0.5, correctedJets[2].pt());
0679         if (correctedJets.size() > 3)
0680           fill("eventLogger_", 6.5, logged_ + 0.5, correctedJets[3].pt());
0681         fill("eventLogger_", 7.5, logged_ + 0.5, wMass);
0682         fill("eventLogger_", 8.5, logged_ + 0.5, topMass);
0683         ++logged_;
0684       }
0685     }
0686   }
0687 }  // namespace TopSingleLepton
0688
0689 TopSingleLeptonDQM::TopSingleLeptonDQM(const edm::ParameterSet& cfg)
0690     : vertexSelect_(nullptr),
0691       beamspot_(""),
0692       beamspotSelect_(nullptr),
0693       MuonStep(nullptr),
0694       ElectronStep(nullptr),
0695       PvStep(nullptr),
0696       METStep(nullptr) {
0697   JetSteps.clear();
0698   CaloJetSteps.clear();
0699   PFJetSteps.clear();
0700   // configure preselection
0701   edm::ParameterSet presel = cfg.getParameter<edm::ParameterSet>("preselection");
0702   if (presel.existsAs<edm::ParameterSet>("trigger")) {
0703     edm::ParameterSet trigger = presel.getParameter<edm::ParameterSet>("trigger");
0704     triggerTable__ = consumes<edm::TriggerResults>(trigger.getParameter<edm::InputTag>("src"));
0705     triggerPaths_ = trigger.getParameter<std::vector<std::string>>("select");
0706   }
0707   if (presel.existsAs<edm::ParameterSet>("beamspot")) {
0708     edm::ParameterSet beamspot = presel.getParameter<edm::ParameterSet>("beamspot");
0709     beamspot_ = beamspot.getParameter<edm::InputTag>("src");
0710     beamspot__ = consumes<reco::BeamSpot>(beamspot.getParameter<edm::InputTag>("src"));
0711     beamspotSelect_ =
0712         std::make_unique<StringCutObjectSelector<reco::BeamSpot>>(beamspot.getParameter<std::string>("select"));
0713   }
0714
0715   // configure the selection
0716   sel_ = cfg.getParameter<std::vector<edm::ParameterSet>>("selection");
0717   setup_ = cfg.getParameter<edm::ParameterSet>("setup");
0718   for (unsigned int i = 0; i < sel_.size(); ++i) {
0719     selectionOrder_.push_back(sel_.at(i).getParameter<std::string>("label"));
0720     selection_[selectionStep(selectionOrder_.back())] =
0721         std::make_pair(sel_.at(i),
0722                        std::make_unique<TopSingleLepton::MonitorEnsemble>(
0723                            selectionStep(selectionOrder_.back()).c_str(), setup_, consumesCollector()));
0724   }
0725   for (std::vector<std::string>::const_iterator selIt = selectionOrder_.begin(); selIt != selectionOrder_.end();
0726        ++selIt) {
0727     std::string key = selectionStep(*selIt), type = objectType(*selIt);
0728     if (selection_.find(key) != selection_.end()) {
0729       if (type == "muons") {
0730         MuonStep = std::make_unique<SelectionStep<reco::PFCandidate>>(selection_[key].first, consumesCollector());
0731       }
0732       if (type == "elecs") {
0733         ElectronStep = std::make_unique<SelectionStep<reco::PFCandidate>>(selection_[key].first, consumesCollector());
0734       }
0735       if (type == "pvs") {
0736         PvStep = std::make_unique<SelectionStep<reco::Vertex>>(selection_[key].first, consumesCollector());
0737       }
0738       if (type == "jets") {
0739         JetSteps.push_back(std::make_unique<SelectionStep<reco::Jet>>(selection_[key].first, consumesCollector()));
0740       }
0741       if (type == "jets/pf") {
0742         PFJetSteps.push_back(std::make_unique<SelectionStep<reco::PFJet>>(selection_[key].first, consumesCollector()));
0743       }
0744       if (type == "jets/calo") {
0745         CaloJetSteps.push_back(
0746             std::make_unique<SelectionStep<reco::CaloJet>>(selection_[key].first, consumesCollector()));
0747       }
0748       if (type == "met") {
0749         METStep = std::make_unique<SelectionStep<reco::MET>>(selection_[key].first, consumesCollector());
0750       }
0751     }
0752   }
0753 }
0754 void TopSingleLeptonDQM::bookHistograms(DQMStore::IBooker& ibooker, edm::Run const&, edm::EventSetup const&) {
0755   for (auto selIt = selection_.begin(); selIt != selection_.end(); ++selIt) {
0756     selIt->second.second->book(ibooker);
0757   }
0758 }
0759 void TopSingleLeptonDQM::analyze(const edm::Event& event, const edm::EventSetup& setup) {
0760   if (!triggerTable__.isUninitialized()) {
0761     edm::Handle<edm::TriggerResults> triggerTable;
0762     if (!event.getByToken(triggerTable__, triggerTable))
0763       return;
0764     if (!accept(event, *triggerTable, triggerPaths_))
0765       return;
0766   }
0767   if (!beamspot__.isUninitialized()) {
0768     edm::Handle<reco::BeamSpot> beamspot;
0769     if (!event.getByToken(beamspot__, beamspot))
0770       return;
0771     if (!(*beamspotSelect_)(*beamspot))
0772       return;
0773   }
0774   //  cout<<" apply selection steps"<<endl;
0775   unsigned int passed = 0;
0776   unsigned int nJetSteps = -1;
0777   unsigned int nPFJetSteps = -1;
0778   unsigned int nCaloJetSteps = -1;
0779   for (std::vector<std::string>::const_iterator selIt = selectionOrder_.begin(); selIt != selectionOrder_.end();
0780        ++selIt) {
0781     std::string key = selectionStep(*selIt), type = objectType(*selIt);
0782     if (selection_.find(key) != selection_.end()) {
0783       if (type == "empty") {
0784         selection_[key].second->fill(event, setup);
0785       }
0786       if (type == "muons" && MuonStep != nullptr) {
0787         if (MuonStep->select(event)) {
0788           ++passed;
0789           //      cout<<"selected event! "<<selection_[key].second<<endl;
0790           selection_[key].second->fill(event, setup);
0791         } else
0792           break;
0793       }
0794       // cout<<" apply selection steps 2"<<endl;
0795       if (type == "elecs" && ElectronStep != nullptr) {
0796         // cout<<"In electrons ..."<<endl;
0797         if (ElectronStep->select(event, "electron")) {
0798           ++passed;
0799           selection_[key].second->fill(event, setup);
0800         } else
0801           break;
0802       }
0803       // cout<<" apply selection steps 3"<<endl;
0804       if (type == "pvs" && PvStep != nullptr) {
0805         if (PvStep->selectVertex(event)) {
0806           ++passed;
0807           selection_[key].second->fill(event, setup);
0808         } else
0809           break;
0810       }
0811       // cout<<" apply selection steps 4"<<endl;
0812       if (type == "jets") {
0813         nJetSteps++;
0814         if (JetSteps[nJetSteps] != nullptr) {
0815           if (JetSteps[nJetSteps]->select(event, setup)) {
0816             ++passed;
0817             selection_[key].second->fill(event, setup);
0818           } else
0819             break;
0820         }
0821       }
0822       if (type == "jets/pf") {
0823         nPFJetSteps++;
0824         if (PFJetSteps[nPFJetSteps] != nullptr) {
0825           if (PFJetSteps[nPFJetSteps]->select(event, setup)) {
0826             ++passed;
0827             selection_[key].second->fill(event, setup);
0828           } else
0829             break;
0830         }
0831       }
0832       if (type == "jets/calo") {
0833         nCaloJetSteps++;
0834         if (CaloJetSteps[nCaloJetSteps] != nullptr) {
0835           if (CaloJetSteps[nCaloJetSteps]->select(event, setup)) {
0836             ++passed;
0837             selection_[key].second->fill(event, setup);
0838           } else
0839             break;
0840         }
0841       }
0842       if (type == "met" && METStep != nullptr) {
0843         if (METStep->select(event)) {
0844           ++passed;
0845           selection_[key].second->fill(event, setup);
0846         } else
0847           break;
0848       }
0849     }
0850   }
0851 }
0852
0853 // Local Variables:
0854 // show-trailing-whitespace: t
0855 // truncate-lines: t
0856 // End: