Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2021-07-30 02:32:59 

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

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!