Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2025-03-23 15:57:39 

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