Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:08:04 

0001 #include "DQM/Physics/src/EwkTauDQM.h"
0002 
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004 
0005 #include "FWCore/ServiceRegistry/interface/Service.h"
0006 #include "DQMServices/Core/interface/DQMStore.h"
0007 
0008 const std::string dqmSeparator = "/";
0009 
0010 std::string dqmDirectoryName(const std::string& dqmRootDirectory, const std::string& dqmSubDirectory) {
0011   //--- concatenate names of dqmRootDirectory and dqmSubDirectory;

0012   //    add "/" separator inbetween if necessary

0013   std::string dirName = dqmRootDirectory;
0014   if (!dirName.empty() && dirName.find_last_of(dqmSeparator) != (dirName.length() - 1))
0015     dirName.append(dqmSeparator);
0016   dirName.append(dqmSubDirectory);
0017   return dirName;
0018 }
0019 
0020 EwkTauDQM::EwkTauDQM(const edm::ParameterSet& cfg) : dqmDirectory_(cfg.getParameter<std::string>("dqmDirectory")) {
0021   maxNumWarnings_ = cfg.exists("maxNumWarnings") ? cfg.getParameter<int>("maxNumWarnings") : 1;
0022 
0023   edm::ParameterSet cfgChannels = cfg.getParameter<edm::ParameterSet>("channels");
0024 
0025   edm::ParameterSet cfgElecTauChannel = cfgChannels.getParameter<edm::ParameterSet>("elecTauChannel");
0026   std::string dqmSubDirectoryElecTauChannel = cfgElecTauChannel.getParameter<std::string>("dqmSubDirectory");
0027   cfgElecTauChannel.addParameter<std::string>("dqmDirectory",
0028                                               dqmDirectoryName(dqmDirectory_, dqmSubDirectoryElecTauChannel));
0029   cfgElecTauChannel.addParameter<int>("maxNumWarnings", maxNumWarnings_);
0030   elecTauHistManager_ = new EwkElecTauHistManager(cfgElecTauChannel);
0031 
0032   edm::ParameterSet cfgMuTauChannel = cfgChannels.getParameter<edm::ParameterSet>("muTauChannel");
0033   std::string dqmSubDirectoryMuTauChannel = cfgMuTauChannel.getParameter<std::string>("dqmSubDirectory");
0034   cfgMuTauChannel.addParameter<std::string>("dqmDirectory",
0035                                             dqmDirectoryName(dqmDirectory_, dqmSubDirectoryMuTauChannel));
0036   cfgMuTauChannel.addParameter<int>("maxNumWarnings", maxNumWarnings_);
0037   muTauHistManager_ = new EwkMuTauHistManager(cfgMuTauChannel);
0038 }
0039 
0040 EwkTauDQM::~EwkTauDQM() {
0041   delete elecTauHistManager_;
0042   delete muTauHistManager_;
0043 }
0044 
0045 void EwkTauDQM::bookHistograms(DQMStore::IBooker& iBooker, edm::Run const&, edm::EventSetup const&) {
0046   elecTauHistManager_->bookHistograms(iBooker);
0047   muTauHistManager_->bookHistograms(iBooker);
0048 }
0049 
0050 void EwkTauDQM::analyze(const edm::Event& evt, const edm::EventSetup& es) {
0051   elecTauHistManager_->fillHistograms(evt, es);
0052   muTauHistManager_->fillHistograms(evt, es);
0053 }
0054 
0055 void EwkTauDQM::dqmEndRun(const edm::Run&, const edm::EventSetup&) {
0056   elecTauHistManager_->finalizeHistograms();
0057   muTauHistManager_->finalizeHistograms();
0058 }
0059 
0060 //-------------------------------------------------------------------------------

0061 // code specific to Z --> e + tau-jet channel

0062 //-------------------------------------------------------------------------------

0063 
0064 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0065 
0066 #include "FWCore/ServiceRegistry/interface/Service.h"
0067 
0068 #include "DataFormats/Common/interface/Handle.h"
0069 #include "DataFormats/Common/interface/View.h"
0070 #include "DataFormats/Common/interface/TriggerResults.h"
0071 #include "FWCore/Common/interface/TriggerNames.h"
0072 #include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
0073 #include "DataFormats/EgammaCandidates/interface/GsfElectronFwd.h"
0074 #include "DataFormats/TauReco/interface/PFTau.h"
0075 #include "DataFormats/TauReco/interface/PFTauFwd.h"
0076 #include "DataFormats/TauReco/interface/PFTauDiscriminator.h"
0077 #include "DataFormats/METReco/interface/CaloMET.h"
0078 #include "DataFormats/METReco/interface/CaloMETFwd.h"
0079 #include "DataFormats/METReco/interface/PFMET.h"
0080 #include "DataFormats/METReco/interface/PFMETFwd.h"
0081 #include "DataFormats/TrackReco/interface/Track.h"
0082 #include "DataFormats/VertexReco/interface/Vertex.h"
0083 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0084 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0085 
0086 #include "TMath.h"
0087 
0088 #include <iostream>
0089 #include <iomanip>
0090 
0091 EwkElecTauHistManager::EwkElecTauHistManager(const edm::ParameterSet& cfg)
0092     : dqmDirectory_(cfg.getParameter<std::string>("dqmDirectory")),
0093       numEventsAnalyzed_(0),
0094       numEventsSelected_(0),
0095       cfgError_(0),
0096       numWarningsTriggerResults_(0),
0097       numWarningsHLTpath_(0),
0098       numWarningsVertex_(0),
0099       numWarningsBeamSpot_(0),
0100       numWarningsElectron_(0),
0101       numWarningsTauJet_(0),
0102       numWarningsTauDiscrByLeadTrackFinding_(0),
0103       numWarningsTauDiscrByLeadTrackPtCut_(0),
0104       numWarningsTauDiscrByTrackIso_(0),
0105       numWarningsTauDiscrByEcalIso_(0),
0106       numWarningsTauDiscrAgainstElectrons_(0),
0107       numWarningsTauDiscrAgainstMuons_(0),
0108       numWarningsCaloMEt_(0),
0109       numWarningsPFMEt_(0) {
0110   triggerResultsSource_ = cfg.getParameter<edm::InputTag>("triggerResultsSource");
0111   vertexSource_ = cfg.getParameter<edm::InputTag>("vertexSource");
0112   beamSpotSource_ = cfg.getParameter<edm::InputTag>("beamSpotSource");
0113   electronSource_ = cfg.getParameter<edm::InputTag>("electronSource");
0114   tauJetSource_ = cfg.getParameter<edm::InputTag>("tauJetSource");
0115   caloMEtSource_ = cfg.getParameter<edm::InputTag>("caloMEtSource");
0116   pfMEtSource_ = cfg.getParameter<edm::InputTag>("pfMEtSource");
0117 
0118   tauDiscrByLeadTrackFinding_ = cfg.getParameter<edm::InputTag>("tauDiscrByLeadTrackFinding");
0119   tauDiscrByLeadTrackPtCut_ = cfg.getParameter<edm::InputTag>("tauDiscrByLeadTrackPtCut");
0120   tauDiscrByTrackIso_ = cfg.getParameter<edm::InputTag>("tauDiscrByTrackIso");
0121   tauDiscrByEcalIso_ = cfg.getParameter<edm::InputTag>("tauDiscrByEcalIso");
0122   tauDiscrAgainstElectrons_ = cfg.getParameter<edm::InputTag>("tauDiscrAgainstElectrons");
0123   tauDiscrAgainstMuons_ = cfg.getParameter<edm::InputTag>("tauDiscrAgainstMuons");
0124 
0125   hltPaths_ = cfg.getParameter<vstring>("hltPaths");
0126 
0127   electronEtaCut_ = cfg.getParameter<double>("electronEtaCut");
0128   electronPtCut_ = cfg.getParameter<double>("electronPtCut");
0129   electronTrackIsoCut_ = cfg.getParameter<double>("electronTrackIsoCut");
0130   electronEcalIsoCut_ = cfg.getParameter<double>("electronEcalIsoCut");
0131   std::string electronIsoMode_string = cfg.getParameter<std::string>("electronIsoMode");
0132   electronIsoMode_ = getIsoMode(electronIsoMode_string, cfgError_);
0133 
0134   tauJetEtaCut_ = cfg.getParameter<double>("tauJetEtaCut");
0135   tauJetPtCut_ = cfg.getParameter<double>("tauJetPtCut");
0136 
0137   visMassCut_ = cfg.getParameter<double>("visMassCut");
0138 
0139   maxNumWarnings_ = cfg.exists("maxNumWarnings") ? cfg.getParameter<int>("maxNumWarnings") : 1;
0140 }
0141 
0142 void EwkElecTauHistManager::bookHistograms(DQMStore::IBooker& iBooker) {
0143   iBooker.setCurrentFolder(dqmDirectory_);
0144   hElectronPt_ = iBooker.book1D("ElectronPt", "P_{T}^{e}", 20, 0., 100.);
0145   hElectronEta_ = iBooker.book1D("ElectronEta", "#eta_{e}", 20, -4.0, +4.0);
0146   hElectronPhi_ = iBooker.book1D("ElectronPhi", "#phi_{e}", 20, -TMath::Pi(), +TMath::Pi());
0147   hElectronTrackIsoPt_ = iBooker.book1D("ElectronTrackIsoPt", "Electron Track Iso.", 20, -0.01, 0.5);
0148   hElectronEcalIsoPt_ = iBooker.book1D("ElectronEcalIsoPt", "Electron Ecal Iso.", 20, -0.01, 0.5);
0149   hTauJetPt_ = iBooker.book1D("TauJetPt", "P_{T}^{#tau-Jet}", 20, 0., 100.);
0150   hTauJetEta_ = iBooker.book1D("TauJetEta", "#eta_{#tau-Jet}", 20, -4.0, +4.0);
0151   hVisMass_ = iBooker.book1D("VisMass", "e + #tau-Jet visible Mass", 20, 20., 120.);
0152   hMtElecPFMEt_ = iBooker.book1D("MtElecPFMEt", "e + E_{T}^{miss} (PF) transverse Mass", 20, 20., 120.);
0153   hElecTauAcoplanarity_ =
0154       iBooker.book1D("ElecTauAcoplanarity", "#Delta #phi_{e #tau-Jet}", 20, -TMath::Pi(), +TMath::Pi());
0155   hElecTauCharge_ = iBooker.book1D("ElecTauCharge", "Q_{e * #tau-Jet}", 5, -2.5, +2.5);
0156   hVertexZ_ = iBooker.book1D("VertexZ", "Event Vertex z-Position", 20, -25., +25.);
0157   hCaloMEtPt_ = iBooker.book1D("CaloMEtPt", "E_{T}^{miss} (Calo)", 20, 0., 100.);
0158   hPFMEtPt_ = iBooker.book1D("PFMEtPt", "E_{T}^{miss} (PF)", 20, 0., 100.);
0159   hCutFlowSummary_ = iBooker.book1D("CutFlowSummary", "Cut-flow Summary", 11, 0.5, 11.5);
0160   hCutFlowSummary_->setBinLabel(kPassedPreselection, "Preselection");
0161   hCutFlowSummary_->setBinLabel(kPassedTrigger, "HLT");
0162   hCutFlowSummary_->setBinLabel(kPassedElectronId, "e ID");
0163   hCutFlowSummary_->setBinLabel(kPassedElectronTrackIso, "e Trk Iso.");
0164   hCutFlowSummary_->setBinLabel(kPassedElectronEcalIso, "e Ecal Iso.");
0165   hCutFlowSummary_->setBinLabel(kPassedTauLeadTrack, "#tau lead. Track");
0166   hCutFlowSummary_->setBinLabel(kPassedTauLeadTrackPt, "#tau lead. Track P_{T}");
0167   hCutFlowSummary_->setBinLabel(kPassedTauDiscrAgainstElectrons, "#tau anti-e Discr.");
0168   hCutFlowSummary_->setBinLabel(kPassedTauDiscrAgainstMuons, "#tau anti-#mu Discr.");
0169   hCutFlowSummary_->setBinLabel(kPassedTauTrackIso, "#tau Track Iso.");
0170   hCutFlowSummary_->setBinLabel(kPassedTauEcalIso, "#tau Ecal Iso.");
0171 }
0172 
0173 void EwkElecTauHistManager::fillHistograms(const edm::Event& evt, const edm::EventSetup& es) {
0174   if (cfgError_)
0175     return;
0176 
0177   //-----------------------------------------------------------------------------

0178   // access event-level information

0179   //-----------------------------------------------------------------------------

0180 
0181   bool readError = false;
0182 
0183   //--- get decision of high-level trigger for the event

0184   edm::Handle<edm::TriggerResults> hltDecision;
0185   readEventData(evt,
0186                 triggerResultsSource_,
0187                 hltDecision,
0188                 numWarningsTriggerResults_,
0189                 maxNumWarnings_,
0190                 readError,
0191                 "Failed to access Trigger results");
0192   if (readError)
0193     return;
0194 
0195   const edm::TriggerNames& triggerNames = evt.triggerNames(*hltDecision);
0196 
0197   bool isTriggered = false;
0198   for (vstring::const_iterator hltPath = hltPaths_.begin(); hltPath != hltPaths_.end(); ++hltPath) {
0199     unsigned int index = triggerNames.triggerIndex(*hltPath);
0200     if (index < triggerNames.size()) {
0201       if (hltDecision->accept(index))
0202         isTriggered = true;
0203     } else {
0204       if (numWarningsHLTpath_ < maxNumWarnings_ || maxNumWarnings_ == -1)
0205         edm::LogWarning("EwkElecTauHistManager") << " Undefined HLT path = " << (*hltPath) << " !!";
0206       ++numWarningsHLTpath_;
0207       continue;
0208     }
0209   }
0210 
0211   //--- get reconstructed primary event vertex of the event

0212   //   (take as "the" primary event vertex the first entry in the collection

0213   //    of vertex objects, corresponding to the vertex associated to the highest

0214   // Pt sum of tracks)

0215   edm::Handle<reco::VertexCollection> vertexCollection;
0216   readEventData(evt,
0217                 vertexSource_,
0218                 vertexCollection,
0219                 numWarningsVertex_,
0220                 maxNumWarnings_,
0221                 readError,
0222                 "Failed to access Vertex collection");
0223   if (readError)
0224     return;
0225 
0226   const reco::Vertex* theEventVertex = (!vertexCollection->empty()) ? &(vertexCollection->at(0)) : nullptr;
0227 
0228   //--- get beam-spot (expected vertex position) for the event

0229   edm::Handle<reco::BeamSpot> beamSpot;
0230   readEventData(
0231       evt, beamSpotSource_, beamSpot, numWarningsBeamSpot_, maxNumWarnings_, readError, "Failed to access Beam-spot");
0232   if (readError)
0233     return;
0234 
0235   //--- get collections of reconstructed electrons from the event

0236   edm::Handle<reco::GsfElectronCollection> electrons;
0237   readEventData(evt,
0238                 electronSource_,
0239                 electrons,
0240                 numWarningsElectron_,
0241                 maxNumWarnings_,
0242                 readError,
0243                 "Failed to access Electron collection");
0244   if (readError)
0245     return;
0246 
0247   const reco::GsfElectron* theElectron = getTheElectron(*electrons, electronEtaCut_, electronPtCut_);
0248 
0249   double theElectronTrackIsoPt = 1.e+3;
0250   double theElectronEcalIsoPt = 1.e+3;
0251   if (theElectron) {
0252     theElectronTrackIsoPt = theElectron->dr03TkSumPt();
0253     theElectronEcalIsoPt = theElectron->dr03EcalRecHitSumEt();
0254 
0255     if (electronIsoMode_ == kRelativeIso && theElectron->pt() > 0.) {
0256       theElectronTrackIsoPt /= theElectron->pt();
0257       theElectronEcalIsoPt /= theElectron->pt();
0258     }
0259   }
0260 
0261   //--- get collections of reconstructed tau-jets from the event

0262   edm::Handle<reco::PFTauCollection> tauJets;
0263   readEventData(evt,
0264                 tauJetSource_,
0265                 tauJets,
0266                 numWarningsTauJet_,
0267                 maxNumWarnings_,
0268                 readError,
0269                 "Failed to access Tau-jet collection");
0270   if (readError)
0271     return;
0272 
0273   //--- get collections of tau-jet discriminators for those tau-jets

0274   edm::Handle<reco::PFTauDiscriminator> tauDiscrByLeadTrackFinding;
0275   readEventData(evt,
0276                 tauDiscrByLeadTrackFinding_,
0277                 tauDiscrByLeadTrackFinding,
0278                 numWarningsTauDiscrByLeadTrackFinding_,
0279                 maxNumWarnings_,
0280                 readError,
0281                 "Failed to access collection of pf. Tau discriminators by "
0282                 "leading Track finding");
0283   edm::Handle<reco::PFTauDiscriminator> tauDiscrByLeadTrackPtCut;
0284   readEventData(evt,
0285                 tauDiscrByLeadTrackPtCut_,
0286                 tauDiscrByLeadTrackPtCut,
0287                 numWarningsTauDiscrByLeadTrackPtCut_,
0288                 maxNumWarnings_,
0289                 readError,
0290                 "Failed to access collection of pf. Tau discriminators by "
0291                 "leading Track Pt cut");
0292   edm::Handle<reco::PFTauDiscriminator> tauDiscrByTrackIso;
0293   readEventData(evt,
0294                 tauDiscrByTrackIso_,
0295                 tauDiscrByTrackIso,
0296                 numWarningsTauDiscrByTrackIso_,
0297                 maxNumWarnings_,
0298                 readError,
0299                 "Failed to access collection of pf. Tau discriminators by "
0300                 "Track isolation");
0301   edm::Handle<reco::PFTauDiscriminator> tauDiscrByEcalIso;
0302   readEventData(evt,
0303                 tauDiscrByTrackIso_,
0304                 tauDiscrByEcalIso,
0305                 numWarningsTauDiscrByEcalIso_,
0306                 maxNumWarnings_,
0307                 readError,
0308                 "Failed to access collection of pf. Tau discriminators by ECAL "
0309                 "isolation");
0310   edm::Handle<reco::PFTauDiscriminator> tauDiscrAgainstElectrons;
0311   readEventData(evt,
0312                 tauDiscrAgainstElectrons_,
0313                 tauDiscrAgainstElectrons,
0314                 numWarningsTauDiscrAgainstElectrons_,
0315                 maxNumWarnings_,
0316                 readError,
0317                 "Failed to access collection of pf. Tau discriminators against "
0318                 "Electrons");
0319   edm::Handle<reco::PFTauDiscriminator> tauDiscrAgainstMuons;
0320   readEventData(evt,
0321                 tauDiscrAgainstMuons_,
0322                 tauDiscrAgainstMuons,
0323                 numWarningsTauDiscrAgainstMuons_,
0324                 maxNumWarnings_,
0325                 readError,
0326                 "Failed to access collection of pf. Tau discriminators against Muons");
0327   if (readError)
0328     return;
0329 
0330   int theTauJetIndex = -1;
0331   const reco::PFTau* theTauJet = getTheTauJet(*tauJets, tauJetEtaCut_, tauJetPtCut_, theTauJetIndex);
0332 
0333   double theTauDiscrByLeadTrackFinding = -1.;
0334   double theTauDiscrByLeadTrackPtCut = -1.;
0335   double theTauDiscrByTrackIso = -1.;
0336   double theTauDiscrByEcalIso = -1.;
0337   double theTauDiscrAgainstElectrons = -1.;
0338   double theTauDiscrAgainstMuons = -1.;
0339   if (theTauJetIndex != -1) {
0340     reco::PFTauRef theTauJetRef(tauJets, theTauJetIndex);
0341     theTauDiscrByLeadTrackFinding = (*tauDiscrByLeadTrackFinding)[theTauJetRef];
0342     theTauDiscrByLeadTrackPtCut = (*tauDiscrByLeadTrackPtCut)[theTauJetRef];
0343     theTauDiscrByTrackIso = (*tauDiscrByTrackIso)[theTauJetRef];
0344     theTauDiscrByEcalIso = (*tauDiscrByEcalIso)[theTauJetRef];
0345     theTauDiscrAgainstElectrons = (*tauDiscrAgainstElectrons)[theTauJetRef];
0346     theTauDiscrAgainstMuons = (*tauDiscrAgainstMuons)[theTauJetRef];
0347   }
0348 
0349   //--- get missing transverse momentum

0350   //    measured by calorimeters/reconstructed by particle-flow algorithm

0351   edm::Handle<reco::CaloMETCollection> caloMEtCollection;
0352   readEventData(evt,
0353                 caloMEtSource_,
0354                 caloMEtCollection,
0355                 numWarningsCaloMEt_,
0356                 maxNumWarnings_,
0357                 readError,
0358                 "Failed to access calo. MET collection");
0359   if (readError)
0360     return;
0361 
0362   const reco::CaloMET& caloMEt = caloMEtCollection->at(0);
0363 
0364   edm::Handle<reco::PFMETCollection> pfMEtCollection;
0365   readEventData(evt,
0366                 pfMEtSource_,
0367                 pfMEtCollection,
0368                 numWarningsPFMEt_,
0369                 maxNumWarnings_,
0370                 readError,
0371                 "Failed to access pf. MET collection");
0372   if (readError)
0373     return;
0374 
0375   const reco::PFMET& pfMEt = pfMEtCollection->at(0);
0376 
0377   if (!(theElectron && theTauJet && theTauJetIndex != -1))
0378     return;
0379 
0380   //-----------------------------------------------------------------------------

0381   // compute EWK tau analysis specific quantities

0382   //-----------------------------------------------------------------------------

0383 
0384   double dPhiElecTau = calcDeltaPhi(theElectron->phi(), theTauJet->phi());
0385 
0386   double mElecTau = (theElectron->p4() + theTauJet->p4()).M();
0387 
0388   // double mtElecCaloMEt = calcMt(theElectron->px(), theElectron->py(),

0389   // caloMEt.px(), caloMEt.py());

0390   double mtElecPFMEt = calcMt(theElectron->px(), theElectron->py(), pfMEt.px(), pfMEt.py());
0391 
0392   // double pZetaCaloMEt = calcPzeta(theElectron->p4(), theTauJet->p4(),

0393   // caloMEt.px(), caloMEt.py());

0394   // double pZetaPFMEt = calcPzeta(theElectron->p4(), theTauJet->p4(),

0395   // pfMEt.px(), pfMEt.py());

0396 
0397   //-----------------------------------------------------------------------------

0398   // apply selection criteria; fill histograms

0399   //-----------------------------------------------------------------------------

0400 
0401   ++numEventsAnalyzed_;
0402 
0403   bool isSelected = false;
0404   bool fullSelect = false;
0405   int cutFlowStatus = -1;
0406 
0407   if (mElecTau > visMassCut_) {
0408     cutFlowStatus = kPassedPreselection;
0409   }
0410   if (cutFlowStatus == kPassedPreselection && (isTriggered || hltPaths_.empty())) {
0411     cutFlowStatus = kPassedTrigger;
0412   }
0413   if (cutFlowStatus == kPassedTrigger && passesElectronId(*theElectron)) {
0414     cutFlowStatus = kPassedElectronId;
0415     hElectronTrackIsoPt_->Fill(theElectronTrackIsoPt);
0416   }
0417   if (cutFlowStatus == kPassedElectronId && theElectronTrackIsoPt < electronTrackIsoCut_) {
0418     cutFlowStatus = kPassedElectronTrackIso;
0419     hElectronEcalIsoPt_->Fill(theElectronEcalIsoPt);
0420   }
0421   if (cutFlowStatus == kPassedElectronTrackIso && theElectronEcalIsoPt < electronEcalIsoCut_) {
0422     cutFlowStatus = kPassedElectronEcalIso;
0423   }
0424   if (cutFlowStatus == kPassedElectronEcalIso && theTauDiscrByLeadTrackFinding > 0.5) {
0425     cutFlowStatus = kPassedTauLeadTrack;
0426     // if ( theTauJet->leadTrack().isAvailable() )

0427     // hTauLeadTrackPt_->Fill(theTauJet->leadTrack()->pt());

0428   }
0429   if (cutFlowStatus == kPassedTauLeadTrack && theTauDiscrByLeadTrackPtCut > 0.5) {
0430     cutFlowStatus = kPassedTauLeadTrackPt;
0431     // hTauTrackIsoPt_->Fill(theTauJet->isolationPFChargedHadrCandsPtSum());

0432   }
0433   if (cutFlowStatus == kPassedTauLeadTrackPt && theTauDiscrAgainstElectrons > 0.5) {
0434     cutFlowStatus = kPassedTauDiscrAgainstElectrons;
0435     // hTauDiscrAgainstMuons_->Fill(theTauDiscrAgainstMuons);

0436   }
0437   if (cutFlowStatus == kPassedTauDiscrAgainstElectrons && theTauDiscrAgainstMuons > 0.5) {
0438     cutFlowStatus = kPassedTauDiscrAgainstMuons;
0439     isSelected = true;
0440   }
0441   if (cutFlowStatus == kPassedTauDiscrAgainstMuons && theTauDiscrByTrackIso > 0.5) {
0442     cutFlowStatus = kPassedTauTrackIso;
0443     // hTauEcalIsoPt_->Fill(theTauJet->isolationPFGammaCandsEtSum());

0444   }
0445   if (cutFlowStatus == kPassedTauTrackIso && theTauDiscrByEcalIso > 0.5) {
0446     cutFlowStatus = kPassedTauEcalIso;
0447     fullSelect = true;
0448     // hTauDiscrAgainstElectrons_->Fill(theTauDiscrAgainstElectrons);

0449   }
0450 
0451   for (int iCut = 1; iCut <= cutFlowStatus; ++iCut) {
0452     hCutFlowSummary_->Fill(iCut);
0453   }
0454 
0455   if (isSelected) {
0456     hElectronPt_->Fill(theElectron->pt());
0457     hElectronEta_->Fill(theElectron->eta());
0458     hElectronPhi_->Fill(theElectron->phi());
0459 
0460     hTauJetPt_->Fill(theTauJet->pt());
0461     hTauJetEta_->Fill(theTauJet->eta());
0462     // hTauJetPhi_->Fill(theTauJet->phi());

0463 
0464     // hTauJetCharge_->Fill(theTauJet->charge());

0465     // if ( theTauJet->signalTracks().isAvailable()    )

0466     // hTauJetNumSignalTracks_->Fill(theTauJet->signalTracks().size());

0467     // if ( theTauJet->isolationTracks().isAvailable() )

0468     // hTauJetNumIsoTracks_->Fill(theTauJet->isolationTracks().size());

0469 
0470     if (fullSelect) {
0471       hVisMass_->Fill(mElecTau);
0472     }
0473     // hMtElecCaloMEt_->Fill(mtElecCaloMEt);

0474     hMtElecPFMEt_->Fill(mtElecPFMEt);
0475     // hPzetaCaloMEt_->Fill(pZetaCaloMEt);

0476     // hPzetaPFMEt_->Fill(pZetaPFMEt);

0477     hElecTauAcoplanarity_->Fill(dPhiElecTau);
0478     hElecTauCharge_->Fill(theElectron->charge() * theTauJet->charge());
0479 
0480     if (theEventVertex) {
0481       // hVertexChi2_->Fill(theEventVertex->normalizedChi2());

0482       hVertexZ_->Fill(theEventVertex->z());
0483       // hVertexD0_->Fill(getVertexD0(*theEventVertex, *beamSpot));

0484     }
0485 
0486     hCaloMEtPt_->Fill(caloMEt.pt());
0487     // hCaloMEtPhi_->Fill(caloMEt.phi());

0488 
0489     hPFMEtPt_->Fill(pfMEt.pt());
0490     // hPFMEtPhi_->Fill(pfMEt.phi());

0491   }
0492 
0493   if (isSelected)
0494     ++numEventsSelected_;
0495 }
0496 
0497 void EwkElecTauHistManager::finalizeHistograms() {
0498   edm::LogInfo("EwkElecTauHistManager") << "Filter-Statistics Summary:" << std::endl
0499                                         << " Events analyzed = " << numEventsAnalyzed_ << std::endl
0500                                         << " Events selected = " << numEventsSelected_;
0501   if (numEventsAnalyzed_ > 0) {
0502     double eff = numEventsSelected_ / (double)numEventsAnalyzed_;
0503     edm::LogInfo("") << "Overall efficiency = " << std::setprecision(4) << eff * 100. << " +/- " << std::setprecision(4)
0504                      << TMath::Sqrt(eff * (1 - eff) / numEventsAnalyzed_) * 100. << ")%";
0505   }
0506 }
0507 
0508 //-------------------------------------------------------------------------------

0509 // code specific to Z --> mu + tau-jet channel

0510 //-------------------------------------------------------------------------------

0511 
0512 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0513 
0514 #include "FWCore/ServiceRegistry/interface/Service.h"
0515 
0516 #include "DataFormats/Common/interface/Handle.h"
0517 #include "DataFormats/Common/interface/View.h"
0518 #include "DataFormats/Common/interface/TriggerResults.h"
0519 #include "FWCore/Common/interface/TriggerNames.h"
0520 #include "DataFormats/MuonReco/interface/Muon.h"
0521 #include "DataFormats/MuonReco/interface/MuonFwd.h"
0522 #include "DataFormats/TauReco/interface/PFTau.h"
0523 #include "DataFormats/TauReco/interface/PFTauFwd.h"
0524 #include "DataFormats/TauReco/interface/PFTauDiscriminator.h"
0525 #include "DataFormats/METReco/interface/CaloMET.h"
0526 #include "DataFormats/METReco/interface/CaloMETFwd.h"
0527 #include "DataFormats/METReco/interface/PFMET.h"
0528 #include "DataFormats/METReco/interface/PFMETFwd.h"
0529 #include "DataFormats/TrackReco/interface/Track.h"
0530 #include "DataFormats/VertexReco/interface/Vertex.h"
0531 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0532 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0533 
0534 #include "TMath.h"
0535 
0536 #include <iostream>
0537 #include <iomanip>
0538 
0539 EwkMuTauHistManager::EwkMuTauHistManager(const edm::ParameterSet& cfg)
0540     : dqmDirectory_(cfg.getParameter<std::string>("dqmDirectory")),
0541       numEventsAnalyzed_(0),
0542       numEventsSelected_(0),
0543       cfgError_(0),
0544       numWarningsTriggerResults_(0),
0545       numWarningsHLTpath_(0),
0546       numWarningsVertex_(0),
0547       numWarningsBeamSpot_(0),
0548       numWarningsMuon_(0),
0549       numWarningsTauJet_(0),
0550       numWarningsTauDiscrByLeadTrackFinding_(0),
0551       numWarningsTauDiscrByLeadTrackPtCut_(0),
0552       numWarningsTauDiscrByTrackIso_(0),
0553       numWarningsTauDiscrByEcalIso_(0),
0554       numWarningsTauDiscrAgainstMuons_(0),
0555       numWarningsCaloMEt_(0),
0556       numWarningsPFMEt_(0) {
0557   triggerResultsSource_ = cfg.getParameter<edm::InputTag>("triggerResultsSource");
0558   vertexSource_ = cfg.getParameter<edm::InputTag>("vertexSource");
0559   beamSpotSource_ = cfg.getParameter<edm::InputTag>("beamSpotSource");
0560   muonSource_ = cfg.getParameter<edm::InputTag>("muonSource");
0561   tauJetSource_ = cfg.getParameter<edm::InputTag>("tauJetSource");
0562   caloMEtSource_ = cfg.getParameter<edm::InputTag>("caloMEtSource");
0563   pfMEtSource_ = cfg.getParameter<edm::InputTag>("pfMEtSource");
0564 
0565   tauDiscrByLeadTrackFinding_ = cfg.getParameter<edm::InputTag>("tauDiscrByLeadTrackFinding");
0566   tauDiscrByLeadTrackPtCut_ = cfg.getParameter<edm::InputTag>("tauDiscrByLeadTrackPtCut");
0567   tauDiscrByTrackIso_ = cfg.getParameter<edm::InputTag>("tauDiscrByTrackIso");
0568   tauDiscrByEcalIso_ = cfg.getParameter<edm::InputTag>("tauDiscrByEcalIso");
0569   tauDiscrAgainstMuons_ = cfg.getParameter<edm::InputTag>("tauDiscrAgainstMuons");
0570 
0571   hltPaths_ = cfg.getParameter<vstring>("hltPaths");
0572 
0573   muonEtaCut_ = cfg.getParameter<double>("muonEtaCut");
0574   muonPtCut_ = cfg.getParameter<double>("muonPtCut");
0575   muonTrackIsoCut_ = cfg.getParameter<double>("muonTrackIsoCut");
0576   muonEcalIsoCut_ = cfg.getParameter<double>("muonEcalIsoCut");
0577   muonCombIsoCut_ = cfg.getParameter<double>("muonCombIsoCut");
0578   std::string muonIsoMode_string = cfg.getParameter<std::string>("muonIsoMode");
0579   muonIsoMode_ = getIsoMode(muonIsoMode_string, cfgError_);
0580 
0581   tauJetEtaCut_ = cfg.getParameter<double>("tauJetEtaCut");
0582   tauJetPtCut_ = cfg.getParameter<double>("tauJetPtCut");
0583 
0584   visMassCut_ = cfg.getParameter<double>("visMassCut");
0585   deltaRCut_ = cfg.getParameter<double>("deltaRCut");
0586 
0587   maxNumWarnings_ = cfg.exists("maxNumWarnings") ? cfg.getParameter<int>("maxNumWarnings") : 1;
0588 }
0589 
0590 void EwkMuTauHistManager::bookHistograms(DQMStore::IBooker& iBooker) {
0591   iBooker.setCurrentFolder(dqmDirectory_);
0592 
0593   hMuonPt_ = iBooker.book1D("MuonPt", "P_{T}^{#mu}", 20, 0., 100.);
0594   hMuonEta_ = iBooker.book1D("MuonEta", "#eta_{#mu}", 20, -4.0, +4.0);
0595   hMuonPhi_ = iBooker.book1D("MuonPhi", "#phi_{#mu}", 20, -TMath::Pi(), +TMath::Pi());
0596   hMuonTrackIsoPt_ = iBooker.book1D("MuonTrackIsoPt", "Muon Track Iso.", 20, -0.01, 10.);
0597   hMuonEcalIsoPt_ = iBooker.book1D("MuonEcalIsoPt", "Muon Ecal Iso.", 20, -0.01, 10.);
0598   hMuonCombIsoPt_ = iBooker.book1D("MuonCombIsoPt", "Muon Comb Iso.", 20, -0.01, 1.);
0599 
0600   hTauJetPt_ = iBooker.book1D("TauJetPt", "P_{T}^{#tau-Jet}", 20, 0., 100.);
0601   hTauJetEta_ = iBooker.book1D("TauJetEta", "#eta_{#tau-Jet}", 20, -4.0, +4.0);
0602   hTauJetPhi_ = iBooker.book1D("TauJetPhi", "#phi_{#tau-Jet}", 20, -TMath::Pi(), +TMath::Pi());
0603   hTauLeadTrackPt_ = iBooker.book1D("TauLeadTrackPt", "P_{T}^{#tau-Jetldg trk}", 20, 0., 50.);
0604   hTauTrackIsoPt_ = iBooker.book1D("TauTrackIsoPt", "Tau Track Iso.", 20, -0.01, 40.);
0605   hTauEcalIsoPt_ = iBooker.book1D("TauEcalIsoPt", "Tau Ecal Iso.", 10, -0.01, 10.);
0606   hTauDiscrAgainstMuons_ = iBooker.book1D("TauDiscrAgainstMuons", "Tau Discr. against Muons", 2, -0.5, +1.5);
0607   hTauJetNumSignalTracks_ = iBooker.book1D("TauJetNumSignalTracks", "Num. Tau signal Cone Tracks", 20, -0.5, +19.5);
0608   hTauJetNumIsoTracks_ = iBooker.book1D("TauJetNumIsoTracks", "Num. Tau isolation Cone Tracks", 20, -0.5, +19.5);
0609 
0610   hVisMass_ = iBooker.book1D("VisMass", "#mu + #tau-Jet visible Mass", 20, 0., 120.);
0611   hVisMassFinal_ = iBooker.book1D("VisMassFinal", "#mu + #tau-Jet visible final Mass", 20, 0., 120.);
0612   hMtMuPFMEt_ = iBooker.book1D("MtMuPFMEt", "#mu + E_{T}^{miss} (PF) transverse Mass", 20, 0., 120.);
0613   hMuTauAcoplanarity_ =
0614       iBooker.book1D("MuTauAcoplanarity", "#Delta #phi_{#mu #tau-Jet}", 20, -TMath::Pi(), +TMath::Pi());
0615   hMuTauDeltaR_ = iBooker.book1D("MuTauDeltaR", "#Delta R_{#mu #tau-Jet}", 20, 0, 5);
0616   hVertexZ_ = iBooker.book1D("VertexZ", "Event Vertex z-Position", 20, -25., +25.);
0617   hCaloMEtPt_ = iBooker.book1D("CaloMEtPt", "E_{T}^{miss} (Calo)", 20, 0., 100.);
0618   hPFMEtPt_ = iBooker.book1D("PFMEtPt", "E_{T}^{miss} (PF)", 20, 0., 100.);
0619   hCutFlowSummary_ = iBooker.book1D("CutFlowSummary", "Cut-flow Summary", 11, 0.5, 11.5);
0620   hCutFlowSummary_->setBinLabel(kPassedPreselection, "Preselection");
0621   hCutFlowSummary_->setBinLabel(kPassedTrigger, "HLT");
0622   hCutFlowSummary_->setBinLabel(kPassedMuonId, "#mu ID");
0623   hCutFlowSummary_->setBinLabel(kPassedMuonTrackIso, "#mu Trk Iso.");
0624   hCutFlowSummary_->setBinLabel(kPassedMuonEcalIso, "#mu Ecal Iso.");
0625   hCutFlowSummary_->setBinLabel(kPassedTauLeadTrack, "#tau lead. Track");
0626   hCutFlowSummary_->setBinLabel(kPassedTauLeadTrackPt, "#tau lead. Track P_{T}");
0627   hCutFlowSummary_->setBinLabel(kPassedTauTrackIso, "#tau Track Iso.");
0628   hCutFlowSummary_->setBinLabel(kPassedTauEcalIso, "#tau Ecal Iso.");
0629   hCutFlowSummary_->setBinLabel(kPassedTauDiscrAgainstMuons, "#tau anti-#mu Discr.");
0630   hCutFlowSummary_->setBinLabel(kPassedDeltaR, "#DeltaR(#mu,#tau) ");
0631 }
0632 
0633 void EwkMuTauHistManager::fillHistograms(const edm::Event& evt, const edm::EventSetup& es) {
0634   if (cfgError_)
0635     return;
0636 
0637   //-----------------------------------------------------------------------------

0638   // access event-level information

0639   //-----------------------------------------------------------------------------

0640 
0641   bool readError = false;
0642 
0643   //--- get decision of high-level trigger for the event

0644   edm::Handle<edm::TriggerResults> hltDecision;
0645   readEventData(evt,
0646                 triggerResultsSource_,
0647                 hltDecision,
0648                 numWarningsTriggerResults_,
0649                 maxNumWarnings_,
0650                 readError,
0651                 "Failed to access Trigger results");
0652   if (readError)
0653     return;
0654 
0655   const edm::TriggerNames& triggerNames = evt.triggerNames(*hltDecision);
0656 
0657   bool isTriggered = false;
0658   for (vstring::const_iterator hltPath = hltPaths_.begin(); hltPath != hltPaths_.end(); ++hltPath) {
0659     unsigned int index = triggerNames.triggerIndex(*hltPath);
0660     if (index < triggerNames.size()) {
0661       if (hltDecision->accept(index))
0662         isTriggered = true;
0663     } else {
0664       if (numWarningsHLTpath_ < maxNumWarnings_ || maxNumWarnings_ == -1)
0665         edm::LogWarning("EwkMuTauHistManager") << " Undefined HLT path = " << (*hltPath) << " !!";
0666       ++numWarningsHLTpath_;
0667       continue;
0668     }
0669   }
0670 
0671   //--- get reconstructed primary event vertex of the event

0672   //   (take as "the" primary event vertex the first entry in the collection

0673   //    of vertex objects, corresponding to the vertex associated to the highest

0674   // Pt sum of tracks)

0675   edm::Handle<reco::VertexCollection> vertexCollection;
0676   readEventData(evt,
0677                 vertexSource_,
0678                 vertexCollection,
0679                 numWarningsVertex_,
0680                 maxNumWarnings_,
0681                 readError,
0682                 "Failed to access Vertex collection");
0683   if (readError)
0684     return;
0685 
0686   const reco::Vertex* theEventVertex = (!vertexCollection->empty()) ? &(vertexCollection->at(0)) : nullptr;
0687 
0688   //--- get beam-spot (expected vertex position) for the event

0689   edm::Handle<reco::BeamSpot> beamSpot;
0690   readEventData(
0691       evt, beamSpotSource_, beamSpot, numWarningsBeamSpot_, maxNumWarnings_, readError, "Failed to access Beam-spot");
0692   if (readError)
0693     return;
0694 
0695   //--- get collections of reconstructed muons from the event

0696   edm::Handle<reco::MuonCollection> muons;
0697   readEventData(
0698       evt, muonSource_, muons, numWarningsMuon_, maxNumWarnings_, readError, "Failed to access Muon collection");
0699   if (readError)
0700     return;
0701 
0702   const reco::Muon* theMuon = getTheMuon(*muons, muonEtaCut_, muonPtCut_);
0703 
0704   double theMuonTrackIsoPt = 1.e+3;
0705   double theMuonEcalIsoPt = 1.e+3;
0706   double theMuonCombIsoPt = 1.e+3;
0707 
0708   if (theMuon) {
0709     theMuonTrackIsoPt = theMuon->isolationR05().sumPt;
0710     // mu.isolationR05().emEt + mu.isolationR05().hadEt +

0711     // mu.isolationR05().sumPt

0712     theMuonEcalIsoPt = theMuon->isolationR05().emEt;
0713 
0714     if (muonIsoMode_ == kRelativeIso && theMuon->pt() > 0.) {
0715       theMuonTrackIsoPt /= theMuon->pt();
0716       theMuonEcalIsoPt /= theMuon->pt();
0717       theMuonCombIsoPt = (theMuon->isolationR05().sumPt + theMuon->isolationR05().emEt) / theMuon->pt();
0718       // std::cout<<"Rel Iso ="<<theMuonCombIsoPt<<std::endl;

0719     }
0720   }
0721 
0722   //--- get collections of reconstructed tau-jets from the event

0723   edm::Handle<reco::PFTauCollection> tauJets;
0724   readEventData(evt,
0725                 tauJetSource_,
0726                 tauJets,
0727                 numWarningsTauJet_,
0728                 maxNumWarnings_,
0729                 readError,
0730                 "Failed to access Tau-jet collection");
0731   if (readError)
0732     return;
0733 
0734   //--- get collections of tau-jet discriminators for those tau-jets

0735   edm::Handle<reco::PFTauDiscriminator> tauDiscrByLeadTrackFinding;
0736   readEventData(evt,
0737                 tauDiscrByLeadTrackFinding_,
0738                 tauDiscrByLeadTrackFinding,
0739                 numWarningsTauDiscrByLeadTrackFinding_,
0740                 maxNumWarnings_,
0741                 readError,
0742                 "Failed to access collection of pf. Tau discriminators by "
0743                 "leading Track finding");
0744   edm::Handle<reco::PFTauDiscriminator> tauDiscrByLeadTrackPtCut;
0745   readEventData(evt,
0746                 tauDiscrByLeadTrackPtCut_,
0747                 tauDiscrByLeadTrackPtCut,
0748                 numWarningsTauDiscrByLeadTrackPtCut_,
0749                 maxNumWarnings_,
0750                 readError,
0751                 "Failed to access collection of pf. Tau discriminators by "
0752                 "leading Track Pt cut");
0753   edm::Handle<reco::PFTauDiscriminator> tauDiscrByTrackIso;
0754   readEventData(evt,
0755                 tauDiscrByTrackIso_,
0756                 tauDiscrByTrackIso,
0757                 numWarningsTauDiscrByTrackIso_,
0758                 maxNumWarnings_,
0759                 readError,
0760                 "Failed to access collection of pf. Tau discriminators by "
0761                 "Track isolation");
0762   edm::Handle<reco::PFTauDiscriminator> tauDiscrByEcalIso;
0763   readEventData(evt,
0764                 tauDiscrByTrackIso_,
0765                 tauDiscrByEcalIso,
0766                 numWarningsTauDiscrByEcalIso_,
0767                 maxNumWarnings_,
0768                 readError,
0769                 "Failed to access collection of pf. Tau discriminators by ECAL "
0770                 "isolation");
0771   edm::Handle<reco::PFTauDiscriminator> tauDiscrAgainstMuons;
0772   readEventData(evt,
0773                 tauDiscrAgainstMuons_,
0774                 tauDiscrAgainstMuons,
0775                 numWarningsTauDiscrAgainstMuons_,
0776                 maxNumWarnings_,
0777                 readError,
0778                 "Failed to access collection of pf. Tau discriminators against Muons");
0779   if (readError)
0780     return;
0781 
0782   int theTauJetIndex = -1;
0783   const reco::PFTau* theTauJet = getTheTauJet(*tauJets, tauJetEtaCut_, tauJetPtCut_, theTauJetIndex);
0784 
0785   double theTauDiscrByLeadTrackFinding = -1.;
0786   double theTauDiscrByLeadTrackPtCut = -1.;
0787   double theTauDiscrByTrackIso = -1.;
0788   double theTauDiscrByEcalIso = -1.;
0789   double theTauDiscrAgainstMuons = -1.;
0790   if (theTauJetIndex != -1) {
0791     reco::PFTauRef theTauJetRef(tauJets, theTauJetIndex);
0792     theTauDiscrByLeadTrackFinding = (*tauDiscrByLeadTrackFinding)[theTauJetRef];
0793     theTauDiscrByLeadTrackPtCut = (*tauDiscrByLeadTrackPtCut)[theTauJetRef];
0794     theTauDiscrByTrackIso = (*tauDiscrByTrackIso)[theTauJetRef];
0795     theTauDiscrByEcalIso = (*tauDiscrByEcalIso)[theTauJetRef];
0796     theTauDiscrAgainstMuons = (*tauDiscrAgainstMuons)[theTauJetRef];
0797   }
0798 
0799   //--- get missing transverse momentum

0800   //    measured by calorimeters/reconstructed by particle-flow algorithm

0801   edm::Handle<reco::CaloMETCollection> caloMEtCollection;
0802   readEventData(evt,
0803                 caloMEtSource_,
0804                 caloMEtCollection,
0805                 numWarningsCaloMEt_,
0806                 maxNumWarnings_,
0807                 readError,
0808                 "Failed to access calo. MET collection");
0809   if (readError)
0810     return;
0811 
0812   const reco::CaloMET& caloMEt = caloMEtCollection->at(0);
0813 
0814   edm::Handle<reco::PFMETCollection> pfMEtCollection;
0815   readEventData(evt,
0816                 pfMEtSource_,
0817                 pfMEtCollection,
0818                 numWarningsPFMEt_,
0819                 maxNumWarnings_,
0820                 readError,
0821                 "Failed to access pf. MET collection");
0822   if (readError)
0823     return;
0824 
0825   const reco::PFMET& pfMEt = pfMEtCollection->at(0);
0826 
0827   if (!(theMuon && theTauJet && theTauJetIndex != -1))
0828     return;
0829 
0830   //-----------------------------------------------------------------------------

0831   // compute EWK tau analysis specific quantities

0832   //-----------------------------------------------------------------------------

0833 
0834   double dPhiMuTau = calcDeltaPhi(theMuon->phi(), theTauJet->phi());
0835   // double dRMuTau = calcDeltaR(theMuon->p4(), theTauJet->p4());

0836   double dRMuTau = fabs(ROOT::Math::VectorUtil::DeltaR(theMuon->p4(), theTauJet->p4()));
0837   double mMuTau = (theMuon->p4() + theTauJet->p4()).M();
0838 
0839   // double mtMuCaloMEt = calcMt(theMuon->px(), theMuon->px(), caloMEt.px(),

0840   // caloMEt.py());

0841   double mtMuPFMEt = calcMt(theMuon->px(), theMuon->px(), pfMEt.px(), pfMEt.py());
0842 
0843   // double pZetaCaloMEt = calcPzeta(theMuon->p4(), theTauJet->p4(),

0844   // caloMEt.px(), caloMEt.py());

0845   // double pZetaPFMEt = calcPzeta(theMuon->p4(), theTauJet->p4(), pfMEt.px(),

0846   // pfMEt.py());

0847 
0848   //-----------------------------------------------------------------------------

0849   // apply selection criteria; fill histograms

0850   //-----------------------------------------------------------------------------

0851 
0852   ++numEventsAnalyzed_;
0853 
0854   bool isSelected = false;
0855   int cutFlowStatus = -1;
0856 
0857   // if ( muonIsoMode_ == kAbsoluteIso){

0858   if (mMuTau > visMassCut_) {
0859     cutFlowStatus = kPassedPreselection;
0860   }
0861   if (cutFlowStatus == kPassedPreselection && (isTriggered || hltPaths_.empty())) {
0862     cutFlowStatus = kPassedTrigger;
0863   }
0864   if (cutFlowStatus == kPassedTrigger && (theMuon->isGlobalMuon() || theMuon->isTrackerMuon())) {
0865     cutFlowStatus = kPassedMuonId;
0866   }
0867 
0868   if (cutFlowStatus == kPassedMuonId && (theTauDiscrByLeadTrackFinding > 0.5) && (theTauJet->eta() < tauJetEtaCut_) &&
0869       (theTauJet->pt() > tauJetPtCut_)) {
0870     cutFlowStatus = kPassedTauLeadTrack;
0871   }
0872   if (cutFlowStatus == kPassedTauLeadTrack && theTauDiscrByLeadTrackPtCut > 0.5) {
0873     cutFlowStatus = kPassedTauLeadTrackPt;
0874     // hTauTrackIsoPt_->Fill(theTauJet->isolationPFChargedHadrCandsPtSum());

0875   }
0876   if (cutFlowStatus == kPassedTauLeadTrackPt && theTauDiscrAgainstMuons > 0.5) {
0877     cutFlowStatus = kPassedTauDiscrAgainstMuons;
0878     // hTauEcalIsoPt_->Fill(theTauJet->isolationPFGammaCandsEtSum());

0879   }
0880   if (cutFlowStatus == kPassedTauDiscrAgainstMuons && dRMuTau > deltaRCut_) {
0881     cutFlowStatus = kPassedDeltaR;
0882     // hTauDiscrAgainstMuons_->Fill(theTauDiscrAgainstMuons);

0883 
0884     hMuonPt_->Fill(theMuon->pt());
0885     hMuonEta_->Fill(theMuon->eta());
0886     hMuonPhi_->Fill(theMuon->phi());
0887 
0888     hTauJetPt_->Fill(theTauJet->pt());
0889     hTauJetEta_->Fill(theTauJet->eta());
0890     hTauJetPhi_->Fill(theTauJet->phi());
0891 
0892     // hTauJetCharge_->Fill(theTauJet->charge());

0893     if (theTauJet->signalTracks().isAvailable())
0894       hTauJetNumSignalTracks_->Fill(theTauJet->signalTracks().size());
0895     if (theTauJet->isolationTracks().isAvailable())
0896       hTauJetNumIsoTracks_->Fill(theTauJet->isolationTracks().size());
0897 
0898     hVisMass_->Fill(mMuTau);
0899     // hMtMuCaloMEt_->Fill(mtMuCaloMEt);

0900     hMtMuPFMEt_->Fill(mtMuPFMEt);
0901     // hPzetaCaloMEt_->Fill(pZetaCaloMEt);

0902     // hPzetaPFMEt_->Fill(pZetaPFMEt);

0903     hMuTauAcoplanarity_->Fill(dPhiMuTau);
0904     hMuTauDeltaR_->Fill(dRMuTau);
0905     // hMuTauCharge_->Fill(theMuon->charge() + theTauJet->charge());

0906 
0907     if (theEventVertex) {
0908       // hVertexChi2_->Fill(theEventVertex->normalizedChi2());

0909       hVertexZ_->Fill(theEventVertex->z());
0910       // hVertexD0_->Fill(getVertexD0(*theEventVertex, *beamSpot));

0911     }
0912 
0913     hCaloMEtPt_->Fill(caloMEt.pt());
0914     // hCaloMEtPhi_->Fill(caloMEt.phi());

0915 
0916     hPFMEtPt_->Fill(pfMEt.pt());
0917     // hPFMEtPhi_->Fill(pfMEt.phi());

0918     hMuonTrackIsoPt_->Fill(theMuonTrackIsoPt);
0919     hMuonEcalIsoPt_->Fill(theMuonEcalIsoPt);
0920     hMuonCombIsoPt_->Fill(theMuonCombIsoPt);
0921     // hMuonCombIsoPt_->Fill((theMuonTrackIsoPt+theMuonEcalIsoPt)/theMuon->pt());

0922 
0923     // std::cout<<"Rel Iso Hist =

0924     // "<<(theMuonTrackIsoPt+theMuonEcalIsoPt)/theMuon->pt()<<std::endl;

0925     hTauEcalIsoPt_->Fill(theTauJet->isolationPFGammaCandsEtSum());
0926     hTauTrackIsoPt_->Fill(theTauJet->isolationPFChargedHadrCandsPtSum());
0927     hTauDiscrAgainstMuons_->Fill(theTauDiscrAgainstMuons);
0928     if (theTauJet->leadTrack().isAvailable())
0929       hTauLeadTrackPt_->Fill(theTauJet->leadTrack()->pt());
0930   }
0931 
0932   if ((cutFlowStatus == kPassedDeltaR) && (((theMuonTrackIsoPt < muonTrackIsoCut_) && (muonIsoMode_ == kAbsoluteIso)) ||
0933                                            ((1 > 0) && (muonIsoMode_ == kRelativeIso)))) {
0934     cutFlowStatus = kPassedMuonTrackIso;
0935     // isSelected = true;

0936   }
0937   if (cutFlowStatus == kPassedMuonTrackIso &&
0938       (((theMuonEcalIsoPt < muonEcalIsoCut_) && (muonIsoMode_ == kAbsoluteIso)) ||
0939        ((theMuonCombIsoPt < muonCombIsoCut_) && (muonIsoMode_ == kRelativeIso)))) {
0940     cutFlowStatus = kPassedMuonEcalIso;
0941     // isSelected = true;

0942   }
0943 
0944   if (cutFlowStatus == kPassedMuonEcalIso && theTauDiscrByTrackIso > 0.5) {
0945     cutFlowStatus = kPassedTauTrackIso;
0946   }
0947 
0948   if (cutFlowStatus == kPassedTauTrackIso && theTauDiscrByEcalIso > 0.5) {
0949     cutFlowStatus = kPassedTauEcalIso;
0950     isSelected = true;
0951   }
0952 
0953   for (int iCut = 1; iCut <= cutFlowStatus; ++iCut) {
0954     hCutFlowSummary_->Fill(iCut);
0955   }
0956 
0957   for (int iCut = 1; iCut <= cutFlowStatus; ++iCut) {
0958     hCutFlowSummary_->Fill(iCut);
0959   }
0960 
0961   //     }

0962 
0963   if (isSelected) {
0964     hVisMassFinal_->Fill(mMuTau);
0965     ++numEventsSelected_;
0966   }
0967 }
0968 
0969 void EwkMuTauHistManager::finalizeHistograms() {
0970   edm::LogInfo("EwkMuTauHistManager") << "Filter-Statistics Summary:" << std::endl
0971                                       << " Events analyzed = " << numEventsAnalyzed_ << std::endl
0972                                       << " Events selected = " << numEventsSelected_;
0973   if (numEventsAnalyzed_ > 0) {
0974     double eff = numEventsSelected_ / (double)numEventsAnalyzed_;
0975     edm::LogInfo("") << "Overall efficiency = " << std::setprecision(4) << eff * 100. << " +/- " << std::setprecision(4)
0976                      << TMath::Sqrt(eff * (1 - eff) / numEventsAnalyzed_) * 100. << ")%";
0977   }
0978 }
0979 
0980 //-------------------------------------------------------------------------------

0981 // common auxiliary functions used by different channels

0982 //-------------------------------------------------------------------------------

0983 
0984 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0985 
0986 #include <TMath.h>
0987 
0988 int getIsoMode(const std::string& isoMode_string, int& error) {
0989   int isoMode_int;
0990   if (isoMode_string == "absoluteIso") {
0991     isoMode_int = kAbsoluteIso;
0992   } else if (isoMode_string == "relativeIso") {
0993     isoMode_int = kRelativeIso;
0994   } else {
0995     edm::LogError("getIsoMode") << " Failed to decode isoMode string = " << isoMode_string << " !!";
0996     isoMode_int = kUndefinedIso;
0997     error = 1;
0998   }
0999   return isoMode_int;
1000 }
1001 
1002 double calcDeltaPhi(double phi1, double phi2) {
1003   double deltaPhi = phi1 - phi2;
1004 
1005   if (deltaPhi < 0.)
1006     deltaPhi = -deltaPhi;
1007 
1008   if (deltaPhi > TMath::Pi())
1009     deltaPhi = 2 * TMath::Pi() - deltaPhi;
1010 
1011   return deltaPhi;
1012 }
1013 
1014 double calcMt(double px1, double py1, double px2, double py2) {
1015   double pt1 = TMath::Sqrt(px1 * px1 + py1 * py1);
1016   double pt2 = TMath::Sqrt(px2 * px2 + py2 * py2);
1017 
1018   double p1Dotp2 = px1 * px2 + py1 * py2;
1019   double cosAlpha = p1Dotp2 / (pt1 * pt2);
1020 
1021   return TMath::Sqrt(2 * pt1 * pt2 * (1 - cosAlpha));
1022 }
1023 
1024 double calcPzeta(const reco::Candidate::LorentzVector& p1,
1025                  const reco::Candidate::LorentzVector& p2,
1026                  double pxMEt,
1027                  double pyMEt) {
1028   double cosPhi1 = cos(p1.phi());
1029   double sinPhi1 = sin(p1.phi());
1030   double cosPhi2 = cos(p2.phi());
1031   double sinPhi2 = sin(p2.phi());
1032   double zetaX = cosPhi1 + cosPhi2;
1033   double zetaY = sinPhi1 + sinPhi2;
1034   double zetaR = TMath::Sqrt(zetaX * zetaX + zetaY * zetaY);
1035   if (zetaR > 0.) {
1036     zetaX /= zetaR;
1037     zetaY /= zetaR;
1038   }
1039 
1040   double pxVis = p1.px() + p2.px();
1041   double pyVis = p1.py() + p2.py();
1042   double pZetaVis = pxVis * zetaX + pyVis * zetaY;
1043 
1044   double px = pxVis + pxMEt;
1045   double py = pyVis + pyMEt;
1046   double pZeta = px * zetaX + py * zetaY;
1047 
1048   return pZeta - 1.5 * pZetaVis;
1049 }
1050 
1051 bool passesElectronPreId(const reco::GsfElectron& electron) {
1052   if ((TMath::Abs(electron.eta()) < 1.479 || TMath::Abs(electron.eta()) > 1.566) &&  // cut ECAL barrel/endcap crack

1053       electron.deltaPhiSuperClusterTrackAtVtx() < 0.8 && electron.deltaEtaSuperClusterTrackAtVtx() < 0.01 &&
1054       electron.sigmaIetaIeta() < 0.03) {
1055     return true;
1056   } else {
1057     return false;
1058   }
1059 }
1060 
1061 bool passesElectronId(const reco::GsfElectron& electron) {
1062   if (passesElectronPreId(electron) && ((TMath::Abs(electron.eta()) > 1.566 &&  // electron reconstructed in ECAL

1063                                                                                 // endcap

1064                                          electron.sigmaEtaEta() < 0.03 && electron.hcalOverEcal() < 0.05 &&
1065                                          TMath::Abs(electron.deltaEtaSuperClusterTrackAtVtx()) < 0.009 &&
1066                                          TMath::Abs(electron.deltaPhiSuperClusterTrackAtVtx()) < 0.7) ||
1067                                         (TMath::Abs(electron.eta()) < 1.479 &&  // electron reconstructed in ECAL

1068                                                                                 // barrel

1069                                          electron.sigmaEtaEta() < 0.01 && electron.hcalOverEcal() < 0.12 &&
1070                                          TMath::Abs(electron.deltaEtaSuperClusterTrackAtVtx()) < 0.007 &&
1071                                          TMath::Abs(electron.deltaPhiSuperClusterTrackAtVtx()) < 0.8))) {
1072     return true;
1073   } else {
1074     return false;
1075   }
1076 }
1077 
1078 const reco::GsfElectron* getTheElectron(const reco::GsfElectronCollection& electrons,
1079                                         double electronEtaCut,
1080                                         double electronPtCut) {
1081   const reco::GsfElectron* theElectron = nullptr;
1082 
1083   for (reco::GsfElectronCollection::const_iterator electron = electrons.begin(); electron != electrons.end();
1084        ++electron) {
1085     if (TMath::Abs(electron->eta()) < electronEtaCut && electron->pt() > electronPtCut &&
1086         passesElectronPreId(*electron)) {
1087       if (theElectron == nullptr || electron->pt() > theElectron->pt())
1088         theElectron = &(*electron);
1089     }
1090   }
1091 
1092   return theElectron;
1093 }
1094 
1095 const reco::Muon* getTheMuon(const reco::MuonCollection& muons, double muonEtaCut, double muonPtCut) {
1096   const reco::Muon* theMuon = nullptr;
1097 
1098   for (reco::MuonCollection::const_iterator muon = muons.begin(); muon != muons.end(); ++muon) {
1099     if (TMath::Abs(muon->eta()) < muonEtaCut && muon->pt() > muonPtCut) {
1100       if (theMuon == nullptr || muon->pt() > theMuon->pt())
1101         theMuon = &(*muon);
1102     }
1103   }
1104 
1105   return theMuon;
1106 }
1107 
1108 const reco::PFTau* getTheTauJet(const reco::PFTauCollection& tauJets,
1109                                 double tauJetEtaCut,
1110                                 double tauJetPtCut,
1111                                 int& theTauJetIndex) {
1112   const reco::PFTau* theTauJet = nullptr;
1113   theTauJetIndex = -1;
1114 
1115   int numTauJets = tauJets.size();
1116   for (int iTauJet = 0; iTauJet < numTauJets; ++iTauJet) {
1117     const reco::PFTau& tauJet = tauJets.at(iTauJet);
1118 
1119     if (fabs(tauJet.eta()) < tauJetEtaCut && tauJet.pt() > tauJetPtCut) {
1120       if (theTauJet == nullptr || tauJet.pt() > theTauJet->pt()) {
1121         theTauJet = &tauJet;
1122         theTauJetIndex = iTauJet;
1123       }
1124     }
1125   }
1126 
1127   return theTauJet;
1128 }
1129 
1130 double getVertexD0(const reco::Vertex& vertex, const reco::BeamSpot& beamSpot) {
1131   double dX = vertex.x() - beamSpot.x0();
1132   double dY = vertex.y() - beamSpot.y0();
1133   return TMath::Sqrt(dX * dX + dY * dY);
1134 }
[ 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!