Back to home page

Project CMSSW displayed by LXR
 
 

    

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

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