Back to home page

Project CMSSW displayed by LXR
 
 

    

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

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