Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-05-10 03:53:58 

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