Back to home page

Project CMSSW displayed by LXR
 
 

    

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

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