Back to home page

Project CMSSW displayed by LXR

 
 

    


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

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