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